# 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)