#-*- R -*-## Script from Fourth Edition of `Modern Applied Statistics with S'# Chapter 13   Survival Analysislibrary(MASS)options(echo=T, width=65, digits=5, height=9999)options(contrasts=c("contr.treatment", "contr.poly"))postscript("ch13.ps", width=8, height=6, pointsize=9)library(survival)# 13.1  Estimators of survivor curvesplot(survfit(Surv(time) ~ ag, data=leuk), lty = 2:3, col = 2:3)legend(80, 0.8, c("ag absent", "ag present"), lty = 2:3, col = 2:3)attach(gehan)Surv(time, cens)plot(log(time) ~ pair)# product-limit estimators with Greenwood's formula for errors:gehan.surv <- survfit(Surv(time, cens) ~ treat, data = gehan,                      conf.type = "log-log")summary(gehan.surv)plot(gehan.surv, conf.int = T, lty = 3:2, log = T,     xlab = "time of remission (weeks)", ylab = "survival")lines(gehan.surv, lty = 3:2, lwd = 2, cex = 2)legend(25, 0.1 , c("control", "6-MP"), lty = 2:3, lwd = 2)detach()survdiff(Surv(time, cens) ~ treat, data = gehan)survdiff(Surv(time) ~ ag, data = leuk)# 13.2  Parametric modelsplot(gehan.surv,  lty = 3:4, col = 2:3, fun = "cloglog",     xlab = "time of remission (weeks)", ylab = "log H(t)")legend(2, 0.5, c("control","6-MP"), lty = 4:3, col = 3:2)survreg(Surv(time) ~ ag*log(wbc), leuk, dist = "exponential")summary(survreg(Surv(time) ~ ag + log(wbc), leuk, dist = "exponential"))summary(survreg(Surv(time) ~ ag + log(wbc), leuk)) # Weibullsummary(survreg(Surv(time) ~ ag + log(wbc), leuk,                 dist="loglogistic"))anova(survreg(Surv(time) ~ log(wbc), data = leuk),   survreg(Surv(time) ~ ag + log(wbc), data = leuk))summary(survreg(Surv(time) ~ strata(ag) + log(wbc), data=leuk))leuk.wei <- survreg(Surv(time) ~ ag + log(wbc), leuk)ntimes <- leuk$time * exp(-leuk.wei$linear.predictors)plot(survfit(Surv(ntimes)), log = T)survreg(Surv(time, cens) ~ factor(pair) + treat, gehan,          dist = "exponential")summary(survreg(Surv(time, cens) ~ treat, gehan, dist = "exponential"))summary(survreg(Surv(time, cens) ~ treat, gehan))plot(survfit(Surv(time, cens) ~ factor(temp), motors),      conf.int = F)motor.wei <- survreg(Surv(time, cens) ~ temp, motors)summary(motor.wei)unlist(predict(motor.wei, data.frame(temp=130), se.fit = T))predict(motor.wei, data.frame(temp=130), type = "quantile",        p = c(0.5, 0.1))t1 <-  predict(motor.wei, data.frame(temp=130),               type = "uquantile", p = 0.5, se = T)exp(c(LL=t1$fit - 2*t1$se, UL=t1$fit + 2*t1$se))t1 <-  predict(motor.wei, data.frame(temp=130),               type = "uquantile", p = 0.1, se = T)exp(c(LL=t1$fit - 2*t1$se, UL=t1$fit + 2*t1$se))# summary(censorReg(censor(time, cens) ~ treat, gehan))# 13.3  Cox proportional hazards modelattach(leuk)leuk.cox <- coxph(Surv(time) ~ ag + log(wbc), data = leuk)summary(leuk.cox)update(leuk.cox, ~ . -ag)(leuk.coxs <- coxph(Surv(time) ~ strata(ag) + log(wbc), data = leuk))(leuk.coxs1 <- update(leuk.coxs, . ~ . + ag:log(wbc)))plot(survfit(Surv(time) ~ ag), lty = 2:3, log = T)lines(survfit(leuk.coxs), lty = 2:3, lwd = 3)legend(80, 0.8, c("ag absent", "ag present"), lty = 2:3)leuk.cox <- coxph(Surv(time) ~ ag, leuk)detach()gehan.cox <- coxph(Surv(time, cens) ~ treat, gehan, method = "exact")summary(gehan.cox)# The next fit is slowcoxph(Surv(time, cens) ~ treat + factor(pair), gehan,        method = "exact")1 - pchisq(45.5 - 16.2, 20)(motor.cox <- coxph(Surv(time, cens) ~ temp, motors))coxph(Surv(time, cens) ~ temp, motors, method = "breslow")coxph(Surv(time, cens) ~ temp, motors, method = "exact")plot( survfit(motor.cox, newdata=data.frame(temp=200),               conf.type = "log-log") )summary( survfit(motor.cox, newdata = data.frame(temp=130)) )# 13.4  Further examples# VA.temp <- as.data.frame(cancer.vet)# dimnames(VA.temp)[[2]] <- c("treat", "cell", "stime",#     "status", "Karn", "diag.time","age","therapy")# attach(VA.temp)# VA <- data.frame(stime, status, treat = factor(treat), age,#     Karn, diag.time, cell = factor(cell), prior = factor(therapy))# detach(VA.temp)(VA.cox <- coxph(Surv(stime, status) ~ treat + age  + Karn +                 diag.time + cell + prior, data = VA))(VA.coxs <- coxph(Surv(stime, status) ~ treat + age + Karn +     diag.time + strata(cell) + prior, data = VA))par(mfrow=c(1,2), pty="s")plot(survfit(VA.coxs), log = T, lty = 1:4, col = 2:5)#legend(locator(1), c("squamous", "small", "adeno", "large"), lty = 1:4, col = 2:5)plot(survfit(VA.coxs), fun = "cloglog", lty = 1:4, col = 2:5)cKarn <- factor(cut(VA$Karn, 5))VA.cox1 <- coxph(Surv(stime, status) ~ strata(cKarn) + cell, data = VA)plot(survfit(VA.cox1), fun="cloglog")VA.cox2 <- coxph(Surv(stime, status) ~ Karn + strata(cell), data = VA)scatter.smooth(VA$Karn, residuals(VA.cox2))VA.wei <- survreg(Surv(stime, status) ~ treat + age + Karn +                  diag.time + cell + prior, data = VA)summary(VA.wei, cor = F)VA.exp <- survreg(Surv(stime, status) ~ Karn + cell,                  data = VA, dist = "exponential")summary(VA.exp, cor = F)cox.zph(VA.coxs)par(mfrow = c(3, 2), pty="m"); plot(cox.zph(VA.coxs))VA2 <- VA ## needed because VA and stepAIC are both in MASSVA2$Karnc <- VA2$Karn - 50VA.coxc <- update(VA.cox, ~ . - Karn + Karnc, data=VA2)VA.cox2 <- stepAIC(VA.coxc, ~ .^2)VA.cox2$anova(VA.cox3 <- update(VA.cox2, ~ treat/Karnc + prior*Karnc   + treat:prior + cell/diag.time))cox.zph(VA.cox3)par(mfrow = c(2, 2))plot(cox.zph(VA.cox3), var = c(1, 3, 7))par(mfrow = c(1, 1))#data(heart) # in package survivalcoxph(Surv(start, stop, event) ~ transplant*    (age + surgery + year), data = heart)(stan <- coxph(Surv(start, stop, event) ~ transplant*year +    age + surgery, data = heart))stan1 <- coxph(Surv(start, stop, event) ~ strata(transplant) +    year + year:transplant + age + surgery, heart)par(mfrow=c(1,2), pty="s")plot(survfit(stan1), conf.int = T, log = T, lty = c(1, 3), col = 2:3)#legend(locator(1), c("before", "after"), lty = c(1, 3), col= 2:3)attach(heart)plot(year[transplant==0], residuals(stan1, collapse = id),     xlab = "year", ylab = "martingale residual")lines(lowess(year[transplant == 0],             residuals(stan1, collapse = id)))par(mfrow = c(1,1), pty = "m")sresid <- resid(stan1, type = "dfbeta", collapse = id)detach()-100 * sresid %*% diag(1/stan1$coef)# Survivor curve for the "average" subjectsummary(survfit(stan))#  follow-up for two yearsstan2 <- data.frame(start = c(0, 183), stop= c(183, 2*365),    event = c(0, 0), year = c(4, 4), age = c(50, 50) - 48,    surgery = c(1, 1), transplant = c(0, 1))summary(survfit(stan, stan2, individual = T,                 conf.type = "log-log"))# Aids analysistime.depend.covar <- function(data) {  id <- row.names(data); n <- length(id)  events <- c(0, 10043, 11139, 12053) # julian days  crit1 <- matrix(events[1:3], n, 3 ,byrow = T)  crit2 <- matrix(events[2:4], n, 3, byrow = T)  diag <- matrix(data$diag,n,3); death <- matrix(data$death,n,3)  incid <- (diag < crit2) & (death >= crit1); incid <- t(incid)  indr <- col(incid)[incid]; indc <- row(incid)[incid]  ind <- cbind(indr, indc); idno <- id[indr]  state <- data$state[indr]; T.categ <- data$T.categ[indr]  age <- data$age[indr]; sex <- data$sex[indr]  late <- indc - 1  start <- t(pmax(crit1 - diag, 0))[incid]  stop <- t(pmin(crit2, death + 0.9) - diag)[incid]  status <- matrix(as.numeric(data$status),n,3)-1 # 0/1  status[death > crit2] <- 0; status <- status[ind]  levels(state) <- c("NSW", "Other", "QLD", "VIC")  levels(T.categ) <- c("hs", "hsid", "id", "het", "haem",                       "blood", "mother", "other")  levels(sex) <- c("F", "M")  data.frame(idno, zid=factor(late), start, stop, status,             state, T.categ, age, sex)}Aids3 <- time.depend.covar(Aids2)attach(Aids3)aids.cox <- coxph(Surv(start, stop, status)     ~ zid + state + T.categ + sex + age, data = Aids3)summary(aids.cox)aids1.cox <- coxph(Surv(start, stop, status)  ~ zid + strata(state) + T.categ + age, data = Aids3)(aids1.surv <- survfit(aids1.cox))plot(aids1.surv, mark.time = F, lty = 1:4, col = 2:5,     xscale = 365.25/12, xlab = "months since diagnosis")#legend(locator(1), levels(state), lty = 1:4, col = 2:5)aids2.cox <- coxph(Surv(start, stop, status)  ~ zid + state + strata(T.categ) + age, data = Aids3)(aids2.surv <- survfit(aids2.cox))par(mfrow = c(1, 2), pty="s")plot(aids2.surv[1:4], mark.time = F, lty = 1:4, col = 2:5,  xscale=365.25/12, xlab="months since diagnosis")#legend(locator(1), levels(T.categ)[1:4], lty = 1:4, col = 2:5)plot(aids2.surv[c(1, 5, 6, 8)], mark.time = F, lty = 1:4, col = 2:5,  xscale=365.25/12, xlab="months since diagnosis")#legend(locator(1), levels(T.categ)[c(1, 5, 6, 8)], lty = 1:4, col = 2:5)par(mfrow=c(1,1), pty="m")cases <- diff(c(0,idno)) != 0aids.res <- residuals(aids.cox, collapse = idno)scatter.smooth(age[cases], aids.res, xlab = "age",  ylab="martingale residual")age2 <- cut(age, c(-1, 15, 30, 40, 50, 60, 100))c.age <- factor(as.numeric(age2), labels = c("0-15", "16-30",  "31-40", "41-50", "51-60", "61+"))table(c.age)c.age <- relevel(c.age, "31-40")summary(coxph(Surv(start, stop, status) ~ zid  + state  + T.categ + age + c.age, data = Aids3))detach()make.aidsp <- function(){  cutoff <- 10043  btime <- pmin(cutoff, Aids2$death) - pmin(cutoff, Aids2$diag)  atime <- pmax(cutoff, Aids2$death) - pmax(cutoff, Aids2$diag)  survtime <- btime + 0.5*atime  status <- as.numeric(Aids2$status)  data.frame(survtime, status = status - 1, state = Aids2$state,    T.categ = Aids2$T.categ, age = Aids2$age, sex = Aids2$sex)}Aidsp <- make.aidsp()aids.wei <- survreg(Surv(survtime + 0.9, status) ~  state    + T.categ + sex + age, data = Aidsp)summary(aids.wei, cor = F)survreg(Surv(survtime + 0.9, status) ~ state + T.categ  + age, data = Aidsp)(aids.ps <- survreg(Surv(survtime + 0.9, status) ~  state   + T.categ + pspline(age, df=6), data = Aidsp))zz <- predict(aids.ps, data.frame(   state = factor(rep("NSW", 83), levels = levels(Aidsp$state)),   T.categ = factor(rep("hs", 83), levels = levels(Aidsp$T.categ)),   age = 0:82), se = T, type = "linear")plot(0:82, exp(zz$fit)/365.25, type = "l", ylim = c(0, 2),   xlab = "age", ylab = "expected lifetime (years)")lines(0:82, exp(zz$fit+1.96*zz$se.fit)/365.25, lty = 3, col = 2)lines(0:82, exp(zz$fit-1.96*zz$se.fit)/365.25, lty = 3, col = 2)rug(Aidsp$age+runif(length(Aidsp$age), -0.5, 0.5), ticksize = 0.015)# End of ch13