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 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))
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])
TipLinear 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
NoteDiscrimination 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.

Censoring-corrected discrimination: the IPCW concordance index

Harrell’s C treats censored subjects by simply excluding incomparable pairs. When censoring is light and roughly uniform, this works well. But when censoring is heavy (more than 30-40% of subjects) or differs across risk groups, the resulting estimate can be biased because the excluded pairs are not missing at random.

The IPCW concordance index of Uno et al. (2011) corrects for this by weighting each concordant pair by the inverse of the censoring probability at the event time. Subjects whose events occur at times when many others have already been censored receive higher weight, compensating for the information lost to censoring. The result is a consistent estimator of the true concordance probability.

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")
gw.concordance_index_ipcw(y, risk)
np.float64(0.595890711822792)

In the lung dataset censoring is moderate (~28%), so the IPCW and Harrell values are close. In datasets with heavier censoring the two can diverge meaningfully. The IPCW version is the safer default for formal model comparisons.

The tau parameter restricts the evaluation to events before a truncation time. This avoids instability from very small censoring probabilities in the tail, where few subjects remain under observation. It also focuses discrimination on a clinically relevant horizon.

gw.concordance_index_ipcw(y, risk, tau=365.0)
np.float64(0.5990691290907791)
TipChoosing tau

When the censoring survival function drops close to zero in the tail, the IPCW weights become very large and the estimate becomes unstable. Setting tau to the end of the clinically relevant follow-up window (for example, one year or five years) avoids this. The default is the largest observed event time, which works well when censoring remains moderate throughout follow-up.

The censoring distribution

The IPCW concordance index above works because it knows how likely each subject was to be censored. That knowledge comes from a CensoringDistribution, which estimates the Kaplan-Meier curve of the censoring process. Greenwood uses it internally for all IPCW metrics, but you can also create one directly when you need the censoring probabilities or IPC weights for a custom analysis.

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

cens = gw.CensoringDistribution(y)
cens
CensoringDistribution (n=228, censored=63)
  60 unique censoring times, tail G = 0.0000

The survival() method evaluates the censoring survival function G(t), which is the probability of not being censored by time t, at any set of query times.

cens.survival([180, 365, 540])
array([0.96644606, 0.69728182, 0.56736787])

The values decrease over follow-up as more subjects are censored. At 180 days about 97% of subjects are still under observation, but by 540 days only 57% remain.

The weights() method returns the IPC weights that correct for censoring. Each uncensored subject receives a weight of 1 / G(t-), where G(t-) is the censoring survival just before their event time. Censored subjects receive weight zero (they are excluded from the reweighted analysis). Subjects whose events occur when censoring is heavy get larger weights because they represent more unobserved individuals.

w = cens.weights()
w[:10]
array([1.35885572, 1.59245456, 0.        , 1.09235351, 3.8210638 ,
       0.        , 1.35885572, 1.41273529, 1.10044502, 1.01015385])

You can restrict weights to a time horizon with tau. Subjects with events after tau receive zero weight, which avoids instability from very small censoring probabilities in the tail.

cens.weights(tau=365.0)[:10]
array([1.35885572, 0.        , 0.        , 1.09235351, 0.        ,
       0.        , 1.35885572, 1.41273529, 1.10044502, 1.01015385])

The full censoring curve is available as a DataFrame via to_frame().

cens.to_frame()
PolarsRows60Columns2
time
f64
survival
f64
0 92 0.995024875622
1 105 0.989948218093
2 173 0.984090536329
3 174 0.978232854565
4 175 0.972375172802
5 177 0.966446055894
6 185 0.96021092005
7 188 0.953935031683
57 965 0.174471482741
58 1010 0.0872357413706
59 1022 0

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

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.21262046811957802
TipEvaluate 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.

Time-dependent discrimination: AUC over time

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.

lp = cox.predict(type="lp")
gw.time_dependent_auc(y, lp, times=[180, 365, 540])
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.

gw.integrated_auc(y, lp, times=[180, 365, 540])
0.6084157733720682

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.

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=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. The same call works for AFT and RoystonParmar models.

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.21443779185001313,
  0.20048397452015598,
  0.21327986359004625,
  0.23470395029714616,
  0.22403878785186804],
 'mean': 0.21738887362184595,
 'std': 0.012800407222447734}

Multi-metric mode

When you want concordance, Brier, and AUC together, pass metrics= (plural) instead of metric=. The model is fit once per fold and scored on every requested metric, so you get consistent fold assignments and avoid redundant refitting.

result_multi = gw.cross_validate(
    gw.CoxPH(), y, lung[["age", "sex"]], k=5,
    metrics=["concordance", "brier", "auc"],
    times=[180, 365, 540], seed=23
)
result_multi
{'k': 5,
 'results': {'concordance': {'scores': [0.6617466174661747,
    0.5448877805486284,
    0.6838235294117647,
    0.5580524344569289,
    0.5926694329183956],
   'mean': 0.6082359589603785,
   'std': 0.06194870045682596},
  'brier': {'scores': [0.21443779185001313,
    0.20048397452015598,
    0.21327986359004625,
    0.23470395029714616,
    0.22403878785186804],
   'mean': 0.21738887362184595,
   'std': 0.012800407222447734},
  'auc': {'scores': [0.7003890238319559,
    0.5241297774787457,
    0.7289478073383704,
    0.5316521748637566,
    0.47168767470244194],
   'mean': 0.5913612916430542,
   'std': 0.11535177220406084}}}

The return value is a keyed dictionary. Each metric has its own scores, mean, and std:

result_multi["results"]["concordance"]["mean"]
0.6082359589603785
result_multi["results"]["brier"]["mean"]
0.21738887362184595

Stratified folds for imbalanced data

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.

WarningWatch the event count

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.

Comparing models

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.

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

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