CoxPH

Semiparametric proportional hazards via the Cox partial likelihood

1 Where it fits

Group: Survival / event-time models

CoxPH estimates log hazard ratios without parameterizing the baseline hazard:

\[ h_i(t \mid x_i) = h_0(t) \exp(x_i'\beta). \]

That means the paved prediction target is relative risk rather than absolute survival. The latent linear index is still useful and is exposed as predict_lin(x).

2 Partial likelihood and ties

For event time \(t_i\) and risk set \(R(t_i)=\{r:t_r\geq t_i\}\), the implementation maximizes

\[ \ell_p(\beta) = \sum_{i:d_i=1} \left[ x_i'\beta -\log\left\{\sum_{r\in R(t_i)}\exp(x_r'\beta)\right\} \right]. \]

When several events share a time, the same full risk-set denominator is subtracted once per event. This is the Breslow partial-likelihood treatment of ties. Newton updates use the analytic score and information, a \(10^{-8}\) ridge in the solve, coordinate-wise step clipping to \([-1,1]\), and backtracking. Convergence requires either \(\|\nabla\ell_p\|_\infty<\tau\) or an accepted partial-likelihood change below \(\tau(1+|\ell_p|)\).

Risk sums use dynamic log rescaling instead of index clipping, so likelihood, score, and Hessian refer to the same Breslow objective. Covariates are centered internally to improve origin invariance; the fitted slopes and original-scale prediction contract are unchanged.

3 Implementation walkthrough

The package owns a cumulative risk-moment kernel and Newton solver.

  1. fit() clears state, validates finite \(X\), positive times, binary events, and at least one event. It copies owned arrays before detaching the pure-Rust fit from the GIL. Slopes start at zero and descending stop-time order is cached once.

  2. At each candidate, the kernel computes centered indices and visits groups of equal stop time in descending order. Adding each group maintains the risk set \(t_r\geq t_i\) before any tied events are evaluated.

  3. Dynamic rescaling maintains the risk mass, first moment, and second moment without clipping. In unscaled notation, an event contributes

    \[ \eta_i-\log S^{(0)},\qquad x_i-S^{(1)}/S^{(0)},\qquad -\{S^{(2)}/S^{(0)}-\bar x\bar x'\} \]

    to likelihood, score, and Hessian. Tied events share the accumulated denominator, giving the Breslow contribution.

  4. The Newton helper subtracts \(10^{-8}\) from the log-likelihood Hessian diagonal and explicitly inverts it. Coordinates of \(H^{-1}g\) are clipped to \([-1,1]\), a tighter bound than in the parametric models. Up to 30 halvings seek a finite partial likelihood no smaller than the current value minus \(10^{-10}\).

  5. The loop stops before a step when the maximum score is below tolerance or after a step when the likelihood change is relatively small. An exhausted budget or failed line search raises. Final covariance reevaluates the Hessian, negates it, and inverts it without the optimization ridge.

  6. Prediction performs only \(X\hat\beta\) and its exponential. No baseline hazard is estimated as a post-fit step, so no absolute-risk or survival prediction is available.

Backtracking uses likelihood-only evaluations, skipping derivative accumulation. A direct Breslow reference and finite-difference derivative tests cover ties and the shared start/stop kernel. No baseline hazard is added by this refactor.

4 Inference and prediction

The covariance is the inverse observed partial-likelihood information. Standard errors, normal \(z\) statistics, and approximate two-sided chi-square-one \(p\)-values are derived from it. There is no robust sandwich, clustering, strata, weights, or bootstrap. The class does not estimate a baseline cumulative hazard, so it cannot return survival probabilities; default prediction is \(\exp(x'\hat\beta)\) and the linear method returns \(x'\hat\beta\).

Reaching max_iterations raises ValueError before covariance or fitted state is stored. A successful summary exposes converged, iterations, termination_reason, and objective, where objective=-\ell_p(\hat\beta). Model-based inference assumes independent subjects and a correctly specified proportional-hazards relative-risk model.

5 Performance and numerical behavior

Sorting costs \(O(n\log n)\) once; derivative evaluation then costs \(O(np^2)\) and the Newton solve adds \(O(p^3)\). Likelihood-only trials avoid second moments. Owned native fitting releases the GIL but does not implement cancellation or signal polling. Prediction still exponentiates the original-scale new-data index, which can overflow; use predict_lin() for extreme values.

6 Python API

Constructor: cm.CoxPH

Use fit(x, time, event). predict_lin(x) returns the log hazard ratio and predict_relative_risk(x) exponentiates it. The default predict(x) is the same relative-risk object. summary() reports coefficients, standard errors, hazard ratios, and optimization diagnostics.

Code 
print(inspect.signature(cm.CoxPH))
()
Code 
cls = cm.CoxPH
display(HTML(html_table(["Public method"], public_methods(cls))))
Public method
CoxPH()
fit(self, /, x, time, event, max_iterations=50, tolerance=1e-08)
predict(self, /, x)
predict_lin(self, /, x)
predict_log_hazard_ratio(self, /, x)
predict_relative_risk(self, /, x)
summary(self, /)

7 Minimal example

rng=np.random.default_rng(33)
x=rng.normal(size=(260,2)); rate=0.05*np.exp(x@np.array([0.5,-0.25])); t_event=rng.exponential(1.0/rate); c=rng.exponential(25,size=260)
time=np.minimum(t_event,c); event=(t_event<=c).astype(float)
model=cm.CoxPH(); model.fit(x,time,event)
fit=model.summary(); print({key: fit[key] for key in ['converged', 'iterations', 'termination_reason', 'objective']})
print(model.predict_lin(x[:5]))
print(model.predict(x[:5]))
{'converged': True, 'iterations': 3, 'termination_reason': 'Relative objective tolerance reached', 'objective': 655.2712229261809}
[ 0.23097646  0.17232858 -0.67524493 -0.14989842  0.51105657]
[1.25982958 1.18806814 0.50903173 0.86079541 1.66705162]

8 summary() contract

The table below is generated by fitting the live class in this repository and then inspecting summary(). Shapes are shown because most values are plain NumPy arrays or scalars.

Code 
rng=np.random.default_rng(133); x=rng.normal(size=(130,2)); rate=0.06*np.exp(x@np.array([0.45,-0.2])); te=rng.exponential(1.0/rate); c=rng.exponential(20,size=130); time=np.minimum(te,c); event=(te<=c).astype(float)
model=cm.CoxPH(); model.fit(x,time,event)
summary = model.summary()
display(HTML(html_table(["summary() key", "shape"], summary_shape_rows(summary))))
summary() key shape
model ()
coef (2,)
hazard_ratio (2,)
se (2,)
z (2,)
p_value (2,)
vcov (2, 2)
log_likelihood ()
converged ()
iterations ()
termination_reason ()
objective ()