Variable selection for prediction models

Change the question and the whole rulebook changes

Every regression in the B2 and B6 families is answering an aetiological question: does this exposure have an effect on this outcome. In that kind of model, “should this variable go in” is a question about causal structure — a confounder must go in, a mediator must stay out, and a collider creates bias if you include it. The criteria live in the cohort chapter and on the DAG page.

A prediction model asks something else entirely: what is this person’s risk? The same decision now runs on completely different criteria:

Three strategies, one dataset

survival::rotterdam is the development cohort (2982 patients, 1713 events) and survival::gbsg is the external validation set (686 patients). There are 8 candidate variables: age, menopausal status, tumour size, grade, node count, PgR, ER and hormonal treatment.

The three approaches you will meet most often in published papers:

1. Full model approach — specify all 8 in advance and keep every one, whatever its p-value.
2. Stepwise regression — start from the full model and delete backwards using AIC.
3. Univariable screening — fit each variable on its own first, and only let the ones with p < 0.05 into the multivariable model.

Stepwise regression’s problem is not that it selects badly, it is that it selects unstably

The model diagnostics page already demonstrated one half of this with 300 simulations: throw 10 pure noise variables at stepAIC() and it keeps 1.81 of them on average, with at least one “significant” random variable surviving in 48.0% of runs. That demonstration is not repeated here.

This page asks a different question: same algorithm, same patients, one new resample — will it give you the same model?

Resample the development cohort 200 times, run the identical backward AIC on each resample, and record how often each variable survives:

Across those 200 resamples, backward AIC produced 25 different models, and the most frequent one accounted for only 22% of the runs. Cut each resample to two hundred patients — much closer to the size of most single-centre studies — and it becomes 69 different models, with the most common one accounting for 12%.

Univariable screening has a fault stepwise regression does not

Univariable screening — fit each variable alone, admit only those with p < 0.05 — is even more common than stepwise regression, because it looks more conservative. It is not.

Look at menopausal status in the table above: its univariable p-value is far below 0.001, so it passes the screen and enters the multivariable model — but backward AIC deletes it. The reason is that it almost duplicates age: once you adjust for age, there is nothing left of it.

The reverse case is the dangerous one: a variable can carry no univariable signal at all and still matter after adjustment. Univariable screening throws it away before it ever gets adjusted, and no number in the final model records that this happened.

So how should you choose?

The current advice (TRIPOD+AI, the methodological literature from Riley and colleagues) comes down to one sentence: fix the candidate variables and the model form before you see any results, then put all of them in.

Run it yourself

library(survival); library(MASS)
data(cancer, package = "survival")

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

cand <- c("age", "meno", "size_mm", "grade3", "nodes", "log_pgr", "log_er", "hormon")
f <- as.formula(paste("Surv(rfs_time, rfs_event) ~", paste(cand, collapse = " + ")))

full <- coxph(f, data = rot)          # full model approach: this one line
step <- stepAIC(full, direction = "backward", trace = 0)

# Univariable screening: note what it is asking -- association without adjustment
sapply(cand, function(v)
  summary(coxph(as.formula(paste("Surv(rfs_time, rfs_event) ~", v)),
                data = rot))$coefficients[1, "Pr(>|z|)"])

# Stability: resample the same data and watch what it picks each time
tally <- setNames(numeric(length(cand)), cand)
for (b in 1:200) {
  d  <- rot[sample(nrow(rot), replace = TRUE), ]
  ff <- f; environment(ff) <- environment()   # without this line the whole loop fails silently
  st <- stepAIC(coxph(ff, data = d), direction = "backward", trace = 0)
  tally[names(coef(st))] <- tally[names(coef(st))] + 1
}
round(100 * tally / 200)

import numpy as np
import pandas as pd
from lifelines import CoxPHFitter

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

# The same derived columns as the R code above. rotterdam ships none of
# them, and without this every line below is a KeyError.
d["rfs_time"]  = d[["rtime", "dtime"]].min(axis=1)
d["rfs_event"] = ((d["recur"] == 1) | (d["death"] == 1)).astype(int)
d["size_mm"]   = d["size"].map({"<=20": 15, "20-50": 35, ">50": 60})
d["grade3"]    = (d["grade"] >= 3).astype(int)
d["log_pgr"]   = np.log1p(d["pgr"]); d["log_er"] = np.log1p(d["er"])

cand = ["age", "meno", "size_mm", "grade3", "nodes", "log_pgr", "log_er", "hormon"]
cols = cand + ["rfs_time", "rfs_event"]

# Full model approach
full = CoxPHFitter().fit(d[cols], "rfs_time", "rfs_event")
print(full.concordance_index_)

# The penalised version (the subject of B5-02) -- not a substitute for variable
# selection but an alternative to it
pen = CoxPHFitter(penalizer=0.1, l1_ratio=1.0).fit(d[cols], "rfs_time", "rfs_event")
print(pen.params_[pen.params_ != 0].index.tolist())

When the paper uses a machine-learning model

“XGBoost predicts AKI”, “random forest predicts sepsis” — papers like these are now common in clinical informatics and critical care, and a medical student usually meets one for the first time at a journal club. This section teaches how to read them, not how to fit them. Swapping the model for a forest or a boosted ensemble is a separate skill tree, but not one item on the acceptance list changes.

In order: variable importance, impurity bias, SHAP, and the one thing that has not changed.

One, variable importance is not a regression coefficient. It has no direction, no confidence interval, and it is not an effect size. A variable ranked first only means it came first under whatever importance definition the package used — permutation importance means that shuffling it degrades prediction the most, while the default in tree packages is often the impurity-based version of point two, and the two can rank the same variables differently. Either way it does not tell you whether risk goes up or down, or by how much, and it cannot be cited as “this factor raises risk by so much”. When a Discussion turns an importance ranking into a list of risk factors, it is treating two entirely different quantities as one.

Two, impurity-based importance is systematically biased towards high-cardinality variables. A tree picks a split at every node, and a variable with more distinct values offers more places to split, so it more easily finds a cut that looks useful on noise alone. Continuous variables and many-levelled categorical ones (admitting specialty, drug codes) therefore rise systematically, and binary variables sink systematically — neither of which has anything to do with whether they are actually useful. Permutation importance, or recomputing on held-out data, reduces this but does not remove it.

Three, SHAP is local attribution, not a causal effect. It answers “for this patient, how far did this variable push the model’s prediction away from the baseline” — a decomposition that takes the fitted model as given. If the association the model learned came from confounding, selection or reverse causation, SHAP will faithfully attribute that association, and will draw it very persuasively. A beautiful SHAP plot is not causal evidence; it is the model describing itself.

Four, the acceptance criteria are unchanged. Discrimination, calibration, the decision curve, and external validation — every one of them, none optional. And the one machine-learning papers most often skip is calibration: tree ensembles shrink their predictions towards the middle by construction, a calibration slope above 1 is a common finding, and a paper that reports only the AUC will never show you that.

Four questions to ask when reading a paper

1. Were the candidate variables prespecified, or chosen from the data? If the Methods do not say “prespecified”, assume they came from the data.
2. Was anything revalidated after selection? Selecting variables spends degrees of freedom; the p-values, confidence intervals and C-index computed afterwards are all optimistic.
3. Will every variable in the model actually be available at deployment? This is the one place a reader can catch data leakage.
4. Are the discarded variables reported? A paper that lists only the final model leaves the reader no way to judge how marginal the deletions were.

Common misuses

Reproducing every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B5-01-variable-selection.R

Watch next

▶What Are Clinical Prediction Rules?

▶AI 臨床研究實戰 EP7｜預測 vs 分類、Data Leakage、Propensity Score

▶Key Steps and Common Pitfalls in Clinical Prediction Model Research

Report a problem on this page