# 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](../reference/Surv.md#greenwood.Surv) object pairing each follow-up time with an event indicator, and it is the target every model on this page is fit against.
``` python
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.
``` python
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](../reference/CoxPH.md#greenwood.CoxPH.cox_zph) returns a [ZPHResult](../reference/ZPHResult.md#greenwood.ZPHResult), an object that bundles the test statistics together. Displaying it directly gives a compact summary of the test.
``` python
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()](../reference/AFT.md#greenwood.AFT.to_frame) on that result. This is the table you will read most often.
``` python
zph.to_frame(format="polars")
```
PolarsRows3Columns4
|
|
term
str
|
chisq
f64
|
df
i64
|
p_value
f64
|
| 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.
> **Warning: 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.
``` python
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.
``` python
cox.residuals("schoenfeld", format="polars")
```
PolarsRows165Columns2
|
|
age
f64
|
sex
f64
|
| 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.
``` python
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")
```
PolarsRows186Columns3
|
|
time
f64
|
cumhaz
f64
|
survival
f64
|
| 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.
``` python
import pandas as pd
newdata = pd.DataFrame({"age": [50, 70], "sex": [1, 2]})
newdata
```
PandasRows2Columns2
|
|
age
i64
|
sex
i64
|
| 0 |
50 |
1 |
| 1 |
70 |
2 |
Given that new data and a set of times, [predict](../reference/AFT.md#greenwood.AFT.predict) returns the estimated survival probability for each subject.
``` python
cox.predict(newdata, type="survival", times=[180, 365, 730], format="polars")
```
PolarsRows3Columns3
|
|
time
f64
|
subject_1
f64
|
subject_2
f64
|
| 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).
``` python
cox.predict(newdata, type="survival", times=[180, 365, 730], ci=True, format="polars")
```
PolarsRows3Columns7
|
|
time
f64
|
subject_1
f64
|
subject_1_lower
f64
|
subject_1_upper
f64
|
subject_2
f64
|
subject_2_lower
f64
|
subject_2_upper
f64
|
| 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)\\.
``` python
cox.predict(newdata, type="survival", times=[365, 730], conditional_after=180, format="polars")
```
PolarsRows2Columns3
|
|
time
f64
|
subject_1
f64
|
subject_2
f64
|
| 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.
``` python
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](../reference/RMSTResult.md#greenwood.RMSTResult.stratified) for a summary, or read its coefficient table as data with [to_frame()](../reference/AFT.md#greenwood.AFT.to_frame).
``` python
stratified.to_frame(format="polars")
```
PolarsRows2Columns7
|
|
term
str
|
estimate
f64
|
std_error
f64
|
statistic
f64
|
p_value
f64
|
conf_low
f64
|
conf_high
f64
|
| 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.
``` python
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"]]
```
PolarsRows2Columns2
|
|
term
str
|
std_error
f64
|
| 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.
``` python
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](prediction-performance.md) extends discrimination and adds the Brier score for calibration, including for external risk scores.
- [Parametric survival models](parametric-models.md) provide an alternative when the proportional hazards assumption is untenable.
- Revisit [Cox regression](cox-regression.md) for the modeling basics if any of this is unfamiliar.