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

# Slide 6: Lowess version
fit <- coxph(Surv(Time, Status) ~ Group + Age, GVHD)
s <- resid(fit, type='scaledsch')[,1]
x <- as.numeric(names(s))
plot(x, s, pch=19, bty='n', las=1, col='gray50', xlab='Time (Days)', ylab='Scaled Schoenfield residuals')
lines(lowess(x, s), col='red', lwd=3)

# Slide 6: cox.zph version
zfit <- cox.zph(fit, transform='identity')
plot(zfit)

# Test (Slide 9)
zfit
cor.test(x,s) # Same as above, except a t-test instead of a z-test

# Slide 7
fit <- coxph(Surv(time/365.25, status) ~ trt + karno + diagtime + age + prior + celltype, data=veteran)
s <- resid(fit, type='scaledsch')[,2]
x <- as.numeric(names(s))
plot(x, s, pch=19, bty='n', las=1, col='gray50', xlab='Time (Years)', ylab='Scaled Schoenfield residuals', xlim=c(0, 1.6))
lines(lowess(x, s), col=pal(2)[1], lwd=3)
abline(h=coef(fit)['karno'], lwd=3, col=pal(2)[2])

# Slide 8
fit <- coxph(Surv(Time, Status) ~ Group + Age, GVHD)
s <- resid(fit, type='scaledsch')[,2]
x <- as.numeric(names(s))
plot(x, s, pch=19, bty='n', las=1, col='gray50', xlab='Time (Days)', ylab='Scaled Schoenfield residuals')
lines(lowess(x, s), col=pal(2)[1], lwd=3)
abline(h=coef(fit), lwd=3, col=pal(2)[2])

# Slide 9: Simulated data
set.seed(3)
n <- 64
lam <- exp(-as.numeric(GVHD$Group))
S <- Surv(rweibull(n, 1.5, 1/lam), sample(GVHD$Status))
fit <- coxph(S ~ Group + Age, GVHD)
s <- resid(fit, type='scaledsch')[,2]
x <- as.numeric(names(s))
cor.test(x,s)
plot(x, s, pch=19, bty='n', las=1, col='gray50', xlab='Time (Days)', ylab='Scaled Schoenfield residuals')
lines(lowess(x, s), col=pal(2)[1], lwd=3)
abline(h=coef(fit), lwd=3, col=pal(2)[2])

# Slide 10-11
fit <- coxph(Surv(Time, Status) ~ Group + Age, GVHD)
cox.zph(fit)
fit <- coxph(Surv(time/365.25, status) ~ trt + karno + diagtime + age + prior + celltype, data=veteran)
cox.zph(fit)

# Effect plot function (for slides 14-21)
effectPlot <- function(fit, t, fun, term, ...) {
  b <- coef(fit)
  y <- lwr <- upr <- numeric(length(t))
  for (i in 1:length(t)) {
    delta <- rep(0, length(b))
    names(delta) <- names(b)
    delta[term] <- 1
    delta[paste0('tt(', term, ')')] <- fun(t[i])
    y[i] <- delta %*% coef(fit)
    v <- delta %*% vcov(fit) %*% delta
    lwr[i] <- y[i] + qnorm(.025)*sqrt(v)
    upr[i] <- y[i] + qnorm(.975)*sqrt(v)
  }
  plot(t, y, ylim=range(c(lwr, upr)), type='l', las=1, bty='n', ...)
  polygon(c(t, rev(t)), c(lwr, rev(upr)), col='gray85', border=NA)
  lines(t, y, lwd=3)
  abline(h=0, col='gray', lwd=2, lty=2)
}

# GVHD Treatment effect (slide 14)
GVHD$Trt <- 1*(GVHD$Group=='MTX+CSP')    # tt doesn't work with factors
fit <- coxph(Surv(Time, Status) ~ Trt + tt(Trt), GVHD, tt=function(x, t, ...) {x*t})
effectPlot(fit, 8:50, as.numeric, 'Trt', xlab='Time (Days)', ylab='Treatment effect')

# Slide 17
summary(fit)

# GVHD Treatment effect, alternate parameterization
f <- function(t) atan((t-20)/10)
fit <- coxph(Surv(Time, Status) ~ Trt + tt(Trt), GVHD, tt=function(x, t, ...) {x*f(t)})
summary(fit)
effectPlot(fit, 8:50, f, 'Trt', xlab='Time (Days)', ylab='Treatment effect')

# Karnofsky effect, linear (slide 18)
fit <- coxph(Surv(time/365.25, status) ~ trt + karno + diagtime + age + prior + celltype + tt(karno), data=veteran,
             tt=function(x, t, ...) {x*t})
effectPlot(fit, seq(0, 1.5, len=99), as.numeric, 'karno', xlab='Time (Days)', ylab='Treatment effect')

# Exponential decay rate (slide 20)
logLik(fit)
lam <- seq(0.1, 2, len=99)
l <- numeric(length(lam))
for (i in 1:length(lam)) {
  fit <- coxph(Surv(time/365.25, status) ~ trt + karno + diagtime + age + prior + celltype + tt(karno), data=veteran,
               tt=function(x, t, ...) {x*exp(-t/lam[i])})
  l[i] <- logLik(fit)
}
par(mar=c(5,6,1,1))
plot(lam, l, type='l', bty='n', las=1, lwd=2,
     xlab=expression(lambda), ylab='')
mtext('Log likelihood', 2, line=4)

# Karnofsky effect: Exponential decay
lamhat <- lam[which.max(l)]
fit <- coxph(Surv(time/365.25, status) ~ trt + karno + diagtime + age + prior + celltype + tt(karno), data=veteran,
             tt=function(x, t, ...) {x*exp(-t/lamhat)})
effectPlot(fit, seq(0, 1.5, len=99), function(t) exp(-t/lamhat), 'karno', xlab='Time (Days)', ylab='Treatment effect')