How to Build a Medical Prediction Model: LASSO, Random Forest and XGBoost, Nomograms and External Validation

Workflow

Workflow: what each step does and which steps must sit inside the resampling loop

Every step that uses the outcome or the data distribution is part of model development and must be redone within each training fold or bootstrap sample during internal validation. Missing one step makes the validation optimistic.

A diagnostic model predicts whether a condition is present now; a prognostic model predicts whether an outcome occurs within a future horizon. They differ in prediction time point and in the variables available. The PROBAST predictor domain asks whether all predictors are available at the time the model is intended to be used, and rates the risk of bias as high when they are not.

Sample size

Sample size: the limits of the EPV rule and the Riley method (pmsampsize)

EPV ≥ 10 (at least 10 events per variable) comes from Peduzzi et al.'s 1996 simulation of regression coefficient bias. It ignores outcome prevalence and expected model performance. van Smeden et al. (2019) showed that EPV is only weakly related to prediction error in new data. Riley et al. (BMJ 2020) proposed a calculation based on target precision; TRIPOD+AI item 10 asks how the study size was arrived at and why it is sufficient, and this method is the calculation currently recommended in the methodological literature.

library(pmsampsize)   # 1.1.3
# Binary outcome: prevalence 0.253, C statistic 0.80 from a published model (prefer an optimism-adjusted value), 30 candidate parameters
pmsampsize(type = "b", cstatistic = 0.80, parameters = 30, prevalence = 0.253)
#> Criteria 1 (shrinkage >= 0.9)          1136
#> Criteria 2 (R2 optimism <= 0.05)        774
#> Criteria 3 (intercept +/- 0.05)         291
#> Final 1136, 288 events, EPP = 9.58

pmsampsize(type = "b", cstatistic = 0.70, parameters = 30, prevalence = 0.253)
#> Final 2732, 692 events, EPP = 23.04

# Time-to-event: Cox-Snell R2 0.10, 15 parameters, event rate 0.08 per person-year, 3-year horizon, mean follow-up 4 years
pmsampsize(type = "s", csrsquared = 0.10, parameters = 15, rate = 0.08,
           timepoint = 3, meanfup = 4)
#> Final 1274, 408 events, EPP = 27.18

Variable selection

LASSO: scaling, choice of lambda and unstable variable sets

cv.glmnet assigns folds at random. Its help page states that "the results of cv.glmnet are random, since the folds are selected at random" and that users can reduce this randomness by running cv.glmnet many times and averaging the error curves. The table below, from one dataset, shows how large the seed effect is.

Instability is worse in high-dimensional expression data (thousands of genes, tens to hundreds of samples) than in the table above: variables far outnumber samples, highly correlated genes can substitute for each other, and LASSO often keeps one gene from a co-expressed group more or less at random. A single "feature gene list" then means little; report bootstrap selection frequencies, or use the elastic net (alpha between 0 and 1) so that correlated genes enter together.

library(glmnet)                                     # version 5.0 in our run
X <- model.matrix(y ~ ., data = dat)[, -1]          # factors become dummy columns; X must be a numeric matrix
y <- dat$y                                          # 0/1

set.seed(2026)
foldid <- sample(rep(1:10, length.out = nrow(X)))   # fix and save the folds; reuse them for reproduction and the supplement
cv <- cv.glmnet(X, y, family = "binomial", foldid = foldid)
# defaults: alpha = 1 (LASSO), standardize = TRUE (internal scaling, coefficients returned on the original scale), type.measure = "deviance"
plot(cv)                                            # left line lambda.min, right line lambda.1se
coef(cv, s = "lambda.min")                          # always pass s; coef(cv) without s returns lambda.1se
coef(cv, s = "lambda.1se")

# Selection stability: one cross-validation per seed for 100 seeds, count how often each variable is selected
sel <- sapply(1:100, function(s) {
  set.seed(s)
  b <- coef(cv.glmnet(X, y, family = "binomial", nfolds = 10), s = "lambda.1se")[-1, 1]
  b != 0
})
sort(rowSums(sel), decreasing = TRUE)               # selection frequency
length(unique(apply(sel, 2, paste, collapse = ""))) # number of distinct variable sets

# Reduce fold randomness: average the CV error curves of 100 fold assignments on one lambda sequence
lam <- glmnet(X, y, family = "binomial")$lambda
cvm <- sapply(1:100, function(s) {
  set.seed(s)
  cv.glmnet(X, y, family = "binomial", lambda = lam,
            foldid = sample(rep(1:10, length.out = nrow(X))))$cvm
})
lam[which.min(rowMeans(cvm))]                       # lambda.min on the averaged curve

Algorithm choice

Choosing between logistic regression, LASSO, random forest, XGBoost and SVM

Christodoulou et al. (2019) reviewed 282 comparisons of logistic regression and machine learning in 71 studies. In the 145 comparisons at low risk of bias, the difference in logit(AUC) was 0.00 (95% CI −0.18 to 0.18); in the 137 at high risk of bias, machine learning was 0.34 higher, which the authors attributed to flawed validation. 79% of studies did not assess calibration.

Order of choice: first fit a logistic or Cox model with splines as the baseline. If the machine learning model does not clearly beat the baseline on both discrimination and calibration in nested cross-validation or external validation, report the baseline. This comparison is the answer when reviewers ask why machine learning was used.

scikit-learn's LogisticRegression applies an L2 penalty with C = 1 by default. When a paper's "logistic regression" baseline is fitted with scikit-learn 1.8 or later, set C = np.inf for an unpenalised model. From 1.8 the penalty argument is deprecated and scheduled for removal in 1.10; L1 is written l1_ratio = 1, and the penalty = "l1" in older tutorials raises a deprecation warning.

In R, predict(rf) without newdata returns out-of-bag (OOB) predictions, which are close to unbiased; predict(rf, newdata = training data) returns predictions from all trees, and the training AUC approaches 1. When reporting "training performance" for a random forest, state which one was used.

Tuning and internal validation

Tuning, nested cross-validation and bootstrap optimism correction

The cross-validation error used to choose hyperparameters is itself optimistic. On simulated data with no real signal, Varma and Simon (2006) found that 18.5% of training sets gave a tuned cross-validation error below 30%; nested cross-validation gave an almost unbiased error.

After nested cross-validation gives the performance estimate, refit the final model once on all data with the same inner tuning procedure. Different outer folds choosing different hyperparameters is expected; it reflects the uncertainty of the pipeline.

In scikit-learn, passing an integer random_state to a cross-validation splitter gives every model the same folds, so folds can be compared one by one; passing an integer to an estimator (such as a random forest) gives every fold the same random state. The documentation recommends integers for splitters and, when comparing algorithms, a RandomState instance or no value for estimators.

rms validate and calibrate default to B = 40, which varies noticeably between seeds; B = 200 or more is common in papers. For stepwise selection, validate(fit, B = 200, bw = TRUE) repeats backward elimination inside each bootstrap sample; LASSO has no such argument, so write the loop yourself as in the code above.

import numpy as np
from sklearn.model_selection import StratifiedKFold, GridSearchCV, cross_validate
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from xgboost import XGBClassifier

inner = StratifiedKFold(n_splits=5, shuffle=True, random_state=1)   # inner loop: tuning
outer = StratifiedKFold(n_splits=5, shuffle=True, random_state=2)   # outer loop: performance estimate
models = {
    # penalty is deprecated from scikit-learn 1.8; L1 is l1_ratio=1. Before 1.8 use penalty="l1"
    "lasso_lr": (Pipeline([("imp", SimpleImputer(strategy="median")),
                           ("sc", StandardScaler()),
                           ("clf", LogisticRegression(l1_ratio=1, solver="saga", max_iter=5000))]),
                 {"clf__C": np.logspace(-3, 1, 20)}),
    "rf": (Pipeline([("imp", SimpleImputer(strategy="median")),
                     ("clf", RandomForestClassifier(n_estimators=500, random_state=0))]),
           {"clf__max_features": ["sqrt", 0.33, 0.5], "clf__min_samples_leaf": [1, 5, 10, 20]}),
    "xgb": (Pipeline([("clf", XGBClassifier(n_estimators=300, eval_metric="logloss"))]),
            {"clf__learning_rate": [0.03, 0.1], "clf__max_depth": [2, 3, 4],
             "clf__min_child_weight": [1, 5]}),
}
for name, (pipe, grid) in models.items():
    search = GridSearchCV(pipe, grid, cv=inner, scoring="neg_log_loss")  # tune on log loss, which also rewards calibration
    res = cross_validate(search, X, y, cv=outer, scoring=["roc_auc", "neg_brier_score"])
    print(name, res["test_roc_auc"].mean().round(3), res["test_roc_auc"].std().round(3),
          (-res["test_neg_brier_score"]).mean().round(4))

library(glmnet); library(pROC)
auc_of <- function(p, y) as.numeric(auc(roc(y, p, direction = "<", quiet = TRUE)))  # explicit direction
fit_lasso <- function(X, y) cv.glmnet(X, y, family = "binomial", nfolds = 10)
p_of <- function(cv, X) predict(cv, X, s = "lambda.1se", type = "response")[, 1]

cv0 <- fit_lasso(X, y)
apparent <- auc_of(p_of(cv0, X), y)                 # apparent AUC: same data for fitting and evaluation
set.seed(2026)
opt <- replicate(200, {
  i <- sample(nrow(X), replace = TRUE)
  cvb <- fit_lasso(X[i, ], y[i])                    # selection and lambda choice are redone in every bootstrap sample
  auc_of(p_of(cvb, X[i, ]), y[i]) - auc_of(p_of(cvb, X), y)
})
c(apparent = apparent, optimism = mean(opt), corrected = apparent - mean(opt))

Class imbalance

Class imbalance: SMOTE is usually unnecessary

van den Goorbergh et al. (JAMIA 2022) compared random undersampling, random oversampling and SMOTE: none improved AUC, and all caused strong overestimation of minority-class probabilities, i.e. worse calibration. A clinical prediction model outputs risks, so poor calibration misleads decisions directly.

from imblearn.pipeline import Pipeline          # the imblearn Pipeline; the sklearn Pipeline does not accept samplers
from imblearn.over_sampling import SMOTE
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold, cross_val_score

pipe = Pipeline([("sc", StandardScaler()),
                 ("smote", SMOTE(random_state=0)),   # applied to each training fold only; test folds keep the original distribution
                 ("clf", LogisticRegression(max_iter=1000))])
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)
print(cross_val_score(pipe, X, y, cv=cv, scoring="roc_auc").mean())
# Predicted probabilities are inflated after resampling; recalibrate on data with the original distribution before reporting risks (TRIPOD+AI item 13)

Fatal errors

Data leakage and validation errors: what goes wrong, consequences and fixes

The figures below come from published methodological studies and official documentation examples and show how much these errors can inflate AUC or accuracy.

A common "intersect several machine learning algorithms to find hub genes" pipeline runs LASSO, SVM-RFE and random forest on one GEO dataset, takes the 3 to 5 genes in the intersection, and draws an ROC curve on the same dataset, often with AUC above 0.9. The intersection is not validation: all three algorithms saw every sample's outcome, so the intersection only reduces the number of genes and the selection bias remains. That ROC curve is the first error in the table. Ambroise and McLachlan's permuted-label experiment shows that under such a pipeline pure noise can reach near-zero error.

There are two acceptable approaches. First, put the whole pipeline (differential expression, selection by the three algorithms, intersection, model fitting) inside every iteration of an outer cross-validation or bootstrap, and report outer performance and each gene's selection frequency. Second, fix genes, coefficients and cut-offs in the development data and only predict and evaluate in another independent cohort (a different centre or period, preferably a different platform). Another series from the same platform and lab in GEO can supplement this, but it resembles the development data and says little about transportability.

In cross-platform validation (for example, development on TCGA RNA-seq and validation on GEO microarrays), the same coefficients are applied to expression values on different scales, so absolute risk scores are not comparable. A common workaround is z-scoring each gene within each cohort; discrimination (C-index, time-dependent AUC) can then be evaluated, while calibration requires re-estimating baseline risk in the new cohort, which should be stated in the paper.

Evaluation

Discrimination, calibration and DCA: what to report and how to read them

TRIPOD+AI item 12e asks for the measures and plots used to evaluate discrimination, calibration and clinical utility, and item 23a asks for performance estimates with confidence intervals. Reporting AUC alone is among the problems reviewers point out most often.

library(pROC); library(rms); library(dcurves)
# p: predicted probabilities from external validation (or outer CV / bootstrap-corrected); y: 0/1 outcome
r <- roc(y, p, direction = "<", quiet = TRUE)    # set direction explicitly on validation data
ci.auc(r)                                        # DeLong 95% CI by default
# roc.test(r1, r2)                               # paired DeLong test for two non-nested models in the same patients

val.prob(p, y)                                   # C, Brier, calibration intercept, slope, Emax, and a smoothed calibration curve

df <- data.frame(y = y, model = p)
dca(y ~ model, data = df, thresholds = seq(0, 0.5, by = 0.01)) |>
  plot(smooth = TRUE)                            # Treat All and Treat None are drawn as references

Nomograms and survival models

rms nomograms, LASSO-Cox, time-dependent ROC and forest plots

A nomogram draws a fitted model as a chart that can be used by hand; it adds nothing to the model's validity. A nomogram needs the same evidence as the model: internal validation, calibration and external validation.

library(rms)                                     # checked against 8.1.1
dd <- datadist(dat); options(datadist = "dd")    # set before fitting, using the same data frame
fit <- lrm(y ~ rcs(age, 4) + sex + scoma + meanbp + pafi + alb + bili + crea,
           data = dat, x = TRUE, y = TRUE)       # x, y = TRUE are needed by validate and calibrate
nom <- nomogram(fit, fun = plogis, funlabel = "Risk of in-hospital death",
                fun.at = c(0.05, 0.1, 0.2, 0.4, 0.6, 0.8))
plot(nom)

set.seed(2026)
validate(fit, B = 200)       # optimism-corrected Dxy, slope, etc.; C = Dxy / 2 + 0.5; B defaults to only 40
cal <- calibrate(fit, B = 200)
plot(cal)                    # Apparent and Bias-corrected lines; report Bias-corrected

library(survival); library(glmnet); library(rms); library(timeROC)
d <- na.omit(subset(pbc, !is.na(trt),
     select = c(time, status, age, albumin, bili, copper, protime, ast, edema, stage)))
d$death <- as.integer(d$status == 2)             # 2 = death; 1 = transplant, treated as censored
d$status <- NULL

# LASSO-Cox: y can be a Surv object (supported since glmnet 4.1)
X <- model.matrix(~ . - time - death, d)[, -1]
set.seed(2026)
cvc <- cv.glmnet(X, Surv(d$time, d$death), family = "cox", type.measure = "C")
coef(cvc, s = "lambda.1se")

# Cox nomogram: surv = TRUE and time.inc equal to the u used later in calibrate
dd <- datadist(d); options(datadist = "dd")
fit <- cph(Surv(time, death) ~ log(bili) + albumin + age + protime + edema,
           data = d, x = TRUE, y = TRUE, surv = TRUE, time.inc = 1826)
cox.zph(fit)                                     # proportional hazards assumption
sv <- Survival(fit)
plot(nomogram(fit, fun = list(function(lp) sv(1826, lp)), funlabel = "5-year survival"))
set.seed(2026)
cal <- calibrate(fit, u = 1826, cmethod = "KM", m = 50, B = 200)   # the default cmethod = "hare" needs the polspline package
plot(cal)

# Time-dependent ROC and C-index
lp <- predict(fit)                               # linear predictor; higher means higher risk
tr <- timeROC(T = d$time, delta = d$death, marker = lp, cause = 1,
              times = c(730, 1826), iid = TRUE)
tr$AUC; confint(tr)
concordance(Surv(time, death) ~ lp, data = d, reverse = TRUE)                   # Harrell's C
concordance(Surv(time, death) ~ lp, data = d, reverse = TRUE, timewt = "n/G2")  # Uno's C

Interpretability

SHAP: how to use it and which conclusions it cannot support

SHAP decomposes a single prediction into feature contributions and is useful for showing what a tree model has learned. The SHAP documentation has an article titled "Be careful when interpreting predictive models in search of causal insights", explaining that SHAP shows the correlations a model exploits.

import shap                                      # 0.52
explainer = shap.TreeExplainer(xgb_model)        # xgb_model: a fitted XGBClassifier
sv = explainer(X_test)                           # compute on the test set or outer folds, not the training set
shap.plots.beeswarm(sv)                          # for binary outcomes the x-axis is in log-odds, not probability
shap.plots.scatter(sv[:, "age"], color=sv)       # dependence plot: non-linearity and interactions

Reporting and peer review

TRIPOD+AI essentials and common reviewer objections

TRIPOD+AI (Collins et al., BMJ 2024;385:e078378) has 27 items, covers regression and machine learning models, and supersedes TRIPOD 2015. PROBAST+AI (Moons et al., BMJ 2025;388:e082505) assesses risk of bias, with 16 signalling questions for model development and 18 for model evaluation, and supersedes PROBAST 2019. Check your manuscript against both before submission.

Chinese-language practice

Prediction model advice from Chinese-language sources: what to follow and what to correct

The points below come from Zhihu columns, "Xiaobai Xue Tongji" articles republished by Mediecogroup, and a bioinformatics tutorial on cnblogs, checked against English methodological literature and software documentation. Templated CSDN tutorials and an article with fabricated reviewer comments were excluded.

Hand it to an agent

Hand this workflow to an agent

Sample size calculation, repeated cross-validation and bootstrap, multi-algorithm comparison and plotting are time-consuming but rule-based steps that a scientific agent can run; the research question, prediction time point and interpretation remain the researcher's decisions.

Example instruction: "Here is my cohort data cohort.csv. The outcome is 30-day death and the prediction time point is 24 hours after admission. First calculate the sample size with pmsampsize, then use spline logistic regression as the baseline and compare LASSO, random forest and XGBoost. Report AUC, calibration intercept and slope, and DCA with 5×5 nested cross-validation, draw a nomogram, and list what the methods section needs under TRIPOD+AI."

The agent works in an isolated cloud computer: it checks variable types, missingness and any variables only available after the prediction time point; runs pmsampsize; places imputation, scaling, selection and tuning in one Pipeline or bootstrap loop; checks LASSO stability across 100 seeds; produces nested cross-validation results, calibration curves, decision curves, a nomogram and SHAP plots for each model; and generates a TRIPOD+AI checklist. scikit-learn, statsmodels, pandas and others are preinstalled, and R packages are installed by the agent in the workspace. The agent reviews results adversarially and checks for leaking steps; the workspace keeps code, random seeds, software versions, logs and results for reproduction.

You still need to check: whether the outcome definition and prediction time point fit the clinical setting; whether every candidate variable is available at the prediction time point; whether the DCA threshold range is reasonable; whether the external validation data are truly independent; and the clinical interpretation of the results.