Prediction performance

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 two complementary ones: the concordance index, which measures how well a model ranks subjects by risk, and the Brier score, which measures how accurate its probability predictions are. This page explains both 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))
y

Surv(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.

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

Discrimination: the concordance index

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

Linear predictors are centered

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.

gw.concordance_index(y, risk)

0.6028530028979714

Discrimination is only part of the story

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.

Accuracy: the Brier score

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

probs = surv.drop("time").to_numpy().T
probs.shape

(228, 3)

With the probabilities in the expected (n_subjects, n_times) shape, we can score them.

gw.brier_score(y, probs, times)

array([0.19160509, 0.23614214, 0.18644928])

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.

gw.integrated_brier_score(y, probs, times)

0.212620468119578

Evaluate on data the model did not see

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.

Cross-validation

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=1)

{'metric': 'concordance',
 'k': 5,
 'scores': [0.5720524017467249,
  0.6159147869674185,
  0.44519704433497537,
  0.6560468140442133,
  0.687603305785124],
 'mean': 0.5953628705756911,
 'std': 0.09448065076706443}

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.

Calibration

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

	bin i64	n i64	predicted f64	observed f64	observed_lower f64	observed_upper f64
PolarsRows5Columns6
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.

Assess calibration out of sample too

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.

Next steps

You can now quantify how well a survival model discriminates and predicts.

Cox model diagnostics covers the assumption checks and residuals that complement these performance metrics.
Parametric survival models and Cox regression produce the predictions these metrics evaluate.
The Quick start ties the whole workflow together.