Sampling the score

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

We sample at the estimate of the coefficients, from coxph:

b <- coef(coxph(y ~ X, ties = "efron"))
b
#>        Xx1        Xx2        Xx3 
#>  0.4530326 -0.3534272  0.3209856

The methods

Each draw r has one entry per observation, and crossprod(X, r) is a draw of the score of the coefficients. We draw at the estimate b, by four methods:

"multinomial": the failures at each event time drawn with replacement, with probabilities proportional to exp(eta). This is the sampling model of Breslow’s likelihood, so the covariance is the Breslow information.
"fisher": the set of failures at each event time drawn from Fisher’s noncentral hypergeometric distribution. This is the sampling model of the exact partial likelihood, whose information is smaller than Breslow’s when ties are heavy.
"bernoulli": each subject at risk fails independently, with probability min(1, d p), where d is the number of failures at the event time and p the subject’s share of the risk. Its mean is zero only while d p < 1: under heavy ties the cap binds, and the draws are biased.
"wild": the score residuals times independent normals. Its covariance is the sandwich, which does not rely on the model being right.

eta_hat <- drop(X %*% b)
methods <- c("multinomial", "fisher", "bernoulli", "wild")
draws <- lapply(methods, function(m) {
  crossprod(X, cox$sample_score(eta_hat, method = m, n_draws = 4000, seed = 1))
})
names(draws) <- methods

Here are the means of the draws, in standard deviations, and their variances relative to the Breslow information:

breslow <- CoxDeviance$new(y, tie_breaking = "breslow")
info_breslow <- crossprod(X, breslow$information(eta_hat) %*% X)
summary_table <- t(sapply(draws, function(d) {
  c(rowMeans(d) / sqrt(diag(info_breslow)), diag(cov(t(d))) / diag(info_breslow))
}))
colnames(summary_table) <- c(paste("mean", colnames(X)), paste("variance", colnames(X)))
knitr::kable(summary_table, digits = 2, caption = paste(
  "Means of the draws of the score, in units of sqrt(Breslow information),",
  "and variances relative to the Breslow information"))

Means of the draws of the score, in units of sqrt(Breslow information), and variances relative to the Breslow information
	mean x1	mean x2	mean x3	variance x1	variance x2	variance x3
multinomial	0.00	-0.01	-0.01	1.02	0.99	0.99
fisher	0.00	0.01	0.00	0.64	0.66	0.68
bernoulli	-0.45	0.25	-0.38	0.56	0.60	0.61
wild	-0.01	-0.02	0.00	0.70	0.70	0.66

The multinomial draws have the Breslow information as covariance, and the Fisher draws a smaller one, the information of the exact likelihood. The Bernoulli draws are biased: at these ties, d p often exceeds 1, so the cap min(1, d p) gives too few failures to the subjects most at risk.

The wild draws estimate the sandwich, crossprod(U) for the score residuals U of the design, the middle of the robust variance solve(I) %*% crossprod(U) %*% solve(I) of coxph(robust = TRUE):

U <- cox$score_residuals(eta_hat, X)
sandwich <- crossprod(U)
fit_robust <- coxph(y ~ X, ties = "efron", id = seq_len(n), robust = TRUE)
I_hat <- solve(fit_robust$naive.var)
rbind(wild = diag(cov(t(draws$wild))), sandwich = diag(sandwich),
      coxph = diag(I_hat %*% fit_robust$var %*% I_hat))
#>                x1       x2       x3
#> wild     224.1504 250.3854 227.6800
#> sandwich 225.5982 238.1139 236.6022
#> coxph    225.5982 238.1139 236.6022

The distributions of the first coefficient’s score:

cols <- c(multinomial = "black", fisher = "firebrick", bernoulli = "darkorange", wild = "steelblue")
dens <- lapply(draws, function(d) density(d[1, ]))
plot(NULL, xlim = c(-1, 1) * 4 * sqrt(info_breslow[1, 1]),
     ylim = c(0, 1.35 * max(sapply(dens, function(d) max(d$y)))),
     xlab = "score of x1", ylab = "density")
for (m in methods) lines(dens[[m]], col = cols[m], lwd = 2)
curve(dnorm(x, sd = sqrt(info_breslow[1, 1])), add = TRUE, lty = 2)
legend("topleft", c(methods, "N(0, Breslow info.)"), cex = 0.8,
       col = c(cols, "black"), lty = c(rep(1, 4), 2), lwd = c(rep(2, 4), 1), bty = "n")

Densities of 4000 draws of the score of x1, by method, with the normal density with the Breslow information as variance.

Choosing the multipliers

wild_score() gives the wild draw for multipliers you choose: the sum over observations of the multiplier times the score residual, G %*% U. With multipliers all one it is the score itself. Other multipliers give other wild bootstraps, here Rademacher signs in place of normals:

all.equal(cox$wild_score(eta_hat, rep(1, n)), -cox$compute(eta_hat)$gradient / 2)
#> [1] TRUE
rademacher <- replicate(4000, {
  drop(crossprod(X, cox$wild_score(eta_hat, sample(c(-1, 1), n, replace = TRUE))))
})
rbind(rademacher = diag(cov(t(rademacher))), sandwich = diag(sandwich))
#>                  x1       x2       x3
#> rademacher 229.9185 249.0698 234.4643
#> sandwich   225.5982 238.1139 236.6022

Reproducibility

The same seed gives the same draws on every platform, and from the Python package. With seed = NULL, the seed comes from R’s generator, so set.seed() makes the draws reproducible.

identical(cox$sample_score(eta_hat, method = "fisher", seed = 7),
          cox$sample_score(eta_hat, method = "fisher", seed = 7))
#> [1] TRUE