source("http://myweb.uiowa.edu/pbreheny/7210/f17/notes/fun.R") # 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 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 Data <- read.delim("http://myweb.uiowa.edu/pbreheny/data/Pike1966.txt") compTable(Surv(Time, Death) ~ Group, Data) dplot(survEst(Surv(Time, Death) ~ 1, Data)) toplegend(legend=c("Exponential", "Weibull"), lwd=3, col=pal(2)) # Pike w/ guarantee time compTable(Surv(Time-100, Death) ~ Group, Data) dplot(survEst(Surv(Time-100, Death) ~ 1, Data)) 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) dplot(survEst(Surv(Time, Status) ~ 1, Data)) toplegend(legend=c("Exponential", "Weibull"), lwd=3, col=pal(2)) dplot(survEst(Surv(pmin(Time, 60), Status) ~ 1, Data)) toplegend(legend=c("Exponential", "Weibull"), lwd=3, col=pal(2))