# PubH 7450
#   Comparing KM estimate with Nelson-Aalen estimator; see Fig 4.1A on p.95
#
leuk<-read.table("/home/merganser/weip/public_html/course/data/leukemia.dat",
                 header=T)
leuk[1:2,]
#  group weeks relapse
#1  6-MP     6      1
#2  6-MP     6      1

library(survival)

km.fit<-survfit(Surv(weeks, relapse) ~ 1, data=leuk, subset=(group=='6-MP'),
                 type='kaplan-meier')
#in Surv(time, event): event=0 for right-censoring; =1 for an event.
# read the manual by ?Surv 

na.fit<-survfit(Surv(weeks, relapse) ~ 1, data=leuk, subset=(group=='6-MP'),
                 type='fleming-harrington')
pdf("ex4.1b.pdf")
# plot K-M estimate with pointwise 95% CI (log-transformed CI):
plot(km.fit, xlab="Time", ylab="Survival")
#draw N-A estimate:
lines(na.fit, lty=3, col="red")
dev.off()


########in survfit, there is a parameter 
########conf.type=c("log","log-log","plain","none")
########what is the default?
km.fit$conf.type
#[1] "log"

na.fit$conf.type
#[1] "log"

##print out CI's:
 km.fit$lower
# [1] 0.7198171 0.6531242 0.6531242 0.5859190 0.5859190 0.5096131 0.4393939
# [8] 0.4393939 0.4393939 0.4393939 0.3370366 0.2487882 0.2487882 0.2487882
#[15] 0.2487882 0.2487882
 km.fit$upper
# [1] 1.0000000 0.9964437 0.9964437 0.9675748 0.9675748 0.9347692 0.8959949
# [8] 0.8959949 0.8959949 0.8959949 0.8582008 0.8073720 0.8073720 0.8073720
#[15] 0.8073720 0.8073720
 na.fit$lower
# [1] 0.7279924 0.6617329 0.6617329 0.5950244 0.5950244 0.5194397 0.4498420
# [8] 0.4498420 0.4498420 0.4498420 0.3489697 0.2616612 0.2616612 0.2616612
#[15] 0.2616612 0.2616612
 na.fit$upper
# [1] 1.0000000 1.0000000 1.0000000 0.9826112 0.9826112 0.9527938 0.9173001
# [8] 0.9173001 0.9173001 0.9173001 0.8885862 0.8491476 0.8491476 0.8491476
#[15] 0.8491476 0.8491476

########compared with linear (or plain) CIs:
km.fit2<-survfit(Surv(weeks, relapse) ~ 1, data=leuk, subset=(group=='6-MP'),
                 type='kaplan-meier', conf.type="plain")
 km.fit2$lower
# [1] 0.7074793 0.6363327 0.6363327 0.5640993 0.5640993 0.4808431 0.4039095
# [8] 0.4039095 0.4039095 0.4039095 0.2864816 0.1843849 0.1843849 0.1843849
#[15] 0.1843849 0.1843849
 km.fit2$upper
# [1] 1.0000000 0.9771127 0.9771127 0.9417830 0.9417830 0.8995491 0.8509924
# [8] 0.8509924 0.8509924 0.8509924 0.7891487 0.7119737 0.7119737 0.7119737
#[15] 0.7119737 0.7119737

########compared with linear (or plain) CIs, log-log-transformed CIs:
km.fit3<-survfit(Surv(weeks, relapse) ~ 1, data=leuk, subset=(group=='6-MP'),
                 type='kaplan-meier', conf.type="log-log")
 km.fit3$lower
# [1] 0.6197180 0.5631466 0.5631466 0.5031995 0.5031995 0.4316102 0.3675109
# [8] 0.3675109 0.3675109 0.3675109 0.2677789 0.1880520 0.1880520 0.1880520
#[15] 0.1880520 0.1880520
 km.fit3$upper
# [1] 0.9515517 0.9228090 0.9228090 0.8893618 0.8893618 0.8490660 0.8049122
# [8] 0.8049122 0.8049122 0.8049122 0.7467907 0.6801426 0.6801426 0.6801426
#[15] 0.6801426 0.6801426