require(survival)
Data <- subset(pbc, !is.na(trt))

# Do-it-yourself Newton-Raphson
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)