#***************************************************************#
# R code for the kidney data in Ch 13  of the Bayes Book.       #
# Data on 863 kidney transplant patients                        #
# Variables:                                                    #
#   Observation number                                          #
#   Time to death or on-study time                              #
#   Death indicator (0=alive, 1=dead)                           #
#   Sex (1=male, 2=female)                                      #
#   Race (1=white, 2=black)                                     #
#   Age in years                                                #
#***************************************************************#


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

kidneydata <- read.table("F:\\MyDocuments\\BAYES\\SurvivalRegression\\Kidney Data\\kidney.txt",header=T)
attach(kidneydata)
kidneydata
### Variables:  ID, time, delta, sex, race, age

kidneyoutput <- read.xls("F:\\MyDocuments\\BAYES\\SurvivalRegression\\Kidney Data\\KidneyOutput.xls",colNames=T,sheet="kidney")
attach(kidneyoutput)
kidneyoutput


timekidney <- c(0, 2, 3, 7, 10, 17, 21, 26, 28, 37, 40, 43, 44, 45, 50, 52, 56, 57, 59, 62, 68, 69, 78, 79, 88, 91, 97, 98, 104, 106, 119, 121, 135, 143, 150, 154, 158, 162, 190, 206, 209, 228, 229, 242, 248, 249, 252, 273, 291, 297, 311, 334, 340, 344, 346, 354, 366, 391, 402, 421, 439, 450, 470, 478, 481, 490, 495, 570, 583, 614, 615, 621, 652, 697, 730, 773, 776, 790, 793, 806, 840, 852, 864, 875, 929, 939, 943, 945, 946, 1001, 1013, 1016, 1105, 1164, 1186, 1191, 1196, 1210, 1275, 1326, 1331, 1357, 1384, 1388, 1418, 1473, 1509, 1734, 1777, 1820, 1835, 1877, 1940, 2034, 2056, 2108, 2171, 2276, 2291, 2301, 2313, 2369, 2414, 2421, 2489, 2557, 2567, 2650, 2795, 3146)

# timekidney <- c(0, unique(sort(time[delta==1])))



female <- sex-1
## if you want female to be baseline use male <- -1*(sex-2)

black <- race-1  
## if you want black to be baseline use white <- -1*(race-2)

mean(delta)
sum(delta)

sum(white & male)
mean(white & male)
sum(white==1 & male==1 & delta==1)

sum(black & male)
mean(black & male)
sum(black==1 & male==1 & delta==1)

sum(white & female)
mean(white & female)
sum(white==1 & female==1 & delta==1)

sum(black & female)
mean(black & female)
sum(black==1 & female==1 & delta==1)

median(age[white==1 & male==1])
sd(age[white==1 & male==1])

median(age[white==0 & male==1])
sd(age[white==0 & male==1])

median(age[white==1 & male==0])
sd(age[white==1 & male==0])

median(age[white==0 & male==0])
sd(age[white==0 & male==0])





#### FIGURE 13.6
# white female, average age: excel column 1 to 130
# white female age 2 sd above the mean:  excel column 131 to 260
wfave <- postmean[1:130]
wf2 <- postmean[131:260]
mean(age[female==1])  ## 42
mean(age[female==1]) + 2*sd(age[female==1])   ## 70

plot(timekidney,wfave,xlab="",ylab="",ylim=c(0,1),type="l",lwd=2)
lines(timekidney,wf2,lty=2,lwd=2)
legend("bottomright",lty=1:2, lwd=c(2,2), c("Age at transplant: 42 years", "Age at transplant: 70 years"))
mtext("t",line=3,side=1,cex=1.2)
mtext("S(t)",line=2.5,side=2,cex=1.2)















#### FIGURE 13.7
# average age, compare WM, BM, WF, BF
wfave <- postmean[1:130]
wmave <- postmean[261:390]
bfave <- postmean[391:520]
bmave <- postmean[521:650]


plot(c(0,timekidney),c(1,bfave),xlab="",ylab="",ylim=c(0.7,1),type="l",lwd=2)
lines(timekidney,bmave,lty=2,lwd=2)
lines(timekidney,wfave,lty=3,lwd=2)
lines(timekidney,wmave,lty=4,lwd=2)
legend("topright",lty=1:4,c("Black female","Black male","White female","White male"),lwd=2)
mtext("t",line=3,side=1,cex=1.2)
mtext("S(t)",line=2.5,side=2,cex=1.2)














#### FIGURE 13.8
## Estimated baseline hazard from Bayes PH model
tk1 <- c(timekidney[2:130],3146)
midp <- (timekidney+tk1)/2

h0 <- postmed[651:780]
#plot(timekidney,h0,type="s",xlab="t",ylab=expression(paste("Estimated  ", h[0](t))),cex=2)
#lines(lowess(midp,h0,f=0.05))
#mtext(expression(h[0](t)),line=2.5,side=2)



par(mfrow=c(2,1))
plot(timekidney[2:130],h0[2:130],type="s",xlab="t",ylab=expression(paste("Estimated  ", h[0](t))))
plot(lowess(midp[2:130],h0[2:130],f=0.1),type="l",xlab="t",ylab=expression(paste("Estimated  ", h[0](t))))
#mtext(expression(h[0](t)),line=2.5,side=2)
















## Frequentist Cox analysis of kidney data
cox1 <- coxph(Surv(time,delta) ~ age+female+black)
cox.zph(cox1)
summary(cox1)
g1 <- survfit(cox1,list(age=42,female=1,black=0))
g2 <- survfit(cox1,list(age=70,female=1,black=0))
plot(g1$time, g1$surv, type="l", lwd=2, xlab="t",ylab="S(t)",ylim=c(0,1))
lines(g2$time, g2$surv,lty=2, lwd=2)
lines(timekidney,wf2,lty=2,lwd=2)
lines(timekidney,wfave,lwd=2)



library(gbm)
h0freq <- basehaz.gbm(time,delta,f.x=cox1$linear.predictors,cumulative=FALSE)