library(survival) source("http://myweb.uiowa.edu/pbreheny/7210/f18/notes/fun.R") Data <- subset(pbc, !is.na(trt)) # Newton-Raphson animation: Setup score <- function(beta) { out <- numeric(length(beta)) for (j in 1:length(beta)) { out[j] <- crossprod(x, d - t*exp(x*beta[j])) } out } tangent_line <- function(beta, next_point=FALSE) { points(beta, score(beta), col='slateblue', pch=19) w <- t*exp(x*beta) intercept <- score(beta) + crossprod(x, w*x) * beta slope <- -crossprod(x, w*x) abline(intercept, slope) if (next_point) { points(-intercept/slope, 0, pch=19) abline(h=0, lwd=2, lty=2, col='gray') cat('Next: ', -intercept/slope, '\n') } } set.seed(7) n <- 100 d <- rbinom(n, 1, 0.5) t <- rexp(n) x <- runif(n, -1, 1) # Newton-Raphson animation bb <- seq(-3, 3, 0.01) plot(bb, score(bb), type='l', lwd=3, col='slateblue', las=1, bty='n', xlab=expression(beta), ylab=expression(u(beta))) b <- -2 tangent_line(-2) tangent_line(-2, next_point=TRUE) plot(bb, score(bb), type='l', lwd=3, col='slateblue', las=1, bty='n', xlab=expression(beta), ylab=expression(u(beta))) tangent_line(-0.758, next_point=TRUE) plot(bb, score(bb), type='l', lwd=3, col='slateblue', las=1, bty='n', xlab=expression(beta), ylab=expression(u(beta))) tangent_line(0.185, next_point=TRUE) # Do-it-yourself Newton-Raphson d <- Data$status != 0 t <- Data$time X <- with(Data, cbind(1, trt, stage, hepato, bili)) b <- rep(0, ncol(X)) for (i in 1:20) { eta <- as.numeric(X%*%b) mu <- t*exp(eta) W <- diag(t*exp(eta)) b <- solve(t(X) %*% W %*% X) %*% t(X) %*% (d-mu) + b } b # Wald CI/test I <- t(X) %*% W %*% X B <- matrix(NA, 4, 4) rownames(B) <- colnames(X)[2:5] colnames(B) <- c("Estimate", "Lower", "Upper", "p") for (j in 2:5) { B[j-1, 1] <- b[j] B[j-1, 2:3] <- b[j] + qnorm(c(.025,.975))*sqrt(solve(I)[j,j]) B[j-1, 4] <- 2*pnorm(-abs(b[j]/sqrt(solve(I)[j,j]))) } B # Slide 18 :Hazard ratio, 5-unit change for bili BB <- B BB[,1:3] <- exp(BB[,1:3] * c(1, 1, 1, 5)) BB # Predictions predline <- function(x, col) { lam <- exp(crossprod(x,b)) tt <- seq(0,4000,len=199) lines(tt, pexp(tt, lam, lower=FALSE), lwd=3, col=col) } fit <- survfit(Surv(time, status!=0) ~ 1, pbc) col <- pal(4) plot(fit, conf.int=FALSE, lwd=3, col='gray', mark.time=FALSE, xmax=4000, xaxt='n', bty='n', las=1, xlab='Time (years)', ylab="Progression-free survival") lab <- seq(0,10,2) axis(1, at=lab*365.25, lab=lab) predline(c(1,0,3,1,1), col[1]) predline(c(1,0,3,1,6), col[2]) predline(c(1,0,1,0,1), col[3]) predline(c(1,0,4,1,8), col[4]) # Information matrix SE.hepato <- sqrt(solve(I)[3,3]) SE.hepato.naive <- sqrt(1/I[3,3]) # Slide 22: Diagnostic plot #1 plot(fit, conf.int=FALSE, lwd=3, col='gray', mark.time=FALSE, xmax=4000, xaxt='n', bty='n', las=1, xlab='Time (years)', ylab="Progression-free survival") lab <- seq(0,10,2) axis(1, at=lab*365.25, lab=lab) tt <- seq(0,4000,len=199) lam <- sum(Data$status!=0)/sum(Data$time) lines(tt, pexp(tt, lam, lower=FALSE), lwd=3, col=pal(2)[2]) # Slide 23: Diagnostic plot #2 y <- -log(summary(fit, time=tt)$surv) plot(tt, y, col='gray', type='l', lwd=3, xaxt='n', bty='n', las=1, xlab='Time (years)', ylab="-log S(t)") axis(1, at=lab*365.25, lab=lab) lines(range(tt), range(tt)*lam, col=pal(2)[2], lwd=3) # Slide 26: Hypothetical diagnostic plot t <- c(rexp(250, 0.5), rexp(250, 4)) fit <- survfit(Surv(t, rep(1, length(t))) ~ 1) tt <- seq(0, max(t), len=99) y <- -log(summary(fit, time=tt)$surv) plot(tt, y, col='gray', type='l', lwd=3, bty='n', las=1, xlab='Time', ylab="-log S(t)") lam <- 1/mean(t) lines(range(tt), range(tt)*lam, col=pal(2)[2], lwd=3)