# Setup/ read in GVHD data
source("http://myweb.uiowa.edu/pbreheny/7210/f15/notes/fun.R")
Data <- read.delim("http://myweb.uiowa.edu/pbreheny/data/gvhd.txt")

# Surv
require(survival)
S <- Surv(Data$Time, Data$Status)
class(S)
head(S)
head(S[,1])
head(S[,2])

# survfit
fit <- survfit(S~Data$Group)
fit
summary(fit, time=60)

# Kaplan-Meier curve
fit <- survfit(S~Data$Group)
plot(fit)
plot(fit, xmax=80)
plot(fit, xmax=80, mark.time=FALSE)
plot(fit, xmax=80, mark.time=FALSE, col=pal(2), lwd=3, bty='n', las=1,
     xlab='Time on study (Days)', ylab='Probability (GVHD-free)')

# Adding an at-risk table (manually)
par(mar=c(5,5,3,5))
plot(fit, xmax=80, mark.time=FALSE, col=pal(2), lwd=3, bty='n', las=1,
     xlab='Time on study (Days)', ylab='Probability (GVHD-free)')
tt <- seq(0, 80, 20)
Y <- rbind(sapply(tt, function(t) sum(Data$Time[Data$Group=="MTX"] >= t)),
           sapply(tt, function(t) sum(Data$Time[Data$Group=="MTX+CSP"] >= t)))
text(tt, 1.05, Y[1,], xpd=TRUE, cex=0.7)
text(tt, 1.1, Y[2,], xpd=TRUE, cex=0.7)
text(90, 1.05, 'MTX', xpd=TRUE, cex=0.7)
text(90, 1.1, 'MTX+CSP', xpd=TRUE, cex=0.7)