import greenwood as gw
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
ySurv(type=right, n=228, events=165)
A model that fits the data in front of it is not necessarily a model that predicts well. To judge a survival model as a predictor, or to compare competing models, you need metrics that account for censoring. Greenwood provides three complementary ones: (1) the concordance index, (measures how well a model ranks subjects by risk), (2) the Brier score (measures how accurate its probability predictions are), and (3) the time-dependent AUC (shows how discrimination changes over the follow-up horizon). This guide page explains all three and shows how to compute them for a fitted model or for any external risk score.
We start from the lung outcomes and a Cox model fit to them. This is the model whose predictions we will judge.
import greenwood as gw
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
ySurv(type=right, n=228, events=165)
The response y records 165 events among 228 subjects. We fit a Cox model to it, then set the model aside: everything below needs only y and the risk scores and survival probabilities that cox produces.
Discrimination is the ability to rank subjects correctly: to give higher risk scores to subjects who go on to have the event sooner. The concordance index, or C-statistic, measures exactly this. It is the probability that, for a randomly chosen comparable pair of subjects, the one who failed first had the higher risk score. A value of 0.5 is no better than a coin flip, and 1.0 is perfect ordering.
You can compute it for any risk score, which makes it useful for evaluating a model on new data or for comparing an external score. Here we use the Cox model’s linear predictor. Asking for type="lp" returns one number per subject, where a larger value means higher predicted risk.
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"]])
risk = cox.predict(type="lp")
risk[:8]array([ 0.3995047 , 0.2972327 , 0.09268872, 0.10973405, 0.16087005,
0.3995047 , -0.21598581, -0.16484982])
The linear predictor returned by predict(type="lp") is always centered around the mean of the training covariates. This ensures consistent, numerically stable predictions across all variable types: continuous, categorical, and mixed. The centering is applied uniformly, so you always get mathematically sound risk scores, even in edge cases like categorical-only models with all positive associations.
These are the scores the index will rank against the observed event order. Because higher risk should correspond to earlier failure, we pass a score where larger means higher risk, such as this linear predictor or its exponential.
A model can rank subjects well yet produce poorly calibrated probabilities, or vice versa. The concordance index says nothing about whether predicted probabilities are numerically accurate. That is what the Brier score adds.
The Brier score is the mean squared error between predicted survival probabilities and actual outcomes, evaluated at a fixed time. Lower is better. Because outcomes are censored, the score uses inverse-probability-of-censoring weighting so that censored subjects contribute correctly rather than being ignored or mishandled.
To compute it you supply predicted survival probabilities for each subject at each evaluation time. Start by asking the Cox model for survival curves at three fixed times.
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"]])
times = [180, 365, 540]
surv = cox.predict(lung[["age", "sex"]], type="survival", times=times, format="polars")
surv.shape(3, 229)
The returned DataFrame has a leading time column and one column per subject, so its shape is 3 rows (the requested times) by 229 columns (228 subjects plus the time column). The Brier score wants the transpose of this: an array with one row per subject and one column per time. We drop the time column and transpose to get there.
With the probabilities in the expected (n_subjects, n_times) shape, we can score them.
The result is one Brier score per time. As a rule of thumb, a useful model should beat the Brier score of a model that ignores covariates, and a score approaching 0.25 at a time where about half the subjects have had the event is close to uninformative.
To summarize accuracy across a range of times in a single number, integrate the Brier score over time.
Computing these metrics on the same data used to fit the model gives optimistic results. For an honest assessment, evaluate on a held-out test set or use cross-validation, fitting on one part of the data and scoring predictions on another.
The concordance index summarises discrimination with a single number across all event times. But a marker’s ability to separate early failures from survivors often changes over follow-up: a biomarker may identify six-month events well but lose its signal at five years. time_dependent_auc() gives you a per-horizon view by computing the area under the ROC curve separately at each requested time.
At each time t, the estimator compares cases (subjects who had an event by t) against controls (subjects still at risk after t), and measures how well the marker separates the two groups. Inverse-probability-of-censoring weights correct for the fact that not all subjects are observed to the end of follow-up. The marker must be on a scale where larger means higher risk, so, the Cox linear predictor is a natural choice.
array([0.64152627, 0.5943264 , 0.6023971 ])
Each value is the AUC at one horizon. Values above 0.5 mean the marker separates early failures from survivors better than chance at that time. Values close to 0.5 signal that the marker has lost predictive power at that horizon.
To reduce these to a single number for model comparison, integrated_auc() computes the time-averaged AUC using the trapezoidal rule over the supplied time range.
The integrated AUC plays the same role as the integrated Brier score: it collapses time-specific discrimination into one comparison-friendly number. Higher is better, and 0.5 is the random baseline.
Rather than splitting the data by hand, cross_validate() runs k-fold cross-validation for you: it fits a fresh copy of the model on each set of training folds and scores its predictions on the held-out fold. Pass an unfitted model, the response, and the covariates.
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
gw.cross_validate(gw.CoxPH(), y, lung[["age", "sex"]], k=5, metric="concordance", seed=23){'metric': 'concordance',
'k': 5,
'scores': [0.6617466174661747,
0.5448877805486284,
0.6838235294117647,
0.5580524344569289,
0.5926694329183956],
'mean': 0.6082359589603785,
'std': 0.06194870045682596}
The result reports the score on each fold plus their mean and standard deviation. The cross-validated concordance is usually a little lower than the value computed on the training data, which is exactly the optimism we are trying to remove. Set metric="brier" with a list of times to cross-validate the integrated Brier score instead, and the same call works for an AFT model.
gw.cross_validate(
gw.CoxPH(), y, lung[["age", "sex"]], k=5,
metric="brier", times=[180, 365, 540], seed=23
){'metric': 'brier',
'k': 5,
'scores': [0.21443779185001324,
0.2004839745201561,
0.21327986359004628,
0.23470395029714647,
0.22403878785186795],
'mean': 0.21738887362184597,
'std': 0.012800407222447782}
By default, cross_validate() uses stratified k-fold splitting. Rather than randomly assigning subjects to folds, it groups subjects by event status first, then draws from each group proportionally. This ensures that each fold contains roughly the same fraction of events as the overall dataset.
This matters whenever events are rare. In a dataset where only 10% of subjects have the event, pure random splitting could easily put all events in one fold and none in another, leaving training sets without enough events to fit the model at all. Stratification prevents that. You can opt out with stratified=False, but there is rarely a good reason to.
cross_validate() warns when the total number of events is less than twice the number of folds. At that point, some test folds will contain only one event, making per-fold scores unreliable:
UserWarning: Only 6 events found for 5-fold cross-validation
(fewer than 2 × k = 10). Some folds may contain too few events
for reliable evaluation. Consider reducing k or collecting more
data with observed events.
If you see this warning, either reduce k= to match the event count (e.g., k=3 with 6 events) or treat the cross-validated scores with caution.
Cross-validation is most useful as a comparison tool. Fit the same folds on two competing models, keeping the same seed, and compare their mean scores directly. The model with higher mean concordance (or lower mean Brier score) generalizes better to unseen subjects.
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cols_simple = lung[["age", "sex"]]
cols_full = lung[["age", "sex", "ph.ecog"]]
simple_cv = gw.cross_validate(gw.CoxPH(), y, cols_simple, seed=1)
full_cv = gw.cross_validate(gw.CoxPH(), y, cols_full, seed=1)
print(f"Simple model C-index: {simple_cv['mean']:.3f} ± {simple_cv['std']:.3f}")
print(f"Full model C-index: {full_cv['mean']:.3f} ± {full_cv['std']:.3f}")Simple model C-index: 0.604 ± 0.077
Full model C-index: 0.639 ± 0.043
The standard deviation across folds gives a rough sense of how stable the estimate is. Large fold-to-fold variation (relative to the difference between models) means the data do not contain enough information to reliably distinguish the two models on this metric alone.
Discrimination and the Brier score both summarize a model in a single number. Calibration asks a more detailed question: at a chosen time, do the predicted survival probabilities match what actually happens? calibration() groups subjects into bins by their predicted survival at a horizon, then compares the mean prediction in each bin against the observed survival, a Kaplan-Meier estimate for that bin.
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"]])
surv_at_year = cox.predict(
lung[["age", "sex"]], type="survival", times=[365.0], format="polars"
)
predicted = surv_at_year.slice(0, 1).drop("time").to_numpy().flatten()
gw.calibration(y, predicted, 365.0, n_bins=5, format="polars")PolarsRows5Columns6 | ||||||
bin i64 |
n i64 |
predicted f64 |
observed f64 |
observed_lower f64 |
observed_upper f64 |
|
|---|---|---|---|---|---|---|
| 0 | 1 | 45 | 0.275643654127 | 0.226515151515 | 0.125630250282 | 0.40841368819 |
| 1 | 2 | 44 | 0.325992610784 | 0.363081617086 | 0.237804025086 | 0.554356725535 |
| 2 | 3 | 45 | 0.388623683685 | 0.408115942029 | 0.284454511367 | 0.585536932909 |
| 3 | 4 | 45 | 0.485377847087 | 0.54000350566 | 0.395496812687 | 0.737310078795 |
| 4 | 5 | 49 | 0.568638173192 | 0.514285348805 | 0.376278476604 | 0.702908713734 |
Each row is a bin, from the lowest predicted survival to the highest. A well-calibrated model has observed close to predicted in every bin, with the observed values inside their confidence limits. Systematic gaps, for example predictions consistently higher than observed, signal miscalibration even when discrimination is good.
Like the other metrics, calibration is optimistic on the training data. For a fair picture, compute it on held-out predictions, for instance the test-fold predictions from a cross-validation split.
You can now quantify how well a survival model discriminates and predicts.