# 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")
require(survival)

# Complimentary log-log
x <- seq(1e-2, 1-1e-2, len=299)
plot(x, log(x), type='l', lwd=3, col='slateblue', bty='n', las=1)
plot(x, log(-log(x)), type='l', lwd=3, col='slateblue', bty='n', las=1)

# By hand (time 9, MTX+CSP group)
fit <- survfit(Surv(Time, Status) ~ Group, Data, conf.type='plain')
summary(fit, time=c(8, 10, 16))
S <- 30/31
s <- sqrt(1/(31*30))
S + qnorm(c(.025,.975))*s*S
fit <- survfit(Surv(Time, Status) ~ Group, Data)
summary(fit, time=c(8, 10, 16))
S*exp(qnorm(c(.025,.975))*s)
fit <- survfit(Surv(Time, Status) ~ Group, Data, conf.type='log-log')
summary(fit, time=c(8, 10, 16))
S^exp(qnorm(c(.025,.975))*s/log(S))

# Plain/linear
fit1 <- survfit(Surv(Time, Status) ~ Group, Data, conf.type='plain')
plot(fit1, xmax=80, bty='n', xlab='Time (days)', ylab="Pr (GVHD-free)", las=1, col=pal(2), lwd=3, mark.time=FALSE)
ciband(fit1, col=pal(2, alpha=0.4))
toplegend(legend=levels(Data$Group), col=pal(2), lwd=3)

# Log transform
fit2 <- survfit(Surv(Time, Status) ~ Group, Data)
plot(fit2, xmax=80, bty='n', xlab='Time (days)', ylab="Pr (GVHD-free)", las=1, col=pal(2), lwd=3, mark.time=FALSE)
ciband(fit2, col=pal(2, alpha=0.4))
toplegend(legend=levels(Data$Group), col=pal(2), lwd=3)

# Log-log transform
fit3 <- survfit(Surv(Time, Status) ~ Group, Data, conf.type='log-log')
plot(fit3, xmax=80, bty='n', xlab='Time (days)', ylab="Pr (GVHD-free)", las=1, col=pal(2), lwd=3, mark.time=FALSE)
ciband(fit3, col=pal(2, alpha=0.4))
toplegend(legend=levels(Data$Group), col=pal(2), lwd=3)

# Median: Illustration
Data <- read.delim("http://myweb.uiowa.edu/pbreheny/data/Pike1966.txt")
fit <- survfit(Surv(Time, Death) ~ Group, subset(Data, Group==1))
plot(fit, bty='n', xlab='Time (days)', ylab="Survival probability", las=1, col=pal(2)[2], lwd=3, mark.time=FALSE, conf.int=FALSE)
ciband(fit, col=rgb(0.5, 0.5, 0.5, alpha=0.3))
lines(fit, col=pal(2)[2], lwd=3, mark.time=FALSE, conf.int=FALSE)
abline(h=0.5, lwd=2, lty=5)

# Median: By hand
fit$time[min(which(fit$lower < 0.5))]
fit$time[max(which(fit$upper > 0.5))+1]
fit