Cox model diagnostics

Fitting a Cox model is only half the work. Before you trust its hazard ratios, you should check that its central assumption holds, look at residuals for influential or poorly fit observations, and understand what the model predicts. This page covers the proportional hazards test, residuals, baseline hazard and survival prediction, stratification, robust standard errors, and the concordance index. Together these turn a fitted model into a defensible one.

We begin with the lung outcomes. The response y is a Surv object pairing each follow-up time with an event indicator, and it is the target every model on this page is fit against.

import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
y

Surv(type=right, n=228, events=165)

The printed response records 165 events among 228 subjects, with a + marking each censored observation. We fit a Cox model to it using age and sex as covariates. This is the model we diagnose throughout the page, so we fit it once here and reuse it below.

cox = gw.CoxPH().fit(y, lung[["age", "sex"]])

Checking the proportional hazards assumption

The Cox model assumes each covariate’s effect on the hazard is constant over time. When this fails, for example if a treatment helps early but not late, the reported hazard ratio is a misleading average. The Grambsch-Therneau test checks the assumption by looking for a trend in the scaled Schoenfeld residuals against time. A small p-value is evidence that the assumption is violated for that covariate.

Calling cox_zph returns a ZPHResult, an object that bundles the test statistics together. Displaying it directly gives a compact summary of the test.

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])

zph = cox.cox_zph()
zph

ZPHResult(transform='identity', age: p=0.7065, sex: p=0.0992, GLOBAL p=0.2425)

For the full per-covariate breakdown, call to_frame() on that result. This is the table you will read most often.

zph.to_frame(format="polars")

	term str	chisq f64	df i64	p_value f64
PolarsRows3Columns4
0	age	0.141747561345	1	0.70654984373
1	sex	2.71840108815	1	0.0991973377972
2	GLOBAL	2.8333913831	2	0.242514035592

The table has one row per covariate plus a GLOBAL row that tests all covariates jointly. Large p-values, as here, are reassuring: they give no evidence against proportional hazards. The test uses a time transform, which defaults to "identity"; "log" is also available.

A non-significant test is not proof

Failing to reject the proportional hazards assumption does not prove it holds, especially in small samples. Complement the test by plotting scaled Schoenfeld residuals or by fitting time-stratified models when you have reason to suspect a time-varying effect.

Residuals

Residuals reveal how individual observations relate to the fitted model. Martingale residuals, one per subject, are the difference between the observed number of events and the number the model expected; large negative values flag subjects who lived much longer than predicted. Schoenfeld residuals, one per event, underlie the proportional hazards test and are useful for spotting time trends in a covariate’s effect.

The martingale residuals come back as a NumPy array with one entry per subject. There are 228 of them, so we look at the first five rather than the whole array.

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])

cox.residuals("martingale")[:5]

array([ 0.00438999, -0.50576203, -3.12981924,  0.53275397, -2.35065745])

The Schoenfeld residuals come back as a DataFrame instead, with one row per event and one column per covariate.

cox.residuals("schoenfeld", format="polars")

	age f64	sex f64
PolarsRows165Columns2
0	0.935464711892	0.727079613965
1	10.0052299319	-0.272302937094
2	17.0052299319	-0.272302937094
3	3.00522993189	-0.272302937094
4	10.1404536933	-0.275789612537
5	12.227706252	-0.278411033642
6	1.22770625199	-0.278411033642
7	5.29450657479	-0.280979164464
162	-2.90825326766	-0.155341657484
163	0.611217430947	-0.196097246718
164	-4.65909168016	-0.171473016943

Baseline hazard and predicted survival

Although the Cox model does not assume a shape for the baseline hazard, it can estimate one after fitting. The baseline cumulative hazard, and the survival curve derived from it, describe a reference subject.

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])

cox.baseline_hazard(format="polars")

	time f64	cumhaz f64	survival f64
PolarsRows186Columns3
0	5	0.00295531012895	0.997049052501
1	11	0.0119063396838	0.988164260305
2	12	0.0149282203866	0.985182653095
3	13	0.0210294286393	0.979190147912
4	15	0.0241081722771	0.976180108421
5	26	0.0272054191316	0.973161315038
6	30	0.0303227228765	0.970132399105
7	31	0.0334599993722	0.967093594809
183	965	2.01305845682	0.133579502131
184	1010	2.01305845682	0.133579502131
185	1022	2.01305845682	0.133579502131

More useful in practice is predicting the survival curve for specific covariate values. First we build a small frame of the subjects we want predictions for: a 50-year-old with sex 1 and a 70-year-old with sex 2. Showing it confirms the covariate values before we predict from them.

import pandas as pd

newdata = pd.DataFrame({"age": [50, 70], "sex": [1, 2]})
newdata

	age i64	sex i64
PandasRows2Columns2
0	50	1
1	70	2

Given that new data and a set of times, predict returns the estimated survival probability for each subject.

cox.predict(newdata, type="survival", times=[180, 365, 730], format="polars")

	time f64	subject_1 f64	subject_2 f64
PolarsRows3Columns3
0	180	0.7297626204	0.767074757103
1	365	0.417781082587	0.479674113531
2	730	0.119648956441	0.167440538104

Each column is a subject from newdata, and each row is a requested time. The type="lp" and type="risk" options instead return the linear predictor and the relative risk.

Confidence bands and conditional survival

Pass ci=True to add a pointwise confidence band around each predicted curve. Every subject then gains _lower and _upper columns, built from the standard error of the cumulative hazard and matching R’s survfit (which combines baseline-hazard and coefficient uncertainty).

cox.predict(newdata, type="survival", times=[180, 365, 730], ci=True, format="polars")

	time f64	subject_1 f64	subject_1_lower f64	subject_1_upper f64	subject_2 f64	subject_2_lower f64	subject_2_upper f64
PolarsRows3Columns7
0	180	0.7297626204	0.65075626184	0.818360903092	0.767074757103	0.696725970364	0.844526697744
1	365	0.417781082587	0.312492619485	0.558544497004	0.479674113531	0.379184415457	0.6067951261
2	730	0.119648956441	0.0574578159532	0.249154489079	0.167440538104	0.0910524175713	0.307914216319

You can also predict conditional on a subject having already survived some time, which is useful for updating a prognosis partway through follow-up. With conditional_after=180, the returned probabilities are \(P(T > t \mid T > 180) = S(t) / S(180)\).

cox.predict(newdata, type="survival", times=[365, 730], conditional_after=180, format="polars")

	time f64	subject_1 f64	subject_2 f64
PolarsRows2Columns3
0	365	0.572489013425	0.625329029654
1	730	0.163955994862	0.218284510803

Stratification

Sometimes a variable violates the proportional hazards assumption but is not of direct interest, such as enrolling hospital. Stratification lets each level of that variable have its own baseline hazard while sharing the coefficients of the covariates you do care about. You pass the stratifying variable to strata. Here we stratify by sex while keeping age and ph.ecog as covariates.

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

stratified = gw.CoxPH().fit(y, lung[["age", "ph.ecog"]], strata=lung["sex"])

As with any fitted model, we can print stratified for a summary, or read its coefficient table as data with to_frame().

stratified.to_frame(format="polars")

	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64
PolarsRows2Columns7
0	age	0.0105662546009	0.00924137389309	1.14336404123	0.252887475744	-0.00754650539718	0.0286790145991
1	ph.ecog	0.462424434358	0.114761097855	4.0294528634	5.59068243623e-05	0.237496815737	0.68735205298

Notice that the stratifying variable, sex, does not appear as a coefficient. Its effect is absorbed into the separate baselines, which is exactly the point of stratification.

Robust and clustered standard errors

When observations are not independent, for example when the same subject contributes several rows or when subjects are grouped into clusters, the model-based standard errors are too small. The robust, or sandwich, variance corrects for this. Set robust=True for the Lin-Wei sandwich estimator, and pass cluster to group correlated observations.

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

robust = gw.CoxPH().fit(y, lung[["age", "sex"]], robust=True)

robust.to_frame(format="polars")[["term", "std_error"]]

	term str	std_error f64
PolarsRows2Columns2
0	age	0.00948922108421
1	sex	0.159919372277

The coefficients are unchanged; only the standard errors, and therefore the p-values and intervals, are adjusted. The model-based standard errors remain available as naive_std_error_ for comparison.

The concordance index

The concordance index, or C-statistic, measures how well the model’s risk ordering agrees with the observed order of events. It is the probability that, for a random comparable pair of subjects, the one who failed first had the higher predicted risk. A value of 0.5 is no better than chance and 1.0 is perfect discrimination.

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])

cox.concordance()

0.6028530028979714

Next steps

You can now validate a Cox model and use it to predict.

Prediction performance extends discrimination and adds the Brier score for calibration, including for external risk scores.
Parametric survival models provide an alternative when the proportional hazards assumption is untenable.
Revisit Cox regression for the modeling basics if any of this is unfamiliar.