Using coxdev

coxdev computes the deviance of the Cox proportional hazards model, its gradient and the diagonal of its Hessian with respect to the linear predictor, the information as a matrix-free operator, the score residuals, the baseline hazard, and draws of the score. This vignette fits a model by Newton’s method, computes robust standard errors and survival curves, and shows what keeping the sorted data in a CoxDeviance object saves. Sampling the score has its own vignette, vignette("sampling", package = "coxdev"). The Python package has the same two, as the notebooks doc/using_coxdev.md and doc/sampling.md.

library(coxdev)
library(survival)

A survival response

Survival times on a coarse grid, so that there are many ties, with delayed entry and censoring.

set.seed(1)
n <- 500
p <- 3
X <- matrix(rnorm(n * p), n, p, dimnames = list(NULL, c("x1", "x2", "x3")))
beta <- c(0.5, -0.5, 0.25)
time <- rexp(n, exp(drop(X %*% beta)))
cens <- rexp(n, 0.5)
stop <- ceiling(pmin(time, cens) * 4) / 4
status <- as.integer(time <= cens)
start <- floor(runif(n, 0, 0.5) * stop * 8) / 8
y <- Surv(start, stop, status)

cox <- CoxDeviance$new(y, tie_breaking = "efron")
cox
#> <CoxDeviance> 500 observations, 344 failures, efron ties, with start times

The response can also be given as vectors, as in the Python package: CoxDeviance$new(event = stop, status = status, start = start).

Deviance, gradient and information

compute() returns the deviance and its derivatives with respect to the linear predictor eta. information() returns the information (half the Hessian of the deviance) as an n by n operator. It is never formed: a product with it costs about as much as a gradient. For a design X, crossprod(X, I %*% X) is the information of the coefficients, so Newton’s method takes a few lines.

b <- rep(0, p)
for (iter in 1:10) {
  eta <- drop(X %*% b)
  res <- cox$compute(eta)
  score <- -drop(crossprod(X, res$gradient)) / 2
  info <- crossprod(X, cox$information(eta) %*% X)
  step <- solve(info, score)
  b <- b + step
  if (max(abs(step)) < 1e-10) break
}
iter
#> [1] 4
b
#>         x1         x2         x3 
#>  0.4530326 -0.3534272  0.3209856

These are the estimates of coxph from survival, with the same standard errors:

fit <- coxph(y ~ X, ties = "efron")
rbind(coxdev = b, coxph = coef(fit))
#>               x1         x2        x3
#> coxdev 0.4530326 -0.3534272 0.3209856
#> coxph  0.4530326 -0.3534272 0.3209856
rbind(coxdev = sqrt(diag(solve(info))), coxph = sqrt(diag(vcov(fit))))
#>                x1         x2         x3
#> coxdev 0.05728271 0.05369269 0.05491812
#> coxph  0.05728271 0.05369269 0.05491812

Robust standard errors and the baseline hazard

The score residuals of the design give the robust (sandwich) variance, solve(I) %*% crossprod(U) %*% solve(I):

eta_hat <- drop(X %*% b)
U <- cox$score_residuals(eta_hat, X)
info_hat <- crossprod(X, cox$information(eta_hat) %*% X)
robust <- solve(info_hat, t(solve(info_hat, crossprod(U))))
fit_robust <- coxph(y ~ X, ties = "efron", id = seq_len(n), robust = TRUE)
rbind(coxdev = sqrt(diag(robust)), coxph = sqrt(diag(vcov(fit_robust))))
#>                x1         x2         x3
#> coxdev 0.04851333 0.04433025 0.04591987
#> coxph  0.04851333 0.04433025 0.04591987

baseline_hazard() gives Efron’s estimate of the cumulative baseline hazard (Breslow’s with tie_breaking = "breslow"), the basehaz(fit, centered = FALSE) of coxph. The survival curve of a subject with covariates x0 is exp(-cumulative_hazard * exp(sum(x0 * b))):

H <- cox$baseline_hazard(eta_hat)
head(H)
#>   time    hazard cumulative_hazard
#> 1 0.25 0.3524866         0.3524866
#> 2 0.50 0.2916936         0.6441801
#> 3 0.75 0.2682672         0.9124474
#> 4 1.00 0.3602220         1.2726694
#> 5 1.25 0.1454810         1.4181504
#> 6 1.50 0.2759841         1.6941344

bh <- basehaz(fit, centered = FALSE)
all.equal(H$cumulative_hazard, bh$hazard[match(H$time, bh$time)])
#> [1] TRUE

plot(NULL, xlim = c(0, max(H$time)), ylim = c(0, 1), xlab = "time", ylab = "survival")
for (x1 in c(-1, 0, 1)) {
  surv <- exp(-H$cumulative_hazard * exp(b[1] * x1))
  lines(c(0, H$time), c(1, surv), type = "s", lty = x1 + 2)
}
legend("topright", paste("x1 =", c(-1, 0, 1)), lty = 1:3, bty = "n")

Estimated survival curves for x1 = -1, 0 and 1, with x2 = x3 = 0.

What keeping state saves

Creating a CoxDeviance sorts the data, once. Each compute() after that is linear in n, and repeating a call with the same linear predictor and weights returns the cached result. A stateless function has to sort on every call. That cost is paid at every step of an optimizer, and at every product with the information.

We time 50 evaluations of the deviance and gradient at 100,000 observations, each at a different linear predictor.

set.seed(2)
N <- 1e5
XN <- matrix(rnorm(N * p), N, p)
etaN <- drop(XN %*% beta)
stopN <- ceiling(rexp(N, exp(etaN)) * 20) / 20
statusN <- rbinom(N, 1, 0.7)
startN <- floor(runif(N, 0, 0.5) * stopN * 40) / 40
yN <- Surv(startN, stopN, statusN)
etas <- lapply(1:50, function(i) etaN + rnorm(N, sd = 0.1))

per_call <- function(expr) {
  t <- system.time(expr)[["elapsed"]]
  sprintf("%.1f ms", 1000 * t / 50)
}
setup <- system.time(coxN <- CoxDeviance$new(yN))[["elapsed"]]
timings <- c(
  "one object, reused" = per_call(for (e in etas) coxN$compute(e)),
  "a new object for each call" = per_call(for (e in etas) CoxDeviance$new(yN)$compute(e)),
  "one object, same eta again (cached)" = per_call(for (e in etas) coxN$compute(etas[[1]])),
  "one information product" = per_call({
    I <- coxN$information(etas[[1]])
    for (e in etas) I %*% XN[, 1]
  })
)
knitr::kable(data.frame(time = timings), caption = sprintf(
  "Time per call, n = %d. Creating the object (the sort) took %.0f ms.",
  N, 1000 * setup))

Time per call, n = 100000. Creating the object (the sort) took 52 ms.
	time
one object, reused	14.9 ms
a new object for each call	58.7 ms
one object, same eta again (cached)	1.5 ms
one information product	3.5 ms

The times depend on the machine, but the pattern is stable. Reusing the object avoids the sort, which costs several times as much as an evaluation. The information operator is what makes the Hessian usable at this size: as a dense matrix it would take 74.5 GB.