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

## PBC

# marginal
col <- pal(2)[2:1]
survdiff(Surv(time, status!=0) ~ hepato, pbc)
fit <- survfit(Surv(time, status!=0) ~ hepato, pbc)
plot(fit, mark.time=FALSE, lwd=3, col=col, bty='n', xscale=365.25, las=1, xlab="Time (years)", ylab="Pr(Survival)")
toplegend(legend=c("No", "Yes"), lwd=3, col=col)

# table
with(pbc, table(hepato, stage))

# strat
par(mfrow=c(2,2), oma=c(0,0,2,0))
for (i in 1:4) {
  fit <- survfit(Surv(time, status!=0) ~ hepato, pbc, subset=(stage==i))
  plot(fit, mark.time=FALSE, lwd=3, col=col, bty='n', xscale=365.25, las=1, xlab="Time (years)", ylab="Pr(Survival)", conf.int=FALSE)
  mtext(paste("Stage", i), line=3)
  if (i == 1) next()
  if (i == 2) toplegend(legend=c("No", "Yes"), lwd=3, col=col)
  Diff <- survdiff(Surv(time, status!=0) ~ hepato, pbc, subset=(stage==i))
  p <- pchisq(Diff$chisq, 1, low=FALSE)
  mtext(paste("p =", format.pval(p, digits=1, eps=1e-3)), line=1)
}

# stratified log-rank
Diff <- survdiff(Surv(time, status != 0) ~ hepato + strata(stage), pbc)
Tab <- cbind(Diff$obs[2,], Diff$exp[2,], Diff$obs[2,] - Diff$exp[2,])
Tab <- rbind(Tab, apply(Tab, 2, sum))
Tab

## GVHD
Data <- read.delim("http://myweb.uiowa.edu/pbreheny/data/gvhd.txt")

# LAF
survdiff(Surv(Time, Status) ~ Group, Data)
survdiff(Surv(Time, Status) ~ Group + strata(LAF), Data)

# Age: Strat
cAge <- cut(Data$Age, breaks=quantile(Data$Age, 0:3/3), include.lowest=TRUE)
a <- levels(cAge)
par(mfrow=c(1,3), oma=c(0,0,5,0))
for (i in 1:3) {
  fit <- survfit(Surv(Time, Status) ~ Group, Data, subset=(cAge==a[i]))
  plot(fit, mark.time=FALSE, xmax=80, lwd=3, col=col, bty='n', las=1, xlab="Time (days)", ylab="Pr(Survival)")
  mtext(paste("Age", a[i]), line=3)
  Diff <- survdiff(Surv(Time, Status) ~ Group, Data, subset=(cAge==a[i]))
  p <- pchisq(Diff$chisq, 1, low=FALSE)
  mtext(paste("p =", format.pval(p, digits=1, eps=1e-3)), line=1)
}
toplegend(legend=levels(Data$Group), lwd=3, col=col)

# Age: Strat log-rank
survdiff(Surv(Time, Status) ~ Group, Data)
survdiff(Surv(Time, Status) ~ Group + strata(cAge), Data)