brier_score()

IPCW (Graf) Brier score of predicted survival probabilities at specified times.

Usage

brier_score(
    surv,
    survival_prob,
    times,
    *,
    data=None,
)

Measures calibration and accuracy of predicted survival probabilities at fixed time points using inverse-probability-of-censoring-weighted (IPCW) averaging. The Brier score is the mean squared difference between predicted and observed outcomes, weighted so censored subjects contribute honestly without bias.

Interpretation:

Ranges from 0 (perfect predictions) to 1 (worst possible).
Lower is better. A Brier score of 0.25 means, on average, predictions are off by \pm 0.5 in terms of squared deviation.
Null model baseline: A model predicting 50% survival at every time has Brier score \approx 0.25. Compare your model to this baseline to assess practical improvement.
Score typically increases with time (harder to predict farther into the future).

Practical use: After fitting a survival model (Cox, parametric, flexible), evaluate calibration at important clinical horizons (e.g., 1-year, 5-year survival). Compute Brier scores at multiple times, then use integrated_brier_score() for a single summary.

Parameters

surv: Surv | Outcome | str: A right-censored Surv response (time-to-event data). An Outcome or a formula response such as 'Surv(time, status == 2)' is also accepted, with its columns read from data.
survival_prob: Any: Predicted survival probabilities, shape (n_subjects, n_times). Each entry is a predicted probability of surviving beyond the corresponding time. Must be between 0 and 1. For example: columns from CoxPH.predict(type="survival", times=...) (excluding the time column), transposed so rows are times and columns are subjects.
times: Any: Evaluation times where Brier scores are computed. 1-D array-like. Must have length equal to the second dimension of survival_prob.
data: Any = None: A data frame (pandas, Polars, PyArrow, DuckDB, a lazy frame, …) holding the columns named by the response and survival_prob. When surv is an Outcome or a formula, rows with a missing value in any column used are dropped first, along with the matching rows of any arrays passed alongside.

Returns

ndarray: Brier score at each time, shape (len(times),). Lower is better.

Details

Graf (IPCW) Brier score: The unbiased Brier score under censoring is

BS(t) = E\left[(S(t) - \hat{S}(t))^2 \cdot \text{weights}\right],

where:

S(t) is the true survival status at time t (1 if alive, 0 if dead)
\hat{S}(t) is the predicted survival probability
\text{weights} are inverse-probability-of-censoring: inverse of the censoring survival function

Mathematically:

BS(t) = \frac{1}{n} \sum_{\text{dead at } t} \frac{(\hat{S}_i(t))^2}{G(t_i^-)} + \frac{1}{n} \sum_{\text{alive at } t} \frac{(1 - \hat{S}_i(t))^2}{G(t)}

where G(u) is the Kaplan-Meier estimate of the censoring distribution (probability of not being censored).

Advantages over MSE: IPCW weighting makes the Brier score unbiased under censoring, unlike naive MSE which would be biased (censored subjects look artificially “correct”).

Time dependence: Brier scores typically increase with time in chronic-disease settings (longer prediction horizons are harder). They may decrease in acute-illness settings (events cluster early, predictions stabilize).

Examples

Fit a Cox model on the lung dataset and evaluate its calibration (how accurately does it predict survival?) at a few horizons.

import greenwood as gw
import numpy as np

# Load data and name the response columns
lung = gw.load_dataset("lung", backend="polars")
death = gw.Outcome.surv(time="time", event="status")
cox = gw.CoxPH().fit(death, covariates=["age", "sex"], data=lung)

# Compute Brier scores at three clinically relevant horizons
times = [180, 365, 540]
surv_pred = cox.predict(lung, type="survival", times=times, format="pandas")
probs = surv_pred.iloc[:, 1:].to_numpy().T
brier = gw.brier_score(death, survival_prob=probs, times=times, data=lung)
brier

array([0.19160509, 0.23614214, 0.18644928])

Brier scores at three time points. Scores typically increase over time. Compare to a null model (all subjects at 50% survival) to assess improvement:

# Compare model calibration against a naive 50% baseline
null_probs = np.full_like(probs, 0.5)
null_brier = gw.brier_score(death, survival_prob=null_probs, times=times, data=lung)
print(f"Null model Brier: {null_brier}")
print(f"Cox model Brier: {brier}")
print(f"Improvement: {null_brier - brier}")

Null model Brier: [0.25 0.25 0.25]
Cox model Brier: [0.19160509 0.23614214 0.18644928]
Improvement: [0.05839491 0.01385786 0.06355072]

Summarize Brier scores across times into a single summary via time-averaged Brier score:

# Summarize calibration as a single integrated Brier score
ibs = gw.integrated_brier_score(death, survival_prob=probs, times=times, data=lung)
print(f"Integrated Brier Score: {ibs:.3f}")

Integrated Brier Score: 0.213