rm(list=ls())

library(rgl) 
#setwd("/Users/thomas/madrid/stat/566/examples")
source("http://www.stat.washington.edu/tsr/s566/labs/y0y1polytopenew-rgl-col.R")

par(mfrow=c(1,1))
#source("http://www.stat.washington.edu/tsr/s566/labs/y0y1polytopenew-rgl.R")
# The p's are stored in the order

# p(y0, x0 | z0), p(y0, x1 | z0), p(y1, x0 | z0), p(y1, x1 | z0), 
# p(y0, x0 | z1), p(y0, x1 | z1), p(y1, x0 | z1), p(y1, x1 | z1)

######################################
#uniform data from one arm to see general structure
unif <- rep(0.25,4)

simp <- do.simplex(phi=30,theta=120,r=1000,main="Rhombic Dodecahedron (Uniform Observed Dist)") # Set up plot
## Set up interactive plot
plot3d(c(1,0,0,1),c(0,1,0,0),c(0,0,1,0), ylab="%HU",xlab="%HE",zlab="%AR",par3d(FOV=1)) 
lines3d( c(1,0,0,1),c(0,1,0,0),c(0,0,1,0),col="red",lty=3)

do.polytope(unif,simp,red="blue", blue="red")
do.polytope.rgl(unif,red="blue", blue="red")


#######################################
## Lipid Data (from lecture)
lipid <- c(158/172,0,14/172,0,52/165,23/165,12/165,78/165)

## Check IV inequalities are obeyed
check.iv.ineqs(lipid)

simp <- do.simplex(phi=30,theta=120,r=1000,main="Lipid Data") # Set up plot
## Set up interactive plot
plot3d(c(1,0,0,1),c(0,1,0,0),c(0,0,1,0), ylab="%HU",xlab="%HE",zlab="%AR",par3d(FOV=1)) 
lines3d( c(1,0,0,1),c(0,1,0,0),c(0,0,1,0),col="red",lty=3)


do.polytope(lipid[1:4],simp,red="green")
do.polytope.rgl(lipid[1:4], red="green")

do.polytope(lipid[5:8],simp,red="purple")
do.polytope.rgl(lipid[5:8], red="purple")


do.polytope(lipid,simp)
do.polytope.rgl(lipid)

###### Finding bounds


## Set up interactive plot
plot3d(c(1,0,0,1),c(0,1,0,0),c(0,0,1,0), ylab="%HU",xlab="%HE",zlab="%AR",par3d(FOV=1)) 
lines3d( c(1,0,0,1),c(0,1,0,0),c(0,0,1,0),col="red",lty=3)

do.polytope(lipid,simp)
do.polytope.rgl(lipid)


## Let's add the upper and lower bound lines

ace.bnds.lipid <- ace.bounds(lipid)
ace.bnds.lipid

do.ace.line(ace.bnds.lipid[1], simp) # lower bound
do.ace.line(ace.bnds.lipid[2], simp) # upper bound

do.ace.line.rgl(ace.bnds.lipid[1]) # lower bound
do.ace.line.rgl(ace.bnds.lipid[2]) # upper bound

#####

######################################
#(2)
# Hirano flu data
flu <- c(0.07127430, 0.02159827, 0.73938085, 0.16774658,
            0.05706522, 0.02105978, 0.63519022, 0.28668478)
flu[1:4]
sum(flu[1:4])

flu[5:8]
sum(flu[5:8])

simp <- do.simplex(phi=30,theta=120,r=1000,main="Flu Data") # Set up plot
plot3d(c(1,0,0,1),c(0,1,0,0),c(0,0,1,0), ylab="%HU",xlab="%HE",zlab="%AR",par3d(FOV=1))  
lines3d( c(1,0,0,1),c(0,1,0,0),c(0,0,1,0),col="red",lty=3)


## Check IV inequalities

check.iv.ineqs(flu)


## Control Arm Polytope

do.polytope(flu[1:4],simp,red="green")
do.polytope.rgl(flu[1:4], red="green")

## Treatment Arm Polytope

do.polytope(flu[5:8],simp,red="purple")
do.polytope.rgl(flu[5:8], red="purple")

## Combined Polytope

do.polytope(flu,simp)
do.polytope.rgl(flu)

### Bounds analysis

## Let's look at the bounds for the ACE of X on Y

ace.bnds.flu <- ace.bounds(flu)
ace.bnds.flu

## Redo plots (but only with the combined polytope)

simp <- do.simplex(phi=30,theta=120,r=1000,main="Flu Data") # Set up plot
plot3d(c(1,0,0,1),c(0,1,0,0),c(0,0,1,0), ylab="%HU",xlab="%HE",zlab="%AR",par3d(FOV=1))  
lines3d( c(1,0,0,1),c(0,1,0,0),c(0,0,1,0),col="red",lty=3)

do.polytope(flu,simp)
do.polytope.rgl(flu)


## Let's add the upper and lower bound lines

do.ace.line.rgl(ace.bnds.flu[1]) # lower bound
do.ace.line.rgl(ace.bnds.flu[2])  # upper bound

# For comparison: add Intent-to-treat line to bottom face of cube
# ACE(Z->Y)
# Note the meaning of HE, HU etc., is now different:
# "HE" means "Y(z=0)=0, Y(z=1)=1"

itt.ace.flu <- do.itt.line(flu,simp)
itt.ace.flu
do.ace.line.rgl(itt.ace.flu)


####
####
# Variation independence  

## create random data:
tst.z0 <- runif(4)
tst.z1 <- runif(4)
tst <- c(tst.z0/sum(tst.z0),tst.z1/sum(tst.z1))

par(mfrow=c(1,2))
simp <- do.simplex(phi=30,theta=120,r=1000,main="Variation independence of margins")
do.polytope(tst,simp,red=rainbow(4)[1])

simp <- do.simplex(phi=-35.26439,theta=135,r=1000,main="Variation independence of margins") 
do.polytope(tst,simp,red=rainbow(4)[1])

# Interactive version
plot3d(c(1,0,0,1),c(0,1,0,0),c(0,0,1,0), ylab="%HU",xlab="%HE",zlab="%AR",par3d(FOV=1))  
lines3d( c(1,0,0,1),c(0,1,0,0),c(0,0,1,0),col="red",lty=3)
do.polytope.rgl(tst)

## Re-run this code a few times to see how 
## the shapes vary, but variation independence 
## remains




###########################
### Multiple levels of Z:
set.seed(511)
par(mfrow=c(1,1))
simp <- do.simplex(phi=30,theta=120,r=1000,main="Simulated Data") # Set up plot
## Set up interactive plot
plot3d(c(1,0,0,1),c(0,1,0,0),c(0,0,1,0), ylab="%HU",xlab="%HE",zlab="%AR",par3d(FOV=1)) 
lines3d( c(1,0,0,1),c(0,1,0,0),c(0,0,1,0),col="red",lty=3)


tst.z0 <- runif(4)
#do.polytope(tst.z0/sum(tst.z0),simp,red=rainbow(4)[1])
do.polytope.rgl(tst.z0/sum(tst.z0),red=rainbow(5)[2])


tst.z1 <- runif(4)
#do.polytope(tst.z1/sum(tst.z1),simp,red=rainbow(4)[2])
do.polytope.rgl(tst.z1/sum(tst.z1),red=rainbow(5)[3])


tst.z2 <- runif(4)
#do.polytope(tst.z2/sum(tst.z2),simp,red=rainbow(4)[3])
do.polytope.rgl(tst.z2/sum(tst.z2),red=rainbow(5)[4])


tst <- c(tst.z0/sum(tst.z0),tst.z1/sum(tst.z1),tst.z2/sum(tst.z2))
do.polytope(tst,simp,red="red")
do.polytope.rgl(tst,red="red")


## This is a bit busy. Here are fresh plots, just showing the intersection:

simp <- do.simplex(phi=30,theta=120,r=1000,main="Simulated Data") # Set up plot
do.polytope(tst,simp,red="red")

plot3d(c(1,0,0,1),c(0,1,0,0),c(0,0,1,0), ylab="%HU",xlab="%HE",zlab="%AR",par3d(FOV=1)) 
lines3d( c(1,0,0,1),c(0,1,0,0),c(0,0,1,0),col="red",lty=3)
do.polytope.rgl(tst,red="red")