## Basic methods
require(survival)
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

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

# Pike
Data <- read.delim("http://myweb.uiowa.edu/pbreheny/data/Pike1966.txt")
compTable(Surv(Time, Death) ~ Group, Data)

compTable(Surv(Time-100, Death) ~ Group, Data)

x <- survEst(Surv(Time, Death) ~ 1, Data)
dplot(x)
mtext('Original')
x <- survEst(Surv(Time-100, Death) ~ 1, Data)
dplot(x)
mtext('Guarantee time = 100 days')
toplegend(legend=c("Exponential", "Weibull"), lwd=3, col=pal(2))

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

# GVHD
Data <- read.delim("http://myweb.uiowa.edu/pbreheny/data/gvhd.txt")
compTable(Surv(Time, Status) ~ Group, Data)
compTable(Surv(pmin(Time, 60), Status) ~ Group, Data)