#
# This file contains R commands that reproduce the analyses 
# presented in Chapter 8
#
# Simple example of a survival curve for uncensored data (Figure 8.1)
#
x<-c(0.9,1.3,1.5,3.5,4.9)
pdf(file="survival2.pdf",height=6,width=8)
plot(x,1-(1:5)/5,xlim=c(0,6),ylim=c(0,1),xlab="t",ylab="S(t)",cex=2,cex.axis=1.5,cex.lab=1.5)
xx<-matrix(x,2,5,byrow=T)
xx<-c(0,c(xx),6)
yy<-1-(0:5)/5
yy<-c(matrix(yy,2,6,byrow=T))
lines(xx,yy)
dev.off()
#
# Load survival package
#
library(survival)
#
# Simple example of a survival curve for a mix of uncensored and censored data (Figure 8.2)
#
pdf(file="survival3.pdf",height=6,width=8)
observed<-c(1,0,1,0,1)
survival.object<-Surv(x,observed)
survival.curve<-survfit(survival.object~1)
plot(survival.curve,conf.int=F,
   xlim=c(0,6),ylim=c(0,1),xlab="t",ylab="S(t)",cex=2,cex.axis=1.5,cex.lab=1.5)
points(x[c(1,3,5)],c(0.8,0.53333,0),cex=2)
dev.off()
#
# Read kidney failure data from file "dialysis.data" and check names of variables
#
dialysis<-read.table(file="../data/dialysis.data",header=T)
names(dialysis) 
#
#     id: patient identifier (1 to 124)
#   days: number of days on treatment
#   dead: indicator 1 if patient died, 0 if still alive 
# method: indicator 0 if treated with CAPD, 1 if treated with APD
#    age: age, in years, at start of treatment
#
# Use variables days and dead to create
# a set of observed or censored survuival times
#
event<-Surv(dialysis$days,dialysis$dead)
#
# Calculate and plot Kaplan-Meier estimates of survival functions in
# the two treatment groups (Figure 8.3)
#
pdf(file="renal_survival.pdf",paper="special",width=6,height=4)
par(mfrow=c(1,2)) # draws two diagrams per page, side-by-side
capd<-dialysis$method==0
km.capd<-survfit(event[capd,]~1)
plot(km.capd,lty=c(1,2),xlab="t (days)",ylab="S(t)")
apd<-dialysis$method==1
km.apd<-survfit(event[apd,]~1)
plot(km.apd,lty=c(1,2),xlab="t (days)",ylab="S(t)")
dev.off()
#
# An illustration of the distinction between unconditional and conditional 
# life-time distributions (Figure 8.4)
#
mu<-78.5
a<-15
b<-(105*a/mu)-a
x<-c(0.1*(0:1050))
y<-dbeta(x/105,a,b)
pdf(file="lifetime_distribution.pdf",height=6,width=8)
plot(x,y,type="l",lwd=2,cex.lab=1.5,cex.axis=1.5,xlab="t (life-time in years)",ylab="relative frequency")
lines(c(60,60),c(0,y[601]),lwd=2,lty=2)
dev.off()
#
# Hazard functions corresponding to unconditional and conditional 
# life-time distributions (Figure 8.5)
#
yc<-c(rep(0,600),y[601:1051])*(sum(y)/sum(y[601:1051]))
S<-1-cumsum(y)/sum(y)
Sc<-1-cumsum(yc)/sum(yc)
pdf(file="lifetime_survival.",height=6,width=8)
plot(x[500:1051],S[500:1051],type="l",lwd=2,cex.lab=1.5,cex.axis=1.5,xlab="t (life-time in years)",ylab="S(t)")
lines(x[500:1051],Sc[500:1051],lty=2,lwd=2)
lines(c(50,105),c(1,1),lty=3)
dev.off()
#
# An increasing hazard function and its corresponding life-time distribution (Figure 8.6)
#
t<-0.1*(0:1000)
alpha<-1/3000000
h1<-alpha*t^3
S1<-exp(-alpha*t^4/4)
f1<-h1*S1
pdf(file="hazard1.pdf",paper="special",width=6,height=4)
par(mfrow=c(1,2))
plot(t,h1,type="l",lwd=2,cex.lab=1.3,cex.axis=1.3,xlab="t",ylab="h(t)")
plot(t,f1,type="l",lwd=2,cex.lab=1.3,cex.axis=1.3,xlab="t",ylab="f(t)")
dev.off()
#
# A non-monotone hazard function and its corresponding life-time distribution (Figure 8.7)
#
t<-0.1*(0:1000)
beta<-1/50
h2<-beta*((t-15)^2/225)
S2<-exp(-beta*((t-15)^3/675))
f2<-h2*S2
pdf(file="hazard2.pdf",paper="special",width=6,height=4)
par(mfrow=c(1,2))
plot(t,h2,type="l",lwd=2,cex.lab=1.3,cex.axis=1.3,xlab="t",ylab="h(t)")
plot(t,f2,type="l",lwd=2,cex.lab=1.3,cex.axis=1.3,xlab="t",ylab="f(t)")
dev.off()
#
# A constant hazard function and its corresponding life-time distribution (Figure 8.8)
#
gamma<-1/60
h3<-rep(gamma,1001)
S3<-exp(-gamma*t)
f3<-h3*S3
pdf(file="hazard3.pdf",paper="special",width=6,height=4)
par(mfrow=c(1,2))
plot(t,h3,type="l",lwd=2,cex.lab=1.3,cex.axis=1.3,xlab="t",ylab="h(t)",ylim=c(0,0.04))
plot(t,f3,type="l",lwd=2,cex.lab=1.3,cex.axis=1.3,xlab="t",ylab="f(t)")
dev.off()
#
# Kaplan-Meier plots of kidney failure data (Figure 8.3) suggested better survival under APD, 
# now fit proportional hazards model to confirm (Section 8.5)
#
fit1<-coxph(event~method,data=dialysis)
summary(fit1)
#
# Difference between CAPD and APD us marginally significant (p-value 0.055)
#
# Now extend model to see if age also affects outcome
#
fit2<-coxph(event~method+age,data=dialysis)
summary(fit2)
#
# Now, both method and age are highly significant (p-values < 0.001)
#