flowchart TD
A[Define the clinical question] --> B[Design study & sample size]
B --> C[Handle missing data]
C --> D[Select predictors]
D --> E[Specify & fit the model]
E --> F[Assess: discrimination & calibration]
F --> G[Validate internally]
G --> H[Report transparently]
F -. revise .-> E
17 Developing Clinical Prediction Models: A Complete Workflow
17.1 Introduction
A clinical prediction model estimates a patient’s risk of an outcome (for example, the chance of a heart attack in the next 10 years) from information you can collect at the bedside. Building one well is not simply a matter of fitting a regression or running a machine learning algorithm. It is a structured process that begins with a clearly defined clinical question and proceeds through study design, data preparation, model specification, and evaluation. At every step, decisions must be guided by clinical knowledge, not just statistical convenience.
This chapter walks you through the complete workflow of developing a clinical prediction model, following the framework described by Smits et al. (2026) and the foundational principles laid out by Steyerberg (2019). We will use the Framingham Heart Study data as a running example to illustrate each step, building a model to predict 10-year cardiovascular disease (CVD) risk.
Figure 17.1 shows the whole journey at a glance. Each step in this chapter corresponds to one box.
17.2 Step 1: Define the Clinical Question
Before touching any data, you must answer four questions:
- Who is the target population? (e.g., adults aged 40–70 without prior CVD)
- What is the outcome? (e.g., first CVD event within 10 years)
- When is the prediction made? (e.g., at a routine primary care visit)
- Why is the prediction needed? (e.g., to guide statin therapy decisions)
Why a clinician should care. The intended use determines everything that follows. A model designed to triage emergency department patients needs different predictors, a different time horizon, and different performance requirements than a model for long-term cardiovascular risk stratification. Getting the question wrong here cannot be fixed later by clever statistics.
TipThe PICOTS Framework
You can adapt the familiar PICOTS framework from clinical research: Population, Intended use (Index model), Comparators (existing tools), Outcome, Timing, Setting.
17.2.1 Common Pitfalls in Problem Definition
- Vague outcomes: “clinical deterioration” is not a well-defined outcome. Define exactly what events count, how they are ascertained, and over what time window.
- Leaky predictors: including information that would not be available at the moment of prediction. If you are predicting ICU admission at the time of hospital arrival, you cannot use ICU vital signs as predictors.
- Wrong population: developing a model in hospitalised patients and hoping to apply it in primary care. The case-mix — the blend of ages, severities and comorbidities that makes up the patient group — will be entirely different, and a model tuned to one blend can behave quite differently in another.
17.3 Step 2: Study Design and Sample Size
17.3.1 Study Design
Clinical prediction models are typically developed from cohort studies (prospective or retrospective) or existing databases (electronic health records, registries). Randomised trials can also be used, though their selected populations may limit generalisability.
Key design requirements:
- Inception cohort: all patients should be enrolled at a similar, well-defined point in their clinical trajectory (the “moment of prediction”).
- Complete follow-up: outcomes must be ascertained for as many patients as possible. Differential loss to follow-up can bias estimates.
- Outcome rate: there must be a sufficient number of events to support the number of candidate predictors.
17.3.2 Sample Size: The Events Per Variable Rule and Beyond
Why a clinician should care. A model built on too few patients will look impressive in the paper that describes it and then fail when used on your own patients. Working out the sample size in advance is how you avoid that disappointment.
The classical rule of thumb requires at least 10 events per variable (EPV) in the model. “Events per variable” means the number of outcome events (here, CVD events) divided by the number of predictors you are estimating. So for a logistic regression with 8 predictors, you would need at least 80 events. With 10% prevalence, that means at least 800 patients. The intuition: each predictor needs a certain amount of “evidence” (events) to be estimated reliably, and it is the rarer of the two groups (the events) that does the heavy lifting.
However, this rule is a rough guide: it asks for the same number of events whether the outcome is common or rare and whether the predictors are strong or weak. Riley and colleagues (Riley et al. 2020) replaced it with a framework that requires the sample size to satisfy three criteria at once when the outcome is binary. Each one guards against a different failure:
- The coefficients should not need much shrinking. A model fitted to too little data exaggerates its predictor effects, so honest predictions require multiplying every coefficient by a shrinkage factor below 1 (discussed later in this chapter). The criterion asks for the expected shrinkage factor to be at least 0.9 — the effects need scaling back by no more than about 10%. This is usually the criterion that drives the answer.
- Apparent performance should not be much better than honest performance. Judge a model on the same patients it was fitted to and it always looks better than it really is; that gap is called optimism. The criterion puts a ceiling on it: the model’s flattering self-assessment and its honest one must differ by no more than 0.05 on a measure of explained variation (the Nagelkerke R-squared, explained in the box below).
- The average risk should be estimated precisely. Every individual prediction starts from the overall event rate in the population and is nudged up or down from there by the patient’s predictors. If that starting point is imprecise, every prediction inherits the imprecision. The criterion asks for the overall risk to be pinned down to within about ±5 percentage points.
(For a continuous outcome the framework adds a fourth criterion, about how much variation the model leaves unexplained, which is why you may see it described as having four parts.)
NoteWhat is a Cox-Snell R-squared, and how do I come up with one?
Criteria 1 and 2 need an advance guess at how well the model will perform, expressed as a Cox-Snell R-squared (\(R^2_{CS}\)). It is not an effect measure for a single predictor, so it is a different kind of quantity from a hazard or risk ratio: an HR says how much one predictor moves the outcome, while \(R^2_{CS}\) summarises how much of the outcome the whole model accounts for.
In ordinary linear regression — where the outcome is a measurement like blood pressure — R-squared has a comfortable meaning. Patients differ from one another, and R-squared is the share of that variation the model accounts for. A yes/no outcome does not work like that: each patient either had the event or did not, so there is no range of values to account for and the familiar definition has nothing to latch onto.
The replacement is built around a different question: how much better does the model predict what actually happened than a model that knows nothing except the overall event rate? A good model should assign a high probability to what really did occur — high risk to the patients who went on to have the event, low risk to those who did not. \(R^2_{CS}\) scores how well it manages that. Zero means it does no better than simply quoting the overall rate to everybody, and higher is better.
One awkward property matters in practice: the maximum is not 1. With a yes/no outcome the scale gets squashed, and how far depends on how common the outcome is. At 15% prevalence, even a hypothetical model that called every patient’s outcome correctly and with total certainty would only score about 0.57. So read 0.10 against a ceiling of 0.57, not against 1 — which makes it a respectable model rather than a hopeless one. The Nagelkerke R-squared simply divides by that ceiling so the number does run from 0 to 1 (\(0.10 / 0.57 \approx 0.18\)). The sample-size formulas use Cox-Snell; Nagelkerke is the friendlier one to report.
You rarely need to guess \(R^2_{CS}\) yourself. Take a C-statistic from a published model in a similar population (almost always reported; 0.70–0.80 is typical) and let pmsampsize convert it via the cstatistic argument, as below. Be conservative — assuming a better model than you will get makes the required sample size look smaller than it is.
pmsampsize is an R package that does these calculations for you; the name is short for “prediction model sample size”. It reports the largest sample size demanded by the three criteria.
Code
library(pmsampsize)
# Sample size for a logistic regression prediction model
# 10 candidate predictors, anticipated outcome prevalence of 15%,
# and an expected C-statistic of 0.75 taken from a published model.
# pmsampsize converts the C-statistic into a Cox-Snell R-squared for us.
ss <- pmsampsize(
type = "b", # binary outcome
cstatistic = 0.75, # anticipated C-statistic
parameters = 10, # number of candidate predictors
prevalence = 0.15
) # anticipated prevalence
print(ss) # one row per criterion, plus the final (largest) requirement
cat("\nMinimum sample size:", ss$sample_size, "\n")
cat("Minimum number of events:", ceiling(ss$events), "\n")
# If you already have a Cox-Snell R-squared in mind, supply it directly
# with csrsquared = (or nagrsquared = for a Nagelkerke R-squared).
# Note: exactly one of cstatistic, csrsquared or nagrsquared must be given;
# the older `rsquared =` argument is no longer accepted on its own.
pmsampsize(type = "b", csrsquared = 0.10, parameters = 10, prevalence = 0.15)Code
# Riley et al. sample size criteria for a binary outcome.
# There is no Python equivalent of pmsampsize, so the three criteria
# are implemented directly here. Unlike the R package, this version
# does not convert a C-statistic into a Cox-Snell R-squared for you.
import numpy as np
def pmsampsize_binary(cs_rsquared, parameters, prevalence,
shrinkage=0.9, delta=0.05, moe=0.05):
"""
Minimum sample size for developing a binary-outcome prediction model.
cs_rsquared : anticipated Cox-Snell R-squared
parameters : number of predictor parameters to be estimated
prevalence : anticipated outcome prevalence
shrinkage : target shrinkage factor for criterion 1 (default 0.9)
delta : acceptable apparent-vs-adjusted R-squared gap (criterion 2)
moe : margin of error for the overall risk (criterion 3)
"""
phi = prevalence
# The maximum attainable Cox-Snell R-squared, which depends only on
# prevalence. Nagelkerke R-squared rescales Cox-Snell by this ceiling.
max_r2 = 1 - np.exp(2 * (phi * np.log(phi) + (1 - phi) * np.log(1 - phi)))
# Criterion 1: expected shrinkage factor of at least `shrinkage`
n1 = parameters / ((shrinkage - 1) * np.log(1 - cs_rsquared / shrinkage))
# Criterion 2: apparent and adjusted Nagelkerke R-squared differ by <= delta
s2 = cs_rsquared / (cs_rsquared + delta * max_r2)
n2 = parameters / ((s2 - 1) * np.log(1 - cs_rsquared / s2))
# Criterion 3: overall risk estimated to within +/- moe
n3 = (1.96 / moe) ** 2 * phi * (1 - phi)
n_min = int(np.ceil(max(n1, n2, n3)))
events = int(np.ceil(n_min * phi))
print(f"Max possible Cox-Snell R-squared at {phi:.0%} prevalence: {max_r2:.3f}")
print(f"Equivalent Nagelkerke R-squared: {cs_rsquared / max_r2:.3f}\n")
print(f"Criterion 1 (shrinkage >= {shrinkage}): n = {int(np.ceil(n1))}")
print(f"Criterion 2 (R-squared gap <= {delta}): n = {int(np.ceil(n2))}")
print(f"Criterion 3 (risk +/- {moe}): n = {int(np.ceil(n3))}")
print(f"\nMinimum sample size: {n_min}")
print(f"Minimum events: {events}")
print(f"Events per parameter: {events / parameters:.2f}")
return n_min
# A C-statistic of 0.75 at 15% prevalence corresponds to a
# Cox-Snell R-squared of roughly 0.10
pmsampsize_binary(cs_rsquared=0.10, parameters=10, prevalence=0.15)What the code shows. This runs a sample-size calculation before any data are collected. You supply the expected outcome prevalence (15%), the number of candidate predictors (10), and an anticipated performance (a C-statistic of 0.75, which pmsampsize converts to a Cox-Snell R-squared of 0.1028). The output gives one row per criterion and a final row that is simply the largest of them: criterion 1 asks for 825 patients, criterion 2 for 327 and criterion 3 for 196, so 825 patients — about 124 events, or 12.4 events per parameter — is the requirement. Criterion 1 is the one that decides the answer here, as it usually is, and that answer sits above the 100 events the old 10-EPV rule would have accepted. The Python implementation uses a rounded Cox-Snell R-squared of 0.10 rather than 0.1028 and so returns 850; the small difference is just the input, not a disagreement between the two.
The point of all this is feasibility. An underpowered study produces a model that looks impressive in development and falls apart in validation, and it is far cheaper to discover that now than after recruitment. If the required sample size is well beyond what you can realistically recruit, the levers are to reduce the number of candidate predictors, find a larger or additional data source, or reconsider the study — not to proceed and hope.
NoteWhat is EPV?
Events per variable (EPV) is the number of outcome events divided by the number of parameters you are estimating. If you have 80 events and 8 predictors, you have an EPV of 10. Low EPV is the single most common cause of overfitting in clinical models: with too few events, the model “memorises” the development sample rather than learning generalisable patterns.
17.4 Step 3: Handle Missing Data
Missing data is the norm in clinical research, not the exception. Laboratory values are missing because they were not ordered. Follow-up data is missing because patients moved. Lifestyle variables are missing because patients did not answer the questionnaire.
Why a clinician should care. How you deal with missing values can change which predictors look important and how trustworthy the model’s risk estimates are. The wrong choice (usually the default one) can quietly bias the whole model.
17.4.1 Missing Data Mechanisms
Understanding why data are missing is essential for choosing the right approach. The three standard patterns differ in one respect only: what determines whether a value ends up missing? It helps to walk through the same variable — a missing HbA1c measurement — under each.
- MCAR (Missing Completely At Random): missingness has nothing to do with the patient at all — not with anything you recorded, and not with the missing value itself. Example: the laboratory analyser broke down one morning, so that day’s samples were lost. Which patients have a missing HbA1c is effectively a coin flip. This is the benign case, and it is rare in practice.
- MAR (Missing At Random): missingness depends on things you did record, but not on the missing value itself. Example: HbA1c is mostly measured in patients already known to have diabetes, so missingness depends on recorded diabetes status. Among patients with the same recorded diabetes status, who has a missing HbA1c is again essentially a coin flip. The name is unhelpful: “at random” here means random once you condition on the observed data, not haphazard.
- MNAR (Missing Not At Random): missingness depends on the unobserved value itself, and no observed variable explains it away. Example: the patients too unwell to attend their appointment are also the patients with the worst glycaemic control — so it is precisely the highest HbA1c values that are missing, and nothing in your dataset tells you that.
The distinction matters because it decides what you can fix. Under MCAR and MAR, the observed data contain enough information to fill in the gaps sensibly, which is exactly what multiple imputation exploits. Under MNAR the information simply is not there: no imputation method can recover it from the data alone, and the honest response is a sensitivity analysis that asks how different the conclusions would be if the missing values were systematically worse (or better) than the observed ones. In practice you cannot prove which mechanism you are facing — MCAR versus MAR is partly testable, but MAR versus MNAR is not — so the mechanism is an assumption you argue for from clinical knowledge of how the data were collected, and then state explicitly.
17.4.2 Why Complete Case Analysis Fails
Deleting rows with any missing data (complete case analysis) is the default in most software, and it is almost always wrong:
- It wastes data: if 10 variables each have 5% missing (independently), only 60% of rows are complete.
- It biases estimates: if missingness is related to the outcome or predictors (MAR or MNAR), the complete cases are a non-random subset.
- It reduces power: smaller effective sample sizes lead to wider confidence intervals and less stable models.
17.4.3 Multiple Imputation
Multiple imputation (MI) is the recommended approach for handling missing data under the MAR assumption. “Imputation” simply means filling in a missing value with an informed guess based on the patient’s other, observed values. “Multiple” means doing this several times so the uncertainty in those guesses is not hidden. The procedure has three steps:
- Impute: create \(m\) complete datasets, each one filling in the missing values with plausible guesses informed by the observed data. The guesses differ slightly between datasets, reflecting our genuine uncertainty.
- Analyse: fit the prediction model separately in each of the \(m\) filled-in datasets.
- Pool: combine the \(m\) sets of results using Rubin’s rules, a standard recipe that widens the confidence intervals to account for the fact that the missing values were guessed, not measured.
Typically, \(m = 20\) or more imputations are recommended.
flowchart LR
A[Dataset with<br/>missing values] --> B[Impute:<br/>m filled-in datasets]
B --> C[Analyse:<br/>fit model in each]
C --> D[Pool:<br/>combine with Rubin's rules]
Code
library(mice)
library(dplyr)
# Simulate clinical data with missing values
set.seed(42)
n <- 500
data <- data.frame(
age = rnorm(n, 55, 10),
sbp = rnorm(n, 130, 20),
cholesterol = rnorm(n, 220, 40),
smoking = rbinom(n, 1, 0.25),
bmi = rnorm(n, 27, 5)
)
data$cvd <- rbinom(
n,
1,
plogis(
-4 +
0.03 * data$age +
0.01 * data$sbp +
0.005 * data$cholesterol +
0.5 * data$smoking
)
)
# Introduce MAR missingness
# Cholesterol more likely missing in younger patients
data$cholesterol[data$age < 50 & runif(n) < 0.3] <- NA
# BMI more likely missing in non-smokers
data$bmi[data$smoking == 0 & runif(n) < 0.2] <- NA
cat("Missing data pattern:\n")
md.pattern(data, rotate.names = TRUE)
# Perform multiple imputation
imp <- mice(data, m = 20, method = "pmm", seed = 42, printFlag = FALSE)
# Fit model in each imputed dataset and pool
fit_imp <- with(
imp,
glm(cvd ~ age + sbp + cholesterol + smoking + bmi, family = binomial)
)
pooled <- pool(fit_imp)
summary(pooled)Code
import numpy as np
import pandas as pd
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer
from sklearn.linear_model import LogisticRegression
from scipy.special import expit
np.random.seed(42)
n = 500
data = pd.DataFrame({
'age': np.random.normal(55, 10, n),
'sbp': np.random.normal(130, 20, n),
'cholesterol': np.random.normal(220, 40, n),
'smoking': np.random.binomial(1, 0.25, n),
'bmi': np.random.normal(27, 5, n)
})
lp = -4 + 0.03*data['age'] + 0.01*data['sbp'] + 0.005*data['cholesterol'] + 0.5*data['smoking']
data['cvd'] = np.random.binomial(1, expit(lp))
# Introduce MAR missingness
mask_chol = (data['age'] < 50) & (np.random.uniform(size=n) < 0.3)
data.loc[mask_chol, 'cholesterol'] = np.nan
mask_bmi = (data['smoking'] == 0) & (np.random.uniform(size=n) < 0.2)
data.loc[mask_bmi, 'bmi'] = np.nan
print("Missing values per column:")
print(data.isnull().sum())
# Simple multiple imputation using IterativeImputer
# Note: for proper MI with pooling, consider the miceforest package
predictors = ['age', 'sbp', 'cholesterol', 'smoking', 'bmi']
X = data[predictors].values
y = data['cvd'].values
# Create multiple imputed datasets and pool coefficients
m = 20
coefs = []
for i in range(m):
imputer = IterativeImputer(random_state=i, max_iter=10, sample_posterior=True)
X_imp = imputer.fit_transform(X)
model = LogisticRegression(max_iter=1000, penalty=None)
model.fit(X_imp, y)
coefs.append(np.concatenate([model.intercept_, model.coef_[0]]))
# Pool results (simplified Rubin's rules)
coefs = np.array(coefs)
pooled_mean = coefs.mean(axis=0)
within_var = coefs.var(axis=0) # simplified
names = ['intercept'] + predictors
print("\nPooled coefficients:")
for name, coef in zip(names, pooled_mean):
print(f" {name}: {coef:.4f}")What the code shows. This demonstrates multiple imputation end-to-end on a simulated cohort where cholesterol and BMI are made deliberately missing in a way that depends on observed variables (MAR). In R, md.pattern() prints a table showing which combinations of variables are missing together — useful for spotting structure in the missingness. mice() then creates 20 filled-in datasets, the model is fitted in each, and pool() combines them with Rubin’s rules so the reported standard errors honestly reflect the extra uncertainty from imputing. The pooled coefficients are what you would report. The key takeaway for a health researcher: you do not get a single “best guess” for each missing value — you get several plausible values, and the spread between them is folded into your confidence intervals. The Python version mimics this with IterativeImputer; note the comment that for fully correct pooling in Python you would reach for a dedicated MI package such as miceforest.
17.5 Step 4: Variable Selection
Variable selection (also called predictor selection) means deciding which patient characteristics go into the model. A “candidate predictor” is any variable you are considering, such as age, blood pressure, or smoking status.
17.5.1 Choosing Candidate Predictors
The best predictor selection strategy is expert knowledge. Decades of cardiovascular research tell us that age, sex, blood pressure, cholesterol, smoking, and diabetes are important predictors of CVD. Starting with established clinical knowledge prevents the discovery of spurious associations (chance patterns that will not hold up in new patients) and produces more generalisable models.
When expert knowledge is insufficient or the domain is novel, data-driven selection can supplement clinical reasoning. However, the approach matters enormously.
17.5.2 The Dangers of Stepwise Selection
Stepwise selection is an automated procedure that adds or removes predictors one at a time based on statistical tests (forward = start empty and add; backward = start full and remove; both = a mix). Despite being among the most commonly used strategies, it is also one of the most harmful:
- It inflates the apparent significance of selected variables (multiple testing without correction).
- It produces unstable models: small changes in the data can lead to entirely different variable selections.
- It biases regression coefficients away from zero (selected variables appear stronger than they are).
- It yields overly optimistic performance estimates.
WarningAvoid Univariable Screening
A particularly damaging practice is pre-selecting variables based on univariable p-values (e.g., keeping only variables with p < 0.20). This ignores confounding, suppression, and the fact that weak individual predictors can be strong in combination. If you must reduce the candidate set, use clinical knowledge, not univariable screening.
17.5.3 Better Approaches
- Full model with all pre-specified predictors: “pre-specified” means you decide the predictor list in advance, before looking at the outcome. If sample size permits, include all clinically plausible predictors and apply shrinkage (see below) to keep the model from over-reacting to noise. This is the recommended approach by Steyerberg (2019).
- Penalised regression: a family of methods that deliberately hold back the model’s effect estimates (“penalise” large coefficients) so it does not overfit. LASSO (L1 penalty) performs automatic variable selection by shrinking some coefficients all the way to zero, effectively dropping those predictors. Elastic net combines L1 and L2 penalties. These are covered in detail in Chapter 6.
- Background knowledge with data-driven fine-tuning: start with a core set of established predictors and consider a small number of additional candidates.
17.6 Step 5: Understand and Combat Overfitting
17.6.1 What Is Overfitting?
Overfitting occurs when a model captures not only the true underlying patterns in the data but also the noise — the random variation specific to the development sample. An overfitted model performs well on the data it was trained on but poorly on new data. Think of a student who memorises the answers to last year’s exam: they ace the practice paper but stumble on this year’s questions.
Why a clinician should care. An overfitted risk score will overstate how confident it is, telling you a patient is at 2% or 90% risk when the truth is closer to the average. Acting on those exaggerated numbers can lead to over- or under-treatment.
Clinical data is especially prone to overfitting because:
- Sample sizes are often small relative to the number of candidate predictors.
- Events may be rare, limiting the information available to estimate parameters.
- Complex interactions are tempting to include but expensive. (An “interaction term” lets one predictor’s effect depend on another — for example, smoking mattering more in older patients. Each one adds a parameter the model must estimate, using up scarce “degrees of freedom”, the budget of information your events can support.)
17.6.2 Detecting Overfitting: Optimism
The optimism of a model is the difference between its apparent performance (how good it looks on the data it was built from) and its expected performance on new data. Every model flatters itself a little on its own data, so the apparent numbers are always a bit too rosy. You can estimate the size of that flattery using internal validation techniques such as bootstrapping — repeatedly re-fitting the model on random resamples of the data (covered in Chapter 18).
If a model shows high optimism, it is overfitted. The solution is shrinkage.
17.6.3 Shrinkage Methods
Shrinkage means deliberately pulling the model’s effect estimates back toward zero so its predictions are less extreme and generalise better.
Uniform shrinkage multiplies all regression coefficients by a single factor \(s\) (a number between 0 and 1) estimated from the data, typically via bootstrapping. Shrinking every coefficient would also drag the overall level of predicted risk downward, so one number — the model’s baseline term, or intercept — is re-estimated afterwards to put the average predicted risk back where it belongs. That keeps the model’s average honest, which is one part of what “calibration” means; the full picture is in Section 18.4. The effect of the whole operation is to pull extreme predictions back toward the average, reducing the impact of noise.
Penalised methods build shrinkage directly into the model-fitting step rather than applying it afterwards:
- Ridge regression (L2 penalty) shrinks all coefficients toward zero but never sets any to exactly zero — every predictor stays in the model, just toned down.
- LASSO (L1 penalty) can shrink coefficients to exactly zero, so it also performs variable selection by dropping weak predictors.
- Elastic net combines both penalties.
How hard to shrink (the penalty strength) is chosen by cross-validation — repeatedly setting aside part of the data, trying different penalty strengths, and keeping the one that predicts the held-out part best.
Code
library(glmnet)
# Simulate Framingham-like data. We deliberately make this a setting where
# overfitting is a real threat: 500 patients, and 20 pure-noise variables
# mixed in among the 6 genuine predictors.
set.seed(42)
n <- 500
age <- rnorm(n, 55, 10)
male <- rbinom(n, 1, 0.5)
sbp <- rnorm(n, 130, 20)
chol <- rnorm(n, 220, 40)
smoking <- rbinom(n, 1, 0.25)
diabetes <- rbinom(n, 1, 0.10)
# True model: only these six variables actually affect risk
lp <- -6 +
0.05 * age +
0.3 * male +
0.01 * sbp +
0.005 * chol +
0.4 * smoking +
0.5 * diabetes
cvd <- rbinom(n, 1, plogis(lp))
# 20 candidate predictors that are unrelated to the outcome
noise <- matrix(rnorm(n * 20), n, 20)
colnames(noise) <- paste0("noise", 1:20)
df <- data.frame(age, male, sbp, chol, smoking, diabetes, noise, cvd)
X <- model.matrix(cvd ~ ., data = df)[, -1]
y <- df$cvd
# Fit the same model three ways
fit_full <- glm(cvd ~ ., data = df, family = binomial)
fit_ridge <- cv.glmnet(X, y, family = "binomial", alpha = 0, nfolds = 10)
fit_lasso <- cv.glmnet(X, y, family = "binomial", alpha = 1, nfolds = 10)
# Compare the six real predictors against their known true values
real <- c("age", "male", "sbp", "chol", "smoking", "diabetes")
pick <- function(fit) {
as.vector(coef(fit, s = "lambda.min"))[match(
real,
c("(Intercept)", colnames(X))
)]
}
comparison <- data.frame(
truth = c(0.05, 0.3, 0.01, 0.005, 0.4, 0.5),
unpenalised = coef(fit_full)[real],
ridge = pick(fit_ridge),
lasso = pick(fit_lasso)
)
print(round(comparison, 4))
# How many of the 20 noise variables did each method keep?
kept_noise <- function(fit) {
b <- as.vector(coef(fit, s = "lambda.min"))[-1]
sum(b[grepl("noise", colnames(X))] != 0)
}
cat(
"\nNoise variables kept (of 20) -- ridge:",
kept_noise(fit_ridge),
"| LASSO:",
kept_noise(fit_lasso),
"\n"
)
# Cross-validation curves
par(mfrow = c(1, 2), mar = c(4.5, 4.4, 4.6, 1.2))
for (m in list(
list(f = fit_ridge, t = "Ridge"),
list(f = fit_lasso, t = "LASSO")
)) {
plot(m$f, ylab = "Prediction error (deviance)")
title(main = paste(m$t, "- cross-validated error"), line = 3.2)
mtext(
"number of predictors kept",
side = 3,
line = 2.0,
cex = 0.7,
col = "grey30"
)
abline(v = -log(m$f$lambda.min), col = "#2166ac", lwd = 2)
legend(
"top",
bty = "n",
cex = 0.7,
lwd = c(2, 1),
lty = c(1, 3),
col = c("#2166ac", "black"),
seg.len = 1.6,
legend = c("lambda.min (lowest error)", "lambda.1se (simpler, within 1 SE)")
)
}Code
import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression, LogisticRegressionCV
from sklearn.preprocessing import StandardScaler
from scipy.special import expit
# Same setup as the R example: 500 patients, 6 real predictors, 20 noise ones.
# The random draws differ from R's, so the numbers will not match exactly.
np.random.seed(42)
n = 500
age = np.random.normal(55, 10, n)
male = np.random.binomial(1, 0.5, n)
sbp = np.random.normal(130, 20, n)
chol = np.random.normal(220, 40, n)
smoking = np.random.binomial(1, 0.25, n)
diabetes = np.random.binomial(1, 0.10, n)
lp = -6 + 0.05*age + 0.3*male + 0.01*sbp + 0.005*chol + 0.4*smoking + 0.5*diabetes
cvd = np.random.binomial(1, expit(lp))
noise = np.random.normal(size=(n, 20))
real_names = ['age', 'male', 'sbp', 'chol', 'smoking', 'diabetes']
names = real_names + [f'noise{i}' for i in range(1, 21)]
X = np.column_stack([age, male, sbp, chol, smoking, diabetes, noise])
# Penalties act on the coefficient scale, so predictors must be standardised
# first (glmnet does this internally; scikit-learn does not). We rescale the
# coefficients afterwards so they are comparable to the true values.
scaler = StandardScaler().fit(X)
Xs = scaler.transform(X)
# Two settings matter here and are easy to get wrong:
# scoring='neg_log_loss' -- the default is accuracy, which a "predict the
# majority class" model maximises, so it would choose maximum shrinkage
# Cs -- the default grid of 10 values is too coarse to locate a good penalty
grid = np.logspace(-3, 3, 60)
model_full = LogisticRegression(penalty=None, max_iter=5000).fit(X, cvd)
model_ridge = LogisticRegressionCV(Cs=grid, penalty='l2', cv=10,
scoring='neg_log_loss', max_iter=5000,
random_state=42).fit(Xs, cvd)
model_lasso = LogisticRegressionCV(Cs=grid, penalty='l1', solver='saga', cv=10,
scoring='neg_log_loss', max_iter=20000,
random_state=42).fit(Xs, cvd)
def original_scale(model):
return model.coef_[0] / scaler.scale_
comparison = pd.DataFrame({
'truth': [0.05, 0.3, 0.01, 0.005, 0.4, 0.5],
'unpenalised': model_full.coef_[0][:6],
'ridge': original_scale(model_ridge)[:6],
'lasso': original_scale(model_lasso)[:6],
}, index=real_names)
print(comparison.round(4))
for label, model in [('ridge', model_ridge), ('lasso', model_lasso)]:
kept = np.sum(model.coef_[0][6:] != 0)
print(f"Noise variables kept (of 20) -- {label}: {kept}")What the code shows. Because these are simulated data we know the true coefficients, so the printed table can be read as a scorecard. The unpenalised column overshoots: diabetes comes out at about 0.88 when the truth is 0.50. That exaggeration is overfitting, and it is exactly what shrinkage exists to correct — the ridge and LASSO columns sit visibly closer to the truth. The noise count then shows the difference between the two penalties in one number: ridge keeps all 20 worthless variables (shrunk, but never removed), while LASSO sets most of them to exactly zero and drops them outright.
How to read the cross-validation plots. These two panels are the workhorse figure of penalised regression, and they repay a slow first read.
- Each red dot is one fitted model. The data are split into 10 parts; the model is fitted on 9 and its prediction error measured on the held-out one, rotating through all 10. The dot is the average of those 10 errors, so it estimates how the model would do on patients it has not seen.
- The y-axis is prediction error. The label says “binomial deviance”, which is just this model family’s way of scoring how far the predicted probabilities were from what actually happened. Lower is better. Do not try to interpret the number itself — 1.29 means nothing in isolation. Only the shape of the curve matters.
- The x-axis is the amount of penalty, and its direction is the part that trips people up.
glmnet(version 5 onwards) plots \(-\log(\lambda)\), the negative log penalty, so left means heavy penalty and strong shrinkage, right means almost no penalty — the far right is essentially the ordinary unpenalised model. Older textbook figures plot \(\log(\lambda)\) instead, which runs the opposite way, so check the axis label before interpreting any such plot you meet. - The grey bars show how much that error estimate wobbles (±1 standard error across the 10 splits) depending on which patients happened to be held out. They are a humility check: where the bars are wide, the exact position of the curve’s lowest point is not to be taken too literally, because a different split of the same patients would have moved it.
- The numbers along the top are how many predictors still have a non-zero coefficient. For ridge this is stuck at 26 the whole way across, because ridge never eliminates anything. For LASSO it climbs from 0 on the left to 26 on the right. That contrast is the L1-versus-L2 difference, made visible.
The curve is U-shaped, and both arms are informative. Move right from the bottom and error rises because the model is too free and starts fitting noise — overfitting. Move left and error rises again because everything has been squashed so hard toward zero that genuine signal is lost — over-shrinking. The bottom of the U is the best available trade-off.
The two vertical lines mark the two conventional choices, and this is where the axis direction matters most:
lambda.min(solid blue) is simply the penalty with the lowest average error — the bottom of the U.lambda.1se(dotted) is the strongest penalty whose error is still within one standard error of that minimum. It therefore sits to the left oflambda.min, because left means more penalty. The logic is that the minimum is itself estimated with noise, so deliberately erring toward a simpler, more heavily shrunk model is often the safer bet — and for LASSO it yields a shorter predictor list, which matters if the model has to be used at the bedside.
Either choice is defensible; what matters is deciding in advance which you will use, rather than fitting both and reporting whichever looks better.
NoteShrinkage in plain terms
A shrinkage factor is a number between 0 and 1 that you multiply coefficients by to make predictions less extreme. An overfitted model produces predictions that are too confident — very high risks too high, very low risks too low. Shrinkage pulls those predictions back toward the average event rate, which is usually closer to the truth for new patients. A shrinkage factor of 0.85 means the development model was about 15% too optimistic and its effects should be scaled back accordingly.
17.7 Step 6: Build the Model — A Framingham Example
Let us now walk through a complete model development process using simulated data based on the Framingham Heart Study. We will predict 10-year CVD risk using established risk factors. This is the “derivation” (or development) cohort — the dataset in which the model is built. Later, the model should be tested in a separate “validation” cohort, ideally patients from a different time or place, to check it still works.
17.7.1 Data Preparation
Code
library(rms)
library(ggplot2)
# Simulate Framingham-like cohort
set.seed(2024)
n <- 2000
framingham <- data.frame(
age = round(runif(n, 30, 74)),
male = rbinom(n, 1, 0.48),
sbp = round(rnorm(n, 130, 18)),
total_chol = round(rnorm(n, 210, 38)),
hdl_chol = round(rnorm(n, 52, 15)),
smoking = rbinom(n, 1, 0.22),
diabetes = rbinom(n, 1, 0.08),
bp_treatment = rbinom(n, 1, 0.15)
)
# Generate outcome based on known risk factors
lp <- with(
framingham,
-7.5 +
0.06 * age +
0.4 * male +
0.012 * sbp +
0.005 * total_chol -
0.02 * hdl_chol +
0.5 * smoking +
0.7 * diabetes +
0.3 * bp_treatment
)
framingham$cvd_10yr <- rbinom(n, 1, plogis(lp))
cat("Event rate:", round(mean(framingham$cvd_10yr) * 100, 1), "%\n")
cat("Number of events:", sum(framingham$cvd_10yr), "\n")
cat("Events per variable:", round(sum(framingham$cvd_10yr) / 8, 1), "\n")
# Set up the rms environment for enhanced model building
dd <- datadist(framingham)
options(datadist = "dd")
# Fit the model using lrm (logistic regression model from rms)
# Allow nonlinearity for continuous variables using restricted cubic splines
fit <- lrm(
cvd_10yr ~ rcs(age, 4) +
male +
rcs(sbp, 3) +
total_chol +
hdl_chol +
smoking +
diabetes +
bp_treatment,
data = framingham,
x = TRUE,
y = TRUE
)
print(fit)
# Check for nonlinearity
anova(fit)
# Visualise partial effects for age and SBP side by side.
# fun = plogis converts the model's log odds onto the risk (probability)
# scale, which is what a clinician actually wants to see.
partial <- function(var, label) {
p <- as.data.frame(Predict(fit, name = var, fun = plogis))
data.frame(
x = p[[var]],
risk = p$yhat,
lo = p$lower,
hi = p$upper,
panel = label,
row.names = NULL
)
}
age_lab <- "Age (years)"
sbp_lab <- "Systolic blood pressure (mmHg)"
curves <- rbind(partial("age", age_lab), partial("sbp", sbp_lab))
# Where the patients actually are, drawn as a rug along the bottom
rugs <- rbind(
data.frame(x = framingham$age, panel = age_lab),
data.frame(x = framingham$sbp, panel = sbp_lab)
)
ggplot(curves, aes(x, risk)) +
geom_ribbon(aes(ymin = lo, ymax = hi), fill = "grey75", alpha = 0.6) +
geom_line(colour = "#2166ac", linewidth = 0.9) +
geom_rug(
data = rugs,
aes(x = x),
inherit.aes = FALSE,
sides = "b",
alpha = 0.05
) +
facet_wrap(~panel, scales = "free_x", strip.position = "bottom") +
scale_y_continuous(labels = function(p) paste0(round(p * 100), "%")) +
labs(
title = "Predicted 10-year CVD risk across the range of each predictor",
subtitle = paste(
"Other predictors held fixed at: female, non-smoker,",
"no diabetes, untreated,\ntotal cholesterol 210,",
"HDL 52. Shaded band = 95% confidence interval."
),
x = NULL,
y = "Predicted 10-year CVD risk"
) +
theme_minimal(base_size = 11) +
theme(
strip.placement = "outside",
panel.grid.minor = element_blank(),
plot.subtitle = element_text(size = 8.5, colour = "grey35")
)Code
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy.special import expit
np.random.seed(2024)
n = 2000
framingham = pd.DataFrame({
'age': np.random.randint(30, 75, n),
'male': np.random.binomial(1, 0.48, n),
'sbp': np.round(np.random.normal(130, 18, n)),
'total_chol': np.round(np.random.normal(210, 38, n)),
'hdl_chol': np.round(np.random.normal(52, 15, n)),
'smoking': np.random.binomial(1, 0.22, n),
'diabetes': np.random.binomial(1, 0.08, n),
'bp_treatment': np.random.binomial(1, 0.15, n)
})
lp = (-7.5 + 0.06 * framingham['age'] + 0.4 * framingham['male'] +
0.012 * framingham['sbp'] + 0.005 * framingham['total_chol'] -
0.02 * framingham['hdl_chol'] + 0.5 * framingham['smoking'] +
0.7 * framingham['diabetes'] + 0.3 * framingham['bp_treatment'])
framingham['cvd_10yr'] = np.random.binomial(1, expit(lp))
print(f"Event rate: {framingham['cvd_10yr'].mean()*100:.1f}%")
print(f"Number of events: {framingham['cvd_10yr'].sum()}")
print(f"Events per variable: {framingham['cvd_10yr'].sum() / 8:.1f}")
# Fit logistic regression
predictors = ['age', 'male', 'sbp', 'total_chol', 'hdl_chol',
'smoking', 'diabetes', 'bp_treatment']
X = sm.add_constant(framingham[predictors])
y = framingham['cvd_10yr']
model = sm.Logit(y, X)
result = model.fit(disp=0)
print("\n", result.summary())
# Calculate predicted probabilities
framingham['pred_prob'] = result.predict(X)
# Summary of predicted risk distribution
print("\nPredicted risk distribution:")
print(framingham['pred_prob'].describe())What the code shows. This builds the actual prediction model on a simulated Framingham-style cohort. The first printed lines report the event rate, the number of events, and the EPV — a sanity check that we have enough events for the number of predictors before trusting anything else. In R, lrm() fits the logistic model and print(fit) reports the coefficients plus the apparent C-statistic; rcs(age, 4) lets age have a flexible (non-straight-line) relationship with risk using restricted cubic splines with 4 knots. The anova(fit) output tests whether that flexibility is actually needed (a significant nonlinear term means a straight line would misfit). The Python version fits the equivalent model with statsmodels and summarises the distribution of predicted risks, so you can see whether the model produces a useful spread of probabilities or bunches everyone near the average.
How to read the partial-effect figure. Each panel answers one question: if I walked a single patient across the whole range of this one predictor and changed nothing else about them, how would their predicted risk move? Three features are worth naming explicitly.
- “Partial” means everything else is pinned. Each curve is drawn for one reference patient — here a female non-smoker without diabetes, untreated, with typical cholesterol — and only the panel’s own predictor varies. Change the reference patient and the curves shift up or down, though their shape stays the same. This is why such plots always carry an “adjusted to” note; it is not fine print, it tells you which patient you are looking at.
- The y-axis is risk, not log odds. Logistic regression reports its results on an internal scale (log odds) that means little at the bedside;
fun = plogisconverts them into plain probabilities, so the numbers are ones you could say out loud to a patient. Both panels share one y-axis, which makes the comparison the real point of the figure: across its observed range, age moves predicted risk from about 2% to about 22%, whereas systolic blood pressure moves it only from about 2% to 7%. Age is doing far more work than blood pressure in this model — a conclusion you cannot draw by glancing at coefficients, because coefficients depend on each predictor’s units. - The shaded band is a 95% confidence interval, and the tick marks along the bottom show where the patients actually are. Read them together. The band flares at the far left and right of each panel precisely because there are few patients out there, so the curve is extrapolating on thin evidence. Treat the ends of these curves with suspicion; that is where a smooth-looking line is least trustworthy.
One trap is worth pointing out, because this figure walks straight into it. The age curve looks bent — flat to about 45, then steepening — which invites the conclusion that age has a genuinely non-straight-line effect and that the splines have earned their keep. They have not. We simulated these data ourselves, and we built age in as a perfectly straight line; the anova(fit) output agrees, reporting the nonlinear term for age as nowhere near significant (p ≈ 0.45).
So where does the bend come from? Logistic regression works internally on a scale that has no ceiling, then squeezes the result into the 0–100% range that a probability has to live in. That squeezing is not uniform: near 0% there is very little room left, so a straight-line increase gets compressed, while higher up it has room to climb. Redraw any straight-line effect as a risk and you get a curve of exactly this shape. The lesson is practical: a bent curve on a risk plot is not evidence of a non-straight-line effect — the anova test is what tells you that. Keep both views, and use each for what it is good for: the risk scale for talking to patients, the anova for deciding whether the extra flexibility is real.
17.7.2 Model Specification Decisions
In the Framingham example above, several important decisions were made:
- Pre-specified predictors: we included the established Framingham risk factors, not variables discovered through data-dredging.
- Nonlinearity: for continuous variables like age and blood pressure, we allowed the relationship with risk to bend rather than follow a straight line, using restricted cubic splines (in R’s
rmspackage). A spline is a smooth curve stitched together from simple pieces; it lets risk rise gently and then steeply with age, for instance, instead of forcing a single fixed slope. The relationship between age and CVD risk is genuinely curved, so a straight line would misfit. - No interactions without prior evidence: we did not fish for interactions (terms where one predictor’s effect depends on another). If clinical knowledge suggests an interaction (e.g., the effect of cholesterol might differ by sex), it should be pre-specified rather than discovered by trial and error.
- Sample size check: with approximately 250 events and 8 predictors (plus a few extra degrees of freedom for the splines), we have approximately 25 events per variable, which is comfortable.
17.8 Step 7: Assess the Final Model
After fitting, we need to examine the model critically. Two questions matter most. Does it discriminate — rank patients who go on to have the event above those who do not? And is it calibrated — when it says 20% risk, do about 20% of such patients actually have the event? A model can discriminate well yet be poorly calibrated, so both must be checked.
TipThis section is a first look; Chapter 18 is the full treatment
Performance assessment is a large enough topic to have the next chapter to itself, and this section deliberately does not duplicate it. What follows is the minimum you should look at before believing your own model, with just enough interpretation to read the output. If you want the definitions, the formulas and the reasoning behind any of these measures, go to Chapter 18, which organises them into five domains (Van Calster et al. 2025) and covers:
- Discrimination and the C-statistic — Section 18.3 and Section 18.3.1, including what the number does not tell you (Section 18.3.2)
- Calibration — Section 18.4, with the O:E ratio (Section 18.4.2), the calibration slope (Section 18.4.3) and how to read a calibration plot (Section 18.4.4)
- Overall performance, including the Brier score and its no-predictor benchmark — Section 18.5
- Internal validation and bootstrap optimism correction — Section 18.8 and Section 18.8.3
- Clinical utility — net benefit and decision curve analysis, which ask the question none of the above do: would using this model actually improve decisions?
Code
library(rms)
# Assuming fit and framingham from above
pred_probs <- predict(fit, type = "fitted")
event_rate <- mean(framingham$cvd_10yr)
top <- ceiling(max(pred_probs) * 10) / 10 # zoom the axes to the data
cat("C-statistic (AUC):", round(fit$stats["C"], 3), "\n")
# Three complementary views of the same model, side by side
par(mfrow = c(1, 3), mar = c(4.6, 4.4, 3.6, 1.2))
# A. Does the model spread patients out, or call everyone average?
hist(
pred_probs,
breaks = 50,
col = "steelblue",
border = "white",
main = "A. Spread of predicted risk",
xlab = "Predicted 10-year risk",
ylab = "Number of patients"
)
abline(v = event_rate, col = "firebrick", lwd = 2, lty = 2)
text(
event_rate,
par("usr")[4] * 0.92,
paste0(" overall rate ", round(event_rate * 100), "%"),
col = "firebrick",
cex = 0.85,
adj = 0
)
# B. Discrimination: do the patients who had events get higher risks?
boxplot(
pred_probs ~ framingham$cvd_10yr,
names = c("No CVD", "CVD"),
main = "B. Do events get higher risks?",
xlab = "Actual outcome",
ylab = "Predicted 10-year risk",
col = c("lightblue", "salmon"),
outcex = 0.4
)
# C. Calibration: when the model says 20%, do 20% have the event?
# We draw this panel by hand rather than using val.prob()'s default plot,
# so the three lines are told apart at a glance and labelled in plain words.
y <- framingham$cvd_10yr
q95 <- quantile(pred_probs, 0.95) # beyond here, few patients: thin evidence
# The straight-line fit: allowed only an intercept and a slope, so it can
# describe "risks uniformly too high/low" and "risks too spread out/bunched",
# and nothing more.
cal_fit <- glm(y ~ lp_hat, family = binomial,
data = data.frame(y = y, lp_hat = qlogis(pred_probs)))
grid <- seq(min(pred_probs), max(pred_probs), length.out = 200)
straight <- plogis(predict(cal_fit, newdata = data.frame(lp_hat = qlogis(grid))))
# The flexible fit: follows the data locally, with no assumed shape at all.
flexible <- lowess(pred_probs, y, iter = 0)
plot(0, 0, type = "n", xlim = c(0, top), ylim = c(-0.06, top), yaxt = "n",
xlab = "Predicted 10-year risk", ylab = "Observed frequency",
main = "C. Calibration: predicted vs observed")
axis(2, at = seq(0, top, 0.2))
# Shade the region where there are hardly any patients
rect(q95, -0.06, top, top, col = rgb(0, 0, 0, 0.05), border = NA)
abline(v = q95, col = "grey45", lty = 3)
text(q95, top * 0.03, " 95% of patients\n to the left",
cex = 0.6, col = "grey35", adj = 0)
abline(0, 1, col = "grey70", lwd = 7) # ideal
lines(grid, straight, col = "#2166ac", lwd = 2.4) # 2-number
lines(flexible$x, flexible$y, col = "#d6604d", lwd = 2.4, lty = 2) # flexible
# Density strip: where the patients actually are
h <- hist(pred_probs, breaks = 60, plot = FALSE)
bw <- diff(h$mids)[1]
rect(h$mids - bw / 2, -0.055,
h$mids + bw / 2, -0.055 + 0.05 * h$counts / max(h$counts),
col = "grey55", border = NA)
mtext("where patients are", side = 1, line = -1.3,
at = top * 0.62, cex = 0.6, col = "grey35")
legend("topleft", bty = "n", cex = 0.7, lwd = c(7, 2.4, 2.4),
lty = c(1, 1, 2), seg.len = 1.8,
col = c("grey70", "#2166ac", "#d6604d"),
legend = c("Perfect calibration",
"Straight-line fit (intercept + slope only)",
"Flexible fit (no assumed shape)"))
# val.prob still gives us the statistics; pl = FALSE suppresses its own plot
val <- val.prob(pred_probs, y, pl = FALSE)
# val.prob returns 18 statistics. These are the six worth reading.
key <- c("C (ROC)", "Brier", "Intercept", "Slope", "Eavg", "Emax")
summary_table <- data.frame(
value = round(val[key], 3),
ideal = c("1.0", "0", "0", "1.0", "0", "0"),
meaning = c(
"discrimination: 0.5 = coin flip, 1.0 = perfect",
"overall accuracy of the probabilities (lower is better)",
"calibration-in-the-large: are risks too high or too low overall?",
"calibration slope: are risks too extreme (<1) or too flat (>1)?",
"average gap between predicted and observed risk",
"worst gap between predicted and observed risk"
)
)
print(summary_table, right = FALSE)
# Internal validation via bootstrapping (optimism-corrected)
set.seed(123)
val_boot <- validate(fit, B = 200)
print(val_boot)
cat("\nApparent C-statistic: ", round(fit$stats["C"], 3), "\n")
cat(
"Optimism-corrected C-statistic:",
round(0.5 * (val_boot["Dxy", "index.corrected"] + 1), 3),
"\n"
)Code
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from sklearn.metrics import roc_auc_score, brier_score_loss
from sklearn.calibration import calibration_curve
# Assuming framingham and result from above
pred_probs = framingham['pred_prob'].values
y_true = framingham['cvd_10yr'].values
# 1. Discrimination
auc = roc_auc_score(y_true, pred_probs)
print(f"C-statistic (AUC): {auc:.3f}")
print(f"Brier score: {brier_score_loss(y_true, pred_probs):.3f}")
# Calibration intercept and slope, fitted on the linear predictor
logit_p = np.log(pred_probs / (1 - pred_probs))
cal = sm.Logit(y_true, sm.add_constant(logit_p)).fit(disp=0)
print(f"Calibration intercept: {cal.params[0]:.3f} (ideal 0)")
print(f"Calibration slope: {cal.params[1]:.3f} (ideal 1)")
# Three complementary views of the same model, side by side
event_rate = y_true.mean()
top = np.ceil(pred_probs.max() * 10) / 10
fig, axes = plt.subplots(1, 3, figsize=(15, 4.6))
# A. Does the model spread patients out, or call everyone average?
axes[0].hist(pred_probs, bins=50, color="steelblue", edgecolor="white")
axes[0].axvline(event_rate, color="firebrick", ls="--", lw=2)
axes[0].annotate(f" overall rate {event_rate:.0%}", (event_rate, 0.92),
xycoords=("data", "axes fraction"), color="firebrick")
axes[0].set(xlabel="Predicted 10-year risk", ylabel="Number of patients",
title="A. Spread of predicted risk")
# B. Discrimination: do the patients who had events get higher risks?
bp = axes[1].boxplot([pred_probs[y_true == 0], pred_probs[y_true == 1]],
tick_labels=["No CVD", "CVD"], patch_artist=True,
flierprops=dict(markersize=2))
for patch, colour in zip(bp['boxes'], ["lightblue", "salmon"]):
patch.set_facecolor(colour)
axes[1].set(xlabel="Actual outcome", ylabel="Predicted 10-year risk",
title="B. Do events get higher risks?")
# C. Calibration: when the model says 20%, do 20% have the event?
prob_true, prob_pred = calibration_curve(y_true, pred_probs, n_bins=10,
strategy='quantile')
axes[2].plot([0, top], [0, top], '--', color='grey', label="Perfect calibration")
axes[2].plot(prob_pred, prob_true, 's-', color="steelblue", label="Model")
axes[2].set(xlim=(0, top), ylim=(0, top), xlabel="Predicted 10-year risk",
ylabel="Observed frequency",
title="C. Calibration: predicted vs observed")
axes[2].legend()
plt.tight_layout()
plt.show()What the code shows. This is a first critical look at the finished model, from three angles that answer different questions. Read them left to right.
Panel A — does the model actually distinguish anyone? This is the distribution of predicted risks across the 2,000 patients, with the dashed red line at the overall event rate (12%). The shape you want is a wide one. Predicted risks here run from under 1% to about 76%, with the middle half of patients between 4% and 16% — so the model is genuinely sorting people, not hedging. The failure mode this panel is designed to catch is a narrow spike sitting on top of the red line: a model that tells every patient “about 12%” is clinically worthless no matter how respectable its other statistics look, because no decision could ever change based on it.
Panel B — discrimination, made concrete. Median predicted risk is 7.5% among patients who did not have a CVD event and 18.5% among those who did. The events do get higher risks, which is what we want. But look at how much the two boxes still overlap: plenty of patients who had events were given lower risks than patients who did not. That overlap is what a C-statistic of 0.75 looks like — and 0.75 is a perfectly respectable clinical model. If you were expecting two cleanly separated boxes, recalibrate that expectation now; you will not see it in real risk prediction, and if you do, suspect a leaked outcome.
Panel C — calibration. Discrimination only asks whether the ranking is right. Calibration asks whether the numbers can be believed: when the model says 20%, do about 20 in 100 such patients actually have the event? Predicted risk runs along the bottom, the proportion who actually had the event up the side, and the thick grey diagonal is the answer you want — every predicted risk matching the observed frequency exactly.
Why there are two lines, not one. This is the part of the figure that confuses people, and the two lines are asking the same question with different amounts of freedom.
- The blue straight-line fit is deliberately restricted. It is allowed exactly two numbers — an intercept and a slope — so it can only describe two kinds of problem: predictions that are uniformly too high or too low (wrong intercept), and predictions that are too spread out or too bunched together (wrong slope). It cannot bend. Those two numbers are precisely the
InterceptandSlopein the table below, so this line is the visual form of that two-number summary. - The red dashed flexible fit is allowed to go wherever the data lead. It assumes no shape whatsoever — it simply follows the local pattern of who actually had events, which is why the technical name for it is nonparametric (“no assumed shape”). It can therefore reveal trouble confined to one stretch of the risk range: a model that is spot-on for low-risk patients and badly off for high-risk ones.
Reading them together is the point. Where the two lines agree, the tidy two-number summary is an honest description of the model’s calibration. Where the flexible line wanders away from the straight one, the miscalibration has a shape that those two numbers are hiding — and you would never know from the table alone.
That is exactly what happens here, and it is also a lesson in not over-reading. Below the dotted vertical line — which holds 95% of the patients — the flexible fit tracks the diagonal within about 1 percentage point: the model is well calibrated where nearly all the patients are. To the right of it, in the shaded strip, the flexible fit peels upward and the gap grows to nearly 12 percentage points. But that whole region contains only about 100 of the 2,000 patients. The grey density strip along the bottom shows the same thing: patients pile up below 0.2 and thin out to almost nothing past 0.4. So the dramatic-looking divergence is built on very few people, and it is the same sparse tail that produces the alarming Emax of 0.117 in the table below. Always check where the patients are before believing the ends of a calibration curve.
ImportantWhat the numbers mean
val.prob() prints 18 statistics. Six are worth your attention, and knowing which to ignore is half the skill.
| Statistic | Here | Ideal | What it means |
|---|---|---|---|
| C (ROC) | 0.750 | 1.0 | Pick one patient who had the event and one who did not, at random. This is the probability the model gave the event patient the higher risk. 0.5 is a coin flip. Full treatment in Section 18.3.1. |
| Brier | 0.092 | 0 | Average squared distance between the predicted probability and what actually happened (0 or 1). Lower is better, but the scale depends on how common the outcome is — see below, and Section 18.5. |
| Intercept | 0.000 | 0 | Calibration-in-the-large: are the predicted risks systematically too high or too low overall? Usually reported as the O:E ratio instead — Section 18.4.2. |
| Slope | 1.000 | 1.0 | Calibration slope: below 1 means predictions are too extreme (high risks too high, low risks too low) — the classic signature of overfitting. See Section 18.4.3. |
| Eavg | 0.005 | 0 | Average absolute gap between predicted and observed risk, in risk units. 0.005 is half a percentage point. |
| Emax | 0.117 | 0 | The worst such gap anywhere on the curve. Almost always driven by the sparse tail, so check panel C before worrying. |
How good is 0.75? This is the question the raw number never answers on its own. As a rough guide: 0.5 is chance, 0.6–0.7 poor, 0.7–0.8 acceptable and where most useful clinical models live, 0.8–0.9 good, above 0.9 excellent — and, in a clinical prediction model built from routine risk factors, above 0.9 is more often a warning sign than a triumph. It usually means something in the predictor set is a consequence of the outcome rather than a cause of it. See Table 15.2 in Chapter 15 for the fuller version, and treat all such bands as context-dependent: 0.75 for ten-year cardiovascular risk is a genuinely useful model, because whether a particular person has a heart attack is substantially a matter of chance that no set of baseline measurements can predict. There is a ceiling, and it is well below 1.
Two traps in that table. First, the Brier score cannot be read on its own. Predicting the overall event rate of 11.75% for every single patient — a model with no predictors at all — would score \(0.1175 \times 0.8825 = 0.104\). Our model manages 0.092. That is an improvement, but a modest one, and it is a useful corrective to the more flattering C-statistic. Always compare a Brier score against this no-predictor benchmark.
Second, and more important: the Intercept of exactly 0 and Slope of exactly 1 are not good news — they are arithmetically guaranteed. We measured calibration on the very data the model was fitted to, and a logistic model is by construction perfectly calibrated in its own development sample. These two numbers only become informative once they are computed on different data, or corrected for optimism as below. Seeing 0 and 1 here tells you nothing except that the arithmetic worked.
Finally, the statistics not in the table — Dxy, R2, D, U, Q, S:z, S:p — mostly restate the above or serve specialist purposes. Two are worth naming: Dxy (Somers’ D) is just the C-statistic rescaled so that 0 rather than 0.5 means no discrimination, via \(D_{xy} = 2 \times (C - 0.5)\), which is why the bootstrap code below converts it back with 0.5 * (Dxy + 1). And R2 is the Nagelkerke R-squared met in the sample-size section — 0.169 here.
Correcting for optimism. validate(fit, B = 200) re-runs the entire modelling process on 200 bootstrap resamples to estimate how much the apparent performance is inflated by having fitted and evaluated on the same patients. For this model the apparent C-statistic of 0.750 corrects down to 0.739, and the calibration slope — pinned at exactly 1 above — corrects to 0.943, meaning the predictions are mildly too extreme and would need shrinking by about 6% to be honest. The optimism is small here because 2,000 patients for 11 parameters is comfortable; in a smaller study the gap would be far wider. The corrected figures, not the apparent ones, are what belongs in a paper. Why this procedure works, and why it is not optional, is set out in Section 18.8 and Section 18.8.3.
None of the above asks the question that ultimately decides whether a model is worth deploying: would using it actually lead to better decisions than not using it? Discrimination and calibration are necessary but not sufficient for that. See the clinical utility and decision curve material in Chapter 18.
NoteOptimism
Optimism is the gap between how well a model appears to perform on the data it was built from and how well it actually performs on new data. Because the model has partly fitted the noise in its own development sample, the apparent C-statistic is always a little too rosy. Bootstrap validation estimates this gap so you can subtract it off and report a realistic number.
17.9 The Complete Workflow: A Summary
The complete prediction model development workflow can be summarised in the following steps:
- Define the clinical question (population, outcome, timing, intended use).
- Design the study with adequate sample size (Riley criteria).
- Prepare the data: handle missing values with multiple imputation.
- Select predictors based on clinical knowledge; avoid stepwise procedures.
- Specify the model: consider nonlinear terms for continuous predictors.
- Fit the model, applying shrinkage to combat overfitting.
- Evaluate performance: discrimination, calibration, clinical utility.
- Validate internally using bootstrapping or cross-validation.
- Report transparently following TRIPOD+AI guidelines (see Chapter 19).
ImportantThe Modelling Pipeline Is Not Linear
In practice, you will iterate. Poor calibration might send you back to model specification. Inadequate discrimination might prompt reconsidering the candidate predictors. Internal validation might reveal severe overfitting, requiring stronger shrinkage. The important thing is that each decision is documented and justified.
17.10 Common Mistakes to Avoid
| Mistake | Why It Is Harmful | Better Approach |
|---|---|---|
| Stepwise variable selection | Inflates significance, unstable | Pre-specify predictors; use penalisation |
| Univariable screening | Ignores confounding and suppression | Include all clinically plausible variables |
| Complete case analysis | Biases estimates, wastes data | Multiple imputation |
| Evaluating only on training data | Overestimates performance | Internal validation (bootstrap) |
| Dichotomising continuous predictors | Loses information, creates arbitrary groups | Model continuous variables continuously |
| Ignoring nonlinearity | Poor fit, biased predictions | Use splines or fractional polynomials |
| Reporting only accuracy or AUC | Incomplete picture | Report calibration, discrimination, clinical utility |
17.11 Exercises
These three exercises work through the same development pipeline this chapter built, mostly on the simulated Framingham cohort of Section 17.7.
TipExercise 1: Sample size calculation
You are developing a model to predict pre-eclampsia (expected prevalence 4%) using 12 candidate predictors, and a published model in a comparable population reported a C-statistic of 0.72.
- Use
pmsampsizeto determine the minimum required sample size. How many pregnancies would you need to observe? - Which of the criteria it evaluates is the binding one here?
- How does the answer compare with what the 10-events-per-variable rule of thumb would have told you?
NoteSolution
Code
# Exercise 1: Sample size calculation
# Pre-eclampsia model: prevalence 4%, 12 candidate predictors,
# an anticipated C-statistic of 0.72 from a published model in a
# comparable population.
library(pmsampsize)
# (a) Minimum sample size ----------------------------------------------------
# type = "b" binary outcome
# cstatistic = anticipated C-statistic; pmsampsize converts this into the
# Cox-Snell R-squared the criteria actually need
# parameters = number of predictor PARAMETERS, not variables (a categorical
# predictor with k levels costs k - 1, a spline costs more)
ss <- pmsampsize(
type = "b",
cstatistic = 0.72,
parameters = 12,
prevalence = 0.04
)
cat("\nMinimum sample size:", ss$sample_size, "pregnancies\n")
cat("Minimum number of events:", ceiling(ss$events), "\n")
cat("Events per parameter:", round(ss$EPP, 2), "\n")
# (b) Which criterion is binding? -------------------------------------------
# pmsampsize's results table has one row per criterion and a "final" row that
# is simply the largest of them. The binding criterion is whichever row the
# final row was taken from.
res <- ss$results_table
print(res)
crit_n <- res[rownames(res) != "Final", "Samp_size"]
binding <- names(which.max(crit_n))
cat("\nSample size demanded by each criterion:\n")
print(crit_n)
cat("\nBinding criterion:", binding, "->", max(crit_n), "pregnancies\n")
# (c) Comparison with the 10-EPV rule of thumb -----------------------------
# The old rule asks only for 10 events per parameter, and says nothing about
# how precisely the model must be estimated or how much it may overfit.
epv_events <- 10 * 12
epv_n <- ceiling(epv_events / 0.04)
cat("\n--- 10 events per variable rule ---\n")
cat("Events required:", epv_events, "\n")
cat("Implied sample size at 4% prevalence:", epv_n, "\n")
cat("\n--- Riley criteria (pmsampsize) ---\n")
cat("Events required:", ceiling(ss$events), "\n")
cat("Sample size:", ss$sample_size, "\n")
cat(
"\nRiley / EPV ratio:",
round(ss$sample_size / epv_n, 2),
"times the EPV recommendation\n"
)
# Note what the low prevalence does. At 4%, events are expensive: every extra
# event costs 25 pregnancies. That is why the required sample size is large
# even though the number of parameters is modest, and it is the argument for
# reducing the candidate predictor list before recruitment rather than after.Code
"""Exercise 1: Sample size calculation.
Pre-eclampsia model: prevalence 4%, 12 candidate predictors, an anticipated
C-statistic of 0.72 from a published model in a comparable population.
There is no Python port of pmsampsize, so both halves of what the R package
does are implemented here: the C-statistic to Cox-Snell R-squared conversion,
and the three Riley criteria. The numbers are checked against pmsampsize at
the bottom.
"""
import numpy as np
import statsmodels.api as sm
from scipy.stats import norm
C_STATISTIC = 0.72
PARAMETERS = 12
PREVALENCE = 0.04
def cstat_to_cs_rsquared(cstatistic, prevalence, n=1_000_000, seed=123456):
"""Convert an anticipated C-statistic into a Cox-Snell R-squared.
This mirrors what pmsampsize does internally. There is no closed form, so
it simulates a large population whose linear predictor separates events
from non-events by exactly the amount implied by the C-statistic, fits a
logistic regression to it, and reads off the Cox-Snell R-squared.
The separation is mu = sqrt(2) * qnorm(C): under two unit-variance normals
one mu apart, the probability that a randomly chosen event scores above a
randomly chosen non-event is exactly the C-statistic.
"""
rng = np.random.default_rng(seed)
mu = np.sqrt(2) * norm.ppf(cstatistic)
n0 = int(prevalence * n)
n1 = int((1 - prevalence) * n)
lp = np.concatenate([rng.normal(0.0, 1.0, n0), rng.normal(mu, 1.0, n1)])
y = np.concatenate([np.zeros(n0), np.ones(n1)])
fit = sm.Logit(y, sm.add_constant(lp)).fit(disp=0)
# Cox-Snell R-squared = 1 - exp(-(null deviance - model deviance) / n).
# statsmodels reports log-likelihoods; deviance = -2 * log-likelihood.
lr_stat = 2 * (fit.llf - fit.llnull)
return 1 - np.exp(-lr_stat / len(y))
def max_cs_rsquared(prevalence):
"""The largest Cox-Snell R-squared attainable at this prevalence.
Cox-Snell cannot reach 1 for a binary outcome; the ceiling depends only on
how common the outcome is. Nagelkerke R-squared is Cox-Snell divided by it.
"""
phi = prevalence
return 1 - (phi**phi * (1 - phi) ** (1 - phi)) ** 2
def riley_sample_size(cs_rsquared, parameters, prevalence,
shrinkage=0.9, delta=0.05, moe=0.05):
"""The three Riley et al. criteria for a binary-outcome model."""
phi = prevalence
max_r2 = max_cs_rsquared(phi)
# Criterion 1: expected shrinkage of at least `shrinkage` (0.9 = at most
# 10% overfitting)
n1 = parameters / ((shrinkage - 1) * np.log(1 - cs_rsquared / shrinkage))
# Criterion 2: apparent and adjusted Nagelkerke R-squared differ by <= delta
s2 = cs_rsquared / (cs_rsquared + delta * max_r2)
n2 = parameters / ((s2 - 1) * np.log(1 - cs_rsquared / s2))
# Criterion 3: the overall risk itself estimated to within +/- moe
n3 = (1.96 / moe) ** 2 * phi * (1 - phi)
return {
"Criterion 1 (shrinkage >= 0.9)": int(np.ceil(n1)),
"Criterion 2 (R-squared gap <= 0.05)": int(np.ceil(n2)),
"Criterion 3 (risk within +/- 0.05)": int(np.ceil(n3)),
}
# (a) Minimum sample size ----------------------------------------------------
r2cs = cstat_to_cs_rsquared(C_STATISTIC, PREVALENCE)
max_r2 = max_cs_rsquared(PREVALENCE)
print(f"C-statistic {C_STATISTIC} at {PREVALENCE:.0%} prevalence")
print(f" Cox-Snell R-squared: {r2cs:.4f}")
print(f" Maximum possible: {max_r2:.4f}")
print(f" Nagelkerke equivalent: {r2cs / max_r2:.4f}\n")
criteria = riley_sample_size(r2cs, PARAMETERS, PREVALENCE)
for name, n in criteria.items():
print(f" {name:38s} n = {n}")
n_min = max(criteria.values())
events = int(np.ceil(n_min * PREVALENCE))
print(f"\nMinimum sample size: {n_min} pregnancies")
print(f"Minimum number of events: {events}")
print(f"Events per parameter: {events / PARAMETERS:.2f}")
# (b) Which criterion is binding? -------------------------------------------
binding = max(criteria, key=criteria.get)
print(f"\nBinding criterion: {binding}")
print("Criterion 1 almost always binds for a binary outcome: it is the one "
"that limits overfitting,\nwhile criterion 3 only asks that the average "
"risk be pinned down, which needs far fewer patients.")
# (c) Comparison with the 10-EPV rule of thumb -----------------------------
epv_events = 10 * PARAMETERS
epv_n = int(np.ceil(epv_events / PREVALENCE))
print("\n--- 10 events per variable rule ---")
print(f"Events required: {epv_events}")
print(f"Implied sample size at {PREVALENCE:.0%} prevalence: {epv_n}")
print("\n--- Riley criteria ---")
print(f"Events required: {events}")
print(f"Sample size: {n_min}")
print(f"\nRiley / EPV ratio: {n_min / epv_n:.2f} times the EPV recommendation")
# Check against pmsampsize -------------------------------------------------
# pmsampsize(type = "b", cstatistic = 0.72, parameters = 12, prevalence = 0.04)
# reports Cox-Snell R-squared 0.0251 and n = 4243 / 825 / 60 for the three
# criteria. The conversion above is simulation-based and Python's random
# numbers are not R's, so agreement to the nearest few patients is what to
# expect, not an exact match.
print("\n--- agreement with pmsampsize (R) ---")
for label, ours, theirs in [
("Cox-Snell R-squared", round(r2cs, 4), 0.0251),
("Criterion 1", criteria["Criterion 1 (shrinkage >= 0.9)"], 4243),
("Criterion 2", criteria["Criterion 2 (R-squared gap <= 0.05)"], 825),
("Criterion 3", criteria["Criterion 3 (risk within +/- 0.05)"], 60),
]:
print(f" {label:22s} python {ours:<10} pmsampsize {theirs}")
TipExercise 2: Missing data simulation
Take the simulated Framingham dataset above and compare the regression coefficients obtained from (a) complete case analysis, (b) single mean imputation, and (c) multiple imputation.
- Which approach yields coefficients closest to the true values used to simulate the data?
- What happens to the standard errors under single imputation, and why is that the more serious problem?
NoteSolution
Code
# Exercise 2: Missing data simulation
# Compare complete case analysis, single mean imputation, and multiple
# imputation on the simulated Framingham cohort from the chapter.
library(mice)
# --- The cohort, exactly as in the chapter ---------------------------------
set.seed(2024)
n <- 2000
framingham <- data.frame(
age = round(runif(n, 30, 74)),
male = rbinom(n, 1, 0.48),
sbp = round(rnorm(n, 130, 18)),
total_chol = round(rnorm(n, 210, 38)),
hdl_chol = round(rnorm(n, 52, 15)),
smoking = rbinom(n, 1, 0.22),
diabetes = rbinom(n, 1, 0.08),
bp_treatment = rbinom(n, 1, 0.15)
)
# The true coefficients. Because we simulated the data we know them, which is
# what makes this comparison possible at all.
truth <- c(
"(Intercept)" = -7.5, age = 0.06, male = 0.4, sbp = 0.012,
total_chol = 0.005, hdl_chol = -0.02, smoking = 0.5,
diabetes = 0.7, bp_treatment = 0.3
)
lp <- with(
framingham,
-7.5 + 0.06 * age + 0.4 * male + 0.012 * sbp + 0.005 * total_chol -
0.02 * hdl_chol + 0.5 * smoking + 0.7 * diabetes + 0.3 * bp_treatment
)
framingham$cvd_10yr <- rbinom(n, 1, plogis(lp))
cat("Cohort:", n, "patients,", sum(framingham$cvd_10yr), "events\n")
model_formula <- cvd_10yr ~ age + male + sbp + total_chol + hdl_chol +
smoking + diabetes + bp_treatment
# --- Make two predictors missing, under MAR -------------------------------
# The mechanism matters more than the amount. Missingness in total_chol is made
# to depend on the OUTCOME and on age -- both recorded, so this is MAR, not
# MNAR. That choice is deliberate: if missingness depended only on the
# predictors already in the model, complete case analysis would still be
# unbiased for the coefficients, and there would be nothing to see.
punch_holes <- function(df, intercept) {
set.seed(7)
p_chol <- plogis(intercept + 1.0 * df$cvd_10yr + 0.03 * (df$age - 52))
p_hdl <- plogis(intercept + 0.2 + 0.9 * df$smoking + 0.02 * (df$sbp - 130))
df$total_chol[runif(nrow(df)) < p_chol] <- NA
df$hdl_chol[runif(nrow(df)) < p_hdl] <- NA
df
}
incomplete <- punch_holes(framingham, intercept = -2.2)
cat("Missing total_chol:", sum(is.na(incomplete$total_chol)), "\n")
cat("Missing hdl_chol: ", sum(is.na(incomplete$hdl_chol)), "\n")
cat(
"Complete rows:", sum(complete.cases(incomplete)),
sprintf("(%.0f%% of the cohort discarded by a complete case analysis)\n",
100 * mean(!complete.cases(incomplete)))
)
# --- The three approaches, plus the full data as a benchmark --------------
fit_full <- glm(model_formula, data = framingham, family = binomial)
fit_cca <- glm(model_formula, data = incomplete, family = binomial)
mean_imputed <- incomplete
for (v in c("total_chol", "hdl_chol")) {
mean_imputed[[v]][is.na(mean_imputed[[v]])] <- mean(mean_imputed[[v]], na.rm = TRUE)
}
fit_mean <- glm(model_formula, data = mean_imputed, family = binomial)
# 20 imputations, predictive mean matching, pooled with Rubin's rules
imp <- mice(incomplete, m = 20, method = "pmm", seed = 42, printFlag = FALSE)
pool_obj <- pool(with(
imp,
glm(cvd_10yr ~ age + male + sbp + total_chol + hdl_chol + smoking +
diabetes + bp_treatment, family = binomial)
))
pooled <- summary(pool_obj)
mi_est <- setNames(pooled$estimate, pooled$term)
mi_se <- setNames(pooled$std.error, pooled$term)
# lambda is the share of the pooled variance that comes from the missing data
# (the between-imputation part). It is precisely what single imputation drops.
mi_lambda <- setNames(pool_obj$pooled$lambda, pool_obj$pooled$term)
# --- (a) Which approach lands closest to the truth? ----------------------
comparison <- data.frame(
truth = truth,
full_data = coef(fit_full),
complete_case = coef(fit_cca),
mean_imputation = coef(fit_mean),
multiple_imputation = mi_est[names(truth)]
)
cat("\n--- Coefficient estimates ---\n")
print(round(comparison, 4))
# Two error measures, and the difference between them is the point. Distance
# from the truth mixes up two things: damage done by the missing data, and the
# sampling noise that was already in this cohort of 2000. Distance from the
# full-data estimates isolates the first.
slopes <- rownames(comparison) != "(Intercept)"
mae <- function(x) mean(abs(x - comparison$truth[slopes]))
mae_full <- function(x) mean(abs(x - comparison$full_data[slopes]))
cat("\n--- Mean absolute error across the 8 slopes ---\n")
err <- data.frame(
vs_truth = sapply(comparison[slopes, -1], mae),
vs_full_data = sapply(comparison[slopes, -1], mae_full)
)
print(round(err, 5))
cat("\nRanked by distance from the full-data estimates:\n")
print(round(sort(err$vs_full_data[-1] |> setNames(rownames(err)[-1])), 5))
# --- (b) What happens to the standard errors? ---------------------------
se_table <- data.frame(
complete_case = summary(fit_cca)$coefficients[, "Std. Error"],
mean_imputation = summary(fit_mean)$coefficients[, "Std. Error"],
multiple_imputation = mi_se[names(truth)]
)
cat("\n--- Standard errors ---\n")
print(round(se_table, 5))
cat("\nMean imputation SE as a percentage of the multiple-imputation SE,\n")
cat("for the two variables that were actually imputed:\n")
for (v in c("total_chol", "hdl_chol")) {
cat(sprintf(
" %-11s %.1f%%\n", v,
100 * se_table[v, "mean_imputation"] / se_table[v, "multiple_imputation"]
))
}
# What single imputation actually discards, as a number. Rubin's rules split the
# pooled variance into a within-imputation part (ordinary sampling uncertainty)
# and a between-imputation part (uncertainty about the guesses themselves).
# Single imputation sets the second to zero by construction.
cat("\nShare of the pooled MI variance that comes from the imputation itself:\n")
for (v in c("total_chol", "hdl_chol")) {
cat(sprintf(" %-11s %.1f%%\n", v, 100 * mi_lambda[[v]]))
}
# Two opposing effects on the mean-imputation standard error, which is why the
# percentages above are close to 100 at this fraction of missing data:
# 1. it ignores the imputation uncertainty just quantified -> SE too small
# 2. it flattens the variable's spread, and less spread in a
# predictor means less information about its coefficient -> SE too large
# The first grows with the fraction missing; the second is the reason the two
# can briefly cancel. Neither makes the SE trustworthy.
cat(sprintf(
"\nSD of total_chol: %.1f observed -> %.1f after mean imputation\n",
sd(incomplete$total_chol, na.rm = TRUE), sd(mean_imputed$total_chol)
))
# --- The same comparison with far more missing data ---------------------
# At 12-16% missing the three approaches disagree modestly. Raise the
# missingness to roughly half and the false precision of single imputation
# becomes impossible to miss.
heavy <- punch_holes(framingham, intercept = -0.4)
heavy_mean <- heavy
for (v in c("total_chol", "hdl_chol")) {
heavy_mean[[v]][is.na(heavy_mean[[v]])] <- mean(heavy_mean[[v]], na.rm = TRUE)
}
imp_h <- mice(heavy, m = 20, method = "pmm", seed = 42, printFlag = FALSE)
pool_h <- pool(with(
imp_h,
glm(cvd_10yr ~ age + male + sbp + total_chol + hdl_chol + smoking +
diabetes + bp_treatment, family = binomial)
))
pooled_h <- summary(pool_h)
lambda_h <- setNames(pool_h$pooled$lambda, pool_h$pooled$term)
cat(sprintf(
"\n--- With %.0f%% of total_chol and %.0f%% of hdl_chol missing ---\n",
100 * mean(is.na(heavy$total_chol)), 100 * mean(is.na(heavy$hdl_chol))
))
s_cca_h <- summary(glm(model_formula, data = heavy, family = binomial))$coefficients
s_mean_h <- summary(glm(model_formula, data = heavy_mean, family = binomial))$coefficients
for (v in c("total_chol", "hdl_chol")) {
cat(sprintf("%-11s truth %+.4f\n", v, truth[[v]]))
cat(sprintf(
" complete case %+.4f (SE %.5f)\n mean imputation %+.4f (SE %.5f)\n MI %+.4f (SE %.5f)\n",
s_cca_h[v, 1], s_cca_h[v, 2], s_mean_h[v, 1], s_mean_h[v, 2],
pooled_h$estimate[pooled_h$term == v],
pooled_h$std.error[pooled_h$term == v]
))
cat(sprintf(
" -> mean imputation's SE is %.0f%% of MI's, and %.0f%% of MI's variance\n for this coefficient now comes from the imputation (was %.0f%%)\n",
100 * s_mean_h[v, 2] / pooled_h$std.error[pooled_h$term == v],
100 * lambda_h[[v]], 100 * mi_lambda[[v]]
))
}
cat("
Conclusions
-----------
(a) Complete case analysis is the clear loser. It is the furthest from the
full-data estimates, and it is biased by construction here: because
missingness depends on the outcome, the retained rows under-represent
patients who had an event, and coefficients such as smoking are distorted
well beyond sampling noise. It also discards a quarter of the cohort, so
its standard errors are the widest of the three.
Mean imputation and multiple imputation give similar point estimates at
this fraction of missing data. Note also that measuring against the true
values alone is misleading: the diabetes coefficient is far from its true
0.7 in every column, including the full-data one, because 2000 patients
and 235 events cannot pin it down. That error is sampling noise, not
missing-data handling.
(b) Filling every gap with the mean asserts that those values were measured
rather than guessed, so nothing in the model widens to reflect the
guessing. Rubin's rules make visible exactly what is being thrown away:
the printed lambda says what share of the pooled uncertainty comes from
the imputation itself, and single imputation sets that share to zero.
That does not always show up as a smaller standard error, and the output
above is a good reminder to check rather than assume. Two effects pull in
opposite directions -- ignoring the imputation uncertainty makes the
standard error too small, while flattening the variable's spread (SD 37.3
to 35.0) makes it too large -- and at 12-16% missing they nearly cancel,
leaving mean imputation within 2% of the multiple-imputation standard
error. Raise the missingness to roughly half and the first effect wins
outright: mean imputation's standard error for total_chol is about two
thirds of the honest one.
Watch the share itself rather than the ratio, because the share is the
part that behaves predictably: it climbs from under a fifth to two thirds
or more as the missingness grows. The fair statement is not that single
imputation always looks more precise, but that its uncertainty is
unaccounted for, and the size of what it ignores grows with the amount you
imputed.
Which method lands nearest the truth in any one dataset is luck; all of
them sit within a standard error of each other. The missing variance
component is systematic.
That is the more serious problem: a biased estimate with an honest
confidence interval announces its own uncertainty, whereas a spuriously
precise one invites a confident claim about a coefficient the data cannot
support. Multiple imputation exists to keep that uncertainty visible.
One caveat: multiple imputation assumes MAR, which holds here by
construction. If the highest cholesterol values were missing precisely
because they were high (MNAR), no method here would recover them, and the
honest response would be a sensitivity analysis.
")Code
"""Exercise 2: Missing data simulation.
Compare complete case analysis, single mean imputation, and multiple
imputation on the simulated Framingham cohort from the chapter.
Python's random numbers are not R's, so the numbers here differ from the R
solution in detail. The pattern -- which approach fails, and how -- is the same.
"""
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy.special import expit
pd.set_option("display.width", 120)
pd.set_option("display.max_columns", 20)
PREDICTORS = ["age", "male", "sbp", "total_chol", "hdl_chol",
"smoking", "diabetes", "bp_treatment"]
# The true coefficients: we simulated the data, so we know them.
TRUTH = pd.Series({"const": -7.5, "age": 0.06, "male": 0.4, "sbp": 0.012,
"total_chol": 0.005, "hdl_chol": -0.02, "smoking": 0.5,
"diabetes": 0.7, "bp_treatment": 0.3})
# --- The cohort, as in the chapter -----------------------------------------
rng = np.random.default_rng(2024)
n = 2000
framingham = pd.DataFrame({
"age": rng.integers(30, 75, n),
"male": rng.binomial(1, 0.48, n),
"sbp": np.round(rng.normal(130, 18, n)),
"total_chol": np.round(rng.normal(210, 38, n)),
"hdl_chol": np.round(rng.normal(52, 15, n)),
"smoking": rng.binomial(1, 0.22, n),
"diabetes": rng.binomial(1, 0.08, n),
"bp_treatment": rng.binomial(1, 0.15, n),
})
lp = (-7.5 + 0.06 * framingham["age"] + 0.4 * framingham["male"]
+ 0.012 * framingham["sbp"] + 0.005 * framingham["total_chol"]
- 0.02 * framingham["hdl_chol"] + 0.5 * framingham["smoking"]
+ 0.7 * framingham["diabetes"] + 0.3 * framingham["bp_treatment"])
framingham["cvd_10yr"] = rng.binomial(1, expit(lp))
print(f"Cohort: {len(framingham)} patients, "
f"{framingham['cvd_10yr'].sum()} events")
def fit_logit(df):
"""Unpenalised logistic regression, returning coefficients and SEs."""
X = sm.add_constant(df[PREDICTORS].astype(float))
res = sm.Logit(df["cvd_10yr"], X).fit(disp=0)
return res.params, res.bse
def punch_holes(df, intercept, seed=7):
"""Delete values of two predictors under a MAR mechanism.
Missingness in total_chol depends on the OUTCOME and on age, both of which
are recorded -- so this is MAR, not MNAR. That is deliberate: if
missingness depended only on predictors already in the model, complete case
analysis would still be unbiased for the coefficients and there would be
nothing to demonstrate.
"""
out = df.copy()
r = np.random.default_rng(seed)
p_chol = expit(intercept + 1.0 * out["cvd_10yr"] + 0.03 * (out["age"] - 52))
p_hdl = expit(intercept + 0.2 + 0.9 * out["smoking"]
+ 0.02 * (out["sbp"] - 130))
out.loc[r.uniform(size=len(out)) < p_chol, "total_chol"] = np.nan
out.loc[r.uniform(size=len(out)) < p_hdl, "hdl_chol"] = np.nan
return out
def mean_impute(df):
out = df.copy()
for v in ["total_chol", "hdl_chol"]:
out[v] = out[v].fillna(out[v].mean())
return out
def impute_once(df, rng, iterations=5):
"""Create ONE filled-in dataset by chained equations.
This is written out rather than delegated to a package because it is short
enough to read, and reading it is the point: each missing value is drawn
from a regression on the other variables, with fresh randomness every time,
which is what makes the m datasets differ from each other.
Two details decide whether it works at all:
* The draw includes both parameter uncertainty (a draw of sigma^2, then of
beta given sigma^2) and residual noise. Using the fitted value alone
would make every dataset identical -- single imputation in disguise.
* The OUTCOME is one of the predictors in the imputation model. Leaving it
out imputes predictors as though they were unrelated to the outcome,
which biases their coefficients towards zero. R's `mice` includes every
column by default, which is why the R solution just hands it the whole
data frame.
In practice use `mice` in R or `miceforest` in Python. Note that
scikit-learn's `IterativeImputer` is not a drop-in substitute: even with
`sample_posterior=True` its imputations are under-dispersed, which
understates the between-imputation variance this exercise is about.
"""
columns = PREDICTORS + ["cvd_10yr"]
work = df[columns].astype(float).copy()
incomplete = [v for v in columns if work[v].isna().any()]
missing = {v: work[v].isna().to_numpy() for v in incomplete}
for v in incomplete: # start each variable off at its observed mean
work.loc[missing[v], v] = work[v].mean()
for _ in range(iterations):
for v in incomplete:
others = [c for c in columns if c != v]
observed = ~missing[v]
X = np.column_stack([np.ones(observed.sum()),
work.loc[observed, others].to_numpy()])
y = work.loc[observed, v].to_numpy()
XtX_inv = np.linalg.pinv(X.T @ X)
beta_hat = XtX_inv @ X.T @ y
residuals = y - X @ beta_hat
dof = max(len(y) - X.shape[1], 1)
# Posterior draws, in the order mice's "norm" method uses them
sigma2 = residuals @ residuals / rng.chisquare(dof)
chol = np.linalg.cholesky(sigma2 * XtX_inv
+ 1e-12 * np.eye(len(beta_hat)))
beta = beta_hat + chol @ rng.standard_normal(len(beta_hat))
X_missing = np.column_stack([
np.ones(missing[v].sum()),
work.loc[missing[v], others].to_numpy(),
])
work.loc[missing[v], v] = (
X_missing @ beta
+ rng.normal(0, np.sqrt(sigma2), missing[v].sum())
)
work["cvd_10yr"] = df["cvd_10yr"].to_numpy() # never imputed
return work
def multiple_impute(df, m=20, seed=42):
"""Impute m times, fit the model in each, pool with Rubin's rules."""
rng = np.random.default_rng(seed)
estimates, variances = [], []
for _ in range(m):
params, ses = fit_logit(impute_once(df, rng))
estimates.append(params)
variances.append(ses**2)
estimates = pd.DataFrame(estimates)
variances = pd.DataFrame(variances)
# Rubin's rules: the pooled variance is the average within-imputation
# variance plus the between-imputation variance, inflated by 1 + 1/m.
within = variances.mean()
between = estimates.var(ddof=1)
total = within + (1 + 1 / m) * between
return estimates.mean(), np.sqrt(total), within, between
incomplete = punch_holes(framingham, intercept=-2.2)
complete_rows = incomplete[PREDICTORS].notna().all(axis=1)
print(f"Missing total_chol: {incomplete['total_chol'].isna().sum()}")
print(f"Missing hdl_chol: {incomplete['hdl_chol'].isna().sum()}")
print(f"Complete rows: {complete_rows.sum()} "
f"({100 * (1 - complete_rows.mean()):.0f}% of the cohort discarded by a "
"complete case analysis)")
# --- The three approaches, plus the full data as a benchmark -------------
b_full, se_full = fit_logit(framingham)
b_cca, se_cca = fit_logit(incomplete[complete_rows])
b_mean, se_mean = fit_logit(mean_impute(incomplete))
b_mi, se_mi, within, between = multiple_impute(incomplete)
comparison = pd.DataFrame({
"truth": TRUTH,
"full_data": b_full,
"complete_case": b_cca,
"mean_imputation": b_mean,
"multiple_imputation": b_mi,
})[["truth", "full_data", "complete_case", "mean_imputation",
"multiple_imputation"]]
print("\n--- Coefficient estimates ---")
print(comparison.round(4))
# --- (a) Which approach lands closest to the truth? ---------------------
# Two error measures, and the difference between them is the point. Distance
# from the truth mixes up two things: damage done by the missing data, and the
# sampling noise already present in this cohort of 2000. Distance from the
# full-data estimates isolates the first.
slopes = comparison.drop(index="const")
errors = pd.DataFrame({
"vs_truth": (slopes.drop(columns="truth")
.sub(slopes["truth"], axis=0).abs().mean()),
"vs_full_data": (slopes.drop(columns=["truth", "full_data"])
.sub(slopes["full_data"], axis=0).abs().mean()),
})
print("\n--- Mean absolute error across the 8 slopes ---")
print(errors.round(5))
print("\nRanked by distance from the full-data estimates:")
print(errors["vs_full_data"].dropna().sort_values().round(5))
# --- (b) What happens to the standard errors? --------------------------
se_table = pd.DataFrame({
"complete_case": se_cca,
"mean_imputation": se_mean,
"multiple_imputation": se_mi,
})
print("\n--- Standard errors ---")
print(se_table.round(5))
print("\nMean imputation SE as a percentage of the multiple-imputation SE,")
print("for the two variables that were actually imputed:")
for v in ["total_chol", "hdl_chol"]:
print(f" {v:11s} {100 * se_mean[v] / se_mi[v]:.1f}%")
# What single imputation actually discards, as a number. Rubin's rules split the
# pooled variance into a within-imputation part (ordinary sampling uncertainty)
# and a between-imputation part (uncertainty about the guesses themselves).
# Single imputation sets the second to zero by construction. This share is what
# mice reports as lambda.
print("\nShare of the pooled MI variance that comes from the imputation itself:")
for v in ["total_chol", "hdl_chol"]:
share = (1 + 1 / 20) * between[v] / (within[v] + (1 + 1 / 20) * between[v])
print(f" {v:11s} {100 * share:.1f}%")
# Two opposing effects act on the mean-imputation standard error, which is why
# the percentages above sit close to 100 at this fraction of missing data:
# 1. it ignores the imputation uncertainty just quantified -> SE too small
# 2. it flattens the variable's spread, and less spread in a
# predictor means less information about its coefficient -> SE too large
# The first grows with the fraction missing; the second is why they can briefly
# cancel, and one of the two variables below may even come out above 100%.
print(f"\nSD of total_chol: {incomplete['total_chol'].std():.1f} observed"
f" -> {mean_impute(incomplete)['total_chol'].std():.1f}"
" after mean imputation")
# --- The same comparison with far more missing data -------------------
heavy = punch_holes(framingham, intercept=-0.4)
heavy_rows = heavy[PREDICTORS].notna().all(axis=1)
b_cca_h, se_cca_h = fit_logit(heavy[heavy_rows])
b_mean_h, se_mean_h = fit_logit(mean_impute(heavy))
b_mi_h, se_mi_h, within_h, between_h = multiple_impute(heavy)
print(f"\n--- With {100 * heavy['total_chol'].isna().mean():.0f}% of total_chol"
f" and {100 * heavy['hdl_chol'].isna().mean():.0f}% of hdl_chol missing ---")
for v in ["total_chol", "hdl_chol"]:
print(f"{v:11s} truth {TRUTH[v]:+.4f}")
print(f" complete case {b_cca_h[v]:+.4f} (SE {se_cca_h[v]:.5f})")
print(f" mean imputation {b_mean_h[v]:+.4f} (SE {se_mean_h[v]:.5f})")
print(f" MI {b_mi_h[v]:+.4f} (SE {se_mi_h[v]:.5f})")
share_now = ((1 + 1 / 20) * between_h[v]
/ (within_h[v] + (1 + 1 / 20) * between_h[v]))
share_was = ((1 + 1 / 20) * between[v]
/ (within[v] + (1 + 1 / 20) * between[v]))
print(f" -> mean imputation's SE is "
f"{100 * se_mean_h[v] / se_mi_h[v]:.0f}% of MI's, and "
f"{100 * share_now:.0f}% of MI's variance")
print(f" for this coefficient now comes from the imputation "
f"(was {100 * share_was:.0f}%)")
print("""
Conclusions
-----------
(a) Complete case analysis is the clear loser. It sits furthest from the
full-data estimates, and it is biased by construction here: because
missingness depends on the outcome, the retained rows under-represent
patients who had an event. It also discards a quarter of the cohort, so
its standard errors are the widest of the three.
Mean imputation and multiple imputation give similar point estimates at
this fraction of missing data. Measuring against the true values alone is
misleading, though: the diabetes coefficient is far from its true 0.7 in
every column, including the full-data one, because 2000 patients and
roughly 200 events cannot pin it down. That error is sampling noise, not
missing-data handling -- which is why the comparison against the full-data
estimates is the informative one.
(b) Filling every gap with the mean asserts those values were measured rather
than guessed, so nothing in the model widens to reflect the guessing.
Rubin's rules make visible exactly what is discarded: the printed share
says how much of the pooled uncertainty comes from the imputation itself,
and single imputation sets that share to zero.
That does not always surface as a smaller standard error, and the output
above is a good reminder to check rather than assume. Two effects pull in
opposite directions -- ignoring the imputation uncertainty makes the
standard error too small, while flattening the variable's spread makes it
too large -- so at this fraction of missing data they nearly cancel, and
with half the values missing one variable still comes out slightly above
the multiple-imputation standard error while the other falls to about two
thirds of it.
Watch the share rather than the ratio, because the share is the part that
behaves predictably: it climbs from a fifth or a quarter to a half or more
as the missingness grows. The fair statement is not that single imputation always
looks more precise, but that its uncertainty is unaccounted for, and the
size of what it ignores grows with the amount you imputed.
Which method lands nearest the truth in any one dataset is luck; all of
them sit within a standard error of each other. The missing variance
component is systematic.
One caveat: multiple imputation assumes MAR, which holds here by
construction. If the highest cholesterol values were missing precisely
because they were high (MNAR), no method here would recover them, and the
honest response would be a sensitivity analysis.
""")
TipExercise 3: Variable selection
Fit the Framingham model using (a) all pre-specified predictors, (b) forward stepwise selection, and (c) LASSO.
- Compare which variables each approach selects. Then add a dozen pure-noise candidates to the list and compare again — selection only has work to do when some candidates do not belong.
- Repeat the selection across 100 bootstrap samples and record how often each variable is chosen. How stable is stepwise selection compared with LASSO?
- What does that instability imply for a paper reporting a single stepwise-selected model?
NoteSolution
Code
# Exercise 3: Variable selection
# Compare a pre-specified model, forward stepwise selection, and LASSO on the
# simulated Framingham cohort, then check how stable each selection is across
# 100 bootstrap samples.
#
# Runtime is about a minute: the bootstrap refits both selection procedures
# 100 times, cross-validating the LASSO penalty inside each one.
library(MASS) # stepAIC
library(glmnet) # cv.glmnet
# --- The cohort, exactly as in the chapter ---------------------------------
set.seed(2024)
n <- 2000
framingham <- data.frame(
age = round(runif(n, 30, 74)),
male = rbinom(n, 1, 0.48),
sbp = round(rnorm(n, 130, 18)),
total_chol = round(rnorm(n, 210, 38)),
hdl_chol = round(rnorm(n, 52, 15)),
smoking = rbinom(n, 1, 0.22),
diabetes = rbinom(n, 1, 0.08),
bp_treatment = rbinom(n, 1, 0.15)
)
lp <- with(
framingham,
-7.5 + 0.06 * age + 0.4 * male + 0.012 * sbp + 0.005 * total_chol -
0.02 * hdl_chol + 0.5 * smoking + 0.7 * diabetes + 0.3 * bp_treatment
)
framingham$cvd_10yr <- rbinom(n, 1, plogis(lp))
REAL <- setdiff(names(framingham), "cvd_10yr")
cat("Events:", sum(framingham$cvd_10yr), "of", n, "\n")
# --- Helpers ---------------------------------------------------------------
forward_stepwise <- function(data, candidates) {
upper <- as.formula(paste("~", paste(candidates, collapse = " + ")))
fit <- stepAIC(
glm(cvd_10yr ~ 1, data = data, family = binomial),
scope = list(lower = ~1, upper = upper),
direction = "forward", trace = 0
)
setdiff(names(coef(fit)), "(Intercept)")
}
lasso_selected <- function(data, candidates, s = "lambda.min", nfolds = 10) {
X <- as.matrix(data[candidates])
cvfit <- cv.glmnet(X, data$cvd_10yr, family = "binomial",
alpha = 1, nfolds = nfolds)
b <- coef(cvfit, s = s)[-1, 1]
names(b)[b != 0]
}
# --- (a) The three approaches on the pre-specified predictors -------------
cat("\n=== Part (a): candidates are the 8 pre-specified predictors ===\n")
fit_all <- glm(cvd_10yr ~ ., data = framingham, family = binomial)
cat("\nPre-specified model (all 8 predictors):\n")
print(round(summary(fit_all)$coefficients[, c(1, 2, 4)], 4))
set.seed(1)
step_a <- forward_stepwise(framingham, REAL)
lasso_min_a <- lasso_selected(framingham, REAL, "lambda.min")
lasso_1se_a <- lasso_selected(framingham, REAL, "lambda.1se")
cat("\nSelected sets:\n")
cat(" pre-specified (8):", paste(REAL, collapse = ", "), "\n")
cat(sprintf(" forward stepwise (%d): %s\n", length(step_a),
paste(sort(step_a), collapse = ", ")))
cat(sprintf(" LASSO lambda.min (%d): %s\n", length(lasso_min_a),
paste(sort(lasso_min_a), collapse = ", ")))
cat(sprintf(" LASSO lambda.1se (%d): %s\n", length(lasso_1se_a),
paste(sort(lasso_1se_a), collapse = ", ")))
cat("
With 8 genuinely predictive variables and 235 events there is nothing for
selection to remove: every predictor is significant, so stepwise keeps all
eight and so does the LASSO at lambda.min. This is the easy case, and the
honest conclusion is that selection added nothing -- the pre-specified model
was already the answer. Note that lambda.1se, the deliberately conservative
choice, drops real predictors.
")
# --- (b) The realistic case: candidates that do not belong ---------------
# Selection only becomes interesting when the candidate list contains
# variables with no relationship to the outcome, which is the usual situation
# when a list is drawn up from "everything we measured".
set.seed(99)
noise <- as.data.frame(matrix(rnorm(n * 12), n, 12))
names(noise) <- paste0("noise", 1:12)
wide <- cbind(framingham, noise)
CANDIDATES <- c(REAL, names(noise))
cat("\n=== Part (b): 8 real predictors + 12 pure-noise candidates ===\n")
set.seed(1)
step_b <- forward_stepwise(wide, CANDIDATES)
lasso_min_b <- lasso_selected(wide, CANDIDATES, "lambda.min")
lasso_1se_b <- lasso_selected(wide, CANDIDATES, "lambda.1se")
report <- function(label, selected) {
cat(sprintf(
" %-18s kept %2d (%d of 8 real, %d of 12 noise)\n",
label, length(selected), sum(selected %in% REAL),
sum(grepl("noise", selected))
))
}
cat("\n")
report("forward stepwise", step_b)
report("LASSO lambda.min", lasso_min_b)
report("LASSO lambda.1se", lasso_1se_b)
cat("\n stepwise kept these noise variables:",
paste(sort(grep("noise", step_b, value = TRUE)), collapse = ", "), "\n")
# --- (c) Stability across 100 bootstrap samples -------------------------
cat("\n=== Part (c): selection frequency across 100 bootstrap samples ===\n")
set.seed(2025)
B <- 100
counts <- list(
stepwise = setNames(numeric(length(CANDIDATES)), CANDIDATES),
lasso_min = setNames(numeric(length(CANDIDATES)), CANDIDATES)
)
sizes <- list(stepwise = numeric(B), lasso_min = numeric(B))
signatures <- list(stepwise = character(B), lasso_min = character(B))
for (b in seq_len(B)) {
boot <- wide[sample(n, n, replace = TRUE), ]
picks <- list(
stepwise = forward_stepwise(boot, CANDIDATES),
# nfolds = 5 inside the loop purely to keep the runtime reasonable
lasso_min = lasso_selected(boot, CANDIDATES, "lambda.min", nfolds = 5)
)
for (m in names(picks)) {
counts[[m]][picks[[m]]] <- counts[[m]][picks[[m]]] + 1
sizes[[m]][b] <- length(picks[[m]])
signatures[[m]][b] <- paste(sort(picks[[m]]), collapse = "|")
}
}
freq <- data.frame(
variable = CANDIDATES,
truth = ifelse(CANDIDATES %in% REAL, "real", "noise"),
stepwise_pct = 100 * counts$stepwise / B,
lasso_pct = 100 * counts$lasso_min / B,
row.names = NULL
)
freq <- freq[order(freq$truth, -freq$stepwise_pct), ]
cat("\nHow often each candidate was selected (%):\n")
print(freq, row.names = FALSE)
cat("\nSummary across the 100 bootstrap samples:\n")
for (m in c("stepwise", "lasso_min")) {
real_pct <- mean(counts[[m]][REAL]) / B * 100
noise_pct <- mean(counts[[m]][grepl("noise", CANDIDATES)]) / B * 100
cat(sprintf(
" %-10s median size %.0f | real predictors kept %.0f%% of the time | noise kept %.0f%% | %d distinct models\n",
m, median(sizes[[m]]), real_pct, noise_pct, length(unique(signatures[[m]]))
))
}
cat("
Conclusions
-----------
(a) When every candidate belongs in the model, selection has no work to do and
both methods keep everything. Selection cannot improve on a well-chosen
pre-specified list; the most it can do is leave it alone.
(b) Add candidates that do not belong and both methods start letting them in.
Forward stepwise admits several noise variables, because a variable enters
on whether it improves AIC in this particular sample, and with 12 chances
some noise variable always looks helpful. The LASSO at lambda.min is no
better here -- lambda.min optimises prediction error, not variable
recovery, so it keeps most of the noise with small coefficients. Only
lambda.1se is clean of noise, and it pays by dropping real predictors.
(c) The bootstrap frequencies are the real finding. Stepwise produced 96
different models in 100 resamples of the same patients, and the LASSO 75.
Real predictors were kept 88% of the time by stepwise and 97% by the
LASSO; noise variables 37% and 75%. Even the strongest real predictors
(age, hdl_chol, diabetes) are selected every time, so the instability is
concentrated exactly where it matters -- the weaker predictors, where you
would actually want the method's advice.
One subtlety worth naming, because it cuts against the obvious reading.
noise10 is selected 75% of the time by stepwise and 99% by the LASSO,
which looks like a reliable predictor. It is not: it happens to correlate
with the outcome in this particular cohort, and every bootstrap sample is
drawn from that same cohort, so the fluke is reproduced. Bootstrap
stability therefore measures robustness to resampling these patients, not
to collecting new ones -- it is a lower bound on the instability a fresh
dataset would reveal.
That is what makes a paper reporting one stepwise-selected model
misleading. The list of retained variables is presented as a finding --
'these are the predictors of risk' -- when a different sample of the same
patients would have produced a different list, and a genuinely new sample
a different one again. The p-values and confidence intervals are also
wrong, because they take no account of the searching that preceded them.
The defensible options are to pre-specify the predictors on clinical
grounds and keep them all, or to use a penalised model and report the
whole procedure rather than the variables that happened to survive it.
")Code
"""Exercise 3: Variable selection.
Compare a pre-specified model, forward stepwise selection, and LASSO on the
simulated Framingham cohort, then check how stable each selection is across
100 bootstrap samples.
Python's random numbers are not R's, so the numbers differ from the R solution
in detail; the pattern is the same. Runtime is a minute or two, because the
bootstrap refits both selection procedures 100 times.
"""
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy.special import expit
from sklearn.linear_model import LogisticRegression, LogisticRegressionCV
from sklearn.preprocessing import StandardScaler
pd.set_option("display.width", 120)
# --- The cohort, as in the chapter -----------------------------------------
rng = np.random.default_rng(2024)
n = 2000
framingham = pd.DataFrame({
"age": rng.integers(30, 75, n),
"male": rng.binomial(1, 0.48, n),
"sbp": np.round(rng.normal(130, 18, n)),
"total_chol": np.round(rng.normal(210, 38, n)),
"hdl_chol": np.round(rng.normal(52, 15, n)),
"smoking": rng.binomial(1, 0.22, n),
"diabetes": rng.binomial(1, 0.08, n),
"bp_treatment": rng.binomial(1, 0.15, n),
})
lp = (-7.5 + 0.06 * framingham["age"] + 0.4 * framingham["male"]
+ 0.012 * framingham["sbp"] + 0.005 * framingham["total_chol"]
- 0.02 * framingham["hdl_chol"] + 0.5 * framingham["smoking"]
+ 0.7 * framingham["diabetes"] + 0.3 * framingham["bp_treatment"])
framingham["cvd_10yr"] = rng.binomial(1, expit(lp))
REAL = [c for c in framingham.columns if c != "cvd_10yr"]
print(f"Events: {framingham['cvd_10yr'].sum()} of {n}")
# --- Helpers ---------------------------------------------------------------
def log_likelihood(X, y):
"""Maximised log-likelihood of an unpenalised logistic fit."""
if X.shape[1] == 0:
p = np.full(len(y), y.mean())
else:
p = (LogisticRegression(penalty=None, max_iter=1000)
.fit(X, y).predict_proba(X)[:, 1])
p = np.clip(p, 1e-12, 1 - 1e-12)
return np.sum(y * np.log(p) + (1 - y) * np.log(1 - p))
def forward_stepwise(df, candidates, outcome="cvd_10yr"):
"""Forward selection on AIC, the equivalent of R's stepAIC(direction='forward').
The candidates are standardised first purely for speed: the likelihood, and
therefore AIC and every selection decision, is unchanged by rescaling a
predictor, but the optimiser converges far faster when age (30-74) and
total cholesterol (around 210) are on a common scale.
"""
y = df[outcome].to_numpy()
Z = pd.DataFrame(
StandardScaler().fit_transform(df[candidates].astype(float)),
columns=candidates,
)
chosen, remaining = [], list(candidates)
current_aic = 2 * 1 - 2 * log_likelihood(np.empty((len(df), 0)), y)
while remaining:
best, best_aic = None, current_aic
for v in remaining:
X = Z[chosen + [v]].to_numpy()
aic = 2 * (X.shape[1] + 1) - 2 * log_likelihood(X, y)
if aic < best_aic:
best, best_aic = v, aic
if best is None: # nothing improves AIC any further
break
chosen.append(best)
remaining.remove(best)
current_aic = best_aic
return chosen
def lasso_selected(df, candidates, outcome="cvd_10yr", cv=10, one_se=False,
seed=1):
"""Variables kept by a cross-validated LASSO.
Two settings matter and are easy to get wrong (both are flagged in the
chapter): predictors must be standardised, because the penalty acts on the
coefficient scale, and scoring must be a probability score -- the default
accuracy is maximised by predicting the majority class, which would choose
maximum shrinkage and drop everything.
"""
X = StandardScaler().fit_transform(df[candidates].astype(float))
y = df[outcome].to_numpy()
grid = np.logspace(-4, 1, 25)
model = LogisticRegressionCV(
Cs=grid, penalty="l1", solver="liblinear", cv=cv,
scoring="neg_log_loss", max_iter=5000, random_state=seed,
).fit(X, y)
if one_se:
# The equivalent of glmnet's lambda.1se: the strongest penalty whose
# mean CV score is within one standard error of the best.
scores = model.scores_[1] # folds x Cs
mean, se = scores.mean(axis=0), scores.std(axis=0) / np.sqrt(cv)
best = mean.argmax()
ok = np.where(mean >= mean[best] - se[best])[0]
chosen_C = model.Cs_[ok.min()] # smallest C = strongest penalty
model = LogisticRegression(C=chosen_C, penalty="l1", solver="liblinear",
max_iter=5000).fit(X, y)
return [c for c, b in zip(candidates, model.coef_[0]) if b != 0]
def describe(label, selected):
n_real = sum(v in REAL for v in selected)
n_noise = sum("noise" in v for v in selected)
print(f" {label:18s} kept {len(selected):2d} "
f"({n_real} of 8 real, {n_noise} of 12 noise)")
# --- (a) The three approaches on the pre-specified predictors -------------
print("\n=== Part (a): candidates are the 8 pre-specified predictors ===")
X_all = sm.add_constant(framingham[REAL].astype(float))
fit_all = sm.Logit(framingham["cvd_10yr"], X_all).fit(disp=0)
print("\nPre-specified model (all 8 predictors):")
print(pd.DataFrame({"coef": fit_all.params, "se": fit_all.bse,
"p": fit_all.pvalues}).round(4))
step_a = forward_stepwise(framingham, REAL)
lasso_min_a = lasso_selected(framingham, REAL)
lasso_1se_a = lasso_selected(framingham, REAL, one_se=True)
print("\nSelected sets:")
print(f" pre-specified (8): {', '.join(REAL)}")
print(f" forward stepwise ({len(step_a)}): {', '.join(sorted(step_a))}")
print(f" LASSO best lambda ({len(lasso_min_a)}): {', '.join(sorted(lasso_min_a))}")
print(f" LASSO 1-SE lambda ({len(lasso_1se_a)}): {', '.join(sorted(lasso_1se_a))}")
print("""
Every one of these 8 variables belongs in the model -- we put them all in the
simulation -- so a perfect selection method would keep all 8. Seven are clearly
significant here. bp_treatment is not (p = 0.34 in this draw, its true effect
being the smallest of the eight), and forward stepwise duly drops it.
That is the first lesson, and it is easy to miss because dropping a variable
feels like tidying up: stepwise did not remove a useless predictor, it removed
a real one that this sample could not resolve. The LASSO at its best lambda
keeps all 8. Selection cannot improve on a well-chosen pre-specified list; it
can only lose parts of it.""")
# --- (b) The realistic case: candidates that do not belong ---------------
# Selection only becomes interesting when the candidate list contains variables
# with no relationship to the outcome, which is the usual situation when the
# list is drawn up from everything that happened to be measured.
noise_rng = np.random.default_rng(99)
wide = framingham.copy()
for i in range(1, 13):
wide[f"noise{i}"] = noise_rng.normal(size=n)
CANDIDATES = REAL + [f"noise{i}" for i in range(1, 13)]
print("\n=== Part (b): 8 real predictors + 12 pure-noise candidates ===\n")
step_b = forward_stepwise(wide, CANDIDATES)
lasso_min_b = lasso_selected(wide, CANDIDATES)
lasso_1se_b = lasso_selected(wide, CANDIDATES, one_se=True)
describe("forward stepwise", step_b)
describe("LASSO best lambda", lasso_min_b)
describe("LASSO 1-SE lambda", lasso_1se_b)
print("\n stepwise kept these noise variables:",
", ".join(sorted(v for v in step_b if "noise" in v)) or "none")
# --- (c) Stability across 100 bootstrap samples -------------------------
print("\n=== Part (c): selection frequency across 100 bootstrap samples ===")
B = 100
boot_rng = np.random.default_rng(2025)
counts = {m: dict.fromkeys(CANDIDATES, 0) for m in ("stepwise", "lasso")}
sizes = {m: [] for m in counts}
signatures = {m: [] for m in counts}
for b in range(B):
idx = boot_rng.integers(0, n, n)
boot = wide.iloc[idx].reset_index(drop=True)
picks = {
"stepwise": forward_stepwise(boot, CANDIDATES),
# cv=5 inside the loop purely to keep the runtime reasonable
"lasso": lasso_selected(boot, CANDIDATES, cv=5),
}
for m, chosen in picks.items():
for v in chosen:
counts[m][v] += 1
sizes[m].append(len(chosen))
signatures[m].append("|".join(sorted(chosen)))
freq = pd.DataFrame({
"variable": CANDIDATES,
"truth": ["real" if v in REAL else "noise" for v in CANDIDATES],
"stepwise_pct": [100 * counts["stepwise"][v] / B for v in CANDIDATES],
"lasso_pct": [100 * counts["lasso"][v] / B for v in CANDIDATES],
}).sort_values(["truth", "stepwise_pct"], ascending=[True, False])
print("\nHow often each candidate was selected (%):")
print(freq.to_string(index=False))
print("\nSummary across the 100 bootstrap samples:")
for m in ("stepwise", "lasso"):
real_pct = np.mean([counts[m][v] for v in REAL]) / B * 100
noise_pct = np.mean([counts[m][v] for v in CANDIDATES
if "noise" in v]) / B * 100
print(f" {m:9s} median size {np.median(sizes[m]):.0f} | real predictors "
f"kept {real_pct:.0f}% of the time | noise kept {noise_pct:.0f}% | "
f"{len(set(signatures[m]))} distinct models")
print("""
Conclusions
-----------
(a) Every one of the 8 pre-specified predictors belongs in the model, so the
only thing selection can do is lose one -- which is what happens: stepwise
drops bp_treatment, the weakest real effect. The best-lambda LASSO keeps
all 8. Selection cannot improve on a well-chosen pre-specified list.
(b) Add candidates that do not belong and every method starts letting them in.
Forward stepwise admits noise variables because a variable enters on
whether it improves AIC in this particular sample, and with 12 chances some
noise variable always looks helpful. The best-lambda LASSO is the worst
offender, keeping 9 of the 12: that lambda minimises prediction error, not
the number of wrong variables, so it prefers to retain noise with small
coefficients. The 1-SE lambda is more parsimonious but not clean either --
it still admits noise, and it still misses bp_treatment.
So none of the three recovers the true model, and the two that look tidiest
are tidy for the wrong reason: they dropped a real predictor along with
some of the noise.
(c) The bootstrap frequencies are the real finding. Stepwise produced 99
different models in 100 resamples of the same patients, and the LASSO 70.
Real predictors were kept 84% of the time by stepwise and 97% by the LASSO;
noise variables 40% and 83%.
The instability sits exactly where you would want the method's advice.
age and hdl_chol are chosen every single time, while bp_treatment -- a
genuine predictor -- is chosen by stepwise in under a third of samples. A
reader of any one such model would conclude that blood-pressure treatment
does not predict risk, and the next sample would tell them otherwise.
One subtlety worth naming, because it cuts against the obvious reading.
Two of the noise variables are selected in about 79% of bootstrap samples
by stepwise and 97% by the LASSO, which looks like reliability. It is not:
they happen to correlate with the outcome in this particular cohort, and
every bootstrap sample is drawn from that same cohort, so the fluke is
faithfully reproduced. Bootstrap stability measures robustness to
resampling these patients, not to collecting new ones -- it is a lower
bound on the instability a fresh dataset would reveal.
That is what makes a paper reporting one stepwise-selected model
misleading. The list of retained variables is presented as a finding --
these are the predictors of risk -- when a different sample of the same
patients would have produced a different list, and a genuinely new sample
a different one again. The p-values and confidence intervals are also
wrong, because they take no account of the searching that preceded them.
The defensible options are to pre-specify the predictors on clinical
grounds and keep them all, or to use a penalised model and report the whole
procedure rather than the variables that happened to survive it.
""")17.12 Summary
- A prediction model starts with a clearly defined clinical question, not with data analysis.
- Sample size should be determined prospectively using the Riley et al. criteria, not just the EPV rule of thumb.
- Missing data should be handled with multiple imputation, not deletion.
- Variable selection should be driven by clinical knowledge. Stepwise methods are harmful.
- Overfitting is the central challenge; combat it with adequate sample size and shrinkage.
- Penalised regression (ridge, LASSO, elastic net) embeds shrinkage into the estimation process.
- The final model should be internally validated before any claims about performance.
17.13 References and Further Reading
- For the prediction modelling workflow, see Smits et al. (2026).
- For handling missing data during model development, see Buuren (2018).
- For penalised regression and modelling strategy, see Harrell (2015) and Hastie et al. (2009).
Buuren, Stef van. 2018. Flexible Imputation of Missing Data. 2nd ed. Chapman; Hall/CRC.
Collins, Gary S, Karel G M Moons, Paula Dhiman, et al. 2024. “TRIPOD+AI Statement: Updated Guidance for Reporting Clinical Prediction Models That Use Regression or Machine Learning Methods.” BMJ 385: e078378. https://doi.org/10.1136/bmj-2023-078378.
Harrell, Frank E. 2015. Regression Modeling Strategies: With Applications to Linear Models, Logistic and Ordinal Regression, and Survival Analysis. 2nd ed. Springer.
Hastie, Trevor, Robert Tibshirani, and Jerome Friedman. 2009. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. 2nd ed. Springer.
Riley, Richard D, Joie Ensor, Kym I E Snell, et al. 2020. “Calculating the Sample Size Required for Developing a Clinical Prediction Model.” BMJ 368: m441.
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.
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.
Wijesuriya, Rushani et al. 2025. “Multiple Imputation for Longitudinal Data: A Tutorial.” Statistics in Medicine, ahead of print. https://doi.org/10.1002/sim.10274.