CoxPH

Semiparametric proportional hazards via the Cox partial likelihood

from _api_doc_utils import *

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'eta). \]

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.

Risk-set exponentials are evaluated after clipping \(x_r'\beta\) to \([-40,40]\) for stability, while the numerator uses the unclipped index. Outside that range the reported criterion, score, and Hessian are therefore no longer exact derivatives of one common clipped objective. Typical well-scaled fits stay inside the clipping region.

3 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\).

The iteration count is stored, but the optimizer does not expose a convergence flag and does not raise merely because it reached the configured maximum. Model-based inference assumes independent subjects and a correctly specified proportional-hazards relative-risk model.

4 Performance and numerical behavior

The implementation rebuilds every risk set by scanning all \(n\) rows for every event and accumulates a dense \(p\times p\) second moment. With \(E\) events, runtime is approximately \(O(Enp^2)\), substantially slower than cumulative-risk-set implementations. The information solve adds \(O(p^3)\) per Newton iteration. Prediction exponentiates the unbounded new-data index without clipping, so extreme values can overflow. There is no baseline-hazard storage cost because no baseline is estimated.

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

print(inspect.signature(cm.CoxPH))

()

cls = cm.CoxPH
display(HTML(html_table(["Public method"], public_methods(cls))))

Public method
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, /)

6 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)
print(model.predict_lin(x[:5]))
print(model.predict(x[:5]))

[ 0.23097646  0.17232858 -0.67524493 -0.14989842  0.51105657]
[1.25982958 1.18806814 0.50903173 0.86079541 1.66705162]

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

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	()
iterations	()