# Setup/ read in GVHD data
source("http://myweb.uiowa.edu/pbreheny/7210/f19/notes/fun.R")
GVHD <- read.delim("https://s3.amazonaws.com/pbreheny-data-sets/gvhd.txt")
library(survival)

# Slide 11: 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)
x <- seq(-2, 2, len=299)
plot(x, exp(x), type='l', lwd=3, col='slateblue', bty='n', las=1)
plot(x, exp(-exp(x)), type='l', lwd=3, col='slateblue', bty='n', las=1)

# Slides 15-16: By hand (time 9, MTX+CSP group)
fit <- survfit(Surv(Time, Status) ~ Group, GVHD, 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.gvhd <- fit <- survfit(Surv(Time, Status) ~ Group, GVHD)
summary(fit, time=c(8, 10, 16))
S*exp(qnorm(c(.025,.975))*s)
fit <- survfit(Surv(Time, Status) ~ Group, GVHD, conf.type='log-log')
summary(fit, time=c(8, 10, 16))
S^exp(qnorm(c(.025,.975))*s/log(S))

# Slide 17: Plain/linear
fit1 <- survfit(Surv(Time, Status) ~ Group, GVHD, conf.type='plain')
Plot(fit1, xmax=80, xlab='Time (days)', ylab="Pr (GVHD-free)")
ciband(fit1)

# Slide 18: Log transform
fit2 <- survfit(Surv(Time, Status) ~ Group, GVHD)
Plot(fit2, xmax=80, xlab='Time (days)', ylab="Pr (GVHD-free)")
ciband(fit2)

# Slide 19: Log-log transform
fit3 <- survfit(Surv(Time, Status) ~ Group, GVHD, conf.type='log-log')
Plot(fit3, xmax=80, xlab='Time (days)', ylab="Pr (GVHD-free)")
ciband(fit3)

# Slide 22: survminer
library(survminer)
ggsurvplot(fit3, xlim=c(0, 60), break.x.by=20, conf.int = TRUE)

# Slide 26: Median illustration
Pike <- read.delim("https://s3.amazonaws.com/pbreheny-data-sets/Pike1966.txt")
fit.pike <- fit <- survfit(Surv(Time, Death) ~ Group, subset(Pike, Group==1))
Plot(fit, xlab='Time (days)')
ciband(fit)
abline(h=0.5, col='gray50', lwd=2, lty=2)

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

# Slide 28: Nelson-Aalen estimator
Plot(fit.gvhd, fun='cumhaz', xmax=80, ylab="Cumulative hazard")
ciband(fit.gvhd, fun=function(x) -log(x))
Plot(fit.pike, fun='cumhaz', ylab="Cumulative hazard")
ciband(fit.pike, fun=function(x) -log(x))