coxdev

A high-performance C++ for computing Cox proportional hazards model deviance, gradients, and Hessian information matrices. Built with Eigen for ease of use in R and Python, this library provides efficient survival analysis computations with support for different tie-breaking methods.

Features

• High Performance: C++ implementation with Eigen linear algebra library
• Comprehensive Support: Handles both Efron and Breslow tie-breaking methods
• Left-Truncated Data: Support for left-truncated survival data
• Efficient Computations: Optimized algorithms for deviance, gradient, and Hessian calculations
• Memory Efficient: Uses linear operators for large-scale computations
• Sampling the Score: Draws of the score from the fitted model (multinomial, Bernoulli or Fisher noncentral hypergeometric draws of the failures at each event time) or from its score residuals (a wild bootstrap), reproducible from a seed on every platform
• Cross-Platform: Works on Linux, macOS, and Windows

Installation

Prerequisites

This package requires the Eigen C++ library headers. The Eigen library is included as a git submodule.

1. Initialize Eigen submodule: The Eigen library is included as a git submodule. Make sure it's initialized:

   git submodule update --init --recursive
   

Standard Installation

pip install .

With Custom Eigen Path

If you have Eigen installed elsewhere, you can specify its location:

env EIGEN_LIBRARY_PATH=/path/to/eigen pip install .

Development Installation

pip install .[dev]
pip install -e . --no-build-isolation

Quick Start

import numpy as np
from coxdev import CoxDeviance

# Generate sample survival data
n_samples = 1000
event_times = np.random.exponential(1.0, n_samples)
status = np.random.binomial(1, 0.7, n_samples)  # 70% events, 30% censored
linear_predictor = np.random.normal(0, 1, n_samples)

# Create CoxDeviance object
coxdev = CoxDeviance(event=event_times, status=status, tie_breaking='efron')

# Compute deviance and related quantities
result = coxdev(linear_predictor)

print(f"Deviance: {result.deviance:.4f}")
print(f"Saturated log-likelihood: {result.loglik_sat:.4f}")
print(f"Gradient norm: {np.linalg.norm(result.gradient):.4f}")

R Package

The R package in R/coxdev runs the same C++ code, with the interface of the Python package. Install it from a clone of the repository:

R CMD INSTALL R/coxdev

or with remotes::install_github("jonathan-taylor/coxdev", subdir = "R/coxdev").

library(coxdev)
set.seed(1)
n <- 1000
X <- matrix(rnorm(n * 3), n, 3)
eta <- drop(X %*% c(0.5, -0.5, 0.25))
event <- rexp(n, exp(eta))
status <- rbinom(n, 1, 0.7)

cox <- CoxDeviance$new(event = event, status = status, tie_breaking = "efron")
# or with a Surv response: Surv(time, status) or Surv(start, stop, status),
# with strata from glmnet::stratifySurv() if set
cox <- CoxDeviance$new(survival::Surv(event, status), tie_breaking = "efron")

res <- cox$compute(eta)          # Python: coxdev(eta)
res$deviance                     # the objective
res$gradient                     # its gradient in eta
res$diag_hessian                 # the diagonal of its Hessian

I <- cox$information(eta)        # a matrix-free operator, never n x n
I %*% X                          # or I$multiply(X); t(X) %*% I also works
info_beta <- crossprod(X, I %*% X)

draws <- cox$sample_score(eta, method = "fisher", n_draws = 2000, seed = 1)
G <- rnorm(n)
w <- cox$wild_score(eta, G)

Python	R
CoxDeviance(event=, status=, start=, strata=, sample_weight=, tie_breaking=)	CoxDeviance$new(event =, status =, start =, strata =, sample_weight =, tie_breaking =)
	CoxDeviance$new(Surv(...), strata =, sample_weight =, tie_breaking =)
coxdev(eta, sample_weight)	cox$compute(eta, sample_weight)
coxdev.information(eta, sample_weight) (a LinearOperator)	cox$information(eta, sample_weight) (a CoxInformation, with %*%, $multiply(), $diagonal(), $to_dense(), dim())
coxdev.sample_score(eta, sample_weight, method, n_draws, seed)	cox$sample_score(eta, sample_weight, method, n_draws, seed)
coxdev.wild_score(eta, multipliers, sample_weight)	cox$wild_score(eta, multipliers, sample_weight)
coxdev.score_residuals(eta, X, sample_weight)	cox$score_residuals(eta, X, sample_weight)
coxdev.baseline_hazard(eta, sample_weight) (a BaselineHazard)	cox$baseline_hazard(eta, sample_weight) (a data frame, with stratum labels as given)

The same seed gives the same draws from R and Python. In R, a seed is a whole number up to 2^53, or a string of digits for any seed below 2^64; with seed = NULL it is drawn from R's generator, so set.seed() makes the draws reproducible. Strata may be integers, a factor or character. %*% with the operator needs R 4.4.0 or later; $multiply() works on any version.

The C++ sources in R/coxdev/src and R/coxdev/inst/include are copies of those in cpp/, so that the R package is self-contained: after editing cpp/, run R/coxdev/tools/sync_cpp.sh (CI checks that the copies agree).

Notebooks

doc/using_coxdev.md (fitting by Newton's method with the information operator, robust standard errors, the baseline hazard and survival curves, and what keeping state saves) and doc/sampling.md (the methods of sample_score compared) are jupytext notebooks, with executed copies in doc/*.ipynb. The R package has the same two as vignettes.

Advanced Usage

Left-Truncated Data

# With start times (left-truncated data)
start_times = np.random.exponential(0.5, n_samples)
coxdev = CoxDeviance(
    event=event_times, 
    status=status, 
    start=start_times,
    tie_breaking='efron'
)

Computing Information Matrix

# Get information matrix as a linear operator
info_matrix = coxdev.information(linear_predictor)

# Matrix-vector multiplication
v = np.random.normal(0, 1, n_samples)
result_vector = info_matrix @ v

# For small problems, you can compute the full matrix
X = np.random.normal(0, 1, (n_samples, 10))
beta = np.random.normal(0, 1, 10)
eta = X @ beta

# Information matrix for coefficients: X^T @ I @ X
I = info_matrix @ X
information_matrix = X.T @ I

Sampling the Score

sample_score draws the score at a linear predictor. Each draw has one entry per observation, on the scale of the score -gradient / 2, so for a design X the vector X.T @ draw is a draw of the score of the coefficients. Ties, start times and strata are handled as for the deviance.

X = np.random.normal(0, 1, (n_samples, 3))
eta = X @ np.array([0.5, -0.5, 0.25])
coxdev = CoxDeviance(event=event_times, status=status, tie_breaking='efron')

# 2000 draws of the score of the coefficients
draws = X.T @ coxdev.sample_score(eta, method='multinomial', n_draws=2000, seed=1)

# their covariance is about the (Breslow) information: the ratios are near 1
I = coxdev.information(eta)
information_matrix = X.T @ (I @ X)
print(np.diag(np.cov(draws)) / np.diag(information_matrix))

The methods are:

method	Draw	Score of the draw	Covariance
'multinomial'	the d failures at each event time drawn with replacement from the risk set, with probabilities proportional to sample_weight * exp(eta)	under Breslow's (1974) likelihood, of which this is the sampling model	the Breslow information
'bernoulli'	each subject at risk fails independently with probability min(1, d * p), a Cox-model analogue of the weird bootstrap (Andersen et al., 1993; Burr, 1994)	Breslow's, with the number of failures drawn	less than the information when ties are heavy
'fisher'	the set of d failures at each event time drawn from Fisher's noncentral hypergeometric distribution (Harkness, 1965; Fog, 2008): independent Bernoullis given d failures, P(A) proportional to the product of exp(eta_j) over the set A	under the exact discrete partial likelihood (Cox, 1972; Kalbfleisch and Prentice, 2002), of which this is the sampling model, centered at its exact mean	the information of the exact likelihood
'wild'	the score residuals of the observations times independent standard normals (Lin, Wei and Ying, 1993; Kline and Santos, 2012)	of the likelihood set by tie_breaking (Breslow, 1974; Efron, 1977)	the sandwich, the sum of the outer products of the score residuals (Lin and Wei, 1989)

Each model-based draw is the score of failures drawn from the model: at each event time, the sum over the failures of x_j minus their expected sum. For 'multinomial' and 'fisher' the draw is from the model whose likelihood gives the score, so every draw has mean exactly zero. For 'bernoulli' this holds only while d * p < 1 for every subject at risk: the cap min(1, d * p) biases it when ties are heavy. The methods differ only at tied event times: with one failure at each event time all three draw the same failure distribution. The model-based methods (all but 'wild') need equal positive weights and do not depend on tie_breaking. 'wild' does not rely on the model: its covariance is the sandwich, which is right even when the proportional hazards model is not.

'fisher' samples and centers with the elementary symmetric polynomials of the risk weights, computed on the log scale by the recursion of Gail, Lubin and Rubinstein (1981), the one used for the exact partial likelihood by R's survival package and SAS PROC PHREG. Each event time with d tied failures among m at risk costs O(m d) once, and O(m) per draw.

wild_score gives the wild draw for given multipliers: with multipliers G, X.T @ coxdev.wild_score(eta, G) is the sum over observations of G[i] times the score residual of observation i. With G all one it is the score.

References

• Andersen, P. K., Borgan, Ø., Gill, R. D. and Keiding, N. (1993). Statistical Models Based on Counting Processes. Springer, New York.
• Breslow, N. E. (1974). Covariance analysis of censored survival data. Biometrics, 30(1), 89–99.
• Burr, D. (1994). A comparison of certain bootstrap confidence intervals in the Cox model. Journal of the American Statistical Association, 89(428), 1290–1302.
• Cox, D. R. (1972). Regression models and life-tables (with discussion). Journal of the Royal Statistical Society, Series B, 34(2), 187–220.
• Efron, B. (1977). The efficiency of Cox's likelihood function for censored data. Journal of the American Statistical Association, 72(359), 557–565.
• Fog, A. (2008). Sampling methods for Wallenius' and Fisher's noncentral hypergeometric distributions. Communications in Statistics: Simulation and Computation, 37(2), 241–257.
• Gail, M. H., Lubin, J. H. and Rubinstein, L. V. (1981). Likelihood calculations for matched case-control studies and survival studies with tied death times. Biometrika, 68(3), 703–707.
• Harkness, W. L. (1965). Properties of the extended hypergeometric distribution. The Annals of Mathematical Statistics, 36(3), 938–945.
• Kalbfleisch, J. D. and Prentice, R. L. (2002). The Statistical Analysis of Failure Time Data, 2nd edition. Wiley, Hoboken, NJ.
• Kline, P. and Santos, A. (2012). A score based approach to wild bootstrap inference. Journal of Econometric Methods, 1(1), 23–41.
• Lin, D. Y. and Wei, L. J. (1989). The robust inference for the Cox proportional hazards model. Journal of the American Statistical Association, 84(408), 1074–1078.
• Lin, D. Y., Wei, L. J. and Ying, Z. (1993). Checking the Cox model with cumulative sums of martingale-based residuals. Biometrika, 80(3), 557–572.

The random numbers come from the C++ code, from the bits of std::mt19937_64, so the same seed gives the same draws on every platform. With seed=None, a fresh seed is drawn.

Score Residuals and the Baseline Hazard

# the score residuals of a design: row i is observation i's contribution
# to the score of the coefficients
U = coxdev.score_residuals(eta, X)
I = X.T @ (coxdev.information(eta) @ X)
robust_cov = np.linalg.solve(I, np.linalg.solve(I, U.T @ U).T)   # the sandwich

# Breslow's (or, with Efron's ties, Efron's) baseline hazard, for eta as given
H = coxdev.baseline_hazard(eta)
H.time, H.hazard, H.cumulative_hazard, H.stratum

The cumulative hazard of observation i at time t is the cumulative_hazard at the last failure time <= t (in its stratum) times exp(eta[i]). Both agree with R's survival (residuals(fit, type = "score"), which is per unit weight, and basehaz(fit, centered = FALSE)).

Schoenfeld Residuals and Expected Survival at Exit

# Schoenfeld residuals: for each failure, x_i less the mean of x over its
# risk set; zero rows for the observations that do not fail
R = coxdev.schoenfeld_residuals(eta, X)

# the expected number of events of each observation over (start, event],
# and the probability of surviving to event, given survival to start
E = coxdev.expected_events(eta)
S = coxdev.survival_at_exit(eta)          # exp(-E)
martingale = status - E

The Schoenfeld residual of a failure is its covariates less their mean over its risk set, with weights sample_weight * exp(eta) (Schoenfeld, 1982). With Efron's ties the mean is the average over the tied failures of their Efron means, as in R's survival. Weighted by sample_weight, they sum to the score. Plotted against time, or against a function of time, they check the proportional hazards assumption (Grambsch and Therneau, 1994; Therneau and Grambsch, 2000, Chapter 6). As for the score residuals, the C++ takes one covariate, an n-vector, at a time, and schoenfeld_residuals loops over the columns of X, dense or sparse.

risk_set_variance(eta, Z) gives, for each failure, the variance of each column of Z over its risk set, the V_z by which Grambsch and Therneau (1994) scale the Schoenfeld residuals: for one covariate, beta + r / V_z has mean about beta(t) at the failure time t. Weighted by sample_weight, the variances sum to the diagonal of the information. The covariance of two columns follows by polarization, (V(y + z) - V(y - z)) / 4.

The expected number of events of observation i is exp(eta[i]) times the increase of the baseline cumulative hazard (Breslow, 1972) over (start[i], event[i]], with Efron's share at a tied failure time for the failures there. status - E is the martingale residual (Barlow and Prentice, 1988; Therneau, Grambsch and Fleming, 1990), and exp(-E) the expected survival probability at exit (Therneau and Grambsch, 2000, Chapter 4). These agree with R's survival: residuals(fit, type = "schoenfeld") (which has the rows of the failures only, in order of failure time), predict(fit, type = "expected"), residuals(fit, type = "martingale") and predict(fit, type = "survival").

References

• Barlow, W. E. and Prentice, R. L. (1988). Residuals for relative risk regression. Biometrika, 75(1), 65–74.
• Breslow, N. E. (1972). Discussion of Professor Cox's paper. Journal of the Royal Statistical Society, Series B, 34(2), 216–217.
• Grambsch, P. M. and Therneau, T. M. (1994). Proportional hazards tests and diagnostics based on weighted residuals. Biometrika, 81(3), 515–526.
• Schoenfeld, D. (1982). Partial residuals for the proportional hazards regression model. Biometrika, 69(1), 239–241.
• Therneau, T. M. and Grambsch, P. M. (2000). Modeling Survival Data: Extending the Cox Model. Springer, New York.
• Therneau, T. M., Grambsch, P. M. and Fleming, T. R. (1990). Martingale-based residuals for survival models. Biometrika, 77(1), 147–160.

The C-index

C = coxdev.c_index(eta)          # Harrell's concordance index

Harrell's concordance index (Harrell et al., 1982; Harrell, Lee and Mark, 1996): for each failure i at time t, every j at risk at t is compared with it, except the other failures at t; the pair has weight w_i * w_j, is concordant if eta[i] > eta[j], and counts half if they are tied. With strata, only pairs within a stratum are compared. It is R's survival::concordance(y ~ -eta, weights = w) and, without strata, glmnet::Cindex. A two-dimensional eta gives one C-index per column.

It takes O(n log n), by the algorithm of the survival package's concordance vignette (Therneau and Atkinson): the subjects in decreasing time, with the weights of the risk set in a Fenwick tree (Fenwick, 1994) over the ranks of the predictor, counting the ordered pairs as Knight (1966) does for Kendall's tau.

References

• Fenwick, P. M. (1994). A new data structure for cumulative frequency tables. Software: Practice and Experience, 24(3), 327–336.
• Harrell, F. E., Califf, R. M., Pryor, D. B., Lee, K. L. and Rosati, R. A. (1982). Evaluating the yield of medical tests. JAMA, 247(18), 2543–2546.
• Harrell, F. E., Lee, K. L. and Mark, D. B. (1996). Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Statistics in Medicine, 15(4), 361–387.
• Knight, W. R. (1966). A computer method for calculating Kendall's tau with ungrouped data. Journal of the American Statistical Association, 61(314), 436–439.
• Therneau, T. M. and Atkinson, E. J. The concordance statistic. Vignette of the R package survival.

Different Tie-Breaking Methods

# Efron's method (default)
coxdev_efron = CoxDeviance(event=event_times, status=status, tie_breaking='efron')

# Breslow's method
coxdev_breslow = CoxDeviance(event=event_times, status=status, tie_breaking='breslow')

API Reference

CoxDeviance

The main class for computing Cox model quantities.

Parameters

• event: Event times (failure times) for each observation
• status: Event indicators (1 for event occurred, 0 for censored)
• start: Start times for left-truncated data (optional)
• strata: Stratum labels (optional): integers, strings or any sortable values; baseline_hazard reports each stratum by its label
• tie_breaking: Method for handling tied event times ('efron' or 'breslow')

Methods

• __call__(linear_predictor, sample_weight=None): Compute deviance and related quantities
• information(linear_predictor, sample_weight=None): Get information matrix as linear operator
• sample_score(linear_predictor, sample_weight=None, method='multinomial', n_draws=None, seed=None): Draws of the score, one entry per observation: a vector if n_draws is None, else an (n, n_draws) array. See Sampling the Score
• wild_score(linear_predictor, multipliers, sample_weight=None): The wild draw of the score for given multipliers
• score_residuals(linear_predictor, X, sample_weight=None): The (n, p) score residuals of a design X
• baseline_hazard(linear_predictor, sample_weight=None): The baseline hazard, a BaselineHazard with time, hazard, cumulative_hazard and stratum
• schoenfeld_residuals(linear_predictor, X, sample_weight=None): The (n, p) Schoenfeld residuals of a design X, zero for the observations that do not fail
• risk_set_variance(linear_predictor, Z, sample_weight=None): The (n, p) risk-set variances of the columns of Z at each failure, zero for the observations that do not fail
• c_index(linear_predictor, sample_weight=None): Harrell's concordance index (C-index) of a linear predictor, or one per column
• expected_events(linear_predictor, sample_weight=None): The expected number of events of each observation over (start, event]
• survival_at_exit(linear_predictor, sample_weight=None): The expected survival probability of each observation at its exit, exp(-expected_events)

CoxDevianceResult

Result object containing computation results.

Attributes

• linear_predictor: The linear predictor values used
• sample_weight: Sample weights used
• loglik_sat: Saturated log-likelihood value
• deviance: Computed deviance value
• gradient: Gradient of deviance with respect to linear predictor
• diag_hessian: Diagonal of Hessian matrix

Performance

The library is optimized for performance:

• C++ Implementation: Core computations in C++ with Eigen
• Memory Efficient: Reuses buffers and uses linear operators
• Vectorized Operations: Leverages Eigen's optimized linear algebra
• Minimal Python Overhead: Heavy computations done in C++

Building from Source

Prerequisites

• C++ compiler with C++17 support
• Eigen library headers
• For Python library: Python >= 3.9 + pybind11

C++ API

The Python package (and the R package in R/coxdev) wrap the C++ classes CoxDeviance and StratifiedCoxDeviance in cpp/include/coxdev.h. After compute_deviance(eta, weights), besides the gradient, the diagonal of the Hessian, products with the information and draws of the score, they give:

• score_residuals(x): for one covariate x (an n-vector), the n-vector of score residuals, each observation's contribution to the score x^T (-gradient / 2), centered within each risk set (with Efron's weights for ties if efron). They are linear in x, with wild_score the adjoint: G^T score_residuals(x) is x^T wild_score(G). The C++ never sees a design: applied to each column of X they give the n x p matrix U, whose crossproduct U^T U is the middle of the robust (sandwich) variance (Lin and Wei, 1989).
• baseline_hazard(): the estimate of the baseline hazard, Breslow's (1972) or, with Efron's ties, Efron's (1977), at each distinct failure time (within each stratum): its increments and cumulative sum, for the linear predictor as given, so that the cumulative hazard of observation i at time t is cumulative_hazard(t) * exp(eta_i).

cpp/tests has tests of these against direct computations:

c++ -std=c++17 -O2 -Ieigen -Icpp/include cpp/src/coxdev.cpp cpp/tests/test_residuals_hazard.cpp -o test_residuals_hazard
./test_residuals_hazard

Testing

Most of the tests compare coxdev with R's survival and glmnet packages, through rpy2. If rpy2, R or those packages are missing, those tests are skipped and the run still passes, so check the summary for skipped tests.

Running the full test suite in a fresh environment

1. Install R (4.x), and in it survival and glmnet 5.0 or later (the tests use glmnet::coxnet.deviance and glmnet:::coxgrad with cox.ties):

   Rscript -e 'install.packages(c("survival", "glmnet"), repos = "https://cloud.r-project.org")'
   Rscript -e 'packageVersion("glmnet")'   # 5.0 or later
   

2. Clone the repository with its Eigen submodule, and make a new virtual environment:

   git clone --recurse-submodules https://github.com/jonathan-taylor/coxdev.git
   cd coxdev
   python -m venv .venv
   source .venv/bin/activate          # on Windows: .venv\Scripts\activate
   python -m pip install --upgrade pip
   

   (In an existing clone without the submodule, run git submodule update --init --recursive.)

3. Build and install coxdev (this needs a C++17 compiler), with the test dependencies and rpy2:

   pip install ".[test]" rpy2
   

4. Check that rpy2 finds R and the two packages:

   python -c "from rpy2.robjects.packages import importr; importr('survival'); importr('glmnet'); print('R is ready')"
   

   • If rpy2 first prints Error importing in API mode ... Trying to import in ABI mode, its compiled interface was built for another version of R. ABI mode works; export RPY2_CFFI_MODE=ABI skips the attempt and the message.
   • If rpy2 cannot find R, set R_HOME: export R_HOME=$(R RHOME).

5. From the top of the repository, run the tests, reporting any skipped:

   python -m pytest tests -rs
   

   It takes a couple of minutes, and should end with every test passed and none skipped.

For a quick check, --test-size N runs a random sample of N of the (many) parametrized cases in each of the large comparison files:

python -m pytest tests -rs --test-size 50

Testing the R package

The R package has its own tests, against survival and glmnet, and tests/test_R_package.py checks that the R and Python packages agree (results, information and draws for the same seed). With the R package installed (R CMD INSTALL R/coxdev, which needs Rcpp, RcppEigen and R6, plus testthat for its tests):

Rscript -e 'testthat::test_local("R/coxdev")'
python -m pytest tests/test_R_package.py -rs

Contributing

1. Fork the repository
2. Create a feature branch
3. Make your changes
4. Add tests for new functionality
5. Ensure all tests pass
6. Submit a pull request

License

This project is licensed under the BSD-3-Clause License - see the LICENSE file for details.

Citation

If you use this library in your research, please cite:

@software{coxdev2024,
  title={coxdev: High-performance Cox proportional hazards deviance computation},
  author={Taylor, Jonathan and Hastie, Trevor and Narasimhan, Balasubramanian},
  year={2024},
  url={https://github.com/jonathan-taylor/coxdev}
}

Acknowledgments

• Built with Eigen for efficient linear algebra
• Uses pybind11 for Python bindings
• Inspired by the R glmnet package for survival analysis