# R code for chapter 12  

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

### p(t), S(t), h(t), and hazard ratio for group 1 and group 2
### FIGURE 12.1


xe <- seq(0.00001,20,0.01)
ye <- dgamma(xe,10,2)
xebar <- seq(0.00001,20,0.01)
yebar <- dgamma(xe,20,2)
he <- ye/(1-pgamma(xe,10,2))
hebar <- yebar/(1-pgamma(xebar,20,2))
hr <- he/hebar

par(mfrow=c(2,2))
plot(xe,ye,type='l',xlab="",ylab="",ylim=c(0,0.32), lwd=2, cex=1.2)
lines(xebar,yebar,lty=2, lwd=2, cex=1.2)
text(5,0.3,"Group 1", lwd=2, cex=1.2)
text(10,0.22,"Group 2", lwd=2, cex=1.2)
mtext("t",line=3,side=1,cex=1.5)
mtext("f(t)",line=2.5,side=2,cex=1.5)

plot(xe,1-pgamma(xe,10,2),xlab="",ylab="",type='l',lwd=2, cex=1.5)
lines(xebar,1-pgamma(xebar,20,2),lty=2,lwd=2, cex=1.2)
text(2,0.28,"Group 1",lwd=2,cex=1.2)
text(16,0.2,"Group 2",lwd=2,cex=1.2)
mtext("t",line=3,side=1,cex=1.5)
mtext("S(t)",line=2.5,side=2,cex=1.5)

plot(xe,he,xlab="",ylab="",type='l',lwd=2, cex=1.5)
lines(xebar,hebar,lty=2,lwd=2, cex=1.2)
text(3,1,"Group 1",lwd=2,cex=1.2)
text(15,0.2,"Group 2",lwd=2,cex=1.2)
mtext("t",line=3,side=1,cex=1.5)
mtext("h(t)",line=2.5,side=2,cex=1.5)

x <- seq(0.5,20,0.01)
ye <- dgamma(x,10,2)
yebar <- dgamma(x,20,2)
he <- ye/(1-pgamma(x,10,2))
hebar <- yebar/(1-pgamma(x,20,2))
hr <- he/hebar

plot(x,hr,xlab="",ylab="",type='l',lwd=2, cex=1.5,ylim=c(0,50))
mtext("t",line=3,side=1,cex=1.5)
mtext(expression(h[1](t)/h[2](t)),line=2.5,side=2,cex=1.2)





### FIGURE 12.2.  HISTOGRAMS OF LEUKEMIA DATA
n1 <- 17
n2 <- 16
y1 <- c(65,156,100,134,16,108,121,4,39,143,56,26,22,1,1,5,65)
y2 <- c(56,65,17,7,16,22,3,4,2,3,8,4,3,30,4,43)
par(mfrow=c(2,1))
hist(y1,main="AG Positive",xlab="",ylab="",xlim=c(0,170),prob=F, ylim=c(0,10))
mtext("Survival time (weeks)",line=3,side=1,cex=1)
mtext("Count",line=2.5,side=2,cex=1)
hist(y2,main="AG Negative",xlab="",ylab="",xlim=c(0,170),prob=F,ylim=c(0,10))
mtext("Survival time (weeks)",line=3,side=1,cex=1)
mtext("Count",line=2.5,side=2,cex=1)







### Example 12.2.2
# Prior on theta1.  Mode=0.02.  95th percentile=0.15
b1 <- seq(0.001,50,0.01)
a1 <- 0.02*b1+1
cbind(a1,b1,qgamma(0.95,a1,b1))   #  Look for 0.15.  Answer is [2637,] 1.52722 26.361 1.500113e-01



# Prior on theta2.  Mode=0.69/20=0.0345.  95th percentile=0.69/5=0.138
b2 <- seq(0.001,50,0.01)
a2 <- 0.0345*b2+1
cbind(a2,b2,qgamma(0.95,a2,b2))   #  Look for 0.138.  Answer is [3796,] 2.309310 37.951    0.1380110




### FIGURE 12.3.  
### Prior and posterior of theta1 and theta2
n1 <- 17
n2 <- 16
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)
a1 <- 1.53
b1 <- 26.4
a2 <- 2.31
b2 <- 37.95

par(mfrow=c(2,1))
t <- seq(0,0.15,0.001)
plot(t,dgamma(t,a1+n1,rate=(b1+sum(t1))),type='l',xlab="",ylab="",main="AG Positive", lwd=3)
lines(t,dgamma(t,a1,rate=b1),lty=2, lwd=3)
mtext(expression(theta[1]),line=2.5,side=1,cex=1.5)

plot(t,dgamma(t,a2+n2,rate=(b2+sum(t2))),type='l',xlab="",ylab="",main="AG Negative", lwd=3)
lines(t,dgamma(t,a2,rate=b2),lty=2, lwd=3)
mtext(expression(theta[2]),line=2.5,side=1,cex=1.5)






### EXAMPLE 12.3.3
### Using R2WinBUGS to fit the 2 parameter exponential model and estimate survival functions
# 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\\SurvivalAnalysis\\Code", 
                 n.chains=1, n.iter=50000, n.thin=1, n.burnin=0)
round(print(leuk.fit),3)
attach.bugs(leuk.fit)





### Figure 12.4 
par(mfrow=c(1,2))

time <- seq(0,170,1)
plot(time,exp(-theta1[1]*time), type="l", xlab="", ylab="", cex=2, lwd=2)
mtext("t",line=3,side=1,cex=1.5)
mtext(expression(S[1](t)),line=2.5,side=2,cex=1.5)
for(i in 2:20){
	lines(time, exp(-theta1[i]*time), lty=1, cex=2, lwd=2)
}
### Plot mean survival curves
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)



### Figure 12.5.  Plot bands for the survival curves
surv1 <- apply(survcurves1,2,quantile,c(0.025,0.50,0.975))
surv2 <- apply(survcurves2,2,quantile,c(0.025,0.50,0.975))
plot(time,surv1[2,], type="l", xlab="", ylab="", lwd=3)
lines(time, surv1[1,], lty=1, lwd=3)
lines(time, surv1[3,], lty=1, lwd=3)
lines(time,surv2[2,],lty=2, lwd=3)
lines(time, surv2[1,], lty=2, lwd=3)
lines(time, surv2[3,], lty=2, lwd=3)
mtext("t",line=3,side=1,cex=1.7,lwd=3)
mtext("S(t)",line=2.5,side=2,cex=1.7,lwd=3)