Estimate a Survival Curve with Cox PH

Fit a GAM-extended Cox proportional hazards model for time-to-event data.

Use the Cox proportional hazards family when your outcome is the time until an event (death, equipment failure, customer churn, or relapse) and some observations may be right-censored (the event had not yet occurred when follow-up ended). The GAM extension allows the log-hazard to be a smooth function of continuous covariates like age, avoiding the need to categorize or assume linearity. Parametric terms like treatment enter as log-hazard ratios, directly interpretable as relative hazard multipliers.

Load data

Load the survival dataset and inspect its columns.

import whittaker as wk

data = wk.load_dataset("survival", as_frame=True)
data.columns.tolist()

['time', 'event', 'age', 'treatment']

The dataset has four columns: time, event, age, and treatment.

data.head()

	time	event	age	treatment
0	5.0	0.0	58.129296	0.0
1	5.0	0.0	70.374621	1.0
2	5.0	0.0	64.905856	0.0
3	5.0	0.0	40.134324	0.0
4	5.0	0.0	43.507483	0.0

The event column is a 0/1 indicator where 1 means the event occurred and 0 means the observation was censored. The CoxPH family reads the event column automatically from the data by name, alongside the time response in the formula.

Fit

Fit a Cox PH GAM with a smooth over age and a linear parametric term for treatment.

model = wk.GAM(
    "time ~ s(age) + treatment",
    family=wk.CoxPH(),
).fit(data)

model.summary()

GAM fit summary
============================================================
Formula:    time ~ s(age) + treatment
Family:     CoxPH(ties='breslow')
Inference:  GCV
Observations: 250
Coefficients: 11

Parametric coefficients:
  Term                       Estimate    Std.Err    z value    p-value
  ------------------------ ---------- ---------- ---------- ----------
  (Intercept)                  0.3934     0.1156      3.402  0.0006686
  treatment                   -0.6502     0.1770     -3.674  0.0002384

Approximate significance of smooth terms:
  Term                        EDF Ref.df     Chi.sq    p-value
  ------------------------ ------ ------ ---------- ----------
  s(age)                     1.09      2     13.872  0.0009724

Total EDF:  3.09
Scale est:  1.000000
Deviance:   1356.2592
Null dev:   1381.1403
Dev. expl:  1.8%
GCV score:  5.561854
AIC:        1362.45
BIC:        1373.34

The summary reports the coefficient for treatment as a log-hazard ratio: exponentiate it to obtain the hazard ratio. The EDF for s(age) indicates how non-linearly age affects the log-hazard (an EDF near 1 implies a roughly log-linear age effect).

Partial effects

wk.partial_effects(model)

Partial effects are on the log-hazard scale. For s(age), an upward slope at a given age means the instantaneous risk of the event increases with age in that region; a downward slope means decreasing risk. The treatment panel shows the linear log-hazard contrast between treatment levels. Centered partial effects are relative to the mean log-hazard in the data, so positive values indicate above-average hazard.

Diagnostics

wk.check(model)

For Cox models, the diagnostic plots check the proportional hazards assumption and the distribution of martingale or deviance residuals. Systematic trends in residuals over time can indicate time-varying effects not captured by the model.