Statistical Tests for Model Comparison

Notation: k folds, per-fold difference d_i = score_A,i − score_B,i, mean d̄, sample variance s².

#Corrected resampled t-test (Nadeau & Bengio 2003)

t = d̄ / sqrt( (1/k + n_test/n_train) · s² ),   df = k − 1

For k-fold CV, n_test/n_train = 1/(k − 1). The ordinary paired t-test uses s²/k and ignores that training sets overlap across folds, so it rejects far too often (inflated type I error). For repeated k-fold (r repeats): use all r·k differences with (1/(r·k) + n_test/n_train).

import numpy as np
from scipy import stats
d = a - b; k = len(d)
t = d.mean() / np.sqrt((1 / k + 1 / (k - 1)) * d.var(ddof=1))
p = 2 * stats.t.sf(abs(t), k - 1)

#5x2cv paired t-test (Dietterich 1998)

Five repetitions of 2-fold CV; statistic uses the first fold's difference over the averaged variance; df = 5. Conservative and well-calibrated; available in mlxtend.evaluate.paired_ttest_5x2cv.

#Wilcoxon signed-rank

Ranks |d_i|, sums ranks of positive differences. Exact distribution for small k. With 5 folds the smallest possible two-sided p-value is 2/32 = 0.0625: it can never reach 0.05 with 5 folds. Use 10 folds or repeated CV if you want this test. scipy.stats.wilcoxon(d).

#Friedman test and post-hoc

For K models over k folds (or datasets): rank models within each fold, average ranks R_j.

χ²_F = 12k / (K(K+1)) · Σ (R_j − (K+1)/2)²,   df = K − 1

Then post-hoc: Nemenyi (all pairs, critical difference q_α sqrt(K(K+1)/(6k))) or pairwise tests vs a control with Holm correction (Demšar 2006). scipy.stats.friedmanchisquare, scikit-posthocs for Nemenyi.

#Bayesian correlated t-test (Corani & Benavoli 2015; Benavoli et al. 2017)

Posterior of the mean difference μ is Student-t:

μ ~ t_{k−1}( d̄, (1/k + ρ/(1−ρ)) · s² ),   ρ = n_test / (n_test + n_train)  (= 1/k for k-fold)

Note ρ/(1−ρ) = n_test/n_train, the same correction as Nadeau-Bengio. With a region of practical equivalence [−r, r]: P(μ > r) = A better, P(|μ| <= r) = equivalent, P(μ < −r) = B better. Decide when one probability exceeds e.g. 0.95. Implemented in compare_models.py; also in the baycomp package.

#DeLong test for correlated ROC-AUCs

AUC = mean of placement values. For each positive i: V10_i = share of negatives it outranks (ties ½); for each negative j: V01_j = share of positives that outrank it. With S10, S01 the covariance matrices of these components across the two models:

Var(AUC_A − AUC_B) = (S10_AA + S10_BB − 2 S10_AB)/m + (S01_AA + S01_BB − 2 S01_AB)/n
z = (AUC_A − AUC_B) / sqrt(Var)

paired_predictions.py uses the midrank (Sun & Xu 2014) form; it matches the brute-force definition to machine precision on the fixture.

#McNemar's test

Only discordant pairs matter: b = A right & B wrong, c = A wrong & B right.

• Exact: two-sided binomial test of b vs Binomial(b + c, 0.5) (use when b + c < ~50).
• Chi-square with continuity correction: (|b − c| − 1)² / (b + c), df 1. statsmodels.stats.contingency_tables.mcnemar([[n00, b], [c, n11]], exact=True).

#Paired bootstrap

Resample test rows with replacement (stratified by class for rare positives), compute both models' metric on the same resample, record the difference; 95% CI from the 2.5/97.5 percentiles. Works for any metric (PR-AUC, log loss, recall@K, MAE). 1,000-2,000 resamples are usually enough for a CI; more for small p-values.

#Choosing among them