Brier score, NRI and IDI

What these three scores add to the C-index

The C-index only looks at ranking. Calibration only looks at the scale. Once both have been checked, a paper usually reports three more numbers:

• Brier score — one overall score that mixes discrimination and calibration together
• NRI (net reclassification improvement) — how many patients the new model moves into the right band
• IDI (integrated discrimination improvement) — how much further apart the mean predictions of the two outcome groups have been pulled

All three are computed here on exactly the same two models and the same patients as the decision curve page: developed on survival::rotterdam (2982 patients, 1713 events), validated on survival::gbsg (686 patients, 299 events), horizon 5 years.

Model A has the higher C-index by 0.020, but it was never recalibrated: it predicts 43.1% on average, where this cohort actually runs at 50.8% — it underpredicts systematically. Every tension on this page grows out of that.

The Brier score: two things inside one number

The Brier score is the mean squared error between predicted probability and observed outcome:

$y_i$ is 0 or 1 and $\hat{p}_i$ is the model’s probability. Lower is better; the range is 0 to 1.

It decomposes into a calibration term and a discrimination term (the calibration-refinement decomposition of DeGroot and Fienberg 1983; Murphy 1973 is the three-term version, which carries an extra uncertainty term that depends only on the event rate). So a Brier score that improves may mean better ranking, or merely a corrected scale — the number alone cannot tell you which. And it has no scale of its own. Is 0.228 good or bad? The only way to answer is against a reference. The reference used here is written into the stats file: hand every patient this cohort’s own Kaplan-Meier 5-year risk of 50.8%. A model that cannot beat that has told the clinic nothing it did not already know from the base rate.

Computing the Brier score directly on censored data is wrong

The awkward part of a 5-year risk model is that some patients were never followed for five years. In this cohort, 280 of 686 were censored before the 5-year mark, or 40.8%. The remaining 406 have a known 5-year status: 285 had the event and 121 did not.

The intuitive fix is to score only the patients whose status is known. That is complete-case scoring, and it is not the same estimate with more noise — it estimates a different quantity: among the patients who were followed the event rate is 70.2%, while this cohort’s 5-year risk is 50.8%.

The correct approach is IPCW (inverse probability of censoring weighting): estimate the censoring distribution $G(t)$ with a reverse Kaplan-Meier, then weight everyone who still carries information:

Patients with the event before $t$ get $1/\hat{G}(T_i^-)$, patients still event-free at $t$ get $1/\hat{G}(t)$, and patients censored earlier get a weight of zero — they are not discarded, they are spoken for by the patients like them who were followed.

Scaled Brier: a reference gives it a scale

Divide the Brier score by the reference Brier score and subtract from one, and you get the scaled Brier score (also called the Brier skill score):

Zero means no better than handing everyone the base rate; one means perfect. This is the number you can actually use to say how much a model adds.

The C-index favours model A; the Brier score favours model B. The two scaled-Brier intervals overlap heavily, but interval overlap is not an argument — the two estimates come from the same patients, so their errors are almost perfectly correlated and what matters is the interval on the difference itself. Resampled inside the same bootstrap, A minus B is −0.004, 95% interval −0.045 to +0.024, which crosses zero. The defensible statement is therefore that no difference in overall accuracy between the two models was detected in this cohort — not that B is better.

Recomputed at every horizon, which also shows where the estimate stops being trustworthy

The Brier score belongs to one horizon; change the horizon and the answer changes with it:

The reclassification table: how it is laid out and how to read it

The NRI and the IDI are both built on a reclassification table. Cut the predicted risks into bands, then look at which band each patient falls into under the old and the new model.

The layout is fixed throughout this page: the old model (B) runs down the rows, one row per old band; the new model (A) runs across the columns. Patients on the diagonal are the ones the two models agree about, cells above the diagonal are patients the new model moved up, and cells below it are patients it moved down.

The event group and the non-event group must be tabulated separately. Moving a patient up is the right call for someone who went on to have the event and the wrong call for someone who did not; a single combined table hides who was helped and who was harmed.

These tables are counts, so they can only hold the 406 patients whose 5-year status is known. The 280 censored earlier do not fit into integer cells. That has consequences, and the last section measures them with IPCW weights.

How the NRI is computed

Add the net movement of the two tables and you have the categorical NRI:

First term: for the event group, moving up is right, so it is up minus down. Second term: for the non-event group, moving down is right, so it is down minus up. Each term runs from −1 to +1, and the sum from −2 to +2.

Why the NRI comes out positive here

Notice something first: apart from the scheme in which nobody moved at all, every scheme has a negative event term and a large positive non-event term. That is not a coincidence. It is the direct consequence of model A underpredicting — it pushes almost every patient down a band.

Take the 15/30 scheme. The patients moved down are 30 who had the event and 38 who did not. The first group was harmed (they were going to have the event, and their risk was lowered); the second was helped. The NRI adds the two proportions directly:

• Event group: 30 out of 285, recorded as −0.105
• Non-event group: 38 out of 121, recorded as +0.314
• Sum: +0.209

The cut-points decide everything

One pair of models, one set of patients. Change the cut-points and the NRI runs from +0.000 to +0.209, and the last scheme is negative.

The net benefit difference is in net true positives per 1000 patients, model A minus model B, defined exactly as on the decision curve page.

Three things are worth taking one at a time:

The first scheme has an NRI of exactly +0.000, because nobody moved. Look back at its two tables: all 285 event patients and all 121 non-event patients land in the top band under both models. Both cut-points sit below this cohort’s risk distribution, so the scheme asked nothing at all — which is what cut-points carried over from another disease and another population, never checked against the local data, look like.

The last scheme reverses the sign of the NRI, and it was chosen after the answer was known. That scheme carries a chosen_to_demonstrate flag in the stats file, and the marks in the two tables above and the asterisk on the figure are rendered from it. It is left on the page rather than tidied away because this choice is exactly what the page is criticising: when a paper reports one scheme, the reader has no way of knowing how many were tried.

Within one scheme, the NRI and the net benefit need not share a sign. At the lower cut of the 50/60 scheme, model A’s net benefit is lower than model B’s by 23.55 per 1000, while at the upper cut it is 3.04 higher.

The IDI: no cut-points, but not free either

The IDI drops cut-points entirely and compares the mean predicted risk of the two outcome groups:

Equivalently, it is the increase in the discrimination slope — the mean prediction among events minus the mean prediction among non-events.

The IDI is 0.0400, bootstrap 95% CI 0.0139 to 0.0563. That interval excludes zero, so it can be stated plainly: model A separates the event and non-event groups further than model B does.

The complete-case tables are biased: the IPCW-weighted comparison

Every reclassification table above holds integer counts, so none of them can hold the 280 patients censored early. Recomputing the same NRI and IDI with IPCW weights shows what that costs:

The IDI behaves the same way: 0.0400 on complete cases against 0.0371 weighted.

The gaps are small, but they point the same way and they are not rounding: apart from the scheme that moved nobody at all, where both are exactly zero, the complete-case version is a little more optimistic every time. After weighting the effective sample sizes are 348.701 events and 337.107 non-events, both larger than the 285 and 121 of the integer tables — the missing patients are back, and they pull the answer down.

In this dataset the bias is too small to change any conclusion. Whether it stays small depends on how heavy the censoring is and on whether censoring is related to the predictions, and both are things a paper has to report rather than assume.

Running it yourself

library(survival)
data(cancer, package = "survival")
# The same cohort harmonisation and model fit as B5-05, spelled out so
# this block runs on its own.


H <- 5 * 365.25

rot <- rotterdam
rot$rfs_time  <- pmin(rot$rtime, rot$dtime)
rot$rfs_event <- as.integer(rot$recur == 1 | rot$death == 1)
rot$size_mm   <- c("<=20" = 15, "20-50" = 35, ">50" = 60)[as.character(rot$size)]
rot$grade3    <- as.integer(rot$grade >= 3)

ext <- gbsg
ext$rfs_time  <- ext$rfstime; ext$rfs_event <- ext$status
ext$size_mm   <- ext$size;    ext$grade3    <- as.integer(ext$grade >= 3)

# The two candidate models, defined exactly as this page's figure script does.
BIG   <- c("age", "meno", "size_mm", "grade3", "nodes", "log_pgr", "log_er", "hormon")
SMALL <- c("age", "size_mm", "nodes", "grade3")
rot$log_pgr <- log1p(rot$pgr); rot$log_er <- log1p(rot$er)
ext$log_pgr <- log1p(ext$pgr); ext$log_er <- log1p(ext$er)
keep <- c("rfs_time", "rfs_event", BIG)
dev  <- rot[complete.cases(rot[, keep]), keep]
ext  <- ext[complete.cases(ext[, keep]), keep]

# Breslow cumulative baseline hazard, written out: basehaz(centered = FALSE)
# does not match predict(type = "lp"), which is centred, and pairing them
# leaves every predicted risk short by a constant factor.
breslow_h0 <- function(time, event, lp, h = H) {
  o <- order(time); time <- time[o]; event <- event[o]; r <- exp(lp[o])
  at_risk <- rev(cumsum(rev(r)))
  sum(1 / at_risk[event == 1 & time <= h])
}
form_of <- function(v) as.formula(paste("Surv(rfs_time, rfs_event) ~",
                                        paste(v, collapse = " + ")))
fit_big   <- coxph(form_of(BIG),   data = dev)
fit_small <- coxph(form_of(SMALL), data = dev)
h0_big   <- breslow_h0(dev$rfs_time, dev$rfs_event,
                       as.numeric(as.matrix(dev[, BIG]) %*% coef(fit_big)))
lp_big   <- as.numeric(as.matrix(ext[, BIG])   %*% coef(fit_big))
lp_small <- as.numeric(as.matrix(ext[, SMALL]) %*% coef(fit_small))

# A: eight predictors, exactly as it left development. B: four predictors,
# then recalibrated to this cohort.
pred_A <- 1 - exp(-h0_big * exp(lp_big))
sl_B   <- unname(coef(coxph(Surv(rfs_time, rfs_event) ~ lp_small, data = ext))[1])
h0_B   <- breslow_h0(ext$rfs_time, ext$rfs_event, sl_B * lp_small)
pred_B <- 1 - exp(-h0_B * exp(sl_B * lp_small))


# 1. Censoring distribution: a Kaplan-Meier on the flipped event indicator
Gfit <- survfit(Surv(rfs_time, 1 - rfs_event) ~ 1, data = ext)
G_at    <- function(u) summary(Gfit, times = u, extend = TRUE)$surv
G_minus <- function(u) sapply(u, function(x)          # left limit, for the cases
  summary(Gfit, times = max(0, x - 1e-8), extend = TRUE)$surv)

# 2. IPCW weights. Patients censored early get zero: not discarded, but
#    represented by the patients like them who were followed.
ipcw <- function(d, t) {
  is_case <- d$rfs_time <= t & d$rfs_event == 1
  is_ctrl <- d$rfs_time >  t
  w <- numeric(nrow(d))
  w[is_case] <- 1 / G_minus(d$rfs_time[is_case])
  w[is_ctrl] <- 1 / G_at(t)
  list(w = w, case = is_case, ctrl = is_ctrl)
}

# 3. The IPCW Brier score and its reference
brier_ipcw <- function(pred, t = H) {
  z <- ipcw(ext, t)
  contrib <- numeric(nrow(ext))
  contrib[z$case] <- (1 - pred[z$case])^2
  contrib[z$ctrl] <- pred[z$ctrl]^2
  mean(z$w * contrib)
}
base <- 1 - summary(survfit(Surv(rfs_time, rfs_event) ~ 1, ext),
                    times = H)$surv
1 - brier_ipcw(pred_A) / brier_ipcw(rep(base, nrow(ext)))   # scaled Brier

# 4. The complete-case version, for contrast. Only the `known` line differs.
known <- ext$rfs_time > H | (ext$rfs_time <= H & ext$rfs_event == 1)
y     <- as.integer(ext$rfs_time <= H & ext$rfs_event == 1)
mean((y[known] - pred_A[known])^2)   # the reference becomes the followed-up event rate,
                                     # which is not this cohort's risk

# 5. The NRI: band, count the moves, add the two proportions
band <- function(p, cuts) cut(p, breaks = c(-Inf, cuts, Inf), labels = FALSE)
nri <- function(cuts) {
  a <- band(pred_A, cuts); b <- band(pred_B, cuts)
  ev <- known & y == 1; ne <- known & y == 0
  (sum(a[ev] > b[ev]) - sum(a[ev] < b[ev])) / sum(ev) +   # events: up minus down
  (sum(a[ne] < b[ne]) - sum(a[ne] > b[ne])) / sum(ne)     # non-events: down minus up
}
sapply(list(c(.10, .20), c(.15, .30), c(.30, .60), c(.50, .60)), nri)
# A new answer for every scheme -- a paper reporting one leaves the reader
# unable to tell how many were tried.

# 6. The IDI: no cut-points, but no absolute levels either
ev <- known & y == 1; ne <- known & y == 0
(mean(pred_A[ev]) - mean(pred_B[ev])) - (mean(pred_A[ne]) - mean(pred_B[ne]))

import numpy as np
from lifelines import KaplanMeierFitter
# The same cohort harmonisation and model fit as B5-05, spelled out so
# this block runs on its own.

import pandas as pd
from lifelines import CoxPHFitter

RD  = "https://vincentarelbundock.github.io/Rdatasets/csv/"
rot = pd.read_csv(RD + "survival/rotterdam.csv")
ext = pd.read_csv(RD + "survival/gbsg.csv")

rot["rfs_time"]  = rot[["rtime", "dtime"]].min(axis=1)
rot["rfs_event"] = ((rot["recur"] == 1) | (rot["death"] == 1)).astype(int)
rot["size_mm"]   = rot["size"].map({"<=20": 15, "20-50": 35, ">50": 60})
rot["grade3"]    = (rot["grade"] >= 3).astype(int)
ext["rfs_time"]  = ext["rfstime"]
ext["rfs_event"] = ext["status"]
ext["size_mm"]   = ext["size"]
ext["grade3"]    = (ext["grade"] >= 3).astype(int)

v    = ["age", "size_mm", "nodes", "grade3"]
beta = (CoxPHFitter()
        .fit(rot[v + ["rfs_time", "rfs_event"]], "rfs_time", "rfs_event")
        .params_[v].to_numpy())

H = 5 * 365.25
G = KaplanMeierFitter().fit(ext["rfs_time"], 1 - ext["rfs_event"])   # censoring

def ipcw(t=H):
    case = (ext["rfs_time"].values <= t) & (ext["rfs_event"].values == 1)
    ctrl = ext["rfs_time"].values > t
    w = np.zeros(len(ext))                       # censored early stay at 0
    w[case] = 1 / G.predict(ext["rfs_time"].values[case] - 1e-8).values
    w[ctrl] = 1 / float(G.predict(t))
    return w, case, ctrl

def brier_ipcw(pred, t=H):
    w, case, ctrl = ipcw(t)
    c = np.zeros(len(ext))
    c[case] = (1 - pred[case]) ** 2
    c[ctrl] = pred[ctrl] ** 2
    return float(np.mean(w * c))

def nri(cuts, pred_new, pred_old, y, sel):
    a, b = np.digitize(pred_new, cuts), np.digitize(pred_old, cuts)
    ev, ne = sel & (y == 1), sel & (y == 0)
    return ((a[ev] > b[ev]).sum() - (a[ev] < b[ev]).sum()) / ev.sum() + \
           ((a[ne] < b[ne]).sum() - (a[ne] > b[ne]).sum()) / ne.sum()

Five things to ask when reading a paper

1. Is there a reference next to the Brier score? Without one it has no scale. At minimum there should be a scaled Brier score, or the score of a model that hands everyone the prevalence.
2. How was censoring handled? A Brier score reported on survival data with no mention of IPCW is usually complete-case scoring, and that can be wrong down to the sign.
3. Where did the NRI cut-points come from? Were they specified in advance, clinically motivated, and is one scheme reported? Several schemes with only the best discussed looks identical on paper to one scheme with several tried.
4. Are the two NRI components reported separately? The sum alone hides an event term that got worse, which is what three of the four schemes on this page do.
5. Is a decision curve or net benefit reported alongside? The NRI is a unitless sum; net benefit has units and is tied to a trade-off someone had to state out loud.

Common misuses

What this page concludes

One pair of models, one set of patients, one 5-year horizon:

• The C-index favours model A
• The Brier and scaled Brier scores favour model B
• The NRI runs from +0.000 to +0.209 across four schemes, and the scheme chosen after the fact gives −0.045
• The net benefit difference changes sign 10 times across the threshold grid, with intervals that generally cross zero

None of these scores is wrong. They answer different questions, and most of those are not the question the clinic is asking. The clinical question is whether deciding by this model beats current practice, and that one is answered by the decision curve — which presupposes that the model’s scale is right in the first place, the business of the calibration page. The Brier score, the NRI and the IDI belong beside those two as supporting description, not in place of them.

Reproducing every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B5-10-brier-nri-idi.R

Report a problem on this page