from _api_doc_utils import *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-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 Implementation walkthrough
The partial likelihood is evaluated by literal event-by-event risk-set scans; there is no delegated survival engine or cumulative-sum optimization.
fit()clears prior state, validates finite \(X\), positive times, and binary events, and initializes all slopes at zero. It does not sort the observations because every risk set is reconstructed from comparisons with the event time.At a candidate \(\beta\), the code computes all raw indices \(\eta=X\beta\) once and all risk scores as \(\exp\{\operatorname{clip}(\eta,-40,40)\}\). It then loops over rows and skips every non-event.
For an event at \(t_i\), an inner scan includes row \(r\) when \(t_r\geq t_i\). It accumulates the denominator, the risk-weighted first moment \(S^{(1)}\), and the full risk-weighted second moment \(S^{(2)}\). The 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. Separate tied events repeat the same denominator, which is exactly the stated Breslow contribution.
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}\).
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.
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.
This direct implementation mirrors the mathematics closely enough to inspect each risk-set moment, but it does redundant work for sorted event times and ties. A production implementation would normally sort once, update cumulative risk-set sums, and compute a Breslow baseline after fitting.
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
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.
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.
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, /) |
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.2712229261812}
[ 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.
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 |
() |