# CoxPH


Cox proportional hazards family for survival analysis.


Usage

``` python
CoxPH(
    status="event",
    ties="breslow",
)
```


[CoxPH](CoxPH.md#whittaker.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](CoxPH.md#whittaker.CoxPH) maximizes the Cox partial likelihood via a custom `irls_update`, so the "response" `y` passed to [GAM.fit()](GAM.md#whittaker.GAM.fit) is the observed survival/censoring time, and the event indicator is supplied separately through [set_data()](CoxPH.md#whittaker.CoxPH.set_data) (populated automatically by [GAM.fit()](GAM.md#whittaker.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()](GAM.md#whittaker.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()](GAM.md#whittaker.GAM.summary) is -2\ell(\beta). After fitting, [baseline_hazard()](CoxPH.md#whittaker.CoxPH.baseline_hazard) and [survival_function()](CoxPH.md#whittaker.CoxPH.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:


``` python
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](#scale_known) | Whether the dispersion parameter is fixed rather than estimated. |
| [ties](#ties) | Tie-handling method used for the partial likelihood. |

------------------------------------------------------------------------


### scale_known


Whether the dispersion parameter is fixed rather than estimated.


`scale_known: bool`


------------------------------------------------------------------------


### ties


Tie-handling method used for the partial likelihood.


`ties: str`


## Methods

| Name | Description |
|----|----|
| [baseline_hazard()](#baseline_hazard) | Breslow estimate of the cumulative baseline hazard H_0(t). |
| [deviance()](#deviance) | Total deviance D = -2\\\ell(\beta) based on the Cox partial log-likelihood. |
| [initialize()](#initialize) | Store the observed survival times and seed the linear predictor at zero. |
| [irls_update()](#irls_update) | Newton step for the Cox partial likelihood, supplying custom working values. |
| [link()](#link) | Identity link, \eta = \mu. |
| [link_derivative()](#link_derivative) | Derivative of the identity link, g'(\mu) = 1. |
| [link_inverse()](#link_inverse) | Inverse identity link, \mu = \eta. |
| [log_likelihood()](#log_likelihood) | Cox partial log-likelihood \ell(\beta) at the current linear predictor. |
| [set_data()](#set_data) | Extract the event/censoring indicator from the fitting data. |
| [simulate()](#simulate) | Simulate survival times by inverting the fitted cumulative baseline hazard. |
| [survival_function()](#survival_function) | Fitted survival probability at the last event time. |
| [unit_deviance()](#unit_deviance) | Per-observation deviance placeholder (not decomposable for a partial likelihood). |
| [variance()](#variance) | Variance function V(\mu) = 1. |

------------------------------------------------------------------------


### baseline_hazard()


Breslow estimate of the cumulative baseline hazard H_0(t).


Usage

``` python
baseline_hazard()
```


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.


Usage

``` python
deviance(
    y,
    mu,
    *,
    weights=None,
)
```


Implements the family-specific version of [Family.deviance](Family.md#whittaker.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()](GAM.md#whittaker.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](CoxPH.md#whittaker.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.


Usage

``` python
initialize(y)
```


Implements the family-specific version of [Family.initialize](Family.md#whittaker.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.


Usage

``` python
irls_update(
    y,
    mu,
    eta,
)
```


Implements the family-specific version of [Family.irls_update](Family.md#whittaker.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](CoxPH.md#whittaker.CoxPH.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()](CoxPH.md#whittaker.CoxPH.baseline_hazard) and [survival_function()](CoxPH.md#whittaker.CoxPH.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()](CoxPH.md#whittaker.CoxPH.set_data) has supplied the event indicator.


------------------------------------------------------------------------


### link()


Identity link, \eta = \mu.


Usage

``` python
link(mu)
```


[CoxPH](CoxPH.md#whittaker.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.


Usage

``` python
link_derivative(mu)
```


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.


Usage

``` python
link_inverse(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.


Usage

``` python
log_likelihood(
    y,
    mu,
    scale,
    *,
    weights=None,
)
```


Implements the family-specific version of [Family.log_likelihood](Family.md#whittaker.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](CoxPH.md#whittaker.CoxPH.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](CoxPH.md#whittaker.CoxPH) has no free dispersion parameter (see `scale_known`).

`weights: NDArray | None = None`  
Unused; [CoxPH](CoxPH.md#whittaker.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.


Usage

``` python
set_data(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()](GAM.md#whittaker.GAM.fit).


#### Raises


`KeyError`  
If the status column is not present in `data`.


------------------------------------------------------------------------


### simulate()


Simulate survival times by inverting the fitted cumulative baseline hazard.


Usage

``` python
simulate(
    mu,
    scale,
    rng,
)
```


Implements the family-specific version of [Family.simulate](Family.md#whittaker.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()](CoxPH.md#whittaker.CoxPH.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](CoxPH.md#whittaker.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.


Usage

``` python
survival_function(eta)
```


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


Usage

``` python
unit_deviance(
    y,
    mu,
)
```


Implements the family-specific version of [Family.unit_deviance](Family.md#whittaker.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](Family.md#whittaker.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.


Usage

``` python
variance(mu)
```


[CoxPH](CoxPH.md#whittaker.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](Family.md#whittaker.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`.