Sta678 Handout #2: Estimating hazard functions 

library(survival)

#####################
## class example data
#####################
time<-c(1,4,2,6.5,3,5)
delta<-c(1,1,0,0,1,1)
cbind(time,delta)

## K-M estimator of the survival function
fitKM<-survfit(Surv(time,delta)~1, type='kaplan', conf.type='log-log')
fitKM
summary(fitKM)
plot(fitKM, xlab='time', ylab='survival probability (K-M)')

## Nelson-Aalen estimator of the survival function
fitNA<-survfit(Surv(time,delta)~1, type='fleming', conf.type='log-log')
summary(fitNA)
plot(fitNA)

## N-A estimator of the cumulative hazard function
fit<-coxph(Surv(time,delta)~1)
basehaz(fit)
# outputs two columns, hazard and time

time1<-c(0,basehaz(fit)[,2])
hzd<-c(0,basehaz(fit)[,1])
plot(time1,hzd,xlab='time', ylab='cumulative hazard',type='s',xlim=c(0,7))
#segments(5,1.25,7,1.25)
n<-dim(basehaz(fit))[1]
segments(basehaz(fit)[n,2], basehaz(fit)[n,1], 
       basehaz(fit)[n,2]+2, basehaz(fit)[n,1])

##########
# 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]

> table(trt)
trt
  1   2 
158 154 

# K-M curves comparing trted and untrted groups

fitKM<-survfit(Surv(time,delta)~trt,subset=(case<=312), type='kaplan',
 conf.type='log-log')
fitKM
summary(fitKM)
plot(fitKM, xlab='time (days)', ylab='survival probability')

# with legend
plot(fitKM, xlab='time (days)', ylab='survival probability',lty=c(1:2))
legend(500,0.3,c("D-penicillamine","placebo"),lty=c(1:2))