# Setup; slide 2
require(survival)
source("http://myweb.uiowa.edu/pbreheny/7210/f15/notes/fun.R")
PBC <- pbc[!is.na(pbc$trt),]
f <- function(x) {pmin(x, 3.5)}
fit1 <- coxph(Surv(time/365.25, status!=0) ~ trt + albumin, PBC)
fit2 <- coxph(Surv(time/365.25, status!=0) ~ trt + stage + hepato + f(albumin) + log(bili), PBC)
set.seed(2)
X <- matrix(rnorm(30*nrow(PBC)), nrow(PBC), 30)
fit3 <- coxph(Surv(time/365.25, status!=0) ~ trt + stage + hepato + f(albumin) + log(bili) + X, PBC)

# Slide 3
sd(fit1$linear.predictors)
sd(fit2$linear.predictors)
sd(fit3$linear.predictors)

# Slide 4
brk <- seq(-7,7,0.5)
at <- log(4^(-(-4:4)))
lab <- paste(c('1/', '')[1*(at>=0)+1], exp(abs(at)))
hist(fit1$linear.predictors, xlim=c(-7,7), breaks=brk, xaxt='n', xlab='Hazard ratio')
axis(1, at=at, labels=lab)
mtext('Model 1')
hist(fit2$linear.predictors, xlim=c(-7,7), breaks=brk, xaxt='n', xlab='Hazard ratio')
axis(1, at=at, labels=lab)
mtext('Model 2')
hist(fit3$linear.predictors, xlim=c(-7,7), breaks=brk, xaxt='n', xlab='Hazard ratio')
axis(1, at=at, labels=lab)
mtext('Model 3')

# Slide 5
survPlot <- function(fit, ...) {
  sfit <- survfit(fit)
  plot(sfit, mark.time=FALSE, conf.int=FALSE, las=1, bty='n', ...)
  for (m in c(-2,-1,1,2)) {
    tmp <- sfit
    tmp$surv <- tmp$surv^exp(m*sd(fit$linear.predictors))
    lines(tmp, mark.time=FALSE, conf.int=FALSE)
  }
}
survPlot(fit1, xlab='Time (years)', ylab='Survival')
mtext('Model 1')
survPlot(fit2, xlab='Time (years)', ylab='Survival')
mtext('Model 2')
survPlot(fit3, xlab='Time (years)', ylab='Survival')
mtext('Model 3')

# R2
summary(fit1)$rsq
summary(fit2)$rsq
summary(fit3)$rsq
LR <- 2*diff(fit1$loglik)
1 - exp(-LR/fit1$n)

# Concordance
survConcordance(Surv(time, status!=0) ~ fit2$linear.predictors, PBC)
summary(fit1)$concordance
summary(fit2)$concordance
summary(fit3)$concordance

# Shrinkage
g <- function(fit) {
  LR <- 2*diff(fit$loglik)
  p <- length(coef(fit))
  (LR-p)/LR
}
g(fit1)
g(fit2)
g(fit3)

# SD of shrunken PI
s1 <- fit1$linear.predictors*g(fit1)
s2 <- fit2$linear.predictors*g(fit2)
s3 <- fit3$linear.predictors*g(fit3)
sd(s1)
sd(s2)
sd(s3)

# Sim
set.seed(1)
beta <- c(rep(log(2), 2), rep(0,28))
eta <- X %*% beta
y <- rexp(length(eta), exp(eta))
fit.sim <- coxph(Surv(y)~X)
g(fit.sim)
par(mar=c(5,5,0.5,0.5))
plot(eta, fit.sim$linear.predictors, pch=19, col='gray60', las=1, bty='n',
     xlab=expression(eta), ylab=expression(hat(eta)))
abline(0, 1, col=pal(2)[1], lwd=3)
abline(0, coef(lm(fit.sim$linear.predictors~0+eta)), col=pal(2)[2], lwd=3)

s <- fit.sim$linear.predictors*g(fit.sim)
plot(eta, s, pch=19, col='gray60', las=1, bty='n',
     xlab=expression(eta), ylab=expression(hat(eta)))
abline(0, 1, col=pal(2)[1], lwd=3)
abline(0, coef(lm(s~0+eta)), col=pal(2)[2], lwd=3)

## Calibrated Rsq ?
summary(fit1)$rsq*g(fit1)
summary(fit2)$rsq*g(fit2)
summary(fit3)$rsq*g(fit3)
0.5 + (summary(fit1)$concordance-0.5)*g(fit1)
0.5 + (summary(fit2)$concordance-0.5)*g(fit2)
0.5 + (summary(fit3)$concordance-0.5)*g(fit3)