Metadata-Version: 2.1
Name: coxdev
Version: 0.3.0
Summary: Library for computing Cox deviance
Keywords: Cox model
Author-Email: Jonathan Taylor <jonathan.taylor@stanford.edu>, Trevor Hastie <hastie@stanford.edu>, Balasubramanian Narasimhan <naras@stanford.edu>
Maintainer-Email: Jonathan Taylor <jonathan.taylor@stanford.edu>
License: BSD-3-Clause
Classifier: Development Status :: 3 - Alpha
Classifier: Environment :: Console
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python
Classifier: Topic :: Scientific/Engineering
Requires-Python: >=3.9
Requires-Dist: numpy
Requires-Dist: scipy
Provides-Extra: test
Requires-Dist: pytest; extra == "test"
Requires-Dist: pandas; extra == "test"
Provides-Extra: dev
Requires-Dist: pytest; extra == "dev"
Requires-Dist: scipy; extra == "dev"
Requires-Dist: ninja; extra == "dev"
Requires-Dist: meson; extra == "dev"
Requires-Dist: meson-python; extra == "dev"
Requires-Dist: pybind11; extra == "dev"
Requires-Dist: setuptools_scm; extra == "dev"
Provides-Extra: doc
Requires-Dist: myst-parser; extra == "doc"
Requires-Dist: jupyter-book; extra == "doc"
Requires-Dist: pandas; extra == "doc"
Requires-Dist: matplotlib; extra == "doc"
Requires-Dist: rpy2; extra == "doc"
Requires-Dist: glmstar>=0.2.0; extra == "doc"
Requires-Dist: sphinx; extra == "doc"
Requires-Dist: sphinx-book-theme; extra == "doc"
Description-Content-Type: text/markdown

# 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:
   ```bash
   git submodule update --init --recursive
   ```


### Standard Installation

```bash
pip install .
```

### With Custom Eigen Path

If you have Eigen installed elsewhere, you can specify its location:
```bash
env EIGEN_LIBRARY_PATH=/path/to/eigen pip install .
```

### Development Installation

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

## Quick Start

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

```bash
R CMD INSTALL R/coxdev
```

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

```r
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](https://jupytext.readthedocs.io)
notebooks, with executed copies in `doc/*.ipynb`. The R package has the
same two as vignettes.

## Advanced Usage

### Left-Truncated Data

```python
# 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

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

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

```python
# 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

```python
# 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

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

```python
# 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](#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:

```bash
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](https://rpy2.github.io/). 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`):

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

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

   ```bash
   pip install ".[test]" rpy2
   ```

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

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

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

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

```bash
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](LICENSE) file for details.

## Citation

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

```bibtex
@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](http://eigen.tuxfamily.org/) for efficient linear algebra
- Uses [pybind11](https://pybind11.readthedocs.io/) for Python bindings
- Inspired by the R `glmnet` package for survival analysis
