### R code from vignette source 'survival' ################################################### ### code chunk number 1: setup ################################################### library(pgfSweave) setCacheDir("cache") options(keep.source=TRUE) ################################################### ### code chunk number 2: initseed ################################################### set.seed(123) ################################################### ### code chunk number 3: loadlibrary ################################################### search() library(survival) ################################################### ### code chunk number 4: surv ################################################### data(leukemia) head(leukemia) help(Surv) mysurv <- Surv(leukemia$time, leukemia$status) head(mysurv) ################################################### ### code chunk number 5: survival ################################################### leuk.km <- survfit(Surv(time, status) ~ x, data=leukemia) plot(leuk.km, lty=1, col=c("darkblue","darkred")) legend("topright", legend=c('Maintain', 'Non-main'), lty=1:2, col=c("darkblue","darkred")) ################################################### ### code chunk number 6: Survival_conf ################################################### leuk.km2 <- survfit(Surv(time, status) ~ x, data=leukemia, conf.type='log-log') summary(leuk.km2) plot(leuk.km2, mark.time=FALSE, conf.int=TRUE, lty=1, col=c("darkblue","darkred")) legend("topright", legend=c('Maintain', 'Non-main'), lty=1, col=c("darkblue","darkred")) ################################################### ### code chunk number 7: survivalTest ################################################### survdiff(Surv(time, status) ~ x, data=leukemia) ################################################### ### code chunk number 8: survivalkm (eval = FALSE) ################################################### ## colon.km <- survfit(Surv(time, status) ~ rx, data=colon) ## plot(colon.km, lty=1, col=c("darkblue", "darkgreen", "darkred")) ## legend("topright", legend=c("Obs", "Lev", "Lev+5FU"), lty=1, col=c("darkblue", "darkgreen", "darkred")) ## ## logrank test ## survdiff(formula = Surv(time, status) ~ rx, data = colon) ################################################### ### code chunk number 9: CoxProp ################################################### leuk.ph <- coxph(Surv(time, status) ~ x, data=leukemia) summary(leuk.ph) ################################################### ### code chunk number 10: survivalcoxph (eval = FALSE) ################################################### ## summary(coxph(Surv(time, status) ~ ., data=colon)) ################################################### ### code chunk number 11: installsurvcomp (eval = FALSE) ################################################### ## source("http://www.bioconductor.org/biocLite.R") ## biocLite("survcomp") ################################################### ### code chunk number 12: loadsurvcomp ################################################### library(survcomp) ################################################### ### code chunk number 13: loadbcdata ################################################### data(breastCancerData) ################################################### ### code chunk number 14: bcataset ################################################### print(vdx7g) print(transbig7g) ################################################### ### code chunk number 15: vdx7gsurvdata ################################################### ## clincal information print(head(phenoData(vdx7g)@data)) dd <- phenoData(vdx7g)@data[ ,c("t.dmfs", "e.dmfs")] colnames(dd) <- c("time", "event") ## append gene expressions ge <- t(exprs(vdx7g)) dd <- cbind(dd, ge) ################################################### ### code chunk number 16: vdxfitcox ################################################### mm <- coxph(Surv(time, event) ~ ., data=data.frame(dd)) print(summary(mm)) predtr <- predict(mm, newdata=data.frame(dd), type="risk") ################################################### ### code chunk number 17: transbigvalidatecox ################################################### dd2 <- t(exprs(transbig7g)) predts <- predict(mm, newdata=data.frame(dd2), type="risk") ################################################### ### code chunk number 18: transbigsurvd ################################################### survdd2 <- phenoData(transbig7g)@data[ ,c("t.dmfs", "e.dmfs")] colnames(survdd2) <- c("time", "event") ################################################### ### code chunk number 19: hr ################################################### perf <- hazard.ratio(x=predts, surv.time=survdd2[ ,"time"], surv.event=survdd2[ ,"event"], na.rm=TRUE) print(perf[1:6]) ################################################### ### code chunk number 20: dindex ################################################### perf <- D.index(x=predts, surv.time=survdd2[ ,"time"], surv.event=survdd2[ ,"event"], na.rm=TRUE) print(perf[1:6]) ################################################### ### code chunk number 21: cindex ################################################### perf <- concordance.index(x=predts, surv.time=survdd2[ ,"time"], surv.event=survdd2[ ,"event"], method="noether", na.rm=TRUE) print(perf[1:5]) ################################################### ### code chunk number 22: tdrocc ################################################### perf <- tdrocc(x=predts, surv.time=survdd2[ ,"time"], surv.event=survdd2[ ,"event"], time=5*365, na.rm=TRUE) plot(x=1-perf$spec, y=perf$sens, type="l", xlab="1 - specificity", ylab="sensitivity", xlim=c(0, 1), ylim=c(0, 1), main="Time-dependent ROC curve\nat 5 years") lines(x=c(0,1), y=c(0,1), lty=3, col="red") ################################################### ### code chunk number 23: sbrier ################################################### ddtr <- cbind("time"=dd[ ,"time"], "event"=dd[ ,"event"], "score"=predtr) ddts <- cbind("time"=survdd2[ ,"time"], "event"=survdd2[ ,"event"], "score"=predts) perf <- sbrier.score2proba(data.tr=data.frame(ddtr), data.ts=data.frame(ddts), method="cox") plot(x=perf$time, y=perf$bsc, xlab="Time (days)", ylab="Brier score", type="l") ## null model ddtr <- cbind("time"=dd[ ,"time"], "event"=dd[ ,"event"], "score"=1) ddts <- cbind("time"=survdd2[ ,"time"], "event"=survdd2[ ,"event"], "score"=1) perfnull <- sbrier.score2proba(data.tr=data.frame(ddtr), data.ts=data.frame(ddts), method="cox") lines(x=perfnull$time, y=perfnull$bsc, col="red", lty=2) ## legend legend("bottomright", title="Integrated Brier score", legend=c(sprintf("Cox model = %.3g", perf$bsc.integrated), sprintf("Null model = %.3g", perfnull$bsc.integrated)), col=c("black", "red"), lty=c(1,2))