###########################################################################
# CHAPTER 04
###########################################################################

############
# Figure 4.2
############

age    <- 55:81
tv     <- c(1, 11) 
h.age  <- h.ageL <- h.ageU <- list(NA, NA)

for(i in 1:length(age)){
  
   for(j in 1:length(tv)){
     
    nd  <- data.frame(age = age[i], SevScore1 = 7, rs2284665 = 1)
    
    res <- haz.surv(outPL, eq = 1, newdata = nd, type = "haz", t.vec = tv[j],  
                   plot = FALSE, n.sim = 3000)
    
    h.age[[j]][i]  <- res$h
    h.ageL[[j]][i] <- res$CIh[1]; h.ageU[[j]][i] <- res$CIh[2]
                         }
}

ylimr <- range(c(unlist(h.ageL), unlist(h.ageU)))
plot(age, h.age[[1]], type = "l", ylim = ylimr, ylab = "hazard")
lines(age, h.ageU[[1]], lty = 2)
lines(age, h.ageL[[1]], lty = 2)

lines(age, h.age[[2]], col = "gray")
lines(age, h.ageU[[2]], lty = 2,col = "gray")
lines(age, h.ageL[[2]], lty = 2,col = "gray")

legend(x = 75, y = 0.06, legend = c("t = 1", "t = 11"), 
       fill = c("black", "gray"), bty = "n")


############
# Figure 4.3
############

rs   <- c(0, 2)
h.rs <- h.rsL <- h.rsU <- list(NA, NA)
tv   <- seq(1, 10, length.out = 20)

for(i in 1:length(rs)){
    nd  <- data.frame(age = 69, SevScore2 = 6, rs2284665 = rs[i])
    res <- haz.surv(outPL, eq = 2, newdata = nd, type = "haz", t.vec = tv, 
                   plot = FALSE, n.sim = 1000)
    h.rs[[i]]  <- res$h 
    h.rsL[[i]] <- res$CIh[, 1] 
    h.rsU[[i]] <- res$CIh[, 2]            
}

ylimr <- range(c(unlist(h.rsL), unlist(h.rsU)))
plot(tv, h.rs[[1]], type = "l", ylim = ylimr, ylab = "hazard", xlab = "time")
lines(tv, h.rsL[[1]], lty = 2)
lines(tv, h.rsU[[1]], lty = 2)

lines(tv, h.rs[[2]], col = "gray")
lines(tv, h.rsL[[2]], col = "gray", lty = 2)
lines(tv, h.rsU[[2]], col = "gray", lty = 2)

legend(x = 6, y = 0.05, legend = c("rs2284665 = 0", "rs2284665 = 2"), 
       fill = c("black", "gray"), bty = "n")


############
# Figure 4.5
############

par(mfrow = c(1, 2), mar = c(0, 1.5, 0, 1.5)) 
size <- 20
x <- y <- seq(0, 12, length = size)
t11 <- rep(x = x,  each = size)
t21 <- rep(x = y, times = size)

nd <- data.frame(t11 = t11, t21 = t21, age = 69, SevScore1 = 7, 
                 SevScore2 = 6, rs2284665 = 1) 
z  <- NA
for(i in 1:dim(nd)[1]) z[i] <- copula.prob(outPL, newdata = nd[i, ])$pr
z  <- matrix(z, nrow = size, byrow = TRUE) 

persp(x = x, y = y, z = z, box = TRUE, theta = 50, phi = 10, ticktype = "detailed",
      xlab = "years (left eye)", ylab = "years (right eye)", zlim = c(0, 1),
      zlab = "joint probabilities")

z  <- NA
for(i in 1:dim(nd)[1]) z[i] <- copula.prob(outPL, joint = FALSE, newdata = nd[i, ])$pr
z  <- matrix(z, nrow = size, byrow = TRUE) 

persp(x = x, y = y, z = z, box = TRUE, theta = 50, phi = 10, ticktype = "detailed",
      xlab = "years (left eye)", ylab = "years (right eye)", zlim = c(0, 1),
      zlab = "joint probabilities")


############
# Figure 4.6
############

ti <- seq(5, 12, length = 20)
nd <- data.frame(t11 = ti, t21 = 5, age = 69, SevScore1 = 7, 
                 SevScore2 = 6, rs2284665 = 1) 

z <- as.data.frame(matrix(NA, dim(nd)[1], 3))
for(i in 1:dim(nd)[1]) z[i, ] <- copula.prob(outPL, newdata = nd[i, ], 
                                             cond = 2, intervals = TRUE, 
                                             n.sim = 2000)

plot(ti, z[, 1], xlab = "years", type = "l", ylim = range(z), 
    ylab = "conditional survival function")
lines(ti, z[, 2], lty = 2)
lines(ti, z[, 3], lty = 2)



###########################################################################
## CHAPTER 05
###########################################################################

############
# Figure 5.4
############

res.BYc <- NA; nd$smoke <- 1

for(i in 1:lr){ 
 nd$county <- i  
res.BYc[i] <- copula.prob(out, y1 = 1, y2 = 1, newdata = nd)$pr*100 
              }
polys.map(NC.polys, res.BYc, scheme = "gray")

####################################################
# extra example to illustrate the calculation of the
# effect of a covariate on the outcome probability

nd <- nd1 <- data.frame(ethnicity = "White", educ = "Tertiary", county = 60,
                 firstbirth = 0, marital = 1, smoke = 0, age = 26)
nd1$smoke <- 1
 
s0 <- marg.mv(out, eq = 1, newdata = nd , n.sim = 2000)
s1 <- marg.mv(out, eq = 1, newdata = nd1, n.sim = 2000)

s1$res["mean"] - s0$res["mean"]

Sim.difm <- s1$sim.mv - s0$sim.mv           # Simulated distribution of differences
quantile(Sim.difm, probs = c(0.025, 0.975)) # if 0 not included then difference is significant


###########################################################################
## CHAPTER 13
###########################################################################

#############
# Figure 13.1
#############

par(mfrow = c(1, 2))
t.grid <- seq(1, 25, by = 3)

SATE(out, trt = "agree", surv.t = t.grid,  
     plot.type = "survival", int.var = c("agree:gender", 0), 
     xlab = "Time in weeks", ylab = "Survival functions")
     
SATE(out, trt = "agree", surv.t = t.grid, 
     plot.type = "sate", int.var = c("agree:gender", 0), 
     xlab = "Time in weeks", ylab = "SATEs")