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.

References

• DeLong, DeLong, and Clarke-Pearson (1988). Comparing the areas under two or more correlated receiver operating characteristic curves: a nonparametric approach
• Nadeau and Bengio (1999). Inference for the Generalization Error, Section 3