Comparing Models: Uncertainty and Significance

Model B scores 0.912 and model A scores 0.907. On 2,000 cases, that is ten additional correct predictions. To assess the evidence, we need to know which cases changed, how the test cases were sampled, and whether either model was selected using those results. Statistical evidence and practical value are separate questions: report the size of the improvement and its uncertainty, then consider what the changed decisions are worth.

What the cross-validation standard error does and does not tell you

A standard error describes how much an estimate would vary across repetitions of the evaluation procedure. A common shortcut divides the standard deviation of the \(k\) fold scores by \(\sqrt{k}\), as if they were independent. Cross-validation fits share data, so covariance between fold scores can affect the variance of their mean. The following simulation compares the shortcut with the observed spread across 300 fresh datasets; it tests this particular setup, not a universal rule.

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score, StratifiedKFold

for n in (200, 600, 2000):
    cv_means, reported_ses = [], []
    for rep in range(300):                      # repeat the whole study on fresh data
        g = np.random.default_rng(rep)
        X = g.normal(size=(n, 8))
        y = (X[:, 0] * 0.9 + X[:, 1] * 0.5 + g.normal(0, 1, n) > 0).astype(int)
        sc = cross_val_score(LogisticRegression(), X, y,
                             cv=StratifiedKFold(10, shuffle=True, random_state=rep))
        cv_means.append(sc.mean())
        reported_ses.append(sc.std(ddof=1) / np.sqrt(10))
    true_sd = np.std(cv_means, ddof=1)          # Monte Carlo estimate of the spread
    print(f"n={n:5d}  empirical sd {true_sd:.4f}   mean reported SE {np.mean(reported_ses):.4f}"
          f"   ratio {np.mean(reported_ses) / true_sd:.2f}")
# n=  200  empirical sd 0.0334   mean reported SE 0.0302   ratio 0.90
# n=  600  empirical sd 0.0168   mean reported SE 0.0173   ratio 1.03
# n= 2000  empirical sd 0.0096   mean reported SE 0.0094   ratio 0.98

The ratios of mean reported SE to empirical standard deviation are 0.90, 1.03, and 0.98. In these runs, the shortcut is about 10% low at \(n=200\) and close to the estimated spread at the larger sizes. Both quantities are Monte Carlo estimates. Agreement of their averages does not establish unbiasedness or correct confidence-interval coverage.

This code does not construct confidence intervals or measure their coverage. To check a nominal 95% interval, specify the target accuracy, construct an interval in each repeated study, and count how often it contains that target. Matching an average standard error is not sufficient: centering, distributional shape, and the variability of the reported standard errors also matter.

Ten-fold CV evaluates fits trained on 90% of the available rows. Its average targets the performance of that training procedure at the smaller sample size, rather than the conditional accuracy of one final model fitted on all rows. Additional training data can improve performance, but neither the direction nor the size of that difference is guaranteed. A fresh test set evaluates the particular final fit; repeated training experiments address variability across fitted models.

Compare on the same test set, and use the pairing

Evaluating both models on the same cases creates paired observations. Let \(D_i=1\) when only A is correct, \(D_i=-1\) when only B is correct, and \(D_i=0\) otherwise. Then the mean of \(D_i\) is their accuracy difference. Positive correlation between the models’ correctness reduces its variance relative to an independent-sample comparison; pairing preserves that relationship instead of assuming it away. The test below assumes independent cases and models fixed before evaluation.

import numpy as np
from scipy import stats

n = 2000
both_right, both_wrong, a_only, b_only = 1500, 300, 120, 80
acc_a = (both_right + a_only) / n
acc_b = (both_right + b_only) / n
print(f"accuracy A {acc_a:.4f}   accuracy B {acc_b:.4f}   difference {acc_a - acc_b:+.4f}")

stat = (abs(a_only - b_only) - 1) ** 2 / (a_only + b_only)  # discordant correctness counts
print(f"McNemar chi2 {stat:.4f}  p {1 - stats.chi2.cdf(stat, 1):.4f}"
      f"   (uses the {a_only + b_only} discordant cases)")

se = np.sqrt(acc_a * (1 - acc_a) / n + acc_b * (1 - acc_b) / n)  # assumes independent samples
print(f"unpaired Wald approximation p {2 * (1 - stats.norm.cdf(abs(acc_a - acc_b) / se)):.4f}")
# accuracy A 0.8100   accuracy B 0.7900   difference +0.0200
# McNemar chi2 7.6050  p 0.0058   (uses the 200 discordant cases)
# unpaired Wald approximation p 0.1137

A is correct on 120 cases that B misses, while B is correct on 80 that A misses. The accuracy difference is \((120-80)/2000=0.02\). McNemar’s null hypothesis is that the two discordant outcomes are equally probable. Its exact conditional version treats A’s wins as a \(\mathrm{Binomial}(200,0.5)\) count; the code uses the continuity-corrected chi-square approximation. The 1,800 concordant cases contribute zero to the paired difference; they are not additional noise and they still determine the denominator used to express the effect as an accuracy gap. The unpaired Wald calculation ignores the dependence and produces a larger p-value here.

For an average per-case loss, compare paired loss differences, using a paired t approximation when its assumptions are suitable or a paired bootstrap. For nonlinear metrics such as F1 or RMSE, resample the same case indices for both models and recompute the complete metrics and their difference. Independent evaluation samples can also be compared with an appropriate independent-sample method, but using the same representative cases often makes the comparison more precise. For grouped or time-dependent observations, account for that dependence too.

How much test data a claim needs

import numpy as np

disagree = 0.10                       # exactly one model is correct on 10% of cases
for n_te in (500, 2000, 10_000, 50_000):
    se_diff = np.sqrt(disagree / n_te)
    print(f"test n={n_te:6d}   approximate gap for 80% power {2.8 * se_diff:.4f}")
# test n=   500   approximate gap for 80% power 0.0396
# test n=  2000   approximate gap for 80% power 0.0198
# test n= 10000   approximate gap for 80% power 0.0089
# test n= 50000   approximate gap for 80% power 0.0040

Power is the probability that a test rejects its null when a specified effect is real. Here 2.8 approximates \(1.96+0.84\), combining a 5% two-sided significance level with 80% power. If \(q=P(|D|=1)\) and the true accuracy gap is \(\Delta=E[D]\), then \(\operatorname{Var}(D)=q-\Delta^2\). For small gaps, the standard error is approximately \(\sqrt{q/n}\). The table uses this normal approximation with fixed \(q=0.10\); it is a planning calculation, not a boundary below which an observed result must be noise.

The disagreement here is disagreement in correctness, not merely different predicted labels. The assumed rate matters: smaller \(q\) can make the same accuracy gap easier to detect, provided the gap is feasible since \(|\Delta|\le q\). For example, the opening 0.005 gap could consist of ten B-only successes and no A-only successes. Under independent cases, its exact two-sided McNemar p-value is \(2(1/2)^{10}\approx0.002\). Thus the two accuracies and the test-set size alone do not determine significance. A p-value measures compatibility with a null model; it is not the probability that the null is true or a measure of practical importance.

AUC and F1 need uncertainty methods suited to their structure. Analytic approaches exist, including DeLong’s covariance method for AUC; a binomial standard-error formula for accuracy cannot simply be reused. In a paired bootstrap, keep each case’s label and both models’ predictions together and draw cases with replacement. Recompute the metric difference in every resample. This estimates test-sampling uncertainty for the fixed fits, not the additional variability from retraining them.

Exercises

1. Does the bootstrap interval actually cover? Build percentile bootstrap intervals for a test-set metric and measure how often they contain the truth.

You should get: coverage close to the nominal level, and a width that shrinks like \(1/\sqrt{n}\).

import numpy as np

true_p = 0.70                                   # the model's real accuracy
for n_te in (200, 1000, 5000):
    widths, hits = [], []
    for rep in range(400):
        g = np.random.default_rng(rep)
        y = g.binomial(1, true_p, n_te)         # 1 = the model got this case right
        boots = [g.choice(y, n_te, replace=True).mean() for _ in range(400)]
        lo, hi = np.percentile(boots, [2.5, 97.5])
        widths.append(hi - lo); hits.append(lo <= true_p <= hi)
    print(f"n={n_te:5d}  mean CI width {np.mean(widths):.4f}"
          f"   coverage {np.mean(hits):.4f}   (nominal 0.95)")
# n=  200  mean CI width 0.1249   coverage 0.9525   (nominal 0.95)
# n= 1000  mean CI width 0.0561   coverage 0.9500   (nominal 0.95)
# n= 5000  mean CI width 0.0252   coverage 0.9550   (nominal 0.95)

2. Twenty models and one baseline. Measure how often at least one of several equally accurate rivals appears significantly different from a baseline in either direction.

You should get: an estimated family-wise false-positive rate near 40% at twenty rivals and 78% at a hundred in this simulation.

import numpy as np

for n_models in (1, 5, 20, 100):
    any_sig, best_gap = [], []
    for rep in range(2000):
        g = np.random.default_rng(rep)
        n = 2000
        base = g.binomial(1, 0.80, n)                # baseline, 80% accurate
        sig, gaps = False, []
        for m in range(n_models):
            rival = g.binomial(1, 0.80, n)           # same population accuracy as the baseline
            d = rival.mean() - base.mean()
            se = np.sqrt((rival.var(ddof=1) + base.var(ddof=1)) / n)
            gaps.append(d)
            if abs(d / se) > 1.96: sig = True
        any_sig.append(sig); best_gap.append(max(gaps))
    print(f"{n_models:3d} rivals: at least one 'significant' {np.mean(any_sig):.4f}"
          f"   best apparent gain {np.mean(best_gap):+.4f}")
#   1 rivals: at least one 'significant' 0.0420   best apparent gain +0.0001
#   5 rivals: at least one 'significant' 0.1840   best apparent gain +0.0103
#  20 rivals: at least one 'significant' 0.4085   best apparent gain +0.0167
# 100 rivals: at least one 'significant' 0.7775   best apparent gain +0.0224

3. What a resampling comparison can establish. Compare rejection frequencies for two fixed-seed configurations of the same random-forest algorithm under naive and corrected tests on repeated train/test splits.

You should get: fewer rejections with the correction. Explain why these rejection counts alone establish neither the null hypothesis nor test power.

import numpy as np
from scipy import stats
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import ShuffleSplit

# Same forest settings, with two fixed seeds.
# This comparison does not prove their population accuracy difference is zero.
n, k, test_frac, n_rep = 400, 10, 0.2, 300
naive_rej, nb_rej = [], []
for rep in range(n_rep):
    g = np.random.default_rng(rep)
    X = g.normal(size=(n, 12))
    y = (X[:, 0] * 1.1 + X[:, 1] * 0.6 + g.normal(0, 1, n) > 0).astype(int)
    d = []
    for tr, te in ShuffleSplit(k, test_size=test_frac, random_state=rep).split(X):
        d.append(RandomForestClassifier(50, random_state=0).fit(X[tr], y[tr]).score(X[te], y[te])
                 - RandomForestClassifier(50, random_state=1).fit(X[tr], y[tr]).score(X[te], y[te]))
    d = np.array(d)
    if d.std(ddof=1) == 0:
        naive_rej.append(False); nb_rej.append(False); continue
    crit = stats.t.ppf(0.975, k - 1)
    t_naive = d.mean() / (d.std(ddof=1) / np.sqrt(k))
    # Nadeau-Bengio: inflate the variance to account for overlapping training sets
    t_nb = d.mean() / np.sqrt(d.var(ddof=1) * (1 / k + test_frac / (1 - test_frac)))
    naive_rej.append(abs(t_naive) > crit); nb_rej.append(abs(t_nb) > crit)

print(f"naive paired t-test rejects at 5%   {np.mean(naive_rej):.4f}")
print(f"Nadeau-Bengio corrected             {np.mean(nb_rej):.4f}")
print(f"variance factor: 1/k + n_te/n_tr = {1/k + test_frac/(1-test_frac):.4f}"
      f"  vs naive 1/k = {1/k:.4f}")
# naive paired t-test rejects at 5%   0.0500
# Nadeau-Bengio corrected             0.0000
# variance factor: 1/k + n_te/n_tr = 0.3500  vs naive 1/k = 0.1000