source("http://myweb.uiowa.edu/pbreheny/7210/f15/notes/fun.R") require(survival) fit <- coxph(Surv(time/365.25, status!=0) ~ trt + stage + hepato + bili, pbc) sfit <- survfit(fit) plot(sfit, conf.int=FALSE, mark.time=FALSE, lwd=3, bty='n', las=1, xlab='Time (years)', ylab='Progression-free survival') ciband(sfit, col=rgb(0.4, 0.4, 0.4, alpha=0.4)) apply(model.matrix(fit), 2, mean) # Subject-specific newdata <- data.frame(trt=0, stage=1:4, hepato=0, bili=1) sfit <- survfit(fit, newdata=newdata) plot(sfit, conf.int=FALSE, col=pal(4), mark.time=FALSE, lwd=3, bty='n', las=1, xlab='Time (years)', ylab='Progression-free survival') toplegend(legend=paste('Stage: ', 1:4), lwd=3, col=pal(4)) sfit # GVHD Data <- read.delim("http://myweb.uiowa.edu/pbreheny/data/gvhd.txt") fit <- coxph(Surv(Time, Status) ~ Group, Data) newdata <- data.frame(Group=c('MTX', 'MTX+CSP')) col1 <- hcl(seq(15, 375, len = 4), l=60, c=150)[c(1, 3)] col2 <- hcl(seq(15, 375, len = 4), l=85, c=50)[c(1, 3)] sfit1 <- survfit(fit, newdata=newdata) sfit2 <- survfit(Surv(Time, Status) ~ Group, Data) plot(sfit2, conf.int=FALSE, col=col1, mark.time=FALSE, lwd=3, bty='n', las=1, xlab='Time (days)', ylab='Progression-free survival', xmax=60) toplegend(legend=newdata$Group, lwd=3, col=col1) lines(sfit1, mark.time=FALSE, col=col2, lwd=3)