Sta678 Handout #3: Comparing two or more survival distributions
                   (logrank and related tests)

library(survival)

##########
# PBC data
##########
pbc<-matrix(scan('http://www-rohan.sdsu.edu/~jjfan/sta678/pbcS.txt'),byrow=T,ncol=20)
> dim(pbc)
[1] 418  20
case<-pbc[,1]
time<-pbc[,2]
status<-pbc[,3]
#status: 0=alive, 1=liver transplant, 2=dead
delta<-(status==2)
trt<-pbc[,4]
sex<-pbc[,6]

## logrank test comparing the trted and untrted groups
help(survdiff)

fit<-survdiff(Surv(time,delta)~trt, subset=(case<=312))
fit

# not significant, let's check the KM curves
fitKM<-survfit(Surv(time,delta)~trt,subset=(case<=312))
plot(fitKM, xlab='time (days)', ylab='survival probability',lty=c(1:2))
legend(500,0.3,c("D-penicillamine","placebo"),lty=c(1:2))

# what if we put more weight on earlier differences
fit<-survdiff(Surv(time,delta)~trt, subset=(case<=312),rho=1)
fit

#########################
# Remission duration data
#########################
data<-matrix(scan('http://www-rohan.sdsu.edu/~jjfan/sta678/remisNH.txt'),byrow=T,ncol=5)
pair<-data[,1]
rem<-data[,2]
timep<-data[,3]
timet<-data[,4]
statust<-data[,5]

## rearrage the remission data into 42 rows (one row per patient)
id<-1:42
pair<-rep(1:21,each=2)
remm<-rep(rem,each=2)
trt<-rep(1:0,21)
time<-as.vector(rbind(timet, timep))
status<-as.vector(rbind(statust,rep(1,21)))

remdata<-data.frame(id,pair,time,status,trt,remm)

# stratified logrank test
fit<-survdiff(Surv(time,status)~trt+strata(pair), data=remdata)
fit
> sqrt(10.7)
[1] 3.271085
#compared to 3.27 on page 22 of the text. The trt effect (adjusted for
#pair) is statistically significant. 

fit<-coxph(Surv(time,status)~trt+strata(pair), data=remdata)
fit$score