David Cox (1972)
Cox introduced a survival regression model with unspecified baseline hazard:
\[ h(t\mid x)=h_0(t)\exp(x'\beta). \]
The baseline \(h_0(t)\) can be left unspecified because the ordered failure times imply a partial likelihood:
\[ L(\beta)=\prod_{i:\delta_i=1} \frac{\exp(x_i'\beta)}{\sum_{j:t_j\ge t_i}\exp(x_j'\beta)}. \]
This estimates relative hazards while avoiding a full parametric survival distribution.
Simulate observed event/censoring times and the risk sets induced by their ordering.
n = 180
x = rng.normal(size=n)
beta_true = 0.8
t_event = rng.exponential(scale=np.exp(-beta_true * x))
t_cens = rng.exponential(scale=2.1, size=n)
time = np.minimum(t_event, t_cens)
event = t_event <= t_cens
event.mean()np.float64(0.6277777777777778)Define the negative partial log-likelihood using the risk set \(R(t_i)=\{j:t_j\ge t_i\}\).
def neg_pl(beta):
xb = beta * x
out = 0.0
for ti, xi, di in zip(time, x, event):
if di:
out += beta * xi - special.logsumexp(xb[time >= ti])
return -out
beta_hat = optimize.minimize_scalar(neg_pl, bracket=(-1, 1)).x
beta_hatnp.float64(0.9973838746212481)Evaluate and plot the profile of the partial likelihood.
grid = np.linspace(-0.2, 1.6, 180)
pll = np.array([-neg_pl(b) for b in grid])
fig, ax = plt.subplots(figsize=(6, 3.4))
ax.plot(grid, pll - pll.max())
ax.axvline(beta_hat, color="#ffcc66", lw=2.5, label=f"partial MLE={beta_hat:.2f}")
ax.axvline(beta_true, color="#66ff99", ls="--", label=f"truth={beta_true:.2f}")
ax.set(xlabel="log hazard ratio", ylabel="relative log partial likelihood")
ax.legend()
plt.show()
lifelines CoxPHFitter, scikit-survival CoxPHSurvivalAnalysis, R survival