Comparing more than two models with stambo
V0.1.6: © Aleksei Tiulpin, PhD, 2025
stambo.compare_models (and stambo.two_sample_test) compare exactly two models/samples. Often you have more than two candidates – several classifiers, several ablations of a model – and want to know which pairs actually differ. This notebook shows stambo.compare_models_pairwise (and its lower-level building block, stambo.pairwise_bootstrap_test), which runs the bootstrap test on every pair among \(N\) models.
Why this needs a correction. If you compare \(N\) models, there are \(N(N-1)/2\) pairs. Even if none of the models are actually different, testing many pairs at \(\alpha=0.05\) means each individual test still has a 5% chance of a false positive – and with enough pairs, the chance that at least one of them looks “significant” purely by chance grows quickly. compare_models_pairwise addresses this with a Holm-Bonferroni correction, applied by default.
Imports
[1]:
import numpy as np
import stambo
from sklearn.neighbors import KNeighborsClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import roc_auc_score
SEED = 2025
stambo.__version__
[1]:
'0.1.6'
Loading the UCI breast cancer dataset and training four models
[2]:
X, y = load_breast_cancer(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.5, random_state=SEED, stratify=y)
scaler = StandardScaler()
scaler.fit(Xtr)
Xtr = scaler.transform(Xtr)
Xte = scaler.transform(Xte)
[3]:
# On purpose, we include one clearly weak model (a decision stump) so that we have
# both large, easy-to-detect differences and small, borderline ones in the same comparison.
models = {
"kNN": KNeighborsClassifier(n_neighbors=3),
"LogReg": LogisticRegression(C=1e-2, random_state=42),
"Stump": DecisionTreeClassifier(max_depth=1, random_state=42),
"RandomForest": RandomForestClassifier(n_estimators=50, max_depth=3, random_state=42),
}
labels = list(models.keys())
preds = []
for name in labels:
model = models[name]
model.fit(Xtr, ytr)
p = model.predict_proba(Xte)[:, 1]
preds.append(p)
print(f"{name:>12s} AUC: {roc_auc_score(yte, p):.4f}")
kNN AUC: 0.9728
LogReg AUC: 0.9888
Stump AUC: 0.8690
RandomForest AUC: 0.9818
Pairwise testing, with and without correction
compare_models_pairwise takes a tuple of prediction arrays (one per model) instead of just two, plus optional labels to name them. By default it applies the Holm-Bonferroni correction (correction="holm"); we also run it once with correction=None so we can see the raw, uncorrected p-values side by side.
[4]:
pred_tuple = tuple(preds)
naive = stambo.compare_models_pairwise(
yte, pred_tuple, metrics=("ROCAUC", "AP"), labels=labels,
n_bootstrap=2000, seed=SEED, correction=None,
)
corrected = stambo.compare_models_pairwise(
yte, pred_tuple, metrics=("ROCAUC", "AP"), labels=labels,
n_bootstrap=2000, seed=SEED, correction="holm",
)
Let’s line up the raw p-value against the Holm-adjusted one for every pair, for both metrics.
[5]:
for stat in naive:
print(f"=== {stat} ===")
for comparison in naive[stat]:
p_raw = naive[stat][comparison]["p_value"]
p_holm = corrected[stat][comparison]["p_value_adjusted"]
flips = " <-- conclusion flips!" if (p_raw < 0.05) != (p_holm < 0.05) else ""
print(f"{comparison:22s} raw p={p_raw:.4f} Holm-adjusted p={p_holm:.4f}{flips}")
print()
=== ROCAUC ===
kNN / LogReg raw p=0.0340 Holm-adjusted p=0.1019 <-- conclusion flips!
kNN / Stump raw p=0.0010 Holm-adjusted p=0.0060
kNN / RandomForest raw p=0.3468 Holm-adjusted p=0.3468
LogReg / Stump raw p=0.0010 Holm-adjusted p=0.0060
LogReg / RandomForest raw p=0.1689 Holm-adjusted p=0.3378
Stump / RandomForest raw p=0.0010 Holm-adjusted p=0.0060
=== AP ===
kNN / LogReg raw p=0.0050 Holm-adjusted p=0.0150
kNN / Stump raw p=0.0010 Holm-adjusted p=0.0060
kNN / RandomForest raw p=0.0590 Holm-adjusted p=0.1179
LogReg / Stump raw p=0.0010 Holm-adjusted p=0.0060
LogReg / RandomForest raw p=0.2649 Holm-adjusted p=0.2649
Stump / RandomForest raw p=0.0010 Holm-adjusted p=0.0060
For ROC-AUC, kNN / LogReg looks significant at the conventional \(\alpha=0.05\) threshold when you only look at the raw p-value (\(p \approx 0.03\)) – but once we account for the fact that we ran 6 comparisons for this metric, the Holm-adjusted p-value rises above 0.05 and the difference is no longer significant. This is exactly the kind of false positive that a multiple-comparison correction is meant to catch: with 6 pairs tested, seeing one raw p-value around 0.03-0.05 by chance
alone is not actually surprising.
The comparisons involving the weak decision stump, on the other hand, stay significant even after correction – that effect is large enough to survive it.
LaTeX report
Just like to_latex for two-model comparisons, pairwise_to_latex turns a compare_models_pairwise report into a copy-paste LaTeX table – one row per comparison, with both the raw and the Holm-adjusted p-value shown (in parentheses) whenever a correction was applied.
[6]:
print(stambo.pairwise_to_latex(corrected, n_digits=3))
% \usepackage{booktabs} <-- do not forget to have this imported.
\begin{tabular}{lllll} \\
\toprule
\textbf{Comparison} & \multicolumn{2}{c}{\textbf{ROCAUC}} & \multicolumn{2}{c}{\textbf{AP}} \\
& \textbf{Diff [CI]} & \textbf{$p$-value} & \textbf{Diff [CI]} & \textbf{$p$-value} \\
\midrule
kNN / LogReg & $0.016$ [$0.001$-$0.036$] & $0.034$ ($0.102$) & $0.023$ [$0.005$-$0.046$] & $0.005$ ($0.015$) \\
kNN / Stump & $-0.104$ [$-0.145$-$-0.066$] & $0.001$ ($0.006$) & $-0.102$ [$-0.144$-$-0.064$] & $0.001$ ($0.006$) \\
kNN / RandomForest & $0.009$ [$-0.009$-$0.030$] & $0.347$ ($0.347$) & $0.018$ [$-0.000$-$0.041$] & $0.059$ ($0.118$) \\
LogReg / Stump & $-0.120$ [$-0.163$-$-0.080$] & $0.001$ ($0.006$) & $-0.125$ [$-0.172$-$-0.082$] & $0.001$ ($0.006$) \\
LogReg / RandomForest & $-0.007$ [$-0.018$-$0.003$] & $0.169$ ($0.338$) & $-0.005$ [$-0.015$-$0.004$] & $0.265$ ($0.265$) \\
Stump / RandomForest & $0.113$ [$0.076$-$0.154$] & $0.001$ ($0.006$) & $0.120$ [$0.080$-$0.167$] & $0.001$ ($0.006$) \\
\bottomrule
\end{tabular}
Pairwise comparisons also support nested / clustered data
Just like two_sample_test and compare_models, both pairwise_bootstrap_test and compare_models_pairwise accept a groups argument to correctly account for repeated / correlated measurements from the same subject – now across all \(N\) models at once, and combined with the multiple-comparison correction.
To see the effect with real models, we reuse stambo.synthetic.generate_non_iid_measurements (the same generator, and the same parameters, as in Classification_non_iid.ipynb) to build a dataset where the test set contains several correlated measurements per subject – e.g. multiple scans or visits from the same patient. We then train three classifiers (a kNN, a logistic regression, and, on purpose, a weak decision stump) and compare all three pairwise, both ignoring and accounting for the
subject structure.
[7]:
from sklearn.model_selection import StratifiedGroupKFold
data, y_niid, subject_ids = stambo.synthetic.generate_non_iid_measurements(
n_data=300,
n_subjects=100,
rho=0.8, # Subject correlation between classes
subj_sigma=2.0, # Strong subject-level effect relative to noise -> more within-subject correlation
noise_sigma=0.5,
gamma=0.8,
mu_cls_1=[2, 2],
mu_cls_2=[2.1, 2.4],
overlap=0.3, # Most subjects appear in both classes -> repeated, correlated measurements
feat_corr=0.5,
seed=SEED,
)
gss = StratifiedGroupKFold(n_splits=3, shuffle=True, random_state=SEED)
niid_train_idx, niid_test_idx = next(gss.split(data, y_niid, groups=subject_ids))
X_niid_train, X_niid_test = data[niid_train_idx], data[niid_test_idx]
y_niid_train, y_niid_test = y_niid[niid_train_idx], y_niid[niid_test_idx]
groups_test = subject_ids[niid_test_idx]
print(f"Test set: {len(y_niid_test)} rows from only {len(np.unique(groups_test))} unique subjects")
Test set: 100 rows from only 28 unique subjects
[8]:
niid_models = {
"kNN": KNeighborsClassifier(n_neighbors=5),
"LogReg": LogisticRegression(),
"Stump": DecisionTreeClassifier(max_depth=1, random_state=42),
}
niid_labels = list(niid_models.keys())
niid_preds = []
for name in niid_labels:
model = niid_models[name]
model.fit(X_niid_train, y_niid_train)
niid_preds.append(model.predict_proba(X_niid_test)[:, 1])
niid_preds = tuple(niid_preds)
We now run compare_models_pairwise four ways – with/without groups, and with/without the Holm correction – to see how much each one matters on its own, and combined.
[9]:
variants = {
"naive, uncorrected": dict(groups=None, correction=None),
"naive, Holm": dict(groups=None, correction="holm"),
"clustered, uncorrected": dict(groups=groups_test, correction=None),
"clustered, Holm": dict(groups=groups_test, correction="holm"),
}
for tag, kwargs in variants.items():
res = stambo.compare_models_pairwise(
y_niid_test, niid_preds, metrics=("ROCAUC",), labels=niid_labels,
seed=SEED, n_bootstrap=2000, **kwargs,
)
print(f"--- {tag} ---")
for comparison, entry in res["ROCAUC"].items():
p = entry["p_value_adjusted"] if entry["p_value_adjusted"] is not None else entry["p_value"]
flag = " significant" if p < 0.05 else ""
print(f"{comparison:18s} p={p:.4f}{flag}")
print()
--- naive, uncorrected ---
kNN / LogReg p=0.0170 significant
kNN / Stump p=0.4048
LogReg / Stump p=0.0010 significant
--- naive, Holm ---
kNN / LogReg p=0.0340 significant
kNN / Stump p=0.4048
LogReg / Stump p=0.0030 significant
--- clustered, uncorrected ---
kNN / LogReg p=0.6857
kNN / Stump p=0.5287
LogReg / Stump p=0.2899
--- clustered, Holm ---
kNN / LogReg p=1.0000
kNN / Stump p=1.0000
LogReg / Stump p=0.8696
Summary: With the naive, row-level bootstrap, two of the three pairwise comparisons look significant – and Holm correction alone does not fully fix that, because the underlying problem is not the multiple comparisons, it’s that treating ~100 correlated rows from only ~28 subjects as 100 independent measurements makes every comparison look far more confident than the data actually supports. Once we resample whole subjects via groups, every p-value grows sharply, and after also correcting
for the 3 comparisons, none of the three models come out as reliably different from each other. Both problems – ignoring the nested/clustered structure, and ignoring the multiple-comparison problem – can each independently make you overconfident, and here, as is common in practice, the clustering is the dominant effect: getting it right matters at least as much as the correction itself.