###############################################
## Author: Joshua M. Tebbs
## Date: 15 November 2021
## Update: 17 November 2023
## STAT 513 course notes: R Code
## Chapter 13
###############################################

# Figure 13.1
# Page 160
library(survival)
csp.mtx = c(3,8,10,12,16,17,22,64,65,77,82,98,155,189,199,247,324,356,378,408,411,
	420,449,490,528,547,691,769,1111,1173,1213,1357)
cens.csp.mtx = c(0,1,1,0,1,1,1,0,1,0,0,0,0,1,0,0,1,0,0,0,1,0,0,1,0,0,1,0,0,1,0,1)
mtx = c(9,11,12,20,20,22,25,25,25,28,28,31,35,35,46,49,104,106,156,218,230,231,316,
	393,395,428,469,602,681,690,1112,1180)
cens.mtx = c(0,1,1,0,1,1,1,0,1,1,1,1,0,0,1,0,0,0,1,1,0,0,0,1,0,0,1,1,0,1,0,1)
agvhd.times = c(csp.mtx,mtx)
cens.agvhd.times = c(cens.csp.mtx,cens.mtx)
treatment = c(rep(1,32),rep(2,32))
# KM estimates for each treatment group
fit = survfit(Surv(agvhd.times,cens.agvhd.times)~treatment)
# Plot KM estimates
plot(fit,xlab="Time to AGVHD diagnosis (in days)",ylab="Survival Probability",lty=c(1,4))
abline(h=0)
# Add legend
legend(900,0.875,lty = c(1,4),c("CSP+MTX","MTX only"))


# Figure 13.2
# Page 162
# Left
par(mar = c(4,4.4,4,2))
t = seq(0,12,0.01)
survivor = exp(-t/2)
plot(t,survivor,type="l",xlab="t",ylab=expression(S[T](t)),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
# Right
t = seq(0,12,0.01)
survivor.1 = exp(-t/2)
plot(t,survivor.1,type="l",xlab="t",ylab=expression(S[T](t)),cex.lab=1.25,lty=1)
lines(t,1-pexp(t,1),lty=4)
abline(h=0)
abline(v=0,lty=2)
legend(6,0.8,lty=c(1,4),c(expression(paste(beta, "=2")),expression(paste(beta, "=1"))))


# Figure 13.3
# Page 163
par(mar = c(4,4.4,4,2))
t = seq(0,1,0.001)
haz.1 = 3*t^2
haz.2 = 1.5*t^(0.5)
haz.3 = 1+t*0
haz.4 = 0.4*t^(-0.4)
plot(t,haz.1,type="l",xlab="t",ylab=expression(lambda[T](t)),lty=1,ylim=c(0,2),xaxp=c(0,1,1),yaxp=c(0,2,2),cex.lab=1.25)
lines(t,haz.2,lty=2)
lines(t,haz.3,lty=4)
lines(t,haz.4,lty=8)
abline(h=0)


# Example 13.4
# Page 165-166
turbo=c(1.6,2.0,2.6,3.0,3.5,3.9,4.5,4.6,4.8,5.0,5.1,5.3,5.4,5.6,5.8,6.0,6.0,6.1,6.3,6.5,
	6.5,6.7,7.0,7.1,7.3,7.3,7.3,7.7,7.7,7.8,7.9,8.0,8.1,8.3,8.4,8.4,8.5,8.7,8.8,9.0)
library(fitdistrplus) # need to install
fitdist(turbo,"weibull") # estimates Weibull model using ML


# Figure 13.4
# Page 166
# Left
t = seq(0,12,0.01)
survivor = exp(-(t/6.92)^(3.87))
plot(t,survivor,type="l",xlab="t",ylab="Estimated survivor function",cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
# Right
t = seq(0,12,0.01)
hazard = (3.87/6.92)*(t/6.92)^(2.87)
plot(t,hazard,type="l",xlab="t",ylab="Estimated hazard function",cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)


# Figure 13.5
# Page 169
t = seq(0,5000,1)
survivor.1 = 1-pexp(t,1/367.2)
plot(t,survivor.1,type="l",xlab="t",ylab="Estimated survivor function",cex.lab=1.25,lty=1)
lines(t,1-pexp(t,1/371.5),lty=2)
lines(t,1-pexp(t,1/914.5),lty=4)
abline(h=0)
abline(v=0,lty=2)
legend(3000,0.8,lty=c(1,2,4),c("Approach 1","Approach 2","Approach 3"))


# Figure 13.6
# Page 173
time = seq(0,10,1)
surv = c(1,0.813,0.681,0.509,0.426,0.417,0.393,0.347,0.325,0.264,0.158)
plot(time,surv,type="o",xlab="t",ylab="Estimated survivor function",cex.lab=1.25,ylim=c(0,1),pch=16,xaxp=c(0,10,10))
abline(h=0)

# Figure 13.7
# Page 175
library(survival)
surv.time = c(4.5,7.5,8.5,11.5,13.5,15.5,16.7,17.5,19.5,21.5)
cens = c(1,1,0,1,0,1,1,0,1,0) # 1 = failure; 0 = censored
fit = survfit(Surv(surv.time,cens)~1,se.fit=F,conf.type="none") # Use '~1' (only one survival curve)
# Left
plot(fit,xlab="Patient time (in years)",ylab="Estimated survivor function",xlim=c(0,25))
abline(h=0)
# Right
fit = survfit(Surv(surv.time,cens)~1,conf.type="plain")
plot(fit,xlab="Patient time (in years)",ylab="Estimated survivor function",xlim=c(0,25))


# Figure 13.8
# Page 176
# Left
t = seq(0,6,0.01)
pdf.1 = dgamma(t,3,2)
pdf.2 = dgamma(t,4,2)
plot(t,pdf.1,type="l",ylab="Probability density functions",xlab="Patient time",lty=1)
lines(t,pdf.2,lty=4)
abline(h=0)
legend(3,0.5,lty=c(1,4),c("Failure distribution","Censoring distribution"))
# Right
par(mar = c(4,4.4,4,2))
plot(t,1-pgamma(t,3,2),type="l",xlab="t",ylab=expression(S[T](t)),cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)


# Figure 13.9 
# Page 177
# Note: If you run this, the simulated data will change (I didn't save the seed)
B = 100 # number of patients
failure.time = rgamma(B,3,2) # T
cens.time = rgamma(B,4,2) # C
delta <- ifelse(failure.time <= cens.time,1,0) # Delta
obs.time = failure.time*delta+cens.time*(1-delta) # X
sim = data.frame(failure.time[1:5],cens.time[1:5],obs.time[1:5],delta[1:5])
sim
library(survival)
fit = survfit(Surv(obs.time,delta)~1,se.fit=F,conf.type="none")
plot(fit,xlab="Patient time (in years)",ylab="Survivor function")
t = seq(0,max(obs.time),0.01)
lines(t,1-pgamma(t,3,2),lty=4)
abline(h=0)
legend(1.5,0.9,lty=c(1,4),c("Kaplan-Meier estimate","True survivor function"))


# Figure 13.10
# Page 180
tongue = c(1,3,3,4,10,13,13,16,16,24,26,27,28,30,30,32,41,51,65,67,70,72,73,77,91,93,96,100,104,157,167,
61,74,79,80,81,87,87,88,89,93,97,101,104,108,109,120,131,150,231,240,400,
1,3,4,5,5,8,12,13,18,23,26,27,30,42,56,62,69,104,104,112,129,181,
8,67,76,104,176,231)
delta = c(rep(1,31),rep(0,21),rep(1,22),rep(0,6))
library(survival)
fit = survfit(Surv(tongue,delta)~1,conf.type="plain")
summary(fit)
plot(fit,xlab="Patient time (in weeks)",ylab="Estimated survivor function",xlim=c(0,208))


# Figure 13.11
# Page 187
csp.mtx = c(3,8,10,12,16,17,22,64,65,77,82,98,155,189,199,247,324,356,378,408,411,
	420,449,490,528,547,691,769,1111,1173,1213,1357)
cens.csp.mtx = c(0,1,1,0,1,1,1,0,1,0,0,0,0,1,0,0,1,0,0,0,1,0,0,1,0,0,1,0,0,1,0,1)
mtx = c(9,11,12,20,20,22,25,25,25,28,28,31,35,35,46,49,104,106,156,218,230,231,316,
	393,395,428,469,602,681,690,1112,1180)
cens.mtx = c(0,1,1,0,1,1,1,0,1,1,1,1,0,0,1,0,0,0,1,1,0,0,0,1,0,0,1,1,0,1,0,1)
agvhd.times = c(csp.mtx,mtx)
cens.agvhd.times = c(cens.csp.mtx,cens.mtx)
treatment = c(rep(1,32),rep(2,32))
# KM estimates for each treatment group
fit = survfit(Surv(agvhd.times,cens.agvhd.times)~treatment)
# Plot KM estimates
plot(fit,xlab="Time to AGVHD diagnosis (in days)",ylab="Survival Probability",lty=c(1,4))
abline(h=0)
# Add legend
legend(900,0.875,lty = c(1,4),c("CSP+MTX","MTX only"))
fit.1 = survfit(Surv(agvhd.times,cens.agvhd.times)~treatment,conf.int=0.95) # Summaries for each treatment
fit.1
fit.2 = survdiff(Surv(agvhd.times,cens.agvhd.times)~treatment) # Logrank test
fit.2