library(survival)
source("http://myweb.uiowa.edu/pbreheny/7210/f18/notes/fun.R")
Data <- subset(pbc, !is.na(trt))

# Newton's method: Setup
score <- function(beta) {
  out <- numeric(length(beta))
  for (j in 1:length(beta)) {
    out[j] <- crossprod(x, d - t*exp(x*beta[j]))
  }
  out
}
tangent_line <- function(beta, next_point=FALSE) {
  points(beta, score(beta), col='slateblue', pch=19)
  w <- t*exp(x*beta)
  intercept <- score(beta) + crossprod(x, w*x) * beta
  slope <- -crossprod(x, w*x)
  abline(intercept, slope)
  if (next_point) {
    points(-intercept/slope, 0, pch=19)
    abline(h=0, lwd=2, lty=2, col='gray')
    cat('Next: ', -intercept/slope, '\n')
  }
}
set.seed(7)
n <- 100
d <- rbinom(n, 1, 0.5)
t <- rexp(n)
x <- runif(n, -1, 1)

# Newton's method animation
bb <- seq(-3, 3, 0.01)
plot(bb, score(bb), type='l', lwd=3, col='slateblue', las=1, bty='n', xlab=expression(beta), ylab=expression(u(beta)))
b <- -2
tangent_line(-2)
tangent_line(-2, next_point=TRUE)
plot(bb, score(bb), type='l', lwd=3, col='slateblue', las=1, bty='n', xlab=expression(beta), ylab=expression(u(beta)))
tangent_line(-0.758, next_point=TRUE)
plot(bb, score(bb), type='l', lwd=3, col='slateblue', las=1, bty='n', xlab=expression(beta), ylab=expression(u(beta)))
tangent_line(0.185, next_point=TRUE)

# Do-it-yourself Newton's method
d <- Data$status != 0
t <- Data$time
X <- with(Data, cbind(1, trt, stage, hepato, bili))
b <- rep(0, ncol(X))
for (i in 1:20) {
  eta <- as.numeric(X%*%b)
  mu <- t*exp(eta)
  W <- diag(t*exp(eta))
  b <- solve(t(X) %*% W %*% X) %*% t(X) %*% (d-mu) + b
}
b

# Wald CI/test
I <- t(X) %*% W %*% X
B <- matrix(NA, 4, 4)
rownames(B) <- colnames(X)[2:5]
colnames(B) <- c("Estimate", "Lower", "Upper", "p")
for (j in 2:5) {
  B[j-1, 1] <- b[j]
  B[j-1, 2:3] <- b[j] + qnorm(c(.025,.975))*sqrt(solve(I)[j,j])
  B[j-1, 4] <- 2*pnorm(-abs(b[j]/sqrt(solve(I)[j,j])))
}
B

# Slide 18 :Hazard ratio, 5-unit change for bili
BB <- B
BB[,1:3] <- exp(BB[,1:3] * c(1, 1, 1, 5))
BB

# Predictions
predline <- function(x, col) {
  lam <- exp(crossprod(x,b))
  tt <- seq(0,4000,len=199)
  lines(tt, pexp(tt, lam, lower=FALSE), lwd=3, col=col)
}
fit <- survfit(Surv(time, status!=0) ~ 1, pbc)
col <- pal(4)
plot(fit, conf.int=FALSE, lwd=3, col='gray', mark.time=FALSE, xmax=4000,
     xaxt='n', bty='n', las=1, xlab='Time (years)', ylab="Progression-free survival")
lab <- seq(0,10,2)
axis(1, at=lab*365.25, lab=lab)
predline(c(1,0,3,1,1), col[1])
predline(c(1,0,3,1,6), col[2])
predline(c(1,0,1,0,1), col[3])
predline(c(1,0,4,1,8), col[4])

# Information matrix
SE.hepato <- sqrt(solve(I)[3,3])
SE.hepato.naive <- sqrt(1/I[3,3])

# Slide 22: Diagnostic plot #1
plot(fit, conf.int=FALSE, lwd=3, col='gray', mark.time=FALSE, xmax=4000,
     xaxt='n', bty='n', las=1, xlab='Time (years)', ylab="Progression-free survival")
lab <- seq(0,10,2)
axis(1, at=lab*365.25, lab=lab)
tt <- seq(0,4000,len=199)
lam <- sum(Data$status!=0)/sum(Data$time)
lines(tt, pexp(tt, lam, lower=FALSE), lwd=3, col=pal(2)[2])

# Slide 23: Diagnostic plot #2
y <- -log(summary(fit, time=tt)$surv)
plot(tt, y, col='gray', type='l', lwd=3, xaxt='n', bty='n', las=1,
     xlab='Time (years)', ylab="-log S(t)")
axis(1, at=lab*365.25, lab=lab)
lines(range(tt), range(tt)*lam, col=pal(2)[2], lwd=3)

# Slide 26: Hypothetical diagnostic plot
t <- c(rexp(250, 0.5), rexp(250, 4))
fit <- survfit(Surv(t, rep(1, length(t))) ~ 1)
tt <- seq(0, max(t), len=99)
y <- -log(summary(fit, time=tt)$surv)
plot(tt, y, col='gray', type='l', lwd=3, bty='n', las=1,
     xlab='Time', ylab="-log S(t)")
lam <- 1/mean(t)
lines(range(tt), range(tt)*lam, col=pal(2)[2], lwd=3)