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

# Slide 4: The Weibull hazard
t <- seq(0, 2, len=299)
hw <- function(t, lam, gam) {lam*gam*(t*lam)^{gam-1}}
plot(t, hw(t, 1, 1), type='l', bty='n', ylab=expression(lambda(t)),
     lwd=3, col=pal(3)[1], ylim=c(0, 2), xaxt='n', yaxt='n')
axis(1, at=0:3, labels=c(0, expression(lambda), expression(2*lambda), expression(3*lambda)))
axis(2, at=0:2, labels=c(0, expression(lambda), expression(2*lambda)), las=1)
text(2, 1.1, expression(gamma==1))
lines(t, hw(t, 1, 0.5), lwd=3, col=pal(3)[2])
text(2, 0.45, expression(gamma==0.5), xpd=1)
lines(t, hw(t, 1, 1.5), lwd=3, col=pal(3)[3])
text(2, 2, expression(gamma==1.5), xpd=1)

# Diagnostic plots: Setup
survEst <- function(form, Data) {
  fit <- survfit(form, Data)
  fit.exp <- survreg(form, Data, dist='exponential')
  fit.wbl <- survreg(form, Data, dist='weibull')
  list(fit=fit, lam.e=exp(-coef(fit.exp)), lam.w=exp(-coef(fit.wbl)), gam=1/fit.wbl$scale)
}
dplot <- function(x, ...) {
  y <- log(-log(x$fit$surv))
  ind <- which(is.finite(y))
  y <- y[ind]
  tt <- x$fit$time[ind]
  plot(log(tt), y, col='gray', type='s', lwd=3, bty='n', las=1,
       xlab='Log(time)', ylab="log(-log S(t))", ...)
  abline(a=log(x$lam.e), b=1, col=pal(2)[1], lwd=3)
  abline(a=x$gam*log(x$lam.w), b=x$gam, col=pal(2)[2], lwd=3)
}
splot <- function(x, ...) {
  plot(x$fit, conf.int=FALSE, bty='n', las=1, mark.time=FALSE, col='gray',
       lwd=3, ...)
  tt <- seq(0, max(x$fit$time), len=99)
  lines(tt, pexp(tt, x$lam.e, low=FALSE), col=pal(2)[1], lwd=3)
  lines(tt, pweibull(tt, x$gam, 1/x$lam.w, low=FALSE), col=pal(2)[2], lwd=3)
}

# Slide 8: Pike rat
Data <- read.delim("https://s3.amazonaws.com/pbreheny-data-sets/Pike1966.txt")
x <- survEst(Surv(Time, Death) ~ 1, Data)
dplot(x)
splot(x, xlab="Time (days)", ylab="Survival")
toplegend(legend=c("Exponential", "Weibull"), lwd=3, col=pal(2))

# Slide 9: PBC
x <- survEst(Surv(time/365.25, status!=0) ~ 1, pbc)
dplot(x)
splot(x, xlab="Time (years)", ylab="Progression-free survival")
toplegend(legend=c("Exponential", "Weibull"), lwd=3, col=pal(2))

# Slide 10: GVHD
Data <- read.delim("https://s3.amazonaws.com/pbreheny-data-sets/gvhd.txt")
x <- survEst(Surv(Time/365.25, Status) ~ 1, Data)
dplot(x)
splot(x, xlab="Time (years)", ylab="P(No GVHD)")
toplegend(legend=c("Exponential", "Weibull"), lwd=3, col=pal(2))