Cox proportional hazards family for survival analysis.
CoxPH(
status="event",
ties="breslow",
)
CoxPH fits a semiparametric proportional hazards model, allowing smooth (via s()) and linear covariate effects on the log hazard while leaving the baseline hazard h0(t) unspecified. Use this family whenever the response is a time-to-event outcome that may be right-censored — for example, time to failure, time to churn, or time to death/relapse in survival data — and the goal is to model how covariates shift the instantaneous risk of the event over time. Rather than a per-observation deviance and log-likelihood in the usual GLM sense, CoxPH maximizes the Cox partial likelihood via a custom irls_update, so the “response” y passed to GAM.fit() is the observed survival/censoring time, and the event indicator is supplied separately through set_data() (populated automatically by GAM.fit() from the column named by status).
Parameters
status: str = "event"-
Name of the column in the data dict containing the event indicator (1 = event observed, 0 = right-censored). This column is looked up automatically from the data passed to GAM.fit().
ties: str = "breslow"-
Tie-handling method for the partial likelihood when multiple observations share the same event time:
"breslow" (default, simpler and faster) or "efron" (more accurate when ties are frequent).
Notes
The hazard is modeled multiplicatively as
h(t \mid x) = h_0(t)\, e^{\eta(x)}, \qquad \eta(x) = X\beta,
where h0(t) is an unspecified baseline hazard and eta is the (possibly smooth) linear predictor, so link and link_inverse are both the identity. There is no closed-form variance function or unit deviance in the usual GLM sense; instead, model fitting maximizes the Cox partial log-likelihood,
\ell(\beta) = \sum_{i:\, \delta_i = 1}
\left[ \eta_i - \log\!\left( \sum_{j \in R(t_i)} e^{\eta_j} \right) \right],
where \delta_i is the event indicator and R(t_i) is the risk set at time t_i (those still under observation just before t_i). The overall “deviance” reported by GAM.summary() is -2\ell(\beta). After fitting, baseline_hazard() and survival_function() expose the Breslow estimate of the cumulative baseline hazard and the implied survival curve.
Examples
Fit a Cox proportional hazards GAM with a smooth effect of age on the hazard:
import numpy as np
import whittaker as wk
rng = np.random.default_rng(0)
n = 300
age = rng.uniform(40, 80, n)
risk = np.exp(0.03 * (age - 60))
time = rng.exponential(1.0 / risk)
censor_time = rng.exponential(2.0, n)
observed_time = np.minimum(time, censor_time)
event = (time <= censor_time).astype(float)
data = {"time": observed_time, "age": age, "event": event}
model = wk.GAM("time ~ s(age)", family=wk.CoxPH(status="event"))
model.fit(data, method="REML")
print(model.summary())
GAM fit summary
============================================================
Formula: time ~ s(age)
Family: CoxPH(ties='breslow')
Inference: REML
Observations: 300
Coefficients: 10
Parametric coefficients:
Term Estimate Std.Err z value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) 0.0787 0.0731 1.077 0.2817
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(age) 1.00 2 23.722 7.06e-06
Total EDF: 2.00
Scale est: 1.000000
Deviance: 1824.6258
Null dev: 1848.4760
Dev. expl: 1.3%
GCV score: 6.164023
AIC: 1828.63
BIC: 1836.04
Attributes
|
Name
|
Description
|
|
scale_known
|
Whether the dispersion parameter is fixed rather than estimated.
|
|
ties
|
Tie-handling method used for the partial likelihood.
|
scale_known
Whether the dispersion parameter is fixed rather than estimated.
ties
Tie-handling method used for the partial likelihood.
Methods
|
Name
|
Description
|
|
baseline_hazard()
|
Breslow estimate of the cumulative baseline hazard H_0(t).
|
|
deviance()
|
Total deviance D = -2\,\ell(\beta) based on the Cox partial log-likelihood.
|
|
initialize()
|
Store the observed survival times and seed the linear predictor at zero.
|
|
irls_update()
|
Newton step for the Cox partial likelihood, supplying custom working values.
|
|
link()
|
Identity link, \eta = \mu.
|
|
link_derivative()
|
Derivative of the identity link, g'(\mu) = 1.
|
|
link_inverse()
|
Inverse identity link, \mu = \eta.
|
|
log_likelihood()
|
Cox partial log-likelihood \ell(\beta) at the current linear predictor.
|
|
set_data()
|
Extract the event/censoring indicator from the fitting data.
|
|
simulate()
|
Simulate survival times by inverting the fitted cumulative baseline hazard.
|
|
survival_function()
|
Fitted survival probability at the last event time.
|
|
unit_deviance()
|
Per-observation deviance placeholder (not decomposable for a partial likelihood).
|
|
variance()
|
Variance function V(\mu) = 1.
|
baseline_hazard()
Breslow estimate of the cumulative baseline hazard H_0(t).
Returns the step-function estimate of H_0(t) = \int_0^t h_0(u)\, du evaluated at each observed event time, computed by _compute_baseline_hazard from the most recent linear predictor seen during fitting.
Returns
tuple of NDArray-
(times, cumhaz): the sorted unique event times and the corresponding cumulative baseline hazard values, each shape (n_events,).
Raises
RuntimeError-
If the model has not yet been fitted, so no baseline hazard exists.
deviance()
Total deviance D = -2\,\ell(\beta) based on the Cox partial log-likelihood.
deviance(
y,
mu,
*,
weights=None,
)
Implements the family-specific version of Family.deviance. Since there is no saturated-model decomposition for the partial likelihood, this simply returns minus twice the partial log-likelihood evaluated at the current linear predictor, which GAM.summary() reports as the model deviance.
Parameters
y: NDArray-
Observed survival/censoring times, shape (n,) (unused; retained for interface compatibility).
mu: NDArray-
Current linear predictor values (eta), shape (n,).
weights: NDArray | None = None-
Unused; CoxPH does not support prior weights.
Returns
float-
-2\ell(\beta), the deviance.
initialize()
Store the observed survival times and seed the linear predictor at zero.
Implements the family-specific version of Family.initialize. Records the sorted order of the observed times y (used throughout this file to walk risk sets from latest to earliest time) and starts P-IRLS from eta = 0 for every observation, i.e. an initial hazard ratio of exp(0) = 1 relative to the baseline.
Parameters
y: NDArray-
Observed survival/censoring times, shape
(n,).
Returns
NDArray-
Array of zeros with the same shape as
y, used as the initial linear predictor.
irls_update()
Newton step for the Cox partial likelihood, supplying custom working values.
irls_update(
y,
mu,
eta,
)
Implements the family-specific version of Family.irls_update. Rather than deriving the working response and weights from link, link_derivative, and variance, this computes the gradient U(\eta) and diagonal Hessian H(\eta) of the partial log-likelihood directly (via _breslow_grad_hess or _efron_grad_hess, depending on ties), and forms a Newton update
z = \eta + H(\eta)^{-1} U(\eta), \qquad W = \max(H(\eta), \epsilon),
so that a weighted least-squares fit of z on W reproduces one Newton-Raphson step on the partial log-likelihood. As a side effect, also refreshes the Breslow estimate of the cumulative baseline hazard via _compute_baseline_hazard, so that baseline_hazard() and survival_function() reflect the current fit.
Parameters
y: NDArray-
Observed survival/censoring times, shape (n,) (unused directly; times were captured by initialize).
mu: NDArray-
Current fitted values, identical to eta for this family, shape (n,).
eta: NDArray-
Current linear predictor values, shape
(n,).
Returns
tuple of NDArray-
(z, W), the pseudo-response and diagonal working-weight vector, each shape (n,).
Raises
RuntimeError-
If called before set_data() has supplied the event indicator.
link()
Identity link, \eta = \mu.
CoxPH does not model a conventional conditional mean; mu is instead identified with the linear predictor eta, so the link function is the identity.
Parameters
mu: NDArray-
Linear predictor values (treated as
mu), shape (n,).
Returns
NDArray-
The input array, unchanged.
link_derivative()
Derivative of the identity link, g'(\mu) = 1.
Since link and link_inverse are both the identity, the derivative is constant 1 everywhere.
Parameters
mu: NDArray-
Conditional mean values, shape
(n,).
Returns
NDArray-
Array of ones with the same shape as
mu.
link_inverse()
Inverse identity link, \mu = \eta.
Since link is the identity, its inverse is also the identity: the linear predictor is returned unchanged as mu.
Parameters
eta: NDArray-
Linear predictor values, shape
(n,).
Returns
NDArray-
The input array, unchanged.
log_likelihood()
Cox partial log-likelihood \ell(\beta) at the current linear predictor.
log_likelihood(
y,
mu,
scale,
*,
weights=None,
)
Implements the family-specific version of Family.log_likelihood. Delegates to _partial_log_likelihood, which sums \eta_i - \log \sum_{j \in R(t_i)} e^{\eta_j} over the observed events, using the Breslow or Efron tie-handling scheme selected by ties.
Parameters
y: NDArray-
Observed survival/censoring times, shape (n,) (unused; retained for interface compatibility).
mu: NDArray-
Current linear predictor values (eta), shape (n,).
scale: float-
Unused; CoxPH has no free dispersion parameter (see scale_known).
weights: NDArray | None = None-
Unused; CoxPH does not support prior weights.
Returns
float-
The partial log-likelihood \ell(\beta).
set_data()
Extract the event/censoring indicator from the fitting data.
Implements the family-specific version of Family.set_data. Looks up the column named by the status constructor argument in data and stores it as the event indicator used by irls_update, deviance, and log_likelihood to identify which observations are actual events (1) versus right-censored (0).
Parameters
data: dict[str, NDArray]-
Mapping of column names to arrays, as passed to GAM.fit().
Raises
KeyError-
If the status column is not present in
data.
simulate()
Simulate survival times by inverting the fitted cumulative baseline hazard.
simulate(
mu,
scale,
rng,
)
Implements the family-specific version of Family.simulate. Draws u \sim \mathrm{Uniform}(0, 1) for each observation and solves for the target cumulative baseline hazard H_0(t) = -\log(u) / e^{\eta} implied by the proportional-hazards model, then looks up the smallest fitted baseline event time whose cumulative baseline hazard (from baseline_hazard()) meets or exceeds that target. Observations whose target exceeds the largest observed cumulative hazard are assigned the last baseline time.
Parameters
mu: NDArray-
Current linear predictor values (eta), shape (n,).
scale: float-
Unused; CoxPH has no free dispersion parameter.
rng: np.random.Generator-
A
numpy.random.Generator instance used to draw the uniform variates.
Returns
NDArray-
Simulated survival times, shape
(n,).
Raises
RuntimeError-
If called before the model has been fitted (so no baseline hazard is available).
survival_function()
Fitted survival probability at the last event time.
S(t \mid x) = \exp(-H_0(t)\, e^{\eta}).
Combines the given linear predictor with the final value of the Breslow cumulative baseline hazard (i.e. H_0 evaluated at the largest observed event time) to give the implied survival probability for each observation at that time, following the proportional-hazards relationship S(t \mid x) = S_0(t)^{\exp(\eta)}.
Parameters
eta: NDArray-
Linear predictor values, shape
(n,).
Returns
NDArray-
Fitted survival probabilities at the last observed event time, shape
(n,).
Raises
RuntimeError-
If the model has not yet been fitted, so no baseline hazard exists.
unit_deviance()
Per-observation deviance placeholder (not decomposable for a partial likelihood).
Implements the family-specific version of Family.unit_deviance. The Cox partial log-likelihood does not decompose into independent per-observation contributions the way an ordinary GLM deviance does, so this returns an array of ones purely to satisfy the Family interface; use deviance or log_likelihood for the actual fit statistics.
Parameters
y: NDArray-
Observed survival/censoring times, shape (n,).
mu: NDArray-
Current linear predictor values, shape
(n,).
Returns
NDArray-
Array of ones with the same shape as
y.
variance()
Variance function V(\mu) = 1.
CoxPH fits via the Cox partial likelihood rather than a GLM variance function, so this returns a constant 1 for every observation; it exists only to satisfy the Family interface and is not used by irls_update, which supplies its own working weights directly from the partial-likelihood Hessian.
Parameters
mu: NDArray-
Conditional mean values, shape
(n,).
Returns
NDArray-
Array of ones with the same shape as
mu.