---
title: "Using coxdev"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Using coxdev}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4)
```

`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`.

```{r}
library(coxdev)
library(survival)
```

## A survival response

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

```{r}
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
```

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.

```{r}
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
b
```

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

```{r}
fit <- coxph(y ~ X, ties = "efron")
rbind(coxdev = b, coxph = coef(fit))
rbind(coxdev = sqrt(diag(solve(info))), coxph = sqrt(diag(vcov(fit))))
```

## Robust standard errors and the baseline hazard

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

```{r}
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))))
```

`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)))`:

```{r}
H <- cox$baseline_hazard(eta_hat)
head(H)
```

```{r, include = FALSE}
# survival::basehaz() fails on some R-devel builds (October 2026, in
# stats::delete.response(), as in survival's own CRAN check): the comparison
# below is then left out. TODO: remove once survival's R-devel checks pass.
basehaz_works <- !inherits(try(basehaz(fit, centered = FALSE), silent = TRUE), "try-error")
```

```{r, eval = basehaz_works}
bh <- basehaz(fit, centered = FALSE)
all.equal(H$cumulative_hazard, bh$hazard[match(H$time, bh$time)])
```

```{r, fig.alt = "Estimated survival curves for x1 = -1, 0 and 1, with x2 = x3 = 0."}
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")
```

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

```{r timing}
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))
```

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
`r format(8 * N^2 / 2^30, digits = 3)` GB.