Machine learning

Lesson 11 — Python

Lesson 11 of 20 Advanced ~120 min

Learning objectives

Split data correctly and avoid leakage
Build a scikit-learn pipeline with preprocessing inside it
Fit and compare logistic regression, SVM and tree ensembles
Choose evaluation metrics that suit imbalanced clinical data
Tune hyperparameters with nested cross-validation
Train a neural network, and know when one is not warranted
Interpret a model well enough to defend it

Prediction is not inference

	Inference	Prediction
Question	Does treatment affect outcome?	Who will respond?
Output	Effect size, CI, p-value	A predicted value or probability
Concern	Confounding, assumptions	Generalisation to new data
Validation	Model diagnostics	Held-out performance
Library	`statsmodels`	`scikit-learn`

Lesson 10 covered inference. This lesson is prediction, and the standards are different: nobody cares whether a coefficient is significant if the model does not generalise.

Machine learning is rarely the answer to a regulatory question

A confirmatory efficacy analysis is a pre-specified statistical test, not a model competition. Machine learning belongs in exploratory work — signal detection, site risk-based monitoring, patient stratification hypotheses, imaging endpoints, digital biomarkers.

If someone proposes a random forest for a primary endpoint, the question to ask is what pre-specified hypothesis it tests.

Setting up

uv add scikit-learn

import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split, StratifiedKFold
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.impute import SimpleImputer

RANDOM_STATE = 42

Fix the seed everywhere. An unseeded model is not reproducible, and a performance difference between two runs then tells you nothing.

Splitting data

X = adsl[["AGE", "SEX", "RACE", "BMIBL", "BASELINE_SCORE", "REGION"]]
y = adsl["RESPONDER"]

X_train, X_test, y_train, y_test = train_test_split(
    X, y,
    test_size=0.25,
    random_state=RANDOM_STATE,
    stratify=y,               # preserve the class balance in both halves
)

stratify=y matters whenever classes are imbalanced. Without it a 10% positive class can end up at 6% in the test set purely by chance, and the metrics become incomparable.

Data leakage is the defining failure of applied ML

Leakage means information from the test set influencing training. The result is excellent held-out performance and a model that fails in production.

The three common forms in clinical data:

1. Preprocessing before splitting.

X_scaled = StandardScaler().fit_transform(X)         # WRONG
X_train, X_test = train_test_split(X_scaled, ...)

The scaler saw the test set’s mean and variance. Fit preprocessing inside a pipeline, on training folds only.

2. Subject-level records split at the row level. If one subject contributes several rows, the same subject appearing in both halves lets the model memorise them.

from sklearn.model_selection import GroupShuffleSplit
splitter = GroupShuffleSplit(test_size=0.25, random_state=RANDOM_STATE)
train_idx, test_idx = next(splitter.split(X, y, groups=df["USUBJID"]))

3. Features that encode the outcome. A DISCONTINUATION_REASON column predicts discontinuation perfectly, because it is recorded afterwards. Any feature only knowable after the outcome must be excluded. Ask of every feature: would this be available at the time of prediction?

Pipelines

Preprocessing belongs inside the model, not before it.

numeric_features = ["AGE", "BMIBL", "BASELINE_SCORE"]
categorical_features = ["SEX", "RACE", "REGION"]

numeric_transformer = Pipeline([
    ("impute", SimpleImputer(strategy="median")),
    ("scale",  StandardScaler()),
])

categorical_transformer = Pipeline([
    ("impute", SimpleImputer(strategy="most_frequent")),
    ("encode", OneHotEncoder(handle_unknown="ignore", drop="first")),
])

preprocessor = ColumnTransformer([
    ("num", numeric_transformer, numeric_features),
    ("cat", categorical_transformer, categorical_features),
])

handle_unknown="ignore" is important: a category present in production but absent from training would otherwise raise at predict time.

Logistic regression

from sklearn.linear_model import LogisticRegression

logreg = Pipeline([
    ("prep", preprocessor),
    ("model", LogisticRegression(
        penalty="l2", C=1.0, max_iter=1000,
        class_weight="balanced",          # for imbalanced outcomes
        random_state=RANDOM_STATE,
    )),
])

logreg.fit(X_train, y_train)
y_pred = logreg.predict(X_test)
y_proba = logreg.predict_proba(X_test)[:, 1]

C is the inverse of regularisation strength — smaller C means stronger penalty. This trips people up because it is the opposite of most conventions.

Recovering the coefficients:

feature_names = logreg.named_steps["prep"].get_feature_names_out()
coefs = pd.DataFrame({
    "feature": feature_names,
    "coefficient": logreg.named_steps["model"].coef_[0],
    "odds_ratio": np.exp(logreg.named_steps["model"].coef_[0]),
}).sort_values("coefficient", key=abs, ascending=False)

These odds ratios are regularised and have no confidence intervals. Do not report them as if they came from statsmodels.

Support vector machines

from sklearn.svm import SVC

svm = Pipeline([
    ("prep", preprocessor),
    ("model", SVC(
        kernel="rbf", C=1.0, gamma="scale",
        class_weight="balanced",
        probability=True,           # needed for predict_proba; slows fitting
        random_state=RANDOM_STATE,
    )),
])
svm.fit(X_train, y_train)

Kernel	Use
`linear`	Many features, interpretable, fast
`rbf`	Non-linear boundaries — the usual default
`poly`	Rarely better than `rbf`, harder to tune

SVMs require scaling — without it, features on larger scales dominate the distance calculation. The pipeline handles this. They also scale poorly beyond roughly 50,000 rows; for larger data use LinearSVC or a gradient-boosted tree.

Tree ensembles

In practice these win on tabular clinical data.

from sklearn.ensemble import RandomForestClassifier, HistGradientBoostingClassifier

rf = Pipeline([
    ("prep", preprocessor),
    ("model", RandomForestClassifier(
        n_estimators=500, max_depth=None,
        min_samples_leaf=5, max_features="sqrt",
        class_weight="balanced_subsample",
        random_state=RANDOM_STATE, n_jobs=-1,
    )),
])

# Handles missing values natively — no imputer needed
hgb = HistGradientBoostingClassifier(
    max_iter=500, learning_rate=0.05,
    max_leaf_nodes=31, early_stopping=True,
    validation_fraction=0.1, random_state=RANDOM_STATE,
)

HistGradientBoostingClassifier is the strongest default for tabular data in scikit-learn — comparable to XGBoost or LightGBM, with no extra dependency, and it handles missing values and categorical features natively.

Evaluation

from sklearn.metrics import (
    accuracy_score, precision_score, recall_score, f1_score,
    roc_auc_score, average_precision_score, confusion_matrix,
    classification_report, brier_score_loss,
)

print(classification_report(y_test, y_pred, digits=3))

tn, fp, fn, tp = confusion_matrix(y_test, y_pred).ravel()

Accuracy is the wrong metric for imbalanced data

With 5% responders, a model that predicts “no response” for everyone is 95% accurate and completely useless.

Metric	Answers	Use when
Sensitivity / recall	Of true positives, how many found?	Missing a case is costly
Specificity	Of true negatives, how many correct?	False alarms are costly
PPV / precision	Of predicted positives, how many real?	Acting on positives is expensive
ROC AUC	Ranking quality across thresholds	Roughly balanced classes
PR AUC	Ranking quality for the positive class	Imbalanced classes
Brier score	Calibration of probabilities	Probabilities will be used directly

ROC AUC is optimistic under heavy imbalance because the large true-negative count inflates it. Report average precision (PR AUC) alongside it.

def evaluate(model, X_test, y_test, threshold: float = 0.5) -> dict:
    """Metrics appropriate for a possibly imbalanced clinical outcome."""
    proba = model.predict_proba(X_test)[:, 1]
    pred = (proba >= threshold).astype(int)
    tn, fp, fn, tp = confusion_matrix(y_test, pred).ravel()

    return {
        "n":            len(y_test),
        "prevalence":   y_test.mean(),
        "sensitivity":  tp / (tp + fn) if (tp + fn) else np.nan,
        "specificity":  tn / (tn + fp) if (tn + fp) else np.nan,
        "ppv":          tp / (tp + fp) if (tp + fp) else np.nan,
        "npv":          tn / (tn + fn) if (tn + fn) else np.nan,
        "roc_auc":      roc_auc_score(y_test, proba),
        "pr_auc":       average_precision_score(y_test, proba),
        "brier":        brier_score_loss(y_test, proba),
    }

Cross-validation

A single split gives a noisy estimate.

from sklearn.model_selection import cross_validate

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=RANDOM_STATE)

scores = cross_validate(
    logreg, X, y, cv=cv,
    scoring=["roc_auc", "average_precision", "recall", "precision"],
    return_train_score=True,
)

pd.DataFrame(scores).agg(["mean", "std"]).round(3)

Comparing train and test scores diagnoses the problem:

Train	Test	Diagnosis	Action
High	High	Working	Ship it
High	Low	Overfitting	Regularise, simplify, more data
Low	Low	Underfitting	More capacity, better features

Hyperparameter tuning

from sklearn.model_selection import GridSearchCV, RandomizedSearchCV

param_grid = {
    "model__C": [0.01, 0.1, 1, 10, 100],
    "model__penalty": ["l1", "l2"],
    "model__solver": ["liblinear", "saga"],
}

search = GridSearchCV(
    logreg, param_grid, cv=cv,
    scoring="average_precision", n_jobs=-1, refit=True,
)
search.fit(X_train, y_train)

search.best_params_
search.best_score_

The model__C syntax addresses a parameter of the step named model inside the pipeline.

Tuning on the test set invalidates it

search.fit(X_train, y_train)
search.score(X_test, y_test)          # honest — test used once

Tuning, looking at the test score, adjusting the grid, and re-scoring makes the test set part of training. The reported performance is then optimistic.

For an unbiased estimate when you tune, use nested cross-validation:

from sklearn.model_selection import cross_val_score

inner = StratifiedKFold(5, shuffle=True, random_state=RANDOM_STATE)
outer = StratifiedKFold(5, shuffle=True, random_state=RANDOM_STATE)

nested = cross_val_score(
    GridSearchCV(logreg, param_grid, cv=inner, scoring="average_precision"),
    X, y, cv=outer, scoring="average_precision",
)
print(f"Nested CV: {nested.mean():.3f} (SD {nested.std():.3f})")

The inner loop tunes; the outer loop evaluates. Expensive, and the only honest answer when the model is tuned.

Neural networks

from sklearn.neural_network import MLPClassifier

mlp = Pipeline([
    ("prep", preprocessor),
    ("model", MLPClassifier(
        hidden_layer_sizes=(64, 32),
        activation="relu", alpha=1e-3,
        learning_rate_init=1e-3, max_iter=500,
        early_stopping=True, validation_fraction=0.1,
        random_state=RANDOM_STATE,
    )),
])

For anything beyond this — images, sequences, transfer learning — use PyTorch:

import torch
import torch.nn as nn

class Net(nn.Module):
    def __init__(self, n_features: int):
        super().__init__()
        self.layers = nn.Sequential(
            nn.Linear(n_features, 64), nn.ReLU(), nn.Dropout(0.3),
            nn.Linear(64, 32),         nn.ReLU(), nn.Dropout(0.3),
            nn.Linear(32, 1),
        )

    def forward(self, x):
        return self.layers(x)          # logits; use BCEWithLogitsLoss


model = Net(n_features=X_train.shape[1])
criterion = nn.BCEWithLogitsLoss()
optimiser = torch.optim.AdamW(model.parameters(), lr=1e-3, weight_decay=1e-2)

BCEWithLogitsLoss combines a sigmoid and binary cross-entropy in a numerically stable way. Applying sigmoid yourself and then BCELoss is a common source of nan losses.

Neural networks rarely beat gradient boosting on tabular data

On tabular clinical data of the size most trials produce — hundreds to tens of thousands of rows — a gradient-boosted tree will usually match or beat a neural network, train in seconds rather than hours, need almost no tuning, and be far easier to explain.

Reach for deep learning when the data is genuinely unstructured: images, waveforms, free text, or sequences. For a table of demographics and lab values, start with HistGradientBoostingClassifier and stop there unless it fails.

Convolutional networks, briefly

For imaging endpoints:

import torch.nn as nn

cnn = nn.Sequential(
    nn.Conv2d(1, 32, 3, padding=1), nn.ReLU(), nn.MaxPool2d(2),
    nn.Conv2d(32, 64, 3, padding=1), nn.ReLU(), nn.MaxPool2d(2),
    nn.Flatten(),
    nn.Linear(64 * 7 * 7, 128), nn.ReLU(), nn.Dropout(0.5),
    nn.Linear(128, 10),
)

In practice, fine-tune a pretrained model rather than training from scratch — a medical imaging dataset of a few thousand images is nowhere near enough to train a network from random initialisation.

Interpretation

from sklearn.inspection import permutation_importance

result = permutation_importance(
    rf, X_test, y_test, n_repeats=20,
    random_state=RANDOM_STATE, scoring="average_precision",
)

importance = pd.DataFrame({
    "feature": X_test.columns,
    "importance": result.importances_mean,
    "sd": result.importances_std,
}).sort_values("importance", ascending=False)

Permutation importance is preferable to a tree’s built-in feature_importances_, which is biased towards high-cardinality features.

# SHAP — per-prediction attribution
import shap

explainer = shap.TreeExplainer(rf.named_steps["model"])
shap_values = explainer.shap_values(X_test_transformed)
shap.summary_plot(shap_values, X_test_transformed)

Importance is not causation

A feature can be important because it is a proxy for something else. In one widely cited example, an asthma indicator appeared protective for pneumonia mortality — because asthmatic patients were admitted to intensive care sooner. The model was right about the correlation and dangerously wrong as a basis for triage.

Feature importance describes the model, not the disease.

Calibration

If you will use predicted probabilities as probabilities, they must be calibrated.

from sklearn.calibration import CalibratedClassifierCV, calibration_curve

prob_true, prob_pred = calibration_curve(y_test, y_proba, n_bins=10)

calibrated = CalibratedClassifierCV(svm, method="isotonic", cv=5)
calibrated.fit(X_train, y_train)

SVMs and random forests are typically poorly calibrated out of the box — a random forest’s “0.9” often does not mean 90%. Logistic regression usually is calibrated, being fitted on log-likelihood.

Reproducibility

import joblib, json, sklearn, sys
from datetime import datetime

joblib.dump(model, "models/responder_v1.joblib")

metadata = {
    "model": "responder_v1",
    "created": datetime.now().isoformat(),
    "sklearn_version": sklearn.__version__,
    "python_version": sys.version.split()[0],
    "random_state": RANDOM_STATE,
    "features": list(X.columns),
    "n_train": len(X_train),
    "n_test": len(X_test),
    "prevalence": float(y.mean()),
    "metrics": evaluate(model, X_test, y_test),
    "hyperparameters": model.get_params(),
}
json.dump(metadata, open("models/responder_v1.json", "w"), indent=2, default=str)

A .joblib file without this metadata is unusable in six months — you will not know what it was trained on, which features it expects, or how well it did.

Pickled models are executable

joblib.load and pickle.load deserialise arbitrary Python objects, and loading a model file from an untrusted source can execute code. Treat model artefacts like binaries: known provenance, checksummed, access-controlled.

Common mistakes

Mistake	Consequence	Fix
Preprocessing before splitting	Leakage; optimistic scores	Preprocess inside a `Pipeline`
Row-level split with repeated subjects	Memorisation	`GroupShuffleSplit` on `USUBJID`
Features knowable only after the outcome	Perfect, useless model	Audit every feature
Accuracy on imbalanced data	Meaningless	PR AUC, sensitivity, specificity
Tuning against the test set	Optimistic estimate	Nested CV
Unseeded randomness	Irreproducible	`random_state` everywhere
Reporting `sklearn` coefficients as odds ratios	Regularised, no CIs	`statsmodels` for inference
Neural network on small tabular data	Slower, worse, unexplainable	Gradient boosting
Uncalibrated probabilities used as risks	Miscalibrated decisions	`CalibratedClassifierCV`
Model saved without metadata	Unusable later	Save a JSON sidecar

Exercise 11.1 — A leakage-free comparison

Build and compare logistic regression, random forest and gradient boosting for a binary clinical outcome. Use cross-validation, report metrics suitable for imbalance, and guard against every leakage form described above.

Show solution

import numpy as np
import pandas as pd
from sklearn.compose import ColumnTransformer
from sklearn.ensemble import RandomForestClassifier, HistGradientBoostingClassifier
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GroupShuffleSplit, StratifiedGroupKFold, cross_validate
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler

RANDOM_STATE = 42

# --- Feature audit: exclude anything known only after the outcome ----------
LEAKY = {
    "DISCONTINUATION_REASON",   # recorded at/after the event
    "TOTAL_DOSE_RECEIVED",      # accumulates over follow-up
    "LAST_VISIT_DATE",          # later for subjects who stayed longer
    "EOSSTT",                   # end-of-study status IS the outcome
}

NUMERIC = ["AGE", "BMIBL", "BASELINE_SCORE", "BASELINE_ALT", "DIAGNOSIS_MONTHS"]
CATEGORICAL = ["SEX", "RACE", "REGION", "PRIOR_THERAPY"]

def audit_features(df: pd.DataFrame, features: list[str]) -> list[str]:
    """Fail loudly if a known-leaky feature was included."""
    leaks = set(features) & LEAKY
    if leaks:
        raise ValueError(
            f"Post-outcome feature(s) in the model: {sorted(leaks)}. "
            "These are not knowable at prediction time."
        )
    missing = set(features) - set(df.columns)
    if missing:
        raise ValueError(f"Features absent from the data: {sorted(missing)}")
    return features


features = audit_features(adsl, NUMERIC + CATEGORICAL)
X = adsl[features]
y = adsl["RESPONDER"].astype(int)
groups = adsl["USUBJID"]              # one row per subject here, but be explicit

# --- Split by SUBJECT, not by row ------------------------------------------
splitter = GroupShuffleSplit(n_splits=1, test_size=0.25, random_state=RANDOM_STATE)
train_idx, test_idx = next(splitter.split(X, y, groups=groups))

X_train, X_test = X.iloc[train_idx], X.iloc[test_idx]
y_train, y_test = y.iloc[train_idx], y.iloc[test_idx]
g_train = groups.iloc[train_idx]

assert not (set(groups.iloc[train_idx]) & set(groups.iloc[test_idx])), \
    "Subject appears in both train and test"

print(f"Train {len(X_train)} (prevalence {y_train.mean():.1%}), "
      f"Test {len(X_test)} (prevalence {y_test.mean():.1%})")

# --- Preprocessing INSIDE the pipeline -------------------------------------
preprocessor = ColumnTransformer([
    ("num", Pipeline([("impute", SimpleImputer(strategy="median")),
                      ("scale", StandardScaler())]), NUMERIC),
    ("cat", Pipeline([("impute", SimpleImputer(strategy="most_frequent")),
                      ("encode", OneHotEncoder(handle_unknown="ignore",
                                               drop="first"))]), CATEGORICAL),
])

models = {
    "Logistic regression": Pipeline([
        ("prep", preprocessor),
        ("model", LogisticRegression(max_iter=1000, class_weight="balanced",
                                     random_state=RANDOM_STATE)),
    ]),
    "Random forest": Pipeline([
        ("prep", preprocessor),
        ("model", RandomForestClassifier(n_estimators=500, min_samples_leaf=5,
                                         class_weight="balanced_subsample",
                                         random_state=RANDOM_STATE, n_jobs=-1)),
    ]),
    "Gradient boosting": Pipeline([
        ("prep", preprocessor),
        ("model", HistGradientBoostingClassifier(max_iter=500, learning_rate=0.05,
                                                 early_stopping=True,
                                                 random_state=RANDOM_STATE)),
    ]),
}

# --- Cross-validate, grouped and stratified --------------------------------
cv = StratifiedGroupKFold(n_splits=5, shuffle=True, random_state=RANDOM_STATE)
SCORING = ["roc_auc", "average_precision", "recall", "precision"]

rows = []
for name, model in models.items():
    cvres = cross_validate(model, X_train, y_train, groups=g_train, cv=cv,
                           scoring=SCORING, return_train_score=True, n_jobs=-1)
    rows.append({
        "model": name,
        "cv_pr_auc":  cvres["test_average_precision"].mean(),
        "cv_pr_sd":   cvres["test_average_precision"].std(),
        "cv_roc_auc": cvres["test_roc_auc"].mean(),
        "cv_recall":  cvres["test_recall"].mean(),
        "train_pr_auc": cvres["train_average_precision"].mean(),
    })

cv_summary = pd.DataFrame(rows).sort_values("cv_pr_auc", ascending=False)
cv_summary["overfit_gap"] = cv_summary["train_pr_auc"] - cv_summary["cv_pr_auc"]
print(cv_summary.round(3).to_string(index=False))

              model  cv_pr_auc  cv_pr_sd  cv_roc_auc  cv_recall  train_pr_auc  overfit_gap
  Gradient boosting      0.612     0.048       0.784      0.588         0.741        0.129
      Random forest      0.594     0.052       0.771      0.612         0.998        0.404
Logistic regression      0.571     0.041       0.762      0.647         0.601        0.030

Final evaluation, on the test set used once:

best_name = cv_summary.iloc[0]["model"]
best = models[best_name].fit(X_train, y_train)

final = evaluate(best, X_test, y_test)
print(f"\n{best_name} — held-out performance")
for k, v in final.items():
    print(f"  {k:<12} {v:.3f}")

What the numbers say

The overfit_gap column is the most informative. The random forest scores 0.998 on training data and 0.594 in cross-validation — it has memorised the training set almost perfectly. It still generalises acceptably, but that gap warns that it is sensitive to the training sample and would benefit from stronger constraints (max_depth, larger min_samples_leaf).

Logistic regression has almost no gap (0.030). It is not memorising, and its cross-validated performance is only 0.04 below the best model — for a 4% difference in PR AUC, a model with interpretable coefficients and stable behaviour is often the better choice in a clinical setting.

Note also that PR AUC (0.61) is far below ROC AUC (0.78) — the signature of an imbalanced outcome. Reporting only ROC AUC would substantially overstate how useful this model is.

The four leakage guards, explicitly

audit_features() raises on any post-outcome feature. The list is maintained by hand and should be reviewed with a clinician — this is a domain question, not a technical one.
GroupShuffleSplit on USUBJID plus an assertion that no subject spans the split.
StratifiedGroupKFold so the cross-validation folds also respect subject boundaries and preserve class balance.
All imputation and scaling inside the Pipeline, so each fold fits its own preprocessing.

One thing this still does not do. The three models were compared on cross-validated scores and then the winner was evaluated on the test set. That selection is itself a form of tuning, so the final test score is slightly optimistic. For a number you would publish, wrap the whole selection procedure in an outer cross-validation loop.

Recap

Prediction and inference are different jobs with different libraries
Leakage is the defining failure: preprocess inside a Pipeline, split by subject, audit features
stratify= and class_weight="balanced" for imbalanced outcomes
Report PR AUC, sensitivity and specificity — not accuracy
Compare train and test scores to diagnose over- and underfitting
Nested cross-validation when the model is tuned
Gradient boosting beats neural networks on tabular clinical data
Feature importance describes the model, not the disease
Save a metadata sidecar with every model

Next: Clinical tables and TLF generation.