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
# Load the lung dataset
lung = gw.load_dataset("lung", backend="polars")
# Name the survival response columns, then bind them to the data to display a summary
death = gw.Outcome.surv(time="time", event="status")
death.bind(lung)Surv(type=right, n=228, events=165)The response death names the follow-up time and event columns, and bound to lung it records 165 events among 228 subjects. We fit a Cox model with the same response, then set the model aside: everything below needs only death, the data, and the risk scores and survival probabilities that cox produces. Each metric reads the response columns from the frame passed as data=.
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.
# Rebuild the dataset, response, and Cox fit for this section
lung = gw.load_dataset("lung", backend="polars")
death = gw.Outcome.surv(time="time", event="status")
cox = gw.CoxPH().fit("Surv(time, status) ~ age + sex", data=lung)
# Extract the linear predictor (risk score) for each subject
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.
# Compute Harrell's concordance index from the risk scores
gw.concordance_index(death, risk=risk, data=lung)0.6028530028979714A 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.
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.
# Rebuild the dataset, response, and Cox fit
lung = gw.load_dataset("lung", backend="polars")
death = gw.Outcome.surv(time="time", event="status")
cox = gw.CoxPH().fit("Surv(time, status) ~ age + sex", data=lung)
# IPCW concordance index corrected for censoring bias
risk = cox.predict(type="lp")
gw.concordance_index_ipcw(death, risk=risk, data=lung)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.
# Restrict evaluation to events before one year
gw.concordance_index_ipcw(death, risk=risk, data=lung, tau=365.0)np.float64(0.5990691290907791)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.
A point estimate of the concordance is not enough for a formal comparison. You also need a standard error and a confidence interval. concordance_index_ci() computes the IPCW concordance and derives the standard error from the influence function decomposition, which avoids the cost of bootstrapping. The confidence interval is constructed on the logit scale (the default), which respects the [0, 1] bounds and provides better coverage when the estimate is near either extreme.
# Rebuild the dataset, response, and Cox fit
lung = gw.load_dataset("lung", backend="polars")
death = gw.Outcome.surv(time="time", event="status")
cox = gw.CoxPH().fit("Surv(time, status) ~ age + sex", data=lung)
# Concordance with standard error and confidence interval
risk = cox.predict(type="lp")
gw.concordance_index_ci(death, risk=risk, data=lung)ConcordanceResult(estimate=0.5959, se=0.0391, CI=[0.5176, 0.6696], conf_level=0.95)The result is a ConcordanceResult with the point estimate, standard error, and confidence limits. You can convert it to a one-row DataFrame with to_frame().
# Convert the result to a one-row DataFrame
gw.concordance_index_ci(death, risk=risk, data=lung).to_frame()PolarsRows1Columns5 | |||||
estimate f64 |
se f64 |
conf_low f64 |
conf_high f64 |
n i64 |
|
|---|---|---|---|---|---|
| 0 | 0.595890711823 | 0.0390599327328 | 0.517606378803 | 0.669580691608 | 228 |
The transform parameter controls the scale on which the CI is constructed. The default "logit" is recommended because it keeps the interval inside [0, 1] even when the estimate is close to the boundary. The "log" and "identity" options are available for comparison.
# Same concordance with an identity-scale confidence interval
gw.concordance_index_ci(death, risk=risk, data=lung, transform="identity")ConcordanceResult(estimate=0.5959, se=0.0391, CI=[0.5193, 0.6724], conf_level=0.95)For competing-risks data, pass cause= and tau= to get the cause-specific concordance with its CI. This uses the incidence-weighted influence function, which accounts for both standard and competing-event pairs.
import numpy as np
# Load competing-risks data
mgus2 = gw.load_dataset("mgus2", backend="polars")
# The first observed endpoint (PCM at `ptime`, or death at `futime`) is the event
cr = gw.Outcome.first_event(
endpoints={"pcm": ("ptime", "pstat"), "death": ("futime", "death")},
)
# Fit a Fine-Gray model for PCM and predict CIF at 240 months
fg = gw.FineGray(cause="pcm").fit(cr, covariates=["age", "sex"], data=mgus2)
cif_at_tau = fg.predict_cumulative_incidence(
mgus2, times=np.array([240.0]), format="pandas"
).drop(columns="time").values.flatten()
# Cause-specific concordance with confidence interval
gw.concordance_index_ci(cr, risk=cif_at_tau, data=mgus2, cause=1, tau=240.0)ConcordanceResult(estimate=0.5780, se=0.0278, CI=[0.5228, 0.6314], conf_level=0.95)When two models are fit on the same subjects, their concordance estimates are correlated. A naive comparison that treats them as independent underestimates uncertainty and inflates the Type I error rate. The concordance_index_compare() function uses the paired influence-function delta method: it computes the difference in concordance indices and derives the standard error from the joint influence function, which properly accounts for the shared subjects.
# Rebuild the dataset and response
lung = gw.load_dataset("lung", backend="polars")
death = gw.Outcome.surv(time="time", event="status")
# Fit two competing Cox models
cox_full = gw.CoxPH().fit("Surv(time, status) ~ age + sex", data=lung)
cox_age = gw.CoxPH().fit("Surv(time, status) ~ age", data=lung)
# Extract linear predictors from each model
lp_full = cox_full.predict(type="lp")
lp_age = cox_age.predict(type="lp")
# Paired comparison of concordance indices
gw.concordance_index_compare(death, risk_a=lp_full, risk_b=lp_age, data=lung)ConcordanceCompareResult(delta=0.0467, se=0.0355, CI=[-0.0228, 0.1162], z=1.318, p=0.1875)The result reports the difference (model A minus model B), its standard error, a confidence interval, a z-statistic, and a two-sided p-value. A small p-value means the two models have detectably different discrimination on this data. The result also has a to_frame() method for tabular output.
# Tabular output of the paired comparison
gw.concordance_index_compare(death, risk_a=lp_full, risk_b=lp_age, data=lung).to_frame()PolarsRows1Columns6 | ||||||
delta f64 |
se f64 |
conf_low f64 |
conf_high f64 |
z f64 |
pvalue f64 |
|
|---|---|---|---|---|---|---|
| 0 | 0.0467304912597 | 0.0354554283048 | -0.0227608712742 | 0.116221853794 | 1.31800667751 | 0.187501412585 |
concordance_index_compare() gives you a formal hypothesis test on a single dataset. Cross- validation with cross_validate() gives you an estimate of out-of-sample performance. They answer different questions. Use the paired test when you want a p-value for the difference on the observed data. Use cross-validation when you want to estimate how much better one model will perform on new subjects.
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.
# Rebuild the dataset and response
lung = gw.load_dataset("lung", backend="polars")
death = gw.Outcome.surv(time="time", event="status")
# Estimate the censoring distribution from the response (it takes a bound `Surv`)
cens = gw.CensoringDistribution(death.bind(lung))
censCensoringDistribution (n=228, censored=63)
60 unique censoring times, tail G = 0.0000The 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.
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.
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.
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().
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 |
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.
# Rebuild the dataset, response, and Cox fit
lung = gw.load_dataset("lung", backend="polars")
death = gw.Outcome.surv(time="time", event="status")
cox = gw.CoxPH().fit("Surv(time, status) ~ age + sex", data=lung)
# Predict survival probabilities at three fixed times
times = [180, 365, 540]
surv = cox.predict(lung, 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.
# Transpose to (n_subjects, n_times) for the Brier score
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.
# IPCW Brier score at each evaluation time
gw.brier_score(death, survival_prob=probs, times=times, data=lung)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.
# Integrated Brier score across the three time horizons
gw.integrated_brier_score(death, survival_prob=probs, times=times, data=lung)0.21262046811957802Computing 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 standard Brier score evaluates predicted survival probabilities. In a competing-risks setting, the quantity of interest is instead the cumulative incidence function (CIF) for a specific cause. brier_score_incidence() extends the IPCW Brier score to this setting: it compares predicted CIF probabilities against observed cause-specific outcomes.
The weighting scheme accounts for competing events. Subjects who experienced any event (including a competing cause) before the evaluation time are weighted by the inverse censoring probability at their event time. Subjects still at risk contribute at the inverse censoring probability at the evaluation time. Subjects censored before the evaluation time receive zero weight.
To illustrate, we describe a competing-risks response for the mgus2 dataset, where subjects may progress to plasma cell malignancy (PCM) or die without progression. Outcome.first_event() names the time and status columns of each endpoint.
import numpy as np
# Load competing-risks data
mgus2 = gw.load_dataset("mgus2", backend="polars")
# The first observed endpoint (PCM at `ptime`, or death at `futime`) is the event
cr = gw.Outcome.first_event(
endpoints={"pcm": ("ptime", "pstat"), "death": ("futime", "death")},
)
# Bind the response to the data to display its summary
cr.bind(mgus2)Surv(type=mright, n=1384, events=975, states=('pcm', 'death'))As a baseline, we use the marginal Aalen-Johansen CIF as a naive model that predicts the same CIF for every subject.
# Fit Aalen-Johansen and extract the marginal PCM incidence curve
aj = gw.AalenJohansen().fit(cr, data=mgus2)
cif_df = aj.to_frame(format="polars")
pcm_cif = cif_df.filter(cif_df["cause"] == "pcm")
# Interpolate the marginal CIF at the evaluation times
times_cr = np.array([60, 120, 240])
marginal = np.array(
[float(np.interp(t, pcm_cif["time"].to_numpy(), pcm_cif["estimate"].to_numpy()))
for t in times_cr]
)
# Tile the marginal CIF across all subjects and compute the Brier score
probs_cr = np.tile(marginal, (mgus2.height, 1))
gw.brier_score_incidence(cr, incidence_prob=probs_cr, times=times_cr, data=mgus2, cause=1)array([0.03293805, 0.05966165, 0.0898196 ])The cause= argument selects which event code to evaluate. Here cause=1 corresponds to progression (the first endpoint listed in endpoints=). To evaluate death predictions instead, pass cause=2.
As with the survival Brier score, you can summarize accuracy over a range of times in a single number with the integrated version.
# Evaluate the integrated Brier score on a finer grid
times_fine = np.linspace(30, 360, 12)
marginal_fine = np.array(
[float(np.interp(t, pcm_cif["time"].to_numpy(), pcm_cif["estimate"].to_numpy()))
for t in times_fine]
)
probs_fine = np.tile(marginal_fine, (mgus2.height, 1))
gw.integrated_brier_score_incidence(
cr, incidence_prob=probs_fine, times=times_fine, data=mgus2, cause=1
)0.07514656777562731A model that uses covariates to produce subject-specific CIF predictions should yield a lower integrated Brier score than this marginal baseline.
The standard concordance index evaluates discrimination for a single event type. In a competing-risks setting, concordance_index_incidence() extends the IPCW concordance to cause-specific cumulative incidence predictions (Wolbers et al., 2014). It compares two types of pairs: standard pairs (a subject who experienced the cause of interest vs. a subject still at risk) and competing-event pairs (a subject who experienced the cause vs. one who experienced a competing event). Each pair type receives its own IPCW weight.
Using the same mgus2 competing-risks data from above, we evaluate the concordance of the marginal Aalen-Johansen CIF at 240 months. Because the marginal model predicts the same value for every subject, its concordance is 0.5 (no discrimination).
# Marginal CIF at tau=240 months: same prediction for everyone
tau = 240.0
marginal_at_tau = float(
np.interp(tau, pcm_cif["time"].to_numpy(), pcm_cif["estimate"].to_numpy())
)
pred_flat = np.full(mgus2.height, marginal_at_tau)
# Concordance for the flat prediction (expected to be 0.5)
gw.concordance_index_incidence(cr, incidence_prob=pred_flat, data=mgus2, cause=1, tau=tau)np.float64(0.49999999999999994)A model that produces subject-specific CIF predictions should achieve a concordance above 0.5.
Beyond discrimination, a model’s predicted CIF values should match the observed incidence on average. calibration_incidence() checks marginal calibration by comparing the mean predicted CIF at each time against the nonparametric Aalen-Johansen reference. The result is the absolute calibration error at each time point.
# Calibration check: mean predicted CIF vs. AJ reference
cal_err = gw.calibration_incidence(
cr, incidence_prob=probs_cr, times=times_cr, data=mgus2, cause=1
)
cal_errarray([6.93889390e-18, 1.38777878e-17, 1.38777878e-17])Because the marginal model is the AJ reference, its calibration error is essentially zero. A subject-specific model will typically show nonzero errors, especially in the tails.
When the model predicts CIFs for all competing causes, accuracy_in_time() evaluates whether the model correctly identifies the most likely outcome at each time horizon. At time \zeta, each subject’s predicted outcome is the argmax of the predicted CIFs plus the implied survival probability. The metric reports the proportion of uncensored subjects whose predicted class matches the observed outcome.
# Interpolate both CIFs at the evaluation times
death_cif = cif_df.filter(cif_df["cause"] == "death")
pcm_vals = np.array(
[float(np.interp(t, pcm_cif["time"].to_numpy(), pcm_cif["estimate"].to_numpy()))
for t in times_cr]
)
death_vals = np.array(
[float(np.interp(t, death_cif["time"].to_numpy(), death_cif["estimate"].to_numpy()))
for t in times_cr]
)
# Stack into (n_subjects, n_causes, n_times) and compute accuracy
probs_3d = np.stack(
[np.tile(pcm_vals, (mgus2.height, 1)), np.tile(death_vals, (mgus2.height, 1))], axis=1
)
gw.accuracy_in_time(cr, incidence_probs=probs_3d, times=times_cr, data=mgus2)array([0.63884786, 0.58021613, 0.83711747])Early times tend to show high accuracy because most subjects have survived and the model predicts survival. The metric becomes most informative at later times, when discriminating between competing events matters.
When evaluating a competing-risks model you typically want concordance, per-time Brier scores, and the integrated Brier score all at once. score_cr() computes them in a single call and returns a tidy summary table. The concordance row includes a confidence interval from the influence function. The Brier rows report point estimates only.
import numpy as np
# Load data and build competing-risks response
mgus2 = gw.load_dataset("mgus2", backend="polars")
# The first observed endpoint (PCM at `ptime`, or death at `futime`) is the event
cr = gw.Outcome.first_event(
endpoints={"pcm": ("ptime", "pstat"), "death": ("futime", "death")},
)
# Fit Fine-Gray for PCM and predict CIF at three horizons
fg = gw.FineGray(cause="pcm").fit(cr, covariates=["age", "sex"], data=mgus2)
eval_times = np.array([120, 240, 360])
cif_pred = fg.predict_cumulative_incidence(
mgus2, times=eval_times, format="pandas"
).drop(columns="time").values.T
# Bundle concordance, Brier, and integrated Brier into one summary
gw.score_cr(cr, incidence_prob=cif_pred, data=mgus2, cause=1, times=eval_times)PolarsRows5Columns5 | |||||
metric str |
estimate f64 |
se f64 |
conf_low f64 |
conf_high f64 |
|
|---|---|---|---|---|---|
| 0 | concordance | 0.664422241225 | 0.0445690503462 | 0.572310071046 | 0.74551736935 |
| 1 | brier(120) | 0.0877376948292 | NaN | NaN | NaN |
| 2 | brier(240) | 0.147587386975 | NaN | NaN | NaN |
| 3 | brier(360) | 0.199861349421 | NaN | NaN | NaN |
| 4 | integrated_brier | 0.14569345455 | NaN | NaN | NaN |
The incidence_prob argument expects a matrix with shape (n_subjects, n_times), where each column holds the predicted CIF at the corresponding evaluation time. The cause argument selects the event code and the tau parameter (defaulting to the last evaluation time) truncates the concordance calculation.
This is the competing-risks counterpart of running concordance_index_ci(), brier_score_incidence(), and integrated_brier_score_incidence() separately. When you only need one of those metrics, calling them individually is fine. The score_cr() function saves repetitive setup when you want the full picture.
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.
# Time-dependent AUC at three horizons using the Cox linear predictor
lp = cox.predict(type="lp")
gw.time_dependent_auc(death, marker=lp, times=[180, 365, 540], data=lung)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.
# Time-averaged AUC across the three horizons
gw.integrated_auc(death, marker=lp, times=[180, 365, 540], data=lung)0.6084157733720682The 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 and a formula naming the response and the covariates.
# Rebuild the dataset
lung = gw.load_dataset("lung", backend="polars")
# Five-fold cross-validated concordance
gw.cross_validate(
gw.CoxPH(), "Surv(time, status) ~ age + sex", data=lung, 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.
# Cross-validated integrated Brier score at three time points
gw.cross_validate(
gw.CoxPH(), "Surv(time, status) ~ age + sex", data=lung, k=5,
metric="brier", times=[180, 365, 540], seed=23
){'metric': 'brier',
'k': 5,
'scores': [0.21443779185001316,
0.200483974520156,
0.2132798635900462,
0.23470395029714639,
0.2240387878518678],
'mean': 0.2173888736218459,
'std': 0.01280040722244777}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. Repeated metric names are evaluated once, preserving the order of their first appearance; each metric returns exactly one score per fold.
# Multi-metric cross-validation: concordance, Brier, and AUC in one pass
result_multi = gw.cross_validate(
gw.CoxPH(), "Surv(time, status) ~ age + sex", data=lung, 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.21443779185001316,
0.200483974520156,
0.2132798635900462,
0.23470395029714639,
0.2240387878518678],
'mean': 0.2173888736218459,
'std': 0.01280040722244777},
'auc': {'scores': [0.7003890238319559,
0.5241297774787457,
0.7289478073383705,
0.5316521748637567,
0.47168767470244194],
'mean': 0.5913612916430542,
'std': 0.11535177220406086}}}The return value is a keyed dictionary. Each metric has its own scores, mean, and std:
0.6082359589603785By 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.
# Rebuild the dataset
lung = gw.load_dataset("lung", backend="polars")
# Define two models: simple and full
formula_simple = "Surv(time, status) ~ age + sex"
formula_full = "Surv(time, status) ~ age + sex + ph.ecog"
# Cross-validate both models with the same seed for matched folds
simple_cv = gw.cross_validate(gw.CoxPH(), formula_simple, data=lung, seed=1)
full_cv = gw.cross_validate(gw.CoxPH(), formula_full, data=lung, seed=1)
# Compare mean concordance and fold-to-fold variability
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.043The 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? The calibration() function 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.
# Rebuild the dataset, response, and Cox fit
lung = gw.load_dataset("lung", backend="polars")
death = gw.Outcome.surv(time="time", event="status")
cox = gw.CoxPH().fit("Surv(time, status) ~ age + sex", data=lung)
# Predict one-year survival for each subject
surv_at_year = cox.predict(
lung, type="survival", times=[365.0], format="polars"
)
# Extract predictions as a flat array and compute calibration by quintile
predicted = surv_at_year.slice(0, 1).drop("time").to_numpy().flatten()
gw.calibration(death, predicted=predicted, time=365.0, data=lung, 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.