18  Performance Assessment and Model Validation

18.1 Introduction

You have developed a clinical prediction model. Its coefficients look reasonable, the predictors make clinical sense, and the apparent performance seems promising. But how good is this model, really? And will it work on patients beyond your development dataset?

This chapter addresses these questions. We cover the five domains of model performance — discrimination, calibration, overall performance, classification, and clinical utility — and the methods for internal and external validation. The structure follows the comprehensive framework proposed by Van Calster et al. (2025) in The Lancet Digital Health, which we strongly recommend as a companion reading.

Two of these terms come up constantly, so it is worth fixing them now in plain language. Discrimination asks whether the model can tell sick patients apart from well ones — can it put the patients who will have the outcome higher up the risk ranking than those who will not? Calibration asks whether the actual numbers are right — if the model says “20% risk,” do about 20 in 100 such patients really have the event? These are two genuinely different questions, and a model can pass one while failing the other (Figure 18.1).

The core message is this: no single metric tells the full story. A model with excellent discrimination (high AUC) can still produce dangerously miscalibrated risk predictions. A well-calibrated model might not discriminate well enough to be useful. Clinical utility must be assessed separately from statistical performance. You need the full picture.

flowchart TB
    M["Your prediction model"]
    M --> D["Discrimination<br/>Does it RANK patients correctly?<br/>(sicker patients get higher risk scores)"]
    M --> C["Calibration<br/>Are the NUMBERS accurate?<br/>(20% risk = 20 in 100 have the event)"]
    D --> N["Both must be checked.<br/>Good ranking does NOT guarantee accurate probabilities."]
    C --> N

Figure 18.1: Discrimination and calibration ask two different questions about the same model. A model can do well on one and badly on the other. Framing follows Steyerberg (2019) and Van Calster et al. (2025).

18.2 The Five Performance Domains

Van Calster et al. (2025) organise model performance assessment into five complementary domains:

Domain	Question	Key Metrics
Discrimination	Can the model rank patients by risk?	C-statistic, AUC
Calibration	Are predicted probabilities accurate?	Calibration plot, O:E ratio, calibration slope
Overall performance	How close are predictions to observed outcomes?	Brier score, Nagelkerke R-squared
Classification	How well does the model classify at a threshold?	Sensitivity, specificity, PPV, NPV
Clinical utility	Does using the model improve clinical decisions?	Net benefit, decision curve analysis

We will examine each in detail.

18.3 Domain 1: Discrimination

18.3.1 The C-Statistic / AUC

Discrimination measures how well a model separates patients who experience the outcome from those who do not. The most common measure is the concordance statistic (C-statistic) — “concordance” simply meaning agreement between the ranking and reality. For yes/no (binary) outcomes the C-statistic is the same number as the area under the ROC curve (AUC) discussed in Chapter 15.

Here is the plain-English interpretation. Pick one patient who had the outcome and one who did not, at random. The C-statistic is the probability that the model gave the higher risk score to the one who actually had the outcome. A value of 0.5 is a coin toss (no better than guessing); 1.0 is a perfect ranking every time. So a C-statistic of 0.78 means that if you compared a patient who died with one who survived, the model gave the higher predicted risk to the patient who died about 78% of the time.

Why a clinician should care. Good discrimination is what lets a model sort patients into “watch closely” versus “reassure.” If the model cannot reliably put higher-risk patients above lower-risk ones, no threshold you pick will separate them well.

R

library(rms)

# Simulate clinical data: predicting 30-day mortality after stroke
set.seed(2024)
n <- 1500

stroke_data <- data.frame(
  age = round(rnorm(n, 72, 12)),
  nihss = round(pmax(0, rnorm(n, 8, 6))), # NIH Stroke Scale
  glucose = round(rnorm(n, 140, 50)),
  afib = rbinom(n, 1, 0.25), # Atrial fibrillation
  thrombolysis = rbinom(n, 1, 0.30)
)

# True model for 30-day mortality
lp <- -5 +
  0.04 * stroke_data$age +
  0.12 * stroke_data$nihss +
  0.003 * stroke_data$glucose +
  0.3 * stroke_data$afib -
  0.5 * stroke_data$thrombolysis

stroke_data$death_30d <- rbinom(n, 1, plogis(lp))
cat("30-day mortality rate:", round(mean(stroke_data$death_30d), 3), "\n")

# Fit prediction model
dd <- datadist(stroke_data)
options(datadist = "dd")

fit <- lrm(
  death_30d ~ age + nihss + glucose + afib + thrombolysis,
  data = stroke_data,
  x = TRUE,
  y = TRUE
)

cat("C-statistic (apparent):", round(fit$stats["C"], 3), "\n")

Python

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from scipy.special import expit

np.random.seed(2024)
n = 1500

stroke_data = pd.DataFrame({
    "age": np.round(np.random.normal(72, 12, n)),
    "nihss": np.round(np.maximum(0, np.random.normal(8, 6, n))),
    "glucose": np.round(np.random.normal(140, 50, n)),
    "afib": np.random.binomial(1, 0.25, n),
    "thrombolysis": np.random.binomial(1, 0.30, n)
})

lp = (-5 + 0.04 * stroke_data["age"] +
      0.12 * stroke_data["nihss"] +
      0.003 * stroke_data["glucose"] +
      0.3 * stroke_data["afib"] -
      0.5 * stroke_data["thrombolysis"])

stroke_data["death_30d"] = np.random.binomial(1, expit(lp))
print(f"30-day mortality rate: {stroke_data['death_30d'].mean():.3f}")

predictors = ["age", "nihss", "glucose", "afib", "thrombolysis"]
X = stroke_data[predictors].values
y = stroke_data["death_30d"].values

model = LogisticRegression(max_iter=5000, random_state=42)
model.fit(X, y)

pred_prob = model.predict_proba(X)[:, 1]
auc = roc_auc_score(y, pred_prob)
print(f"C-statistic (apparent): {auc:.3f}")

What the code shows. This simulates a stroke cohort, fits a model for 30-day mortality, and prints the C-statistic. The number to read is that final C-statistic: 0.5 means the model is no better than a coin toss at ranking patients, while 1.0 means it perfectly ranks every patient who died above every patient who survived. In clinical practice values around 0.7–0.8 are common and often useful. This is the apparent C-statistic — measured on the same data the model was built from — so treat it as an optimistic upper bound until it has been validated. Note that discrimination only asks whether the ranking is right; it says nothing about whether the predicted probabilities themselves are accurate, which is what calibration (next section) checks.

18.3.2 Limitations of the C-Statistic

The C-statistic has important limitations that Van Calster et al. (2025) highlight:

1. Insensitive to calibration: A model can rank patients perfectly (C = 1.0) yet attach completely wrong probabilities to them. Good ranking does not mean the numbers are right.
2. Context-blind: It averages the model’s ranking ability over every possible risk threshold, weighting them all equally — including thresholds so extreme that no clinician would ever act on them. The thresholds you would actually use in your clinic get no more say in the number than the ones you would never touch.
3. Hard to improve: In many clinical areas the highest achievable C-statistic is naturally limited. A genuinely useful new biomarker might barely move the C-statistic even though it improves predictions.
4. Misleading comparisons: A C-statistic of 0.75 in one population cannot be compared directly with 0.75 in another if the patient mix (“case-mix”) differs. An easier mix of patients inflates the number.

18.4 Domain 2: Calibration

18.4.1 Why Calibration Matters More Than You Think

Calibration assesses whether predicted probabilities match observed outcomes. If your model says a patient has a 20% risk, do approximately 20 out of 100 such patients actually experience the outcome?

This matters enormously in clinical practice. When a clinician tells a patient “your model-predicted 10-year cardiovascular risk is 15%,” the patient and clinician rely on that number being accurate — not just on the patient’s relative ranking. A model can rank patients perfectly and still tell every one of them the wrong number. If treatment is offered above a fixed risk threshold (say, statins above 10%), miscalibrated risks send the wrong patients for treatment. Risk-based treatment decisions require calibrated predictions.

18.4.2 Calibration-in-the-Large: The O:E Ratio

The coarsest calibration check asks whether the model gets the overall risk level right — on average, does it predict about as many events as actually happened? This is called calibration-in-the-large (“in-the-large” meaning across the whole cohort, not patient by patient). The simplest way to measure it is the observed-to-expected (O:E) ratio:

\[ \text{O:E ratio} = \frac{\text{Observed event rate}}{\text{Mean predicted probability}} \]

An O:E ratio of 1 means the model gets the overall event rate exactly right. An O:E ratio of 1.5 means more events happened than the model expected — it underestimates risk by 50%. A ratio below 1 means it overestimates.

18.4.3 Calibration Slope

Calibration-in-the-large only checks the average. The calibration slope goes further and checks whether the spread of the predictions is right. To get it, you take the model’s own predicted risks and ask how steeply the actual outcomes track them — fitting a fresh, deliberately simple model with the original predictions as its only input. The slope that comes back tells you how stretched or compressed the predictions are:

• A slope of 1 is ideal — the predictions are spread out by exactly the right amount.
• A slope below 1 means the predictions are too extreme: the high risks are too high and the low risks too low. This is the classic fingerprint of overfitting (the model learned quirks of the training data that do not generalise).
• A slope above 1 means the opposite — predictions are too bunched towards the middle.

18.4.4 Calibration Plots

The most informative way to assess calibration is visually, through a calibration plot. This plots predicted probabilities (x-axis) against the proportion of patients who actually had the outcome (y-axis). A perfectly calibrated model follows the 45-degree diagonal: predicted and observed agree everywhere. Points above the line mean the model underestimated risk for those patients; points below mean it overestimated.

Two types are common:

• Grouped calibration plot: Patients are divided into groups (e.g., deciles) by predicted probability, and the observed proportion within each group is plotted.
• Smoothed calibration plot: Instead of grouping, a smooth curve is drawn through the predicted-versus-observed data (by loess or a spline — both are ways of letting a curve follow the local pattern of the data without assuming any particular shape in advance).

R

library(rms)

# Use the stroke model from above
pred_prob <- predict(fit, type = "fitted")

# Grouped calibration plot.
# Note: val.prob() takes no `main` argument, so the title is added
# afterwards with title() rather than passed into the call.
cal <- val.prob(
  pred_prob,
  stroke_data$death_30d,
  m = 100,
  cex = 0.5
)
title(main = "Calibration Plot: 30-Day Stroke Mortality Model")

# Calculate calibration metrics manually
obs_rate <- mean(stroke_data$death_30d)
mean_pred <- mean(pred_prob)

cat("Observed event rate:", round(obs_rate, 3), "\n")
cat("Mean predicted probability:", round(mean_pred, 3), "\n")
cat("O:E ratio:", round(obs_rate / mean_pred, 3), "\n")

# Calibration slope
cal_model <- glm(stroke_data$death_30d ~ qlogis(pred_prob), family = binomial)
cat("Calibration slope:", round(coef(cal_model)[2], 3), "\n")
cat("Calibration intercept:", round(coef(cal_model)[1], 3), "\n")

# Brier score, its no-predictor benchmark, and the scaled version
brier <- mean((pred_prob - stroke_data$death_30d)^2)
brier_max <- obs_rate * (1 - obs_rate)   # score of a no-predictor model
brier_scaled <- 1 - brier / brier_max

cat("Brier score:", round(brier, 4), "\n")
cat("Brier score, no-predictor benchmark:", round(brier_max, 4), "\n")
cat("Scaled Brier score:", round(brier_scaled, 3), "\n")

Python

import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy.special import logit
from sklearn.calibration import calibration_curve
from sklearn.metrics import brier_score_loss

# Use the stroke model predictions from above
pred_prob = model.predict_proba(X)[:, 1]
y_true = stroke_data["death_30d"].values

# Grouped calibration curve
prob_true, prob_pred = calibration_curve(y_true, pred_prob,
                                         n_bins=10, strategy="quantile")

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Left panel: calibration plot
axes[0].plot(prob_pred, prob_true, "o-", color="steelblue", lw=2,
             markersize=8, label="Model")
axes[0].plot([0, 1], [0, 1], "--", color="grey", label="Perfect calibration")
axes[0].set_xlabel("Mean Predicted Probability")
axes[0].set_ylabel("Observed Proportion")
axes[0].set_title("Calibration Plot (Decile Groups)")
axes[0].legend()
axes[0].set_xlim(0, max(prob_pred) * 1.1)
axes[0].set_ylim(0, max(prob_true) * 1.1)

# Right panel: distribution of predicted probabilities
axes[1].hist(pred_prob[y_true == 0], bins=30, alpha=0.5, color="steelblue",
             label="No event", density=True)
axes[1].hist(pred_prob[y_true == 1], bins=30, alpha=0.5, color="coral",
             label="Event", density=True)
axes[1].set_xlabel("Predicted Probability")
axes[1].set_ylabel("Density")
axes[1].set_title("Distribution of Predicted Probabilities")
axes[1].legend()

plt.tight_layout()
plt.show()

# Calibration metrics
obs_rate = y_true.mean()
mean_pred = pred_prob.mean()
oe_ratio = obs_rate / mean_pred
brier = brier_score_loss(y_true, pred_prob)

# Calibration slope and intercept: regress the outcome on the model's
# own predictions (on the log-odds scale), exactly as the R tab does.
lp = logit(np.clip(pred_prob, 1e-6, 1 - 1e-6))
cal_fit = sm.GLM(y_true, sm.add_constant(lp),
                 family=sm.families.Binomial()).fit()
cal_intercept, cal_slope = cal_fit.params

# Brier score, its no-predictor benchmark, and the scaled version
brier_max = obs_rate * (1 - obs_rate)
brier_scaled = 1 - brier / brier_max

print(f"\nCalibration metrics:")
print(f"  Observed event rate: {obs_rate:.3f}")
print(f"  Mean predicted probability: {mean_pred:.3f}")
print(f"  O:E ratio: {oe_ratio:.3f}")
print(f"  Calibration slope: {cal_slope:.3f}")
print(f"  Calibration intercept: {cal_intercept:.3f}")
print(f"  Brier score: {brier:.4f}")
print(f"  Brier score, no-predictor benchmark: {brier_max:.4f}")
print(f"  Scaled Brier score: {brier_scaled:.3f}")

What the code shows. This produces and quantifies a calibration plot, which answers the question: when the model says “20% risk,” do about 20 in 100 such patients actually have the event? In R, val.prob() draws the plot and reports calibration statistics in one step. In Python the left panel plots observed proportion against mean predicted probability for groups of patients (here deciles), and the right panel shows how predicted risks are distributed among those with and without the event. Read the plot by comparing the model’s curve to the 45-degree diagonal: points on the diagonal mean perfect calibration, points above mean the model underestimates risk, points below mean it overestimates. The printed metrics summarise the same thing numerically — the O:E ratio for overall calibration-in-the-large, the calibration slope for whether predictions are too extreme, and the Brier score for overall accuracy.

NoteCalibration slope, in one line

The calibration slope measures how steeply the actual outcomes track the model’s own predicted risks. A slope of 1 is ideal. A slope below 1 (the usual sign of overfitting) means the predictions are too spread out — the high risks are too high and the low risks too low — and is exactly the kind of problem that shrinkage in the previous chapter is designed to prevent.

18.4.5 Interpreting Miscalibration

Common patterns and their clinical meaning:

• O:E > 1 (underestimation): The model systematically underestimates risk. Patients may not receive treatments they need.
• O:E < 1 (overestimation): The model systematically overestimates risk. Patients may receive unnecessary treatments.
• Calibration slope < 1: Predictions are too extreme. Patients at high predicted risk actually have lower risk than predicted, and patients at low predicted risk have higher risk. This is the hallmark of overfitting.
• Calibration slope > 1: Predictions are too moderate. This is less common and may occur after excessive shrinkage.

18.5 Domain 3: Overall Performance

The Brier score rolls discrimination and calibration into a single number. It is just the average squared gap between what the model predicted and what actually happened:

\[ \text{Brier score} = \frac{1}{N} \sum_{i=1}^{N} (p_i - y_i)^2 \]

where \(p_i\) is the predicted probability and \(y_i\) is the observed outcome (0 or 1). In words: for each patient, take the predicted risk minus the actual outcome, square it, and average over everyone. Lower is better. The Brier score ranges from 0 (perfect) to 1 (worst). As a reference point, a model that just predicts the overall event rate for everyone scores the prevalence times (1 minus the prevalence) — any useful model should beat that.

A raw Brier score is therefore hard to read on its own, because how good 0.19 is depends entirely on how common the outcome is. The scaled Brier score (or Brier skill score) fixes that by expressing the score as an improvement over the do-nothing reference model:

\[ \text{Scaled Brier} = 1 - \frac{\text{Brier score}}{\text{Brier}_{\text{max}}} \qquad \text{where} \qquad \text{Brier}_{\text{max}} = \bar{y}\,(1 - \bar{y}) \]

The symbol \(\bar{y}\) (“y-bar”) is just the observed event rate — the proportion of patients in your data who actually had the outcome. So if 33% of the stroke cohort died within 30 days, \(\bar{y} = 0.33\) and \(\text{Brier}_{\text{max}} = 0.33 \times 0.67 = 0.22\). That is the score you would get from a model with no predictors at all, one that simply announces “33%” for every patient.

Read the scaled Brier score as a percentage of the way from useless to perfect: 0 means the model is no better than quoting the average risk to everyone, and 1 means perfect prediction. A scaled Brier of 0.12 means the model closed 12% of the gap between the do-nothing model and perfection — which is typical, and a useful antidote to a flattering-looking C-statistic. Note that both the raw and scaled Brier scores are now printed by the calibration code in Section 18.4.4, so you can see all three numbers — Brier, benchmark, and scaled — side by side.

18.6 Domain 4: Classification

Classification metrics (sensitivity, specificity, PPV, NPV) were covered in Chapter 15. The key point for model validation is that classification performance depends entirely on the chosen threshold. When reporting classification metrics, always specify the threshold and justify its choice based on the clinical context.

18.7 Domain 5: Clinical Utility

18.7.1 Net Benefit and Decision Curve Analysis

A model might have good discrimination and calibration, yet using it may not actually improve clinical decisions. Decision curve analysis (DCA), proposed by Vickers and Elkin (2006), addresses exactly this. It asks a practical question: if you acted on this model, would patients be better off than if you simply treated everyone or treated no-one? It answers using net benefit — a single score that counts the patients correctly treated and then subtracts a penalty for those treated unnecessarily.

Why a clinician should care. Discrimination and calibration are statistical properties; net benefit speaks the language of decisions. It is the metric that tells you whether the model is worth putting in front of a patient at all.

The net benefit at a given threshold probability \(p_t\) (the risk level above which you would act — treat, refer, monitor) is:

\[ \text{Net benefit} = \frac{TP}{N} - \frac{FP}{N} \times \frac{p_t}{1 - p_t} \]

Here \(TP\) is the number of true positives (correctly flagged patients who had the outcome) and \(FP\) the false positives (flagged patients who did not). The term \(\frac{p_t}{1 - p_t}\) is the exchange rate between false and true positives, set by the threshold you chose.

Where the threshold comes from

This is the part of decision curve analysis that people skip, and it is the part that matters most — because the threshold is not something the data can tell you. It is a clinical judgement that you have to make and defend.

Here is the logic. Suppose you decide to treat any patient whose predicted risk is 10% or higher. By drawing the line exactly there, you have said something quite specific: that a patient at 10% risk is a borderline case, someone for whom treating and not treating feel equally acceptable. And if treating a 10%-risk patient is a coin flip, then you are implicitly saying that one missed event is about nine times worse than one unnecessary treatment — because for every 10 such patients you treat, roughly 1 truly needs it and 9 do not. That ratio, 9 to 1, is exactly \(\frac{1 - p_t}{p_t}\), and its reciprocal \(\frac{p_t}{1 - p_t}\) is the weight the formula puts on each false positive.

So the threshold encodes your view of how the two kinds of mistake compare. Different thresholds mean different clinical worlds:

Threshold \(p_t\)	Exchange rate \(\frac{p_t}{1-p_t}\)	You are saying…	Typical situation
2%	1 : 49	I would treat 49 people unnecessarily to prevent one event	The action is cheap and safe, the outcome catastrophic (e.g. a follow-up scan to rule out a serious cancer)
10%	1 : 9	I would accept 9 unnecessary treatments per event prevented	Statins for cardiovascular prevention — a cheap, well-tolerated daily tablet
30%	1 : 2.3	I would accept only about 2 unnecessary treatments per event prevented	A treatment with real side-effects, or one that commits the patient to months of monitoring
50%	1 : 1	An unnecessary treatment is as bad as a missed event	Major surgery, or a toxic therapy where the harm of over-treatment is comparable to the harm of the disease

Read down that table and notice which way it runs: a low threshold means you are willing to over-treat a lot, and a high threshold means you are not. The more burdensome, risky, or expensive the action you would take, the higher the threshold, because each unnecessary intervention costs the patient more.

Two practical consequences. First, you rarely commit to one number — which is exactly why a decision curve plots a range of thresholds, so that you and your colleagues can find your own on the x-axis. Second, the range worth plotting is the range clinicians would actually entertain. For a decision where nobody would ever treat at 60% risk, the curve beyond 60% is irrelevant, and a model that only wins out there has not shown anything useful. Choose the range from the clinical problem, state it in your paper, and say who decided it — ideally clinicians and patients, not the analyst alone.

R

library(dcurves)

# Add the model's predicted risks to the dataset. The column is deliberately
# named `model_risk` because dcurves uses the column name as the curve label,
# so the figure legend then says where that curve came from.
stroke_data$model_risk <- predict(fit, type = "fitted")

# Decision curve analysis over the thresholds a clinician might plausibly use.
# `label` replaces the bare column name in the legend with a readable one, so
# the figure says what the curve actually is.
dca_result <- dca(
  death_30d ~ model_risk,
  data = stroke_data,
  thresholds = seq(0, 0.5, by = 0.01),
  label = list(model_risk = "Stroke mortality model")
)

plot(dca_result, smooth = TRUE, show_ggplot_code = FALSE) +
  ggplot2::ggtitle("Decision Curve Analysis:\n30-Day Stroke Mortality Model")

Python

import numpy as np
import matplotlib.pyplot as plt

def net_benefit(y_true, y_pred, threshold):
    """Calculate net benefit at a given threshold."""
    n = len(y_true)
    pred_pos = (y_pred >= threshold).astype(int)
    tp = np.sum((pred_pos == 1) & (y_true == 1))
    fp = np.sum((pred_pos == 1) & (y_true == 0))
    nb = tp / n - fp / n * (threshold / (1 - threshold))
    return nb

thresholds = np.arange(0.01, 0.51, 0.01)
# Named `model_risk` so it is obvious these risks come from the model
model_risk = model.predict_proba(X)[:, 1]
y_true = stroke_data["death_30d"].values
n_total = len(y_true)

# Net benefit for model
nb_model = [net_benefit(y_true, model_risk, t) for t in thresholds]

# Net benefit for "treat all" strategy
nb_all = [(y_true.mean() - (1 - y_true.mean()) * t / (1 - t))
          for t in thresholds]

# Net benefit for "treat none" is always 0

plt.figure(figsize=(9, 6))
plt.plot(thresholds, nb_model, color="steelblue", lw=2,
         label="Prediction Model")
plt.plot(thresholds, nb_all, color="coral", lw=2, ls="--",
         label="Treat All")
plt.axhline(y=0, color="black", lw=1, label="Treat None")
plt.xlabel("Threshold Probability")
plt.ylabel("Net Benefit")
plt.title("Decision Curve Analysis: 30-Day Stroke Mortality Model")
plt.legend()
plt.ylim(-0.05, max(max(nb_model), max(nb_all)) * 1.1)
plt.tight_layout()
plt.show()

What the code shows. This computes net benefit across a range of decision thresholds and plots three curves: using the model, treating everyone, and treating no-one. In R, dca() from the dcurves package does the calculation; the Python version writes out the net-benefit formula explicitly in the net_benefit() function so you can see the arithmetic. The result you care about is whether the model’s curve sits above both the “treat all” and “treat none” lines over the thresholds clinicians would actually use — if it does, the model improves decisions there; if it dips below either default, the model adds nothing at that threshold and you may as well follow the simpler default strategy.

NoteNet benefit in clinician’s terms

Net benefit puts true positives and false positives on a common scale so they can be added up. The threshold probability you pick reflects how you trade them off: a 10% threshold means you would accept treating 9 patients unnecessarily to catch 1 who truly needs it. Net benefit then counts the true positives the model finds and subtracts the false positives, weighted by that exchange rate. Unlike the C-statistic, it directly reflects the clinical consequences of acting on the model.

18.7.2 Interpreting the Decision Curve

The decision curve shows net benefit (y-axis) across a range of threshold probabilities (x-axis). Three strategies are compared:

• Treat none (horizontal line at 0): no patients receive the intervention.
• Treat all (declining curve): all patients receive the intervention regardless of risk.
• Model-guided (the model’s curve): only patients above the threshold receive the intervention.

The model is clinically useful at thresholds where its net benefit exceeds both the “treat all” and “treat none” lines. If the model curve is always below one of the default strategies, the model adds no value.

Reading the curve at a specific threshold. Pick your threshold on the x-axis, read the three curves vertically above it, and the highest curve is the strategy you should adopt at that threshold. Working through our stroke example (Section 18.7.1.1 explains where these thresholds come from):

• At 10%: you have said you would accept about 9 unnecessary interventions per death prevented. With 30-day mortality running at roughly a third of this cohort, treating everyone is still a reasonable strategy at such a permissive threshold, so “treat all” is competitive here and the model may add little.
• At 30%: you would accept only about 2 unnecessary interventions per death prevented — appropriate for an intervention with real burden. “Treat all” now carries a heavy false-positive penalty and its curve has fallen steeply, so this is typically where a model earns its keep: it lets you concentrate the intervention on the patients most likely to need it.
• At 40%: you are close to saying an unnecessary intervention is about as bad as a missed death (the exchange rate is roughly 1 to 1.5). Only patients the model puts at high risk clear this bar. Few patients get treated, so the absolute net benefit is small — but if the model’s curve is still above zero, acting on it beats doing nothing for that narrow group.

The unit of the y-axis. Net benefit is measured in “true positives per patient assessed,” which sounds abstract but has a concrete reading: multiply by 100 and you get the number of events correctly identified per 100 patients, after the unnecessary treatments have been paid for. In our stroke example at a 30% threshold, the model scores about 0.118 against about 0.050 for treating everyone. The difference, 0.068, means that using the model identifies roughly 7 more true high-risk patients per 100 assessed than treating everyone does, at no extra cost in over-treatment. That is the sentence to put in a paper — far more informative than reporting the AUC alone.

For contrast, run the same comparison at a 10% threshold and the two strategies come out level (about 0.261 each). At a permissive threshold where you are content to over-treat, a blanket policy does just as well as the model, and the model’s value only appears once the threshold rises. This is why the range of thresholds you plot decides whether a model looks useful or useless — and why that range has to be justified clinically rather than chosen to flatter the model.

A warning. Net benefit always falls as the threshold rises, for every strategy, because fewer and fewer patients qualify. Do not read a declining curve as a deteriorating model. What matters is only the ranking of the three curves at the thresholds you care about, never the height of any one of them on its own.

18.8 Internal Validation

18.8.1 When Do You Do This?

Before the methods, the practical question: at what point in a project does internal validation happen?

The answer is: once, at the end of model development, after every modelling decision has been made — and before you report a single performance number. Concretely, you have finished choosing which predictors to include, how to handle continuous variables, whether to apply shrinkage, and how to deal with missing data. Nothing about the model is still up for negotiation. That is the moment to run internal validation.

Two mistakes are worth naming explicitly, because both are common:

• Do not treat it as an optional final polish once you are already happy with the model. The apparent performance that made you happy is the very number internal validation exists to correct. If you only run it to confirm a figure you have already decided to believe, you have inverted the point of the exercise.
• Do not use it to choose between candidate models. If you run bootstrap validation on five model variants and keep the best-looking one, that selection is itself a modelling decision made on the same data, and the optimism creeps straight back in. Selection must happen inside the validation loop — which is exactly what step 2 below means by “including all modelling decisions” — or not on this data at all.

So the sequence for a development study is: build the model → internally validate it → report the optimism-corrected performance → then, as a separate study on separate data, externally validate it (Section 18.9). Internal validation answers “how well will this model do on patients like the ones I already have?” It does not, and cannot, tell you whether the model travels somewhere new.

18.8.2 Why Internal Validation Is Necessary

As discussed in Chapter 17, apparent performance — performance measured on the very same data the model was trained on — is optimistically biased. The model has, to some degree, memorised that data, so it looks better there than it will on new patients. This gap is called optimism. Internal validation estimates that optimism and subtracts it off, giving a more honest figure. It does not replace external validation, but it is a necessary first step.

18.8.3 Bootstrap Optimism Correction

The recommended approach (Steyerberg 2019; Van Calster et al. 2025) is bootstrap optimism correction. “Bootstrapping” means repeatedly drawing random samples of patients with replacement from your dataset (so some patients appear more than once, others not at all), refitting the model each time, and watching how much the results wobble. Here that wobble is used to measure optimism:

1. Draw a bootstrap sample (with replacement) from the original data.
2. Develop the model in the bootstrap sample (including all modelling decisions).
3. Measure performance in the bootstrap sample (apparent bootstrap performance).
4. Apply that same bootstrap model to the original data (test performance).
5. Optimism = apparent bootstrap performance minus test performance.
6. Repeat steps 1–5 many times (e.g., 200+).
7. Average the optimism across repetitions.
8. Corrected performance = original apparent performance minus average optimism.

The diagram below shows the loop (Figure 18.2).

flowchart TB
    A["Apparent performance<br/>(model scored on its own training data)"]
    B["Draw a bootstrap sample<br/>(resample patients with replacement)"]
    B --> C["Refit the model on the bootstrap sample"]
    C --> D["Score on the bootstrap sample<br/>(optimistic)"]
    C --> E["Score the same model on the original data<br/>(realistic)"]
    D --> F["Optimism = bootstrap score − original score"]
    E --> F
    F -. "repeat ~200×" .-> B
    F --> G["Average the optimism<br/>over all repeats"]
    A --> H["Corrected performance =<br/>apparent − average optimism"]
    G --> H

Figure 18.2: The bootstrap optimism-correction loop. The optimism measured over many resamples is subtracted from the apparent performance to give an honest, corrected estimate. Based on the optimism-correction procedure of Harrell (2015) and Steyerberg (2019).

R

library(rms)

# Using the stroke model
fit <- lrm(
  death_30d ~ age + nihss + glucose + afib + thrombolysis,
  data = stroke_data,
  x = TRUE,
  y = TRUE
)

# Bootstrap validation with 200 resamples
set.seed(42)
val <- validate(fit, B = 200)

knitr::kable(
  as.data.frame(unclass(val)),
  digits = 3,
  caption = "Bootstrap internal validation (200 resamples). `index.corrected` is the optimism-corrected estimate; for `Dxy`, the C-statistic is C = 0.5 + Dxy/2."
)

# Extract optimism-corrected C-statistic
dxy_corrected <- val["Dxy", "index.corrected"]
c_corrected <- (dxy_corrected + 1) / 2

cat("\nApparent C-statistic:", round(fit$stats["C"], 3), "\n")
cat("Optimism:", round(val["Dxy", "optimism"] / 2, 3), "\n")
cat("Optimism-corrected C-statistic:", round(c_corrected, 3), "\n")

# Calibration slope from bootstrap
cat(
  "\nCalibration slope (optimism-corrected):",
  val["Slope", "index.corrected"],
  "\n"
)

Python

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, brier_score_loss

np.random.seed(42)

predictors = ["age", "nihss", "glucose", "afib", "thrombolysis"]
X = stroke_data[predictors].values
y = stroke_data["death_30d"].values
n = len(y)

# Apparent performance
full_model = LogisticRegression(max_iter=5000, random_state=42)
full_model.fit(X, y)
pred_full = full_model.predict_proba(X)[:, 1]
apparent_auc = roc_auc_score(y, pred_full)
apparent_brier = brier_score_loss(y, pred_full)

# Bootstrap optimism correction
n_boot = 200
optimism_auc = []
optimism_brier = []

for b in range(n_boot):
    boot_idx = np.random.choice(n, n, replace=True)
    X_boot, y_boot = X[boot_idx], y[boot_idx]

    boot_model = LogisticRegression(max_iter=5000, random_state=42)
    boot_model.fit(X_boot, y_boot)

    # Apparent performance on bootstrap sample
    pred_boot = boot_model.predict_proba(X_boot)[:, 1]
    auc_boot = roc_auc_score(y_boot, pred_boot)
    brier_boot = brier_score_loss(y_boot, pred_boot)

    # Test performance on original data
    pred_orig = boot_model.predict_proba(X)[:, 1]
    auc_orig = roc_auc_score(y, pred_orig)
    brier_orig = brier_score_loss(y, pred_orig)

    optimism_auc.append(auc_boot - auc_orig)
    optimism_brier.append(brier_boot - brier_orig)

mean_opt_auc = np.mean(optimism_auc)
mean_opt_brier = np.mean(optimism_brier)
corrected_auc = apparent_auc - mean_opt_auc
corrected_brier = apparent_brier - mean_opt_brier

print(f"Apparent AUC:              {apparent_auc:.3f}")
print(f"Mean optimism (AUC):       {mean_opt_auc:.4f}")
print(f"Corrected AUC:             {corrected_auc:.3f}")
print(f"\nApparent Brier:            {apparent_brier:.4f}")
print(f"Mean optimism (Brier):     {mean_opt_brier:.4f}")
print(f"Corrected Brier:           {corrected_brier:.4f}")

What the code shows. This implements bootstrap optimism correction, the recommended way to get an honest internal estimate of performance. In R, validate(fit, B = 200) does everything: it reports the apparent statistics, the average optimism, and the optimism-corrected versions. (Note that rms works in Somers’ Dxy, which relates to the C-statistic by C = (Dxy + 1) / 2 — that is why the code converts it.) The Python version spells out the loop: for each of 200 bootstrap resamples it refits the model, measures performance on the resample, then measures the same model on the original data; the gap between the two is the optimism. The numbers to compare are the apparent versus corrected AUC. A large drop (say from 0.85 to 0.72) signals serious overfitting and tells you the apparent figure should never have been reported on its own; a small drop means the model is reasonably stable.

18.8.4 Cross-Validation

Cross-validation is an alternative way to estimate performance on new data. The patients are split into \(k\) roughly equal groups (“folds”; typically 10). Each fold is held out in turn as a test set while the model is built on the other folds, and the held-out results are averaged. Because every patient is tested on a model that did not see them, the average reflects performance on unseen data.

Figure 18.3 shows the idea with \(k = 5\). Each row is one round. The shaded block is the fold held back for testing in that round; the rest of the row is used for fitting. Read down the shaded diagonal and you will see that every patient is held out exactly once, so each patient contributes one prediction made by a model that never saw them.

flowchart TB
    S["Full cohort, split into 5 equal folds"]
    S --> R1["<b>Round 1</b>&nbsp;&nbsp; [TEST] · train · train · train · train"]
    S --> R2["<b>Round 2</b>&nbsp;&nbsp; train · [TEST] · train · train · train"]
    S --> R3["<b>Round 3</b>&nbsp;&nbsp; train · train · [TEST] · train · train"]
    S --> R4["<b>Round 4</b>&nbsp;&nbsp; train · train · train · [TEST] · train"]
    S --> R5["<b>Round 5</b>&nbsp;&nbsp; train · train · train · train · [TEST]"]
    R1 --> A["Average the 5 held-out results<br/>= estimated performance on new patients"]
    R2 --> A
    R3 --> A
    R4 --> A
    R5 --> A
    style S fill:#eef3fb,stroke:#4a6fa5
    style A fill:#e8f1ea,stroke:#4a7a55
    style R1 fill:#fbf3e6,stroke:#c79a3b
    style R2 fill:#fbf3e6,stroke:#c79a3b
    style R3 fill:#fbf3e6,stroke:#c79a3b
    style R4 fill:#fbf3e6,stroke:#c79a3b
    style R5 fill:#fbf3e6,stroke:#c79a3b

Figure 18.3: Five-fold cross-validation. The cohort is split into five equal folds. In each of the five rounds a different fold is held out for testing (shaded) while the model is fitted on the remaining four. Every patient is tested exactly once, on a model that did not see them, and the five held-out results are averaged.

Two practical points. The choice of \(k\) trades off bias against computation: \(k = 10\) is the usual default, while \(k = 5\) is cheaper and \(k = N\) (leave-one-out) is the most expensive and often the noisiest. And as with the bootstrap, every modelling decision must happen inside the loop — if you select predictors on the full dataset and only then cross-validate, the folds have already been contaminated and the estimate goes right back to being optimistic.

Cross-validation is simpler to reason about, but it has drawbacks for clinical prediction models: it does not hand you a single final model, and it gives no clean correction for the calibration slope. For that reason bootstrap optimism correction is generally preferred here.

18.9 External Validation

18.9.1 Types of External Validation

It helps to be clear on the distinction. Internal validation reuses your own development data (via bootstrapping or cross-validation) to estimate how optimistic the apparent performance is — it asks “how well will this model do on patients like mine?” External validation applies the finished model to a genuinely separate dataset — different patients, ideally a different place or time — and asks “does it still work out there?” That is the real test, and the further the new setting is from the original, the more stringent it is. Van Calster et al. (2025) and Smits et al. (2026) distinguish three rungs on this ladder (Figure 18.4):

• Temporal validation: Same setting, different time period. Example: model developed on 2015–2019 data, validated on 2020–2022 data. This captures drift in clinical practice and patient populations over time.
• Geographical validation: Different hospital or region, same time period. This tests whether the model travels to other settings (its “transportability”).
• Domain validation: A different clinical context entirely (e.g., a model developed in a tertiary referral centre, validated in primary care). This is the most demanding test.

A useful refinement for multi-centre data is internal-external cross-validation (IECV), where each centre is held out in turn: the model is built on all the other centres and tested on the one left out. This reuses every centre as both a development and a validation set, and shows how well the model is likely to transport to a new centre.

flowchart TB
    A["Apparent performance<br/>(same data the model was built on - optimistic)"]
    A --> B["Internal validation<br/>(bootstrap / cross-validation:<br/>corrects for optimism)"]
    B --> C["Temporal external validation<br/>(same place, later time)"]
    C --> D["Geographical external validation<br/>(different hospital or region)"]
    D --> E["Domain external validation<br/>(different clinical setting entirely)<br/>STRONGEST EVIDENCE"]

Figure 18.4: The validation ladder. Each rung tests the model on data further from where it was built; the higher you climb, the stronger the evidence that the model will work in practice. Based on the validation hierarchy of Steyerberg (2019).

18.9.2 What to Report in External Validation

Following Van Calster et al. (2025), an external validation study should report at minimum:

1. C-statistic with 95% confidence interval
2. Calibration plot (both grouped and smoothed)
3. O:E ratio with 95% confidence interval
4. Calibration slope
5. Net benefit (decision curve) at clinically relevant thresholds
6. Distribution of predicted probabilities (to understand the range and spread)
7. Performance by subgroups (age, sex, comorbidities, etc.)

18.10 Model Updating

18.10.1 When Performance Degrades

External validation often reveals that a model does not perform as well in a new setting as it did at home. This does not necessarily mean the model is useless. More often it is simply mis-tuned for the new population. Model updating adapts an existing model to that population — a cheaper and more stable option than building a brand-new model from scratch, which would need a large dataset of its own.

18.10.2 Updating Strategies

From least to most complex:

1. Recalibration-in-the-large: Adjust only the model’s baseline term (the intercept), which shifts every predicted risk up or down by the same amount while leaving the ranking of patients untouched. Corrects for differences in baseline risk (e.g., a higher mortality rate in the new population).

2. Logistic recalibration: Adjust the intercept and apply a correction to the calibration slope. This corrects for both level and spread of predictions.

3. Model revision: Re-estimate some or all coefficients, or add new predictors. This requires larger sample sizes in the new population.

R

library(rms)
library(pROC)

set.seed(99)
n_ext <- 1600

# The new setting differs in TWO ways, and it is worth being clear about both.
#
# 1. Case-mix: these patients are older, more severe, less often thrombolysed.
# 2. The outcome relationship itself is different: baseline risk is HIGHER
#    (the +0.6 shift) and the predictors matter LESS than they did at home
#    (the 0.70 multiplier).
#
# Point 2 is the one that actually causes miscalibration. Case-mix alone does
# not: a correctly specified model transports perfectly to a sicker population,
# because the higher predicted risks are then genuinely correct. Only a change
# in the underlying relationship makes the predicted numbers wrong.

ext_data <- data.frame(
  age = round(rnorm(n_ext, 78, 10)), # Older
  nihss = round(pmax(0, rnorm(n_ext, 11, 7))), # Higher severity
  glucose = round(rnorm(n_ext, 155, 55)),
  afib = rbinom(n_ext, 1, 0.35),
  thrombolysis = rbinom(n_ext, 1, 0.20) # Less thrombolysis
)

# What the development model's equation would say for these patients
lp_dev_equation <- -5 +
  0.04 * ext_data$age +
  0.12 * ext_data$nihss +
  0.003 * ext_data$glucose +
  0.3 * ext_data$afib -
  0.5 * ext_data$thrombolysis

# The TRUTH in this new setting: higher baseline, weaker predictor effects
ext_data$death_30d <- rbinom(n_ext, 1, plogis(0.6 + 0.70 * lp_dev_equation))

# Apply the original, unmodified model to the external data
ext_data$pred_original <- predict(fit, newdata = ext_data, type = "fitted")

# Helper: the three calibration numbers plus the Brier score
cal_stats <- function(y, p) {
  co <- coef(glm(y ~ qlogis(p), family = binomial))
  c(
    oe = mean(y) / mean(p), intercept = co[[1]],
    slope = co[[2]], brier = mean((p - y)^2)
  )
}

before <- cal_stats(ext_data$death_30d, ext_data$pred_original)

cat("External validation of the ORIGINAL model:\n")
cat("  Observed mortality:   ", round(mean(ext_data$death_30d), 3), "\n")
cat("  Mean predicted risk:  ", round(mean(ext_data$pred_original), 3), "\n")
cat("  O:E ratio:            ", round(before[["oe"]], 3), "\n")
cat("  Calibration slope:    ", round(before[["slope"]], 3), "\n")
cat("  Calibration intercept:", round(before[["intercept"]], 3), "\n")
cat("  C-statistic:          ", round(
  as.numeric(auc(roc(ext_data$death_30d, ext_data$pred_original,
    quiet = TRUE
  ))), 3
), "\n")

The update itself must be fitted on one set of patients and judged on another. If you fit the recalibration and then evaluate it on the very same rows, the O:E ratio comes back as exactly 1 no matter what — that is forced by the arithmetic of fitting an intercept, not evidence that anything improved. So we split the external cohort in half:

# Split the external cohort: half to fit the update, half to judge it on
set.seed(7)
in_update <- sample(rep(c(TRUE, FALSE), length.out = n_ext))
update_set <- ext_data[in_update, ]
test_set <- ext_data[!in_update, ]

# Logistic recalibration: re-estimate an intercept and a slope for the
# model's own predictions, using ONLY the update half
recal <- glm(death_30d ~ qlogis(pred_original),
  data = update_set, family = binomial
)

cat("Recalibration fitted on the update half:\n")
cat(
  "  new intercept:", round(coef(recal)[[1]], 3),
  "   new slope:", round(coef(recal)[[2]], 3), "\n\n"
)

# Apply the update to the held-out half
test_set$pred_updated <- predict(recal, newdata = test_set, type = "response")

after_orig <- cal_stats(test_set$death_30d, test_set$pred_original)
after_upd <- cal_stats(test_set$death_30d, test_set$pred_updated)

comparison <- data.frame(
  Model = c("Original model", "After recalibration"),
  `O:E ratio` = round(c(after_orig[["oe"]], after_upd[["oe"]]), 3),
  `Calib. slope` = round(c(after_orig[["slope"]], after_upd[["slope"]]), 3),
  `Calib. intercept` = round(
    c(after_orig[["intercept"]], after_upd[["intercept"]]), 3
  ),
  Brier = round(c(after_orig[["brier"]], after_upd[["brier"]]), 4),
  `C-statistic` = round(c(
    as.numeric(auc(roc(test_set$death_30d, test_set$pred_original,
      quiet = TRUE
    ))),
    as.numeric(auc(roc(test_set$death_30d, test_set$pred_updated,
      quiet = TRUE
    )))
  ), 3),
  check.names = FALSE
)

knitr::kable(
  comparison,
  caption = "Performance in the held-out half of the external cohort, before and after logistic recalibration."
)

# Calibration plots before and after. Assigning the result keeps val.prob
# from dumping its full 18-statistic table to the console.
par(mfrow = c(1, 2))
invisible(val.prob(test_set$pred_original, test_set$death_30d,
  m = 60, cex = 0.5
))
title(main = "A. Original model in the new setting")
invisible(val.prob(test_set$pred_updated, test_set$death_30d,
  m = 60, cex = 0.5
))
title(main = "B. After logistic recalibration")
par(mfrow = c(1, 1))

Python

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy.special import expit, logit
from sklearn.calibration import calibration_curve
from sklearn.metrics import roc_auc_score

np.random.seed(99)
n_ext = 1600

# The new setting differs in TWO ways:
#   1. Case-mix: older, more severe, less often thrombolysed.
#   2. The outcome relationship itself: baseline risk is HIGHER (the +0.6
#      shift) and the predictors matter LESS (the 0.70 multiplier).
# Only the second causes miscalibration. A correctly specified model
# transports perfectly to a merely sicker population.
ext_data = {
    "age": np.round(np.random.normal(78, 10, n_ext)),
    "nihss": np.round(np.maximum(0, np.random.normal(11, 7, n_ext))),
    "glucose": np.round(np.random.normal(155, 55, n_ext)),
    "afib": np.random.binomial(1, 0.35, n_ext),
    "thrombolysis": np.random.binomial(1, 0.20, n_ext)
}

# What the development model's equation would say for these patients
lp_dev_equation = (-5 + 0.04 * ext_data["age"] +
                   0.12 * ext_data["nihss"] +
                   0.003 * ext_data["glucose"] +
                   0.3 * ext_data["afib"] -
                   0.5 * ext_data["thrombolysis"])

# The truth in the new setting: higher baseline, weaker predictor effects
y_ext = np.random.binomial(1, expit(0.6 + 0.70 * lp_dev_equation))

# Apply the original, unmodified model
X_ext = np.column_stack([ext_data[k] for k in predictors])
pred_original = full_model.predict_proba(X_ext)[:, 1]


def cal_stats(y, p):
    """O:E ratio, calibration intercept and slope, and Brier score."""
    lp = logit(np.clip(p, 1e-6, 1 - 1e-6))
    res = sm.GLM(y, sm.add_constant(lp),
                 family=sm.families.Binomial()).fit()
    return {
        "oe": y.mean() / p.mean(),
        "intercept": res.params[0],
        "slope": res.params[1],
        "brier": np.mean((p - y) ** 2),
    }


before = cal_stats(y_ext, pred_original)
print("External validation of the ORIGINAL model:")
print(f"  Observed mortality:    {y_ext.mean():.3f}")
print(f"  Mean predicted risk:   {pred_original.mean():.3f}")
print(f"  O:E ratio:             {before['oe']:.3f}")
print(f"  Calibration slope:     {before['slope']:.3f}")
print(f"  Calibration intercept: {before['intercept']:.3f}")
print(f"  C-statistic:           {roc_auc_score(y_ext, pred_original):.3f}")

# ---- Fit the update on one half, judge it on the other ----
# Evaluating a recalibration on the same rows that fitted it forces the
# O:E ratio to exactly 1, which proves nothing.
rng = np.random.default_rng(7)
in_update = rng.permutation(np.arange(n_ext) % 2 == 0)
upd, tst = in_update, ~in_update

lp_original = logit(np.clip(pred_original, 1e-6, 1 - 1e-6))
recal = sm.GLM(y_ext[upd], sm.add_constant(lp_original[upd]),
               family=sm.families.Binomial()).fit()
print(f"\nRecalibration fitted on the update half:")
print(f"  new intercept: {recal.params[0]:.3f}"
      f"   new slope: {recal.params[1]:.3f}")

pred_updated_tst = recal.predict(sm.add_constant(lp_original[tst]))

after_orig = cal_stats(y_ext[tst], pred_original[tst])
after_upd = cal_stats(y_ext[tst], pred_updated_tst)

comparison = pd.DataFrame({
    "Model": ["Original model", "After recalibration"],
    "O:E ratio": [after_orig["oe"], after_upd["oe"]],
    "Calib. slope": [after_orig["slope"], after_upd["slope"]],
    "Calib. intercept": [after_orig["intercept"], after_upd["intercept"]],
    "Brier": [after_orig["brier"], after_upd["brier"]],
    "C-statistic": [roc_auc_score(y_ext[tst], pred_original[tst]),
                    roc_auc_score(y_ext[tst], pred_updated_tst)],
}).round(3)
print("\nHeld-out half of the external cohort:")
print(comparison.to_string(index=False))

# ---- Calibration plots before and after ----
fig, axes = plt.subplots(1, 2, figsize=(12, 5.5), sharex=True, sharey=True)
panels = [
    (pred_original[tst], "A. Original model in the new setting"),
    (pred_updated_tst, "B. After logistic recalibration"),
]
for ax, (p, ttl) in zip(axes, panels):
    obs, exp = calibration_curve(y_ext[tst], p, n_bins=10,
                                 strategy="quantile")
    ax.plot([0, 1], [0, 1], "--", color="grey", label="Perfect calibration")
    ax.plot(exp, obs, "o-", color="steelblue", lw=2, label="Observed")
    ax.set_xlabel("Predicted probability")
    ax.set_title(ttl)
    ax.set_xlim(0, 1)
    ax.set_ylim(0, 1)
    ax.legend(loc="upper left", fontsize=9)
axes[0].set_ylabel("Observed proportion")
plt.tight_layout()
plt.show()

What the code shows. This applies the stroke model to a new setting where it is genuinely mis-tuned, then repairs it and checks the repair on patients that played no part in fitting it.

Read the first block of output as the external validation report. The model’s mean predicted risk comes out near 0.5 while observed mortality is over 0.6, giving an O:E ratio of roughly 1.3 — the model is underestimating risk by about 30%, so it would systematically under-treat this population. The calibration slope of roughly 0.7 says the predictions are also too spread out: the model’s high risks are too high and its low risks too low, relative to what actually happens here. (The exact figures differ a little between the R and Python tabs, since each draws its own random sample.)

The comparison table then shows what a light-touch update buys you, measured on the held-out half:

• The O:E ratio moves from about 1.29 to about 0.99 and the calibration intercept from about 0.62 to about 0.07 — the systematic underestimation is essentially gone.
• The calibration slope moves substantially towards 1 (from roughly 0.7). It does not land exactly on 1.0, and it may overshoot slightly in one language and undershoot in the other, because we are judging the update on patients it never saw. That wobble is honest. A version reporting a perfect slope of exactly 1.000 would be measuring itself on its own training rows.
• The Brier score improves (about 0.240 to 0.215), confirming the predictions are genuinely closer to the observed outcomes.
• The C-statistic does not move at all (0.663 in both rows).

That last point is the one to take away, and the calibration plots make it visible: recalibration slides and stretches the predicted risks, but it never reorders patients. Since the C-statistic depends only on the ranking, it cannot change. So a poor C-statistic in a new setting is not something recalibration can fix — but poor calibration usually is. The two problems have different remedies, which is precisely why the chapter insists on measuring them separately.

In the calibration plots, Panel A shows the original model’s curve sitting well above the diagonal at every risk level — the visual signature of consistent underestimation — while in Panel B the curve tracks the diagonal closely. The clinical message: a model that miscalibrates in a new setting is often not broken, just mis-tuned, and an update needing only an intercept and a slope can rescue it without collecting enough data to build a model from scratch.

ImportantCase-mix differences alone do not cause miscalibration

This is worth stating plainly, because it is a common misconception. If the new population is simply older and sicker but the relationships between predictors and outcome are unchanged, a correctly specified model transports perfectly — it will predict higher risks for these patients, and those higher risks will be correct. Miscalibration appears only when the underlying relationship differs: a different baseline risk, or predictors that carry different weight. In the code above, the + 0.6 shift and the 0.70 multiplier are what actually break the model; the older, sicker covariate distribution on its own would not have.

18.11 Recommended Reporting: The Minimum Set

Based on Van Calster et al. (2025) and Smits et al. (2026), every prediction model evaluation should include at minimum:

1. AUROC (C-statistic) with confidence interval — for discrimination
2. Calibration plot (smoothed and/or grouped) — for calibration
3. Decision curve — for clinical utility
4. Risk distribution plot — to show the spread of predicted probabilities

These four visualisations together give a comprehensive picture. A model may look excellent on one plot and problematic on another — which is exactly why you need all four.

R

library(pROC)

pred_prob <- predict(fit, type = "fitted")
y_obs <- stroke_data$death_30d

par(mfrow = c(2, 2))

# Panel 1: ROC curve
roc_obj <- roc(y_obs, pred_prob, quiet = TRUE)
plot(
  roc_obj,
  col = "steelblue",
  lwd = 2,
  main = "A. Discrimination (ROC Curve)",
  print.auc = TRUE,
  print.auc.y = 0.4
)

# Panel 2: Calibration plot.
# val.prob() accepts no `main` argument, so the title is added afterwards.
# Assigning the result also stops it printing its 18-statistic table.
invisible(val.prob(pred_prob, y_obs, m = 100, cex = 0.5))
title(main = "B. Calibration Plot")

# Panel 3: Risk distribution.
# Scale the y-axis to the data rather than hard-coding a limit, otherwise the
# histogram is squashed into the bottom of the panel.
brk <- seq(0, max(pred_prob) + 0.05, length.out = 31)
h0 <- hist(pred_prob[y_obs == 0], breaks = brk, plot = FALSE)
h1 <- hist(pred_prob[y_obs == 1], breaks = brk, plot = FALSE)
y_top <- max(h0$density, h1$density) * 1.25

hist(
  pred_prob[y_obs == 0],
  breaks = brk,
  col = rgb(0.4, 0.6, 0.8, 0.5),
  main = "C. Risk Distribution",
  xlab = "Predicted Probability",
  freq = FALSE,
  ylim = c(0, y_top)
)
hist(
  pred_prob[y_obs == 1],
  breaks = brk,
  col = rgb(0.9, 0.4, 0.4, 0.5),
  add = TRUE,
  freq = FALSE
)
legend(
  "topright",
  legend = c("No event", "Event"),
  fill = c(rgb(0.4, 0.6, 0.8, 0.5), rgb(0.9, 0.4, 0.4, 0.5))
)

# Panel 4: Simplified decision curve
thresholds <- seq(0.01, 0.5, by = 0.01)
nb_model <- nb_all <- numeric(length(thresholds))
n_total <- length(y_obs)

for (i in seq_along(thresholds)) {
  t <- thresholds[i]
  pred_pos <- as.numeric(pred_prob >= t)
  tp <- sum(pred_pos == 1 & y_obs == 1)
  fp <- sum(pred_pos == 1 & y_obs == 0)
  nb_model[i] <- tp / n_total - fp / n_total * (t / (1 - t))
  nb_all[i] <- mean(y_obs) - (1 - mean(y_obs)) * (t / (1 - t))
}

plot(
  thresholds,
  nb_model,
  type = "l",
  lwd = 2,
  col = "steelblue",
  main = "D. Decision Curve Analysis",
  xlab = "Threshold Probability",
  ylab = "Net Benefit",
  ylim = c(-0.05, max(nb_model, nb_all) * 1.1)
)
lines(thresholds, nb_all, lwd = 2, col = "coral", lty = 2)
abline(h = 0, col = "black")
legend(
  "topright",
  legend = c("Model", "Treat All", "Treat None"),
  col = c("steelblue", "coral", "black"),
  lty = c(1, 2, 1),
  lwd = 2,
  cex = 0.8
)

par(mfrow = c(1, 1))

Python

import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import roc_curve, auc
from sklearn.calibration import calibration_curve

pred_prob = full_model.predict_proba(X)[:, 1]
y_obs = stroke_data["death_30d"].values

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Panel A: ROC curve
fpr, tpr, _ = roc_curve(y_obs, pred_prob)
roc_auc = auc(fpr, tpr)
axes[0, 0].plot(fpr, tpr, color="steelblue", lw=2,
                label=f"AUC = {roc_auc:.3f}")
axes[0, 0].plot([0, 1], [0, 1], "--", color="grey")
axes[0, 0].set_xlabel("1 - Specificity")
axes[0, 0].set_ylabel("Sensitivity")
axes[0, 0].set_title("A. Discrimination (ROC Curve)")
axes[0, 0].legend()

# Panel B: Calibration plot
prob_true, prob_pred = calibration_curve(y_obs, pred_prob, n_bins=10,
                                          strategy="quantile")
axes[0, 1].plot(prob_pred, prob_true, "o-", color="steelblue", lw=2)
axes[0, 1].plot([0, 1], [0, 1], "--", color="grey")
axes[0, 1].set_xlabel("Mean Predicted Probability")
axes[0, 1].set_ylabel("Observed Proportion")
axes[0, 1].set_title("B. Calibration Plot")

# Panel C: Risk distribution
axes[1, 0].hist(pred_prob[y_obs == 0], bins=30, alpha=0.5,
                color="steelblue", label="No event", density=True)
axes[1, 0].hist(pred_prob[y_obs == 1], bins=30, alpha=0.5,
                color="coral", label="Event", density=True)
axes[1, 0].set_xlabel("Predicted Probability")
axes[1, 0].set_ylabel("Density")
axes[1, 0].set_title("C. Risk Distribution")
axes[1, 0].legend()

# Panel D: Decision curve
thresholds = np.arange(0.01, 0.51, 0.01)
nb_model = []
nb_all = []
n_total = len(y_obs)

for t in thresholds:
    pred_pos = (pred_prob >= t).astype(int)
    tp = np.sum((pred_pos == 1) & (y_obs == 1))
    fp = np.sum((pred_pos == 1) & (y_obs == 0))
    nb_model.append(tp / n_total - fp / n_total * (t / (1 - t)))
    nb_all.append(y_obs.mean() - (1 - y_obs.mean()) * (t / (1 - t)))

axes[1, 1].plot(thresholds, nb_model, color="steelblue", lw=2,
                label="Model")
axes[1, 1].plot(thresholds, nb_all, color="coral", lw=2, ls="--",
                label="Treat All")
axes[1, 1].axhline(y=0, color="black", lw=1, label="Treat None")
axes[1, 1].set_xlabel("Threshold Probability")
axes[1, 1].set_ylabel("Net Benefit")
axes[1, 1].set_title("D. Decision Curve Analysis")
axes[1, 1].legend(fontsize=9)

plt.tight_layout()
plt.show()

What the code shows. This assembles the four plots that, together, give the complete picture of a model’s performance and that you should report as a set. Panel A (ROC curve) shows discrimination; B (calibration plot) shows whether predicted risks are accurate; C (risk distribution) shows whether the model produces a useful spread of predictions across patients with and without the event; and D (decision curve) shows clinical utility. The reason to view all four at once is that a model can look excellent on one and alarming on another — strong discrimination with poor calibration, for instance — and only the combined display reveals that. When you read someone else’s model paper, check whether they reported all four; a paper that shows only the AUC has told you only a quarter of the story.

18.12 Exercises

TipExercise 1: Calibration Assessment

Using the stroke mortality model developed in this chapter:

1. Create a calibration plot using deciles of predicted risk. Is the model well-calibrated?
2. Calculate the O:E ratio. What does it tell you?
3. Calculate the calibration slope. Is there evidence of overfitting?
4. Apply bootstrap optimism correction. How much does the C-statistic decrease? How much does the calibration slope decrease?

NoteSolution

R

# Chapter 11, Exercise 1: Calibration Assessment
# Using the stroke mortality model from the chapter

library(rms)

# ---- Simulate the stroke data (same as chapter) ----
set.seed(2024)
n <- 1500

stroke_data <- data.frame(
  age = round(rnorm(n, 72, 12)),
  nihss = round(pmax(0, rnorm(n, 8, 6))),
  glucose = round(rnorm(n, 140, 50)),
  afib = rbinom(n, 1, 0.25),
  thrombolysis = rbinom(n, 1, 0.30)
)

lp <- -5 + 0.04 * stroke_data$age +
  0.12 * stroke_data$nihss +
  0.003 * stroke_data$glucose +
  0.3 * stroke_data$afib -
  0.5 * stroke_data$thrombolysis

stroke_data$death_30d <- rbinom(n, 1, plogis(lp))

# Fit the prediction model
dd <- datadist(stroke_data)
options(datadist = "dd")

fit <- lrm(death_30d ~ age + nihss + glucose + afib + thrombolysis,
           data = stroke_data, x = TRUE, y = TRUE)

pred_prob <- predict(fit, type = "fitted")

# ---- (a) Calibration plot using deciles of predicted risk ----
cat("=== Part (a): Calibration Plot ===\n")
# val.prob() takes no `main` argument, so the title is added with title()
invisible(val.prob(pred_prob, stroke_data$death_30d, m = 150, cex = 0.5))
title(main = "Calibration Plot (Deciles): 30-Day Stroke Mortality")
# The model appears well-calibrated since the points cluster near the diagonal.
# This is expected because we are evaluating apparent performance on the
# training data.

# ---- (b) O:E ratio ----
cat("\n=== Part (b): O:E Ratio ===\n")
obs_rate <- mean(stroke_data$death_30d)
mean_pred <- mean(pred_prob)
oe_ratio <- obs_rate / mean_pred

cat("Observed event rate:", round(obs_rate, 3), "\n")
cat("Mean predicted probability:", round(mean_pred, 3), "\n")
cat("O:E ratio:", round(oe_ratio, 3), "\n")
# An O:E ratio near 1.0 indicates good calibration-in-the-large.
# The model's average predictions match the observed event rate.

# ---- (c) Calibration slope ----
cat("\n=== Part (c): Calibration Slope ===\n")
cal_model <- glm(stroke_data$death_30d ~ qlogis(pred_prob), family = binomial)
cal_slope <- coef(cal_model)[2]
cal_intercept <- coef(cal_model)[1]

cat("Calibration slope:", round(cal_slope, 3), "\n")
cat("Calibration intercept:", round(cal_intercept, 3), "\n")
# A calibration slope of 1.0 indicates no overfitting.
# The apparent slope is typically close to 1 on the training data.
# Values < 1 on new data would suggest overfitting (predictions too extreme).

# ---- (d) Bootstrap optimism correction ----
cat("\n=== Part (d): Bootstrap Optimism Correction ===\n")
set.seed(42)
val <- validate(fit, B = 200)

cat("Bootstrap Validation Results:\n")
print(val)

# Extract optimism-corrected C-statistic
dxy_corrected <- val["Dxy", "index.corrected"]
c_apparent <- fit$stats["C"]
c_corrected <- (dxy_corrected + 1) / 2

cat("\nApparent C-statistic:", round(c_apparent, 4), "\n")
cat("Optimism-corrected C-statistic:", round(c_corrected, 4), "\n")
cat("C-statistic decrease:", round(c_apparent - c_corrected, 4), "\n")

# Calibration slope from bootstrap
slope_corrected <- val["Slope", "index.corrected"]
cat("\nApparent calibration slope: 1.000\n")
cat("Optimism-corrected calibration slope:", round(slope_corrected, 4), "\n")
cat("Slope decrease:", round(1 - slope_corrected, 4), "\n")

# The C-statistic decreases slightly after optimism correction, reflecting
# the mild optimism in apparent performance. The calibration slope also
# decreases below 1, indicating some overfitting that inflates apparent
# performance.

Python

# Chapter 11, Exercise 1: Calibration Assessment
# Using the stroke mortality model from the chapter

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, brier_score_loss
from sklearn.calibration import calibration_curve
from scipy.special import expit, logit

# ---- Simulate the stroke data (same as chapter) ----
np.random.seed(2024)
n = 1500

stroke_data = pd.DataFrame({
    "age": np.round(np.random.normal(72, 12, n)),
    "nihss": np.round(np.maximum(0, np.random.normal(8, 6, n))),
    "glucose": np.round(np.random.normal(140, 50, n)),
    "afib": np.random.binomial(1, 0.25, n),
    "thrombolysis": np.random.binomial(1, 0.30, n)
})

lp = (-5 + 0.04 * stroke_data["age"] +
      0.12 * stroke_data["nihss"] +
      0.003 * stroke_data["glucose"] +
      0.3 * stroke_data["afib"] -
      0.5 * stroke_data["thrombolysis"])

stroke_data["death_30d"] = np.random.binomial(1, expit(lp))

predictors = ["age", "nihss", "glucose", "afib", "thrombolysis"]
X = stroke_data[predictors].values
y = stroke_data["death_30d"].values

model = LogisticRegression(max_iter=5000, random_state=42)
model.fit(X, y)
pred_prob = model.predict_proba(X)[:, 1]

# ---- (a) Calibration plot using deciles ----
print("=== Part (a): Calibration Plot ===")
prob_true, prob_pred = calibration_curve(y, pred_prob, n_bins=10, strategy="quantile")

fig, ax = plt.subplots(figsize=(7, 6))
ax.plot(prob_pred, prob_true, "o-", color="steelblue", lw=2, markersize=8, label="Model")
ax.plot([0, 1], [0, 1], "--", color="grey", label="Perfect calibration")
ax.set_xlabel("Mean Predicted Probability")
ax.set_ylabel("Observed Proportion")
ax.set_title("Calibration Plot (Deciles): 30-Day Stroke Mortality")
ax.legend()
plt.tight_layout()
plt.savefig("ch11_ex1_calibration.png", dpi=150)
plt.show()
# The points cluster near the diagonal, suggesting good apparent calibration.

# ---- (b) O:E ratio ----
print("\n=== Part (b): O:E Ratio ===")
obs_rate = y.mean()
mean_pred = pred_prob.mean()
oe_ratio = obs_rate / mean_pred

print(f"Observed event rate: {obs_rate:.3f}")
print(f"Mean predicted probability: {mean_pred:.3f}")
print(f"O:E ratio: {oe_ratio:.3f}")
# An O:E ratio near 1.0 indicates good calibration-in-the-large.

# ---- (c) Calibration slope ----
print("\n=== Part (c): Calibration Slope ===")
lp_pred = logit(np.clip(pred_prob, 1e-8, 1 - 1e-8))
cal_model = LogisticRegression(max_iter=5000, penalty=None)
cal_model.fit(lp_pred.reshape(-1, 1), y)

print(f"Calibration slope: {cal_model.coef_[0][0]:.3f}")
print(f"Calibration intercept: {cal_model.intercept_[0]:.3f}")
# A slope near 1.0 on training data is expected. Values < 1 on new data
# indicate overfitting.

# ---- (d) Bootstrap optimism correction ----
print("\n=== Part (d): Bootstrap Optimism Correction ===")
np.random.seed(42)

# Apparent performance
apparent_auc = roc_auc_score(y, pred_prob)
apparent_brier = brier_score_loss(y, pred_prob)

n_boot = 200
optimism_auc = []
optimism_brier = []

for b in range(n_boot):
    boot_idx = np.random.choice(n, n, replace=True)
    X_boot, y_boot = X[boot_idx], y[boot_idx]

    boot_model = LogisticRegression(max_iter=5000, random_state=42)
    boot_model.fit(X_boot, y_boot)

    # Apparent on bootstrap sample
    pred_boot = boot_model.predict_proba(X_boot)[:, 1]
    auc_boot = roc_auc_score(y_boot, pred_boot)

    # Test on original data
    pred_orig = boot_model.predict_proba(X)[:, 1]
    auc_orig = roc_auc_score(y, pred_orig)

    optimism_auc.append(auc_boot - auc_orig)

mean_opt_auc = np.mean(optimism_auc)
corrected_auc = apparent_auc - mean_opt_auc

print(f"Apparent C-statistic:          {apparent_auc:.4f}")
print(f"Mean optimism (AUC):           {mean_opt_auc:.4f}")
print(f"Optimism-corrected C-statistic:{corrected_auc:.4f}")
print(f"C-statistic decrease:          {mean_opt_auc:.4f}")

# The C-statistic decreases slightly after correction, reflecting mild
# optimism in the apparent performance. With 1500 observations and 5
# predictors, overfitting is modest.

TipExercise 2: External Validation Simulation

R

# Create three external validation populations that differ from development:
# Population A: Same demographics, 3 years later (temporal)
# Population B: Different hospital, younger patients (geographical)
# Population C: Primary care setting with lower severity (domain)

# For each population:
# a. Calculate C-statistic, O:E ratio, calibration slope
# b. Create calibration plots
# c. Determine which population shows worst calibration and explain why
# d. Perform logistic recalibration and show the improvement

Python

# Same exercise as R tab - create three external populations
# and assess model performance in each

NoteSolution

R

# Chapter 11, Exercise 2: External Validation Simulation
# Create three external validation populations and assess model performance

library(rms)
library(pROC)

# ---- Simulate development data and fit model (same as chapter) ----
set.seed(2024)
n <- 1500

stroke_data <- data.frame(
  age = round(rnorm(n, 72, 12)),
  nihss = round(pmax(0, rnorm(n, 8, 6))),
  glucose = round(rnorm(n, 140, 50)),
  afib = rbinom(n, 1, 0.25),
  thrombolysis = rbinom(n, 1, 0.30)
)

lp <- -5 + 0.04 * stroke_data$age +
  0.12 * stroke_data$nihss +
  0.003 * stroke_data$glucose +
  0.3 * stroke_data$afib -
  0.5 * stroke_data$thrombolysis

stroke_data$death_30d <- rbinom(n, 1, plogis(lp))

dd <- datadist(stroke_data)
options(datadist = "dd")

fit <- lrm(death_30d ~ age + nihss + glucose + afib + thrombolysis,
           data = stroke_data, x = TRUE, y = TRUE)

# ---- Helper function: evaluate model performance ----
evaluate_ext <- function(ext_data, fit, label) {
  ext_data$pred <- predict(fit, newdata = ext_data, type = "fitted")
  obs_rate <- mean(ext_data$death_30d)
  mean_pred <- mean(ext_data$pred)
  oe <- obs_rate / mean_pred

  # C-statistic
  roc_obj <- roc(ext_data$death_30d, ext_data$pred, quiet = TRUE)
  c_stat <- auc(roc_obj)


  # Calibration slope
  lp_ext <- qlogis(ext_data$pred)
  cal_model <- glm(death_30d ~ lp_ext, data = ext_data, family = binomial)
  cal_slope <- coef(cal_model)[2]
  cal_int <- coef(cal_model)[1]

  cat(sprintf("\n=== %s ===\n", label))
  cat(sprintf("  N = %d, Observed mortality: %.3f\n", nrow(ext_data), obs_rate))
  cat(sprintf("  Mean predicted: %.3f\n", mean_pred))
  cat(sprintf("  C-statistic: %.3f\n", c_stat))
  cat(sprintf("  O:E ratio: %.3f\n", oe))
  cat(sprintf("  Calibration slope: %.3f\n", cal_slope))
  cat(sprintf("  Calibration intercept: %.3f\n", cal_int))

  return(ext_data)
}

# ---- Population A: Temporal validation (3 years later) ----
# Same demographics, slightly different practice patterns
set.seed(101)
n_a <- 800
pop_a <- data.frame(
  age = round(rnorm(n_a, 73, 12)),        # Similar age
  nihss = round(pmax(0, rnorm(n_a, 8, 6))),
  glucose = round(rnorm(n_a, 138, 48)),
  afib = rbinom(n_a, 1, 0.27),
  thrombolysis = rbinom(n_a, 1, 0.40)     # More thrombolysis over time
)

lp_a <- -5 + 0.04 * pop_a$age + 0.12 * pop_a$nihss +
  0.003 * pop_a$glucose + 0.3 * pop_a$afib - 0.5 * pop_a$thrombolysis
pop_a$death_30d <- rbinom(n_a, 1, plogis(lp_a))

pop_a <- evaluate_ext(pop_a, fit, "Population A: Temporal (3 years later)")

# ---- Population B: Geographical (younger patients) ----
set.seed(202)
n_b <- 600
pop_b <- data.frame(
  age = round(rnorm(n_b, 62, 10)),        # Younger
  nihss = round(pmax(0, rnorm(n_b, 6, 5))),  # Lower severity
  glucose = round(rnorm(n_b, 130, 45)),
  afib = rbinom(n_b, 1, 0.15),            # Less AF
  thrombolysis = rbinom(n_b, 1, 0.35)
)

lp_b <- -5 + 0.04 * pop_b$age + 0.12 * pop_b$nihss +
  0.003 * pop_b$glucose + 0.3 * pop_b$afib - 0.5 * pop_b$thrombolysis
pop_b$death_30d <- rbinom(n_b, 1, plogis(lp_b))

pop_b <- evaluate_ext(pop_b, fit, "Population B: Geographical (younger patients)")

# ---- Population C: Domain (primary care, lower severity) ----
set.seed(303)
n_c <- 500
pop_c <- data.frame(
  age = round(rnorm(n_c, 68, 14)),
  nihss = round(pmax(0, rnorm(n_c, 3, 3))),  # Much lower severity
  glucose = round(rnorm(n_c, 120, 35)),
  afib = rbinom(n_c, 1, 0.20),
  thrombolysis = rbinom(n_c, 1, 0.10)        # Rarely used in primary care
)

# Different outcome model: lower baseline risk in primary care
lp_c <- -6 + 0.03 * pop_c$age + 0.10 * pop_c$nihss +
  0.002 * pop_c$glucose + 0.2 * pop_c$afib - 0.3 * pop_c$thrombolysis
pop_c$death_30d <- rbinom(n_c, 1, plogis(lp_c))

pop_c <- evaluate_ext(pop_c, fit, "Population C: Domain (primary care)")

# ---- (b) Calibration plots ----
# val.prob() takes no `main` argument, so each title is added with title()
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
invisible(val.prob(pop_a$pred, pop_a$death_30d, m = 80, cex = 0.5))
title(main = "A: Temporal")
invisible(val.prob(pop_b$pred, pop_b$death_30d, m = 60, cex = 0.5))
title(main = "B: Geographical")
invisible(val.prob(pop_c$pred, pop_c$death_30d, m = 50, cex = 0.5))
title(main = "C: Domain")
par(mfrow = c(1, 1))

# ---- (c) Which population shows worst calibration? ----
cat("\n=== Part (c): Worst Calibration ===\n")
cat("Population C (primary care / domain validation) shows the worst calibration.\n")
cat("This is because the outcome model differs from the development setting:\n")
cat("  - Different baseline risk (intercept)\n")
cat("  - Different predictor-outcome relationships (coefficients)\n")
cat("  - Different case-mix (lower severity patients)\n")
cat("Domain validation is the most stringent test of transportability.\n")

# ---- (d) Logistic recalibration ----
cat("\n=== Part (d): Logistic Recalibration ===\n")

recalibrate <- function(ext_data, label) {
  lp_ext <- qlogis(ext_data$pred)
  cal_fit <- glm(death_30d ~ lp_ext, data = ext_data, family = binomial)
  ext_data$pred_recal <- predict(cal_fit, type = "response")

  oe_before <- mean(ext_data$death_30d) / mean(ext_data$pred)
  oe_after  <- mean(ext_data$death_30d) / mean(ext_data$pred_recal)

  cat(sprintf("\n%s:\n", label))
  cat(sprintf("  O:E before recalibration: %.3f\n", oe_before))
  cat(sprintf("  O:E after recalibration:  %.3f\n", oe_after))
}

recalibrate(pop_a, "Population A")
recalibrate(pop_b, "Population B")
recalibrate(pop_c, "Population C")

cat("\nLogistic recalibration corrects for differences in baseline risk\n")
cat("and prediction spread, bringing O:E ratios closer to 1.0.\n")

Python

# Chapter 11, Exercise 2: External Validation Simulation
# Create three external validation populations and assess model performance

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, brier_score_loss
from sklearn.calibration import calibration_curve
from scipy.special import expit, logit

# ---- Simulate development data and fit model (same as chapter) ----
np.random.seed(2024)
n = 1500

stroke_data = pd.DataFrame({
    "age": np.round(np.random.normal(72, 12, n)),
    "nihss": np.round(np.maximum(0, np.random.normal(8, 6, n))),
    "glucose": np.round(np.random.normal(140, 50, n)),
    "afib": np.random.binomial(1, 0.25, n),
    "thrombolysis": np.random.binomial(1, 0.30, n)
})

lp = (-5 + 0.04 * stroke_data["age"] +
      0.12 * stroke_data["nihss"] +
      0.003 * stroke_data["glucose"] +
      0.3 * stroke_data["afib"] -
      0.5 * stroke_data["thrombolysis"])

stroke_data["death_30d"] = np.random.binomial(1, expit(lp))

predictors = ["age", "nihss", "glucose", "afib", "thrombolysis"]
X = stroke_data[predictors].values
y = stroke_data["death_30d"].values

model = LogisticRegression(max_iter=5000, random_state=42)
model.fit(X, y)


# ---- Helper function ----
def evaluate_ext(ext_df, model, predictors, label):
    """Evaluate model on external data and return predictions."""
    X_ext = ext_df[predictors].values
    y_ext = ext_df["death_30d"].values
    pred = model.predict_proba(X_ext)[:, 1]

    obs_rate = y_ext.mean()
    mean_pred = pred.mean()
    oe = obs_rate / mean_pred
    auc = roc_auc_score(y_ext, pred)

    # Calibration slope
    lp_ext = logit(np.clip(pred, 1e-8, 1 - 1e-8))
    cal = LogisticRegression(max_iter=5000, penalty=None)
    cal.fit(lp_ext.reshape(-1, 1), y_ext)

    print(f"\n=== {label} ===")
    print(f"  N = {len(y_ext)}, Observed mortality: {obs_rate:.3f}")
    print(f"  Mean predicted: {mean_pred:.3f}")
    print(f"  C-statistic: {auc:.3f}")
    print(f"  O:E ratio: {oe:.3f}")
    print(f"  Calibration slope: {cal.coef_[0][0]:.3f}")
    print(f"  Calibration intercept: {cal.intercept_[0]:.3f}")

    return pred


# ---- Population A: Temporal validation (3 years later) ----
np.random.seed(101)
n_a = 800
pop_a = pd.DataFrame({
    "age": np.round(np.random.normal(73, 12, n_a)),
    "nihss": np.round(np.maximum(0, np.random.normal(8, 6, n_a))),
    "glucose": np.round(np.random.normal(138, 48, n_a)),
    "afib": np.random.binomial(1, 0.27, n_a),
    "thrombolysis": np.random.binomial(1, 0.40, n_a)  # More thrombolysis
})
lp_a = (-5 + 0.04 * pop_a["age"] + 0.12 * pop_a["nihss"] +
        0.003 * pop_a["glucose"] + 0.3 * pop_a["afib"] -
        0.5 * pop_a["thrombolysis"])
pop_a["death_30d"] = np.random.binomial(1, expit(lp_a))
pred_a = evaluate_ext(pop_a, model, predictors, "Population A: Temporal")

# ---- Population B: Geographical (younger patients) ----
np.random.seed(202)
n_b = 600
pop_b = pd.DataFrame({
    "age": np.round(np.random.normal(62, 10, n_b)),
    "nihss": np.round(np.maximum(0, np.random.normal(6, 5, n_b))),
    "glucose": np.round(np.random.normal(130, 45, n_b)),
    "afib": np.random.binomial(1, 0.15, n_b),
    "thrombolysis": np.random.binomial(1, 0.35, n_b)
})
lp_b = (-5 + 0.04 * pop_b["age"] + 0.12 * pop_b["nihss"] +
        0.003 * pop_b["glucose"] + 0.3 * pop_b["afib"] -
        0.5 * pop_b["thrombolysis"])
pop_b["death_30d"] = np.random.binomial(1, expit(lp_b))
pred_b = evaluate_ext(pop_b, model, predictors, "Population B: Geographical")

# ---- Population C: Domain (primary care, lower severity) ----
np.random.seed(303)
n_c = 500
pop_c = pd.DataFrame({
    "age": np.round(np.random.normal(68, 14, n_c)),
    "nihss": np.round(np.maximum(0, np.random.normal(3, 3, n_c))),
    "glucose": np.round(np.random.normal(120, 35, n_c)),
    "afib": np.random.binomial(1, 0.20, n_c),
    "thrombolysis": np.random.binomial(1, 0.10, n_c)
})
# Different outcome model in primary care
lp_c = (-6 + 0.03 * pop_c["age"] + 0.10 * pop_c["nihss"] +
        0.002 * pop_c["glucose"] + 0.2 * pop_c["afib"] -
        0.3 * pop_c["thrombolysis"])
pop_c["death_30d"] = np.random.binomial(1, expit(lp_c))
pred_c = evaluate_ext(pop_c, model, predictors, "Population C: Domain")

# ---- (b) Calibration plots ----
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
pops = [(pop_a, pred_a, "A: Temporal"), (pop_b, pred_b, "B: Geographical"),
        (pop_c, pred_c, "C: Domain")]

for ax, (pop, pred, title) in zip(axes, pops):
    y_ext = pop["death_30d"].values
    prob_true, prob_pred = calibration_curve(y_ext, pred, n_bins=10,
                                             strategy="quantile")
    ax.plot(prob_pred, prob_true, "o-", color="steelblue", lw=2)
    ax.plot([0, 1], [0, 1], "--", color="grey")
    ax.set_xlabel("Predicted Probability")
    ax.set_ylabel("Observed Proportion")
    ax.set_title(title)

plt.tight_layout()
plt.savefig("ch11_ex2_calibration_plots.png", dpi=150)
plt.show()

# ---- (c) Worst calibration ----
print("\n=== Part (c): Worst Calibration ===")
print("Population C (primary care / domain validation) shows worst calibration.")
print("The outcome model differs from the development setting:")
print("  - Different baseline risk and predictor-outcome relationships")
print("  - Much lower severity patients (NIHSS ~ 3 vs 8)")
print("  - Domain validation is the most stringent test of transportability.")

# ---- (d) Logistic recalibration ----
print("\n=== Part (d): Logistic Recalibration ===")

for pop, pred, label in pops:
    y_ext = pop["death_30d"].values
    lp_ext = logit(np.clip(pred, 1e-8, 1 - 1e-8))
    recal = LogisticRegression(max_iter=5000, penalty=None)
    recal.fit(lp_ext.reshape(-1, 1), y_ext)
    pred_recal = recal.predict_proba(lp_ext.reshape(-1, 1))[:, 1]

    oe_before = y_ext.mean() / pred.mean()
    oe_after = y_ext.mean() / pred_recal.mean()

    print(f"\n{label}:")
    print(f"  O:E before: {oe_before:.3f}")
    print(f"  O:E after:  {oe_after:.3f}")

print("\nLogistic recalibration adjusts intercept and slope, bringing")
print("O:E ratios closer to 1.0 in each external population.")

TipExercise 3: Decision Curve Interpretation

Consider a model for predicting preeclampsia in pregnant women. The model has an AUC of 0.82 and good calibration. You plot the decision curve.

1. At what range of threshold probabilities is the model useful?
2. A colleague argues that the model should not be used because the AUC is “only” 0.82. Using the decision curve, construct a counter-argument.
3. How would the decision curve change if the follow-up action (closer monitoring) were very low-cost versus very high-cost?

NoteSolution

R

# Chapter 11, Exercise 3: Decision Curve Interpretation
# Conceptual exercise about a preeclampsia prediction model

# This exercise is primarily conceptual. We simulate a plausible preeclampsia
# model and generate a decision curve to support the discussion.

set.seed(42)
n <- 2000

# Simulate data for preeclampsia prediction
preeclampsia_data <- data.frame(
  age = round(rnorm(n, 30, 5)),
  bmi = round(rnorm(n, 27, 5), 1),
  nulliparous = rbinom(n, 1, 0.4),
  history_pe = rbinom(n, 1, 0.05),
  map = round(rnorm(n, 85, 10))  # mean arterial pressure
)

lp <- -5.5 + 0.03 * preeclampsia_data$age +
  0.04 * preeclampsia_data$bmi +
  0.5 * preeclampsia_data$nulliparous +
  1.5 * preeclampsia_data$history_pe +
  0.02 * preeclampsia_data$map

preeclampsia_data$pe <- rbinom(n, 1, plogis(lp))
cat("Preeclampsia rate:", mean(preeclampsia_data$pe), "\n")

# Fit model
model <- glm(pe ~ age + bmi + nulliparous + history_pe + map,
             data = preeclampsia_data, family = binomial)

pred_prob <- predict(model, type = "response")
y_obs <- preeclampsia_data$pe

# ---- Generate decision curve ----
thresholds <- seq(0.01, 0.50, by = 0.01)
n_total <- length(y_obs)
nb_model <- nb_all <- numeric(length(thresholds))

for (i in seq_along(thresholds)) {
  t <- thresholds[i]
  pred_pos <- as.numeric(pred_prob >= t)
  tp <- sum(pred_pos == 1 & y_obs == 1)
  fp <- sum(pred_pos == 1 & y_obs == 0)
  nb_model[i] <- tp / n_total - fp / n_total * (t / (1 - t))
  nb_all[i] <- mean(y_obs) - (1 - mean(y_obs)) * (t / (1 - t))
}

plot(thresholds, nb_model, type = "l", lwd = 2, col = "steelblue",
     main = "Decision Curve: Preeclampsia Prediction Model\n(AUC ~ 0.82)",
     xlab = "Threshold Probability", ylab = "Net Benefit",
     ylim = c(-0.05, max(c(nb_model, nb_all)) * 1.1))
lines(thresholds, nb_all, lwd = 2, col = "coral", lty = 2)
abline(h = 0, col = "black")
legend("topright", legend = c("Model", "Treat All", "Treat None"),
       col = c("steelblue", "coral", "black"),
       lty = c(1, 2, 1), lwd = 2)

# ---- Part (a): Range of useful thresholds ----
cat("\n=== Part (a): Range of Useful Thresholds ===\n")
cat("The model is useful at thresholds where its net benefit curve\n")
cat("exceeds BOTH the 'treat all' and 'treat none' lines.\n")
cat("From the decision curve, the model appears useful approximately\n")
cat("in the range of threshold probabilities from about 2% to 30-40%.\n")
cat("Below ~2%, 'treat all' has similar or higher net benefit.\n")
cat("Above ~30-40%, the model's net benefit approaches zero.\n")

# ---- Part (b): Counter-argument to 'AUC is only 0.82' ----
cat("\n=== Part (b): Counter-argument ===\n")
cat("An AUC of 0.82 is not 'only' -- it is strong discrimination.\n")
cat("More importantly, the AUC does not directly tell you whether\n")
cat("using the model improves clinical decisions.\n\n")
cat("The decision curve shows that across the clinically relevant\n")
cat("threshold range (e.g., 5-20%), the model provides positive\n")
cat("net benefit above both default strategies.\n\n")
cat("This means that using the model to guide closer monitoring\n")
cat("decisions would identify more true preeclampsia cases per\n")
cat("'unnecessary' monitoring than either monitoring everyone or\n")
cat("monitoring no one. Clinical utility is what matters for\n")
cat("patient care -- not the AUC alone.\n")

# ---- Part (c): Effect of action cost on the decision curve ----
cat("\n=== Part (c): Effect of Action Cost ===\n")
cat("The threshold probability reflects the implicit cost-benefit\n")
cat("ratio of the action (closer monitoring).\n\n")
cat("LOW-COST action (e.g., closer monitoring, extra visits):\n")
cat("  - Clinicians would use a LOW threshold (e.g., 3-5%)\n")
cat("  - They accept many false positives to avoid missing cases\n")
cat("  - The relevant region of the decision curve shifts LEFT\n")
cat("  - 'Treat all' remains competitive at low thresholds\n\n")
cat("HIGH-COST action (e.g., prophylactic medication with side effects):\n")
cat("  - Clinicians would use a HIGHER threshold (e.g., 15-25%)\n")
cat("  - They want more certainty before acting\n")
cat("  - The relevant region shifts RIGHT\n")
cat("  - The model has more value here because it avoids unnecessary\n")
cat("    treatment of low-risk patients\n\n")
cat("In summary: the decision curve itself does not change, but the\n")
cat("clinically relevant REGION changes depending on the cost of action.\n")

Python

# Chapter 11, Exercise 3: Decision Curve Interpretation
# Conceptual exercise about a preeclampsia prediction model

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from scipy.special import expit

# Simulate a plausible preeclampsia model to generate a decision curve
np.random.seed(42)
n = 2000

age = np.round(np.random.normal(30, 5, n))
bmi = np.round(np.random.normal(27, 5, n), 1)
nulliparous = np.random.binomial(1, 0.4, n)
history_pe = np.random.binomial(1, 0.05, n)
map_val = np.round(np.random.normal(85, 10, n))

lp = (-5.5 + 0.03 * age + 0.04 * bmi + 0.5 * nulliparous +
      1.5 * history_pe + 0.02 * map_val)
pe = np.random.binomial(1, expit(lp))
print(f"Preeclampsia rate: {pe.mean():.3f}")

X = np.column_stack([age, bmi, nulliparous, history_pe, map_val])
model = LogisticRegression(max_iter=5000, random_state=42)
model.fit(X, pe)
pred_prob = model.predict_proba(X)[:, 1]
print(f"AUC: {roc_auc_score(pe, pred_prob):.3f}")

# ---- Generate decision curve ----
thresholds = np.arange(0.01, 0.51, 0.01)
n_total = len(pe)


def net_benefit(y_true, y_pred, threshold):
    pred_pos = (y_pred >= threshold).astype(int)
    tp = np.sum((pred_pos == 1) & (y_true == 1))
    fp = np.sum((pred_pos == 1) & (y_true == 0))
    return tp / n_total - fp / n_total * (threshold / (1 - threshold))


nb_model = [net_benefit(pe, pred_prob, t) for t in thresholds]
nb_all = [pe.mean() - (1 - pe.mean()) * t / (1 - t) for t in thresholds]

plt.figure(figsize=(9, 6))
plt.plot(thresholds, nb_model, color="steelblue", lw=2, label="Prediction Model")
plt.plot(thresholds, nb_all, color="coral", lw=2, ls="--", label="Treat All")
plt.axhline(y=0, color="black", lw=1, label="Treat None")
plt.xlabel("Threshold Probability")
plt.ylabel("Net Benefit")
plt.title("Decision Curve: Preeclampsia Prediction Model (AUC ~ 0.82)")
plt.legend()
plt.ylim(-0.05, max(max(nb_model), max(nb_all)) * 1.1)
plt.tight_layout()
plt.savefig("ch11_ex3_decision_curve.png", dpi=150)
plt.show()

# ---- Part (a): Range of useful thresholds ----
print("\n=== Part (a): Range of Useful Thresholds ===")
print("The model is useful at thresholds where its net benefit curve")
print("exceeds BOTH 'treat all' and 'treat none'.")
print("From the decision curve, the model is useful approximately")
print("from about 2% to 30-40% threshold probability.")
print("Below ~2%, 'treat all' has similar or higher net benefit.")
print("Above ~30-40%, the model's net benefit approaches zero.")

# ---- Part (b): Counter-argument to 'AUC is only 0.82' ----
print("\n=== Part (b): Counter-argument ===")
print("An AUC of 0.82 is strong discrimination, not 'only'.")
print("More importantly, AUC does not directly measure clinical utility.")
print("")
print("The decision curve shows the model provides positive net benefit")
print("above both default strategies across the clinically relevant")
print("threshold range (e.g., 5-20%). This means using the model to")
print("guide monitoring decisions identifies more true preeclampsia")
print("cases per 'unnecessary' monitoring episode than either monitoring")
print("everyone or no one. Clinical utility -- not AUC -- determines")
print("whether a model should be used in practice.")

# ---- Part (c): Effect of action cost ----
print("\n=== Part (c): Effect of Action Cost ===")
print("The threshold probability reflects the cost-benefit ratio of acting.")
print("")
print("LOW-COST action (closer monitoring, extra visits):")
print("  - Clinicians use a LOW threshold (e.g., 3-5%)")
print("  - Many false positives are tolerable")
print("  - The relevant region of the decision curve is on the LEFT")
print("  - 'Treat all' remains competitive at low thresholds")
print("")
print("HIGH-COST action (prophylactic medication with side effects):")
print("  - Clinicians use a HIGHER threshold (e.g., 15-25%)")
print("  - More certainty is needed before acting")
print("  - The relevant region shifts RIGHT")
print("  - The model adds more value by avoiding unnecessary treatment")
print("")
print("The curve itself does not change, but the clinically relevant")
print("region changes depending on the cost of the follow-up action.")

TipExercise 4: Complete Evaluation

Choose a clinical prediction model from the literature (e.g., QRISK3, Wells score, APACHE II). Using published validation data:

1. Report the C-statistic with confidence interval.
2. Describe the calibration (if a calibration plot is available).
3. Is there a decision curve analysis? If not, what threshold range would be clinically relevant?
4. Has the model been externally validated? In what populations?
5. Based on your assessment, would you recommend implementing this model in your setting?

NoteSolution

R

# Chapter 11, Exercise 4: Complete Evaluation of a Published Model
# Example: QRISK3 Cardiovascular Risk Prediction Model
#
# This exercise is primarily a literature review exercise. Below we provide
# a structured answer based on published validation data for QRISK3.

# ---- Part (a): C-statistic with confidence interval ----
cat("=== Part (a): C-statistic ===\n")
cat("QRISK3 (Hippisley-Cox et al., 2017, BMJ) was developed to predict\n")
cat("10-year cardiovascular disease risk.\n\n")
cat("Development cohort:\n")
cat("  C-statistic (women): 0.880 (95% CI: 0.878-0.882)\n")
cat("  C-statistic (men):   0.858 (95% CI: 0.856-0.860)\n\n")
cat("Validation cohort (separate 25% held out):\n")
cat("  C-statistic (women): 0.879 (95% CI: 0.876-0.882)\n")
cat("  C-statistic (men):   0.858 (95% CI: 0.855-0.861)\n")
cat("The discrimination is excellent and stable between development\n")
cat("and validation sets.\n")

# ---- Part (b): Calibration ----
cat("\n=== Part (b): Calibration ===\n")
cat("QRISK3 provides calibration plots in the original publication.\n")
cat("In the validation cohort:\n")
cat("  - The calibration was generally good, with predicted risks\n")
cat("    closely matching observed event rates across deciles.\n")
cat("  - Some overestimation of risk was observed at higher risk\n")
cat("    levels, particularly in older age groups.\n")
cat("  - External validations in different populations (e.g., outside\n")
cat("    the UK) have shown variable calibration, often with\n")
cat("    overestimation in lower-risk populations.\n")

# ---- Part (c): Decision curve analysis ----
cat("\n=== Part (c): Decision Curve Analysis ===\n")
cat("The original QRISK3 paper does not include decision curve analysis.\n\n")
cat("Clinically relevant threshold range for statin initiation:\n")
cat("  - UK NICE guidelines: 10% 10-year CVD risk threshold\n")
cat("  - US ACC/AHA guidelines: 7.5% threshold\n")
cat("  - A decision curve would be most informative in the 5-20%\n")
cat("    threshold range, where the decision to start statins is\n")
cat("    most uncertain.\n")
cat("  - Below 5%, most clinicians would not start statins\n")
cat("  - Above 20%, most clinicians would start statins regardless\n")

# ---- Part (d): External validation ----
cat("\n=== Part (d): External Validation ===\n")
cat("QRISK3 has been validated in multiple populations:\n")
cat("  - Internal-external: UK CPRD data (separate time period)\n")
cat("  - Geographic: Various European populations\n")
cat("  - Ethnic subgroups: South Asian, Black, Chinese populations\n")
cat("  - Temporal: Validated across different calendar periods\n\n")
cat("Key findings from external validations:\n")
cat("  - Discrimination generally remains good (C > 0.80)\n")
cat("  - Calibration can deteriorate in non-UK populations\n")
cat("  - May overestimate risk in populations with lower baseline\n")
cat("    cardiovascular event rates\n")

# ---- Part (e): Recommendation for implementation ----
cat("\n=== Part (e): Recommendation ===\n")
cat("Recommendation depends on the clinical setting:\n\n")
cat("FOR a UK primary care setting:\n")
cat("  - YES, recommend implementation. QRISK3 was developed and\n")
cat("    validated in UK primary care, has excellent discrimination,\n")
cat("    good calibration, and is already integrated into UK\n")
cat("    clinical guidelines.\n\n")
cat("FOR a non-UK setting:\n")
cat("  - CONDITIONAL recommendation. Would first require:\n")
cat("    1. External validation in the local population\n")
cat("    2. Assessment of calibration (likely needs recalibration)\n")
cat("    3. Decision curve analysis at locally relevant thresholds\n")
cat("    4. Comparison with locally developed risk scores\n")
cat("  - If calibration is poor, logistic recalibration should be\n")
cat("    considered before implementation.\n")

Python

# Chapter 11, Exercise 4: Complete Evaluation of a Published Model
# Example: QRISK3 Cardiovascular Risk Prediction Model
#
# This exercise is primarily a literature review exercise. Below we provide
# a structured answer based on published validation data for QRISK3.

# ---- Part (a): C-statistic with confidence interval ----
print("=== Part (a): C-statistic ===")
print("QRISK3 (Hippisley-Cox et al., 2017, BMJ) was developed to predict")
print("10-year cardiovascular disease risk.")
print()
print("Development cohort:")
print("  C-statistic (women): 0.880 (95% CI: 0.878-0.882)")
print("  C-statistic (men):   0.858 (95% CI: 0.856-0.860)")
print()
print("Validation cohort (separate 25% held out):")
print("  C-statistic (women): 0.879 (95% CI: 0.876-0.882)")
print("  C-statistic (men):   0.858 (95% CI: 0.855-0.861)")
print("The discrimination is excellent and stable between development")
print("and validation sets.")

# ---- Part (b): Calibration ----
print("\n=== Part (b): Calibration ===")
print("QRISK3 provides calibration plots in the original publication.")
print("In the validation cohort:")
print("  - Calibration was generally good, with predicted risks closely")
print("    matching observed event rates across deciles.")
print("  - Some overestimation at higher risk levels, particularly in")
print("    older age groups.")
print("  - External validations outside the UK have shown variable")
print("    calibration, often with overestimation in lower-risk populations.")

# ---- Part (c): Decision curve analysis ----
print("\n=== Part (c): Decision Curve Analysis ===")
print("The original QRISK3 paper does not include decision curve analysis.")
print()
print("Clinically relevant threshold range for statin initiation:")
print("  - UK NICE guidelines: 10% 10-year CVD risk threshold")
print("  - US ACC/AHA guidelines: 7.5% threshold")
print("  - A decision curve would be most informative in the 5-20%")
print("    threshold range, where the statin decision is most uncertain.")
print("  - Below 5%, most clinicians would not start statins")
print("  - Above 20%, most clinicians would start statins regardless")

# ---- Part (d): External validation ----
print("\n=== Part (d): External Validation ===")
print("QRISK3 has been validated in multiple populations:")
print("  - Internal-external: UK CPRD data (separate time period)")
print("  - Geographic: Various European populations")
print("  - Ethnic subgroups: South Asian, Black, Chinese populations")
print("  - Temporal: Validated across different calendar periods")
print()
print("Key findings from external validations:")
print("  - Discrimination generally remains good (C > 0.80)")
print("  - Calibration can deteriorate in non-UK populations")
print("  - May overestimate risk in populations with lower baseline")
print("    cardiovascular event rates")

# ---- Part (e): Recommendation ----
print("\n=== Part (e): Recommendation ===")
print("Recommendation depends on the clinical setting:")
print()
print("FOR a UK primary care setting:")
print("  YES, recommend implementation. QRISK3 was developed and")
print("  validated in UK primary care, has excellent discrimination,")
print("  good calibration, and is integrated into UK clinical guidelines.")
print()
print("FOR a non-UK setting:")
print("  CONDITIONAL recommendation. Would first require:")
print("    1. External validation in the local population")
print("    2. Assessment of calibration (likely needs recalibration)")
print("    3. Decision curve analysis at locally relevant thresholds")
print("    4. Comparison with locally developed risk scores")
print("  If calibration is poor, logistic recalibration should be")
print("  considered before implementation.")

18.13 References and Further Reading

• For model validation and calibration, see Smits et al. (2026), Van Calster et al. (2016), and Steyerberg et al. (2010).
• For decision curve analysis, see Vickers and Elkin (2006).
• For scoring rules, see Brier (1950) (formalisation of the Brier score).

Brier, Glenn W. 1950. “Verification of Forecasts Expressed in Terms of Probability.” Monthly Weather Review 78 (1): 1–3. https://doi.org/10.1175/1520-0493(1950)078<0001:VOFEIT>2.0.CO;2. Formalised the Brier score and its decomposition.

Harrell, Frank E. 2015. Regression Modeling Strategies: With Applications to Linear Models, Logistic and Ordinal Regression, and Survival Analysis. 2nd ed. Springer.

Smits, Luc J M, Sander M J van Kuijk, and Laure Wynants. 2026. Improving Health Care with Clinical Prediction Models: From Idea to Impact. Maastricht University Press. https://doi.org/10.26481/mup.2603. CC BY 4.0

Steyerberg, Ewout W. 2019. Clinical Prediction Models: A Practical Approach to Development, Validation, and Updating. 2nd ed. Springer.

Steyerberg, Ewout W, Andrew J Vickers, Nancy R Cook, et al. 2010. “Assessing the Performance of Prediction Models: A Framework for Traditional and Novel Measures.” Epidemiology 21 (1): 128–38. https://doi.org/10.1097/EDE.0b013e3181c30fb2. Application of the Brier score and other performance measures to clinical prediction models.

Van Calster, Ben, Gary S Collins, Andrew J Vickers, et al. 2025. “Evaluation of Performance Measures in Predictive Artificial Intelligence Models to Support Medical Decisions: Overview and Guidance.” The Lancet Digital Health 7: e100916. https://doi.org/10.1016/j.landig.2025.100916.

Van Calster, Ben, Daan Nieboer, Yvonne Vergouwe, Bavo De Cock, Michael J Pencina, and Ewout W Steyerberg. 2016. “A Calibration Hierarchy for Risk Models Was Defined: From Utopia to Empirical Data.” Journal of Clinical Epidemiology 74: 167–76. https://doi.org/10.1016/j.jclinepi.2015.12.005. An excellent treatment of different calibration measures and their interpretation.

Vickers, Andrew J, and Elena B Elkin. 2006. “Decision Curve Analysis: A Novel Method for Evaluating Prediction Models.” Medical Decision Making 26 (6): 565–74.