#***************************************************************#
# R code for the leukemia analysis in Ch 13.                    #
#***************************************************************#


library(survival)
library(xlsReadWrite)
library(R2WinBUGS)

## R2WinBUGS to fit the 2 group exponential model 
# Group 1 is AG+
# Group 2 is AG-
library(R2WinBUGS)
n <- c(17,16)
a1 <- 0.0001
b1 <- 0.0001
a2 <- 0.0001
b2 <- 0.0001
t1 <- c(65,156,100,134,16,108,121,4,39,143,56,26,22,1,1,5,65)
t2 <- c(56,65,17,7,16,22,3,4,2,3,8,4,3,30,4,43)
leukemiadata <- list("n","t1","t2","a1","b1","a2","b2")
parameters <- c("theta1","theta2","median1","median2","relmedian","S","Sdiff") 
inits <- list( list(theta1=0.05, theta2=0.02))

leuk.fit <- bugs(leukemiadata, inits, parameters,"leukemia.txt", 
                 working.directory = "H:\\MyDocuments\\BAYES\\SurvivalRegression\\Leukemia", 
                 n.chains=1, n.iter=50000, n.thin=1, n.burnin=0)
round(print(leuk.fit),3)
attach.bugs(leuk.fit)



## Plot mean survival curves from the Exp model
time <- seq(0,170,1)
survcurves1 <- exp(-theta1%*%t(time))
survmean1 <- apply(survcurves1,2,mean)
survcurves2 <- exp(-theta2%*%t(time))
survmean2 <- apply(survcurves2,2,mean)
plot(time,survmean1, type="l", xlab="", ylab="", lwd=3)
lines(time, survmean2, lty=2, lwd=3)
legend("topright", c("AG Positive","AG Negative"), lty=1:2, lwd=2, cex=1,bty="n")
mtext("t",line=3,side=1,cex=1.5)
mtext("S(t)",line=2.5,side=2,cex=1.5)


## Frequentist Cox analysis of leukemia data
n <- 33
y <- c(65,156,100,134,16,108,121,4,39,143,56,26,22,1,1,5,65,56,65,17,7,16,22,3,4,2,3,8,4,3,30,4,43)
g <- c(rep(1,17),rep(0,16))
status <- rep(1,33)
cox <- coxph(Surv(y, status) ~ g)
cox.zph(cox)
summary(cox)
agpos <- survfit(cox,list(g=1))
agneg <- survfit(cox,list(g=0))
plot(agpos$time, agpos$surv, type="l")
lines(agneg$time, agneg$surv,lty=1)




## Plot estimated (mean) survival curves from (i) Bayes Exp model, (ii) Bayes PH model, and (iii) Freq PH model
## FIGURE 13.5

timecox <- c(0,1,2,3,4,5,7,8,16,17,22,26,30,39,43,56,65,100,108,121,134,143,156)  

# Estimates from WinBUGS code 
bayesagneg <- c(0.995,0.905,0.8615,0.7319,0.5623,0.5219,0.4829,0.4422,0.3694,0.3311,0.2617,0.2294,0.1979,0.1692,0.1409,0.09246,0.03141,0.02115,
              0.01319,0.007724,0.004056,0.001636,3.62E-04)

bayesagpos <-  c(0.9985,0.9713,0.9575,0.9133,0.8465,0.8285,0.8102,0.79,0.7501,0.7268,0.6788,0.6533,0.6255,0.5972,0.5652,0.4966,0.3377,0.288,
                 0.2346,0.1828,0.1321,0.07805,0.02725)

# Figure 13.4
plot(time,survmean1, type="l", xlab="", ylab="", lwd=2)
lines(time, survmean2, lty=1, lwd=2)
lines(timecox,bayesagpos,lty=2,lwd=2)
lines(timecox,bayesagneg,lty=2,lwd=2)
lines(agneg$time, agneg$surv,lty=1)
lines(agpos$time, agpos$surv,lty=1)
mtext("t",line=3,side=1,cex=1.2)
mtext("S(t)",line=2.5,side=2,cex=1.2)
legend("topright",lty=c(1,1,2), c("Bayes Exponential","PH Frequentist","PH Bayes"),lwd=c(2,1,2))
text(125,0.3,"AG+")
text(17,0.2,"AG-")








####  Figure 13.5. Leukemia data.  Using output from SAS
betapostmean <- -1.32
hoAGneg <- c(0,0.088,0.2818,0.0537,0.0315,0.0504,0.029,0.037,0.0749,0.0571,0.2481,0.2481)
hoAGpos <- exp(betapostmean)*hoAGneg
timeho <- c(0,0,3.5,4.5,7.5,16.5,24,41,60.6,82.5,127.5,160)
plot(timeho,hoAGneg,type="s",xlab="t",ylab="Estimated Hazard",cex=2)
lines(timeho,hoAGpos,type="s", lty=2)
legend("top",lty=1:2,c("AG-","AG+"),bty="n")