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

# Lognormal AFT: Same as lm if no censoring
fit <- survreg(Surv(time)~trt + stage + hepato + bili, pbc, dist='lognormal')
coef(fit)
fit.lm <- lm(log(time) ~ trt + stage + hepato + bili, pbc)
coef(fit.lm)

# Slide 5: Lognormal hazard
tt <- exp(seq(-2,3, len=99))
plot(tt, dlnorm(tt, 1)/plnorm(tt, 1, low=FALSE), type='l', bty='n', las=1,
     xlab='t', ylab=expression(lambda(t)), lwd=3)

# The arbitrary distribution for slides 7-10
t <- seq(0, 1, len=99)
h <- function(tt) {exp(-3*(tt-0.5)*dnorm(tt-0.5, sd=0.2))}
H <- function(tt) {
  val <- numeric(length(tt))
  for (i in 1:length(tt)) val[i] <- integrate(h, 0, tt[i])$value
  val
}
S <- function(tt) exp(-H(tt))
col <- c('gray', pal(2))

# Slide 7: Survivial
plot(t, S(t), type='l', ylim=c(0,1), bty='n', lwd=3, col=col[1], ylab='S(t)', las=1)
lines(t, S(t*3/2), lwd=3, col=col[2])
lines(t, S(t*2/3), lwd=3, col=col[3])
toplegend(legend=c('Baseline', expression(e^eta[i]==2), expression(e^eta[i]==1/2)), lwd=3, col=col)

# Slide 8: Hazard
plot(t, h(t), type='l', ylim=c(0,3), bty='n', lwd=3, col=col[1],
     ylab=expression(lambda(t)), las=1)
lines(t, h(1.5*t)*1.5, lwd=3, col=col[2])
lines(t, h(t/1.5)/1.5, lwd=3, col=col[3])
toplegend(legend=c('Baseline', expression(e^eta[i]==3/2), expression(e^eta[i]==2/3)), lwd=3, col=col)

# Slide 9: PH vs AFT (hazard)
plot(t, h(t), type='l', bty='n', ylim=c(0,2), lwd=3, col=col[1],
     ylab=expression(lambda(t)), las=1)
lines(t, h(t)/1.5, col=col[2], lwd=3)
lines(t, h(t/1.5)/1.5, col=col[3], lwd=3)
toplegend(legend=c('Baseline', 'PH', 'AFT'), lwd=3, col=col)

# Slide 10: PH vs AFT (log-hazard)
t <- seq(0.1, 1.5, len=99)
y <- log(t)
eta <- log(2/3)
h0 <- function(y) h(exp(y))*exp(y)
plot(y, log(h0(y)), type='l', bty='n', ylim=c(-2.5,0.5), lwd=3, col=col[1],
     ylab=expression(log(lambda(y))), las=1)
lines(y, log(h0(y)*exp(eta)), col=col[2], lwd=3)
lines(y, log(h0(y+eta)), col=col[3], lwd=3)
toplegend(legend=c('Baseline', 'PH', 'AFT'), lwd=3, col=col)

# Slides 15-16: Do-it-yourself Weibull (with step halving)
Data <- subset(pbc, !is.na(trt))
d <- Data$status != 0
y <- log(Data$time)
X <- with(Data, cbind(1, trt, stage, hepato, bili))
b <- c(mean(y), rep(0, ncol(X)-1))
s <- 1
for (i in 1:20) {
  w <- (y - as.numeric(X%*%b))/s
  mu <- exp(w)
  W <- diag(exp(w))
  b <- solve(t(X) %*% W %*% X) %*% t(X) %*% (mu-d)*s + b
  ss <- as.numeric(-crossprod(y - as.numeric(X%*%b), d-mu)/sum(d))
  s <- .5*s + .5*ss
}
b
s

# Same as
fit <- survreg(Surv(time, status!=0) ~ trt + stage + hepato + bili, pbc)
coef(fit)
fit$scale

# Exponential (compare to results from 10-03.R)
fit <- survreg(Surv(time, status!=0) ~ trt + stage + hepato + bili, dist='exponential', pbc)
coef(fit)
vcov(fit)
fit$scale

# Slide 24: Comparison
fit.e <- survreg(Surv(time, status!=0) ~ trt + stage + hepato + bili, dist='exponential', pbc)
fit.w <- survreg(Surv(time, status!=0) ~ trt + stage + hepato + bili, dist='weibull', pbc)
cbind(coef(fit.e), coef(fit.w), -coef(fit.e), -coef(fit.w)/fit.w$scale)