brier_score_incidence()

IPCW Brier score for a cause-specific cumulative incidence function at specified times.

Usage

brier_score_incidence(
    surv,
    incidence_prob,
    times,
    *,
    data=None,
    cause,
)

Extends the Graf (1999) IPCW Brier score to competing risks following Kretowska (2018). Measures how accurately predicted cumulative incidence probabilities for a specific cause match the observed cause-specific outcomes, correcting for censoring bias.

This is suitable for use after fitting a competing-risks model (FineGray, AalenJohansen, MultiState) that produces predicted CIF curves, evaluate whether the predicted cause-specific probabilities are well-calibrated.

Interpretation:

Ranges from 0 (perfect predictions) to 1 (worst possible).
Lower is better. Compare to a null model that predicts the marginal CIF at each time.
Score typically increases with time (longer horizons are harder to predict).

Parameters

surv: Surv | Outcome | str: A multi-state Surv response (from gw.Surv() with a categorical event). The status column encodes 0 = censored, 1 = first cause, 2 = second cause, etc. An Outcome or a formula response such as 'Surv(time, status == 2)' is also accepted, with its columns read from data.
incidence_prob: Any: Predicted cumulative incidence probabilities for the cause of interest, shape (n_subjects, n_times). Each entry is a predicted probability that cause cause has occurred by the corresponding time. Values should be in [0, 1].
times: Any: Evaluation times where Brier scores are computed. 1-D array-like. Must have length equal to the second dimension of incidence_prob.
cause: int: The cause of interest, as an integer event code from the Surv response. For example, if the response’s states are ("relapse", "death"), then cause=1 evaluates predictions for relapse and cause=2 evaluates predictions for death.
data: Any = None: A data frame (pandas, Polars, PyArrow, DuckDB, a lazy frame, …) holding the columns named by the response and incidence_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

The cause-specific IPCW Brier score at time t for cause k is:

BS_k(t) = \frac{1}{n} \sum_{i=1}^{n} \hat{\omega}_i(t) \bigl(\mathbb{1}(T_i \le t,\, \Delta_i = k) - \hat{F}_k(t \mid \mathbf{x}_i)\bigr)^2

where \hat{F}_k(t \mid \mathbf{x}_i) is the predicted CIF for cause k, and the IPCW weights are:

\hat{\omega}_i(t) = \frac{\mathbb{1}(T_i \le t,\, \Delta_i > 0)}{\hat{G}(T_i)} + \frac{\mathbb{1}(T_i > t)}{\hat{G}(t)}

Subjects censored before t receive zero weight. Subjects who experienced any event (including competing causes) before t are weighted by the inverse censoring probability at their event time. Subjects still at risk after t are weighted by the inverse censoring probability at t.

Examples

Name the competing endpoints in the mgus2 dataset (progression to PCM vs. death) and evaluate a marginal Aalen-Johansen CIF as a naive baseline model:

import greenwood as gw
import numpy as np

# The first of progression (PCM) or death, each with its own time and status columns
mgus2 = gw.load_dataset("mgus2", backend="polars")
first = gw.Outcome.first_event(
    endpoints={"pcm": ("ptime", "pstat"), "death": ("futime", "death")}
)

# Fit Aalen-Johansen for the marginal CIF
aj = gw.AalenJohansen().fit(first, data=mgus2)
cif_df = aj.to_frame(format="polars")
pcm_cif = cif_df.filter(cif_df["cause"] == "pcm")

# Interpolate the marginal CIF at evaluation times
times = 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]
)

# Naive model: same marginal CIF for every subject
probs = np.tile(marginal, (mgus2.height, 1))
bs = gw.brier_score_incidence(
    first, incidence_prob=probs, times=times, cause=1, data=mgus2
)
bs

array([0.03293805, 0.05966165, 0.0898196 ])