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

# Basic methods
fit <- coxph(Surv(time, status!=0) ~ trt + stage + hepato + bili, pbc)
coef(fit)
vcov(fit)
head(X <- model.matrix(fit))
dim(pbc)
dim(X)

# Summary
summary(fit)

# Anova
fit0 <- coxph(Surv(time, status!=0) ~ trt + stage + hepato + bili, pbc)
fit1 <- coxph(Surv(time, status!=0) ~ trt + factor(stage) + hepato + bili, pbc)
anova(fit0, fit1)

# logLik
logLik(fit0)
logLik(fit1)

# AIC
AIC(fit0)
AIC(fit1)

# BIC
BIC(fit0)
BIC(fit1)
p <- length(coef(fit1))
-2*logLik(fit1) + log(fit1$nevent)*p   # BIC is this
-2*logLik(fit1) + log(fit1$n)*p        # Not this
w <- exp(-c(BIC(fit0), BIC(fit1))/2)
w/sum(w)                               # Posterior model probabilities
b <- c(BIC(fit0), BIC(fit1))
b <- b - min(b)
w <- exp(-b/2)
w/sum(w)                               # Posterior model probabilities

# Predict
newData <- data.frame(trt=0, stage=2, hepato=1, bili=1)
predict(fit, newData)
XX <- as.matrix(newData)
XX %*% coef(fit)
m <- apply(model.matrix(fit), 2, mean)
(XX-m) %*% coef(fit)

# Predictions are invariant
Data <- pbc
Data$bili <- 5 + Data$bili/3
newFit <- coxph(Surv(time, status!=0) ~ trt + stage + hepato + bili, Data)
newData$bili <- 5 + newData$bili/3
predict(newFit, newData)

# Functions to make comparison table and diagnostic plot
compTable <- function(form, Data) {
  B <- matrix(NA, 3, 3, dimnames=list(c('Cox', 'Weibull', 'Exponential'), c('Est', 'Lower', 'Upper')))
  fit <- coxph(form, Data)
  B[1,] <- c(coef(fit)[1], confint(fit, 1))
  fit <- survreg(form, Data)
  B[2,] <- -c(coef(fit)[2], rev(confint(fit, 2)))/fit$scale
  fit <- survreg(form, Data, dist='exponential')
  B[3,] <- -c(coef(fit)[2], rev(confint(fit, 2)))
  B
}
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)
}

# PBC
compTable(Surv(time, status!=0) ~ stage + trt + hepato + bili, pbc)
dplot(survEst(Surv(time, status!=0) ~ 1, pbc))
toplegend(legend=c("Exponential", "Weibull"), lwd=3, col=pal(2))

# Pike
compTable(Surv(Time, Death) ~ Group, Pike)
dplot(survEst(Surv(Time, Death) ~ 1, Pike))
toplegend(legend=c("Exponential", "Weibull"), lwd=3, col=pal(2))

# Pike w/ guarantee time
compTable(Surv(Time-100, Death) ~ Group, Pike)
dplot(survEst(Surv(Time-100, Death) ~ 1, Pike))
toplegend(legend=c("Exponential", "Weibull"), lwd=3, col=pal(2))

logLik(survreg(Surv(Time, Death) ~ Group, Pike))
logLik(survreg(Surv(Time-100, Death) ~ Group, Pike))

# GVHD
compTable(Surv(Time, Status) ~ Group, GVHD)
compTable(Surv(pmin(Time, 60), Status) ~ Group, GVHD)
dplot(survEst(Surv(Time, Status) ~ 1, GVHD))
toplegend(legend=c("Exponential", "Weibull"), lwd=3, col=pal(2))
dplot(survEst(Surv(pmin(Time, 60), Status) ~ 1, GVHD))
toplegend(legend=c("Exponential", "Weibull"), lwd=3, col=pal(2))