source("http://myweb.uiowa.edu/pbreheny/7210/f19/notes/fun.R")
library(survival)

# Baseline (slide 19)
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 (slide 21)
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 (slide 22)
Data <- read.delim("https://s3.amazonaws.com/pbreheny-data-sets/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)