(How) should multiple comparison correction be applied for hypothesis testing of cross validation performance?

Question

I Would like to compare if results from two models are significantly different.

The models have been trained on the same samples in a K-fold cross validation setup, so both will spit out K performance scores, for which we can test is the mean of the performance is significantly different with a related T-test.

At first, setting the significance level alpha at 0.05, my results are not significant:

>>> from scipy.stats import ttest_rel
>>> np.random.seed(12)
>>> n_splits = 5
>>> performance_model_A = np.random.normal(0.5, 0.2, n_splits)
>>> performance_model_B = np.random.normal(0.6, 0.1, n_splits)
>>> _, pval = ttest_rel(performance_model_A, performance_model_B)

At 5-fold cross validation the p-value is 0.180

However, if you add repeats to the cross-validation set up, I end up with a p-value lower than alpha:

>>> n_repeats = 10
>>> performance_model_A = np.random.normal(0.5, 0.2, n_splits * n_repeats)
>>> performance_model_B = np.random.normal(0.6, 0.1, n_splits * n_repeats)
>>> _, pval = ttest_rel(performance_model_A, performance_model_B)
>>> print(f'At {n_splits}-fold, {n_repeats}-repeat cross validation the p-value is {pval:.3f}.')

At 5-fold, 10-repeat cross validation the p-value is 0.002.

Now, I think we should perform multiple comparison correction because I used the same sample multiple times (e.g. using Bonferroni correction for n_repeats and n_splits, reducing alpha to 0.01 and 0.001 respectively.

But I am reluctant, since all I am comparing is the performance of the models on a given sample of data, so all results are further evidence of differences between the models (which is what we are testing). However, this leads to the bizarre conclusion that -with sufficient repeats- any model is significantly different from another, so then I might not be performing the right test for this application.

Alternatively, I could just strike it in the middle: correct for n_repeats and not for n_splits, since this way I correct for time a sample is in the test data partition more than once. But would not be backed by any strong statistical insight.

I have not been able to find a definitive advice online/best practice/papers on multiple comparison testing for comparison of cross-validation results, and any help is greatly appreciated.

---

Thanks, @Firebug for the reference to Bouckaert and Frank's corrected repeated k-fold cv test.

I could not find a Python implementation for this method, so finally I would like to share the following Python solution if that is allowed:

def corr_rep_kfold_cv_test(a:list, b:list, n_splits:int, n_samples:int) -> tuple[float, float]:
    """
        Implementation of Bouckaert and Franks (2004) corrected repeated k-fold cv test.
    """
    r  = len(a) // n_splits                             # number of r-times repeats as integer
    k  = n_splits                                       # number of k-folds as integer
    n1 = n_samples // n_splits * (n_splits - 1)         # number of instances used for training
    n2 = n_samples // n_splits                          # number of instances used for testing
    x  = np.subtract(a, b)                              # observed differences
    m  = x.mean()                                       # mean estimate
    s  = np.sum((x - m) ** 2) / (k * r - 1)             # variance estimate
    t_stat = m / np.sqrt((1 / (k * r) + n2 / n1) * s)   # corrected test statistic
    p_val = stats.t.sf(np.abs(t_stat), r) * k           # p-value
    return t_stat, p_val

Answer 1

See the linked answer to this question: Comparing two models using Repeated K-fold Cross Validation

Basically, in [1], authors propose a correction for interdependence to the t-statistic in repeated K-fold cross-validation (see section 3.3):

$$t = \frac{ \frac{1}{k \cdot r} \sum_{i=1}^k \sum_{k=1}^r x_{ij} }{ \sqrt{\left(\frac{1}{k\cdot r}+\frac{n_2}{n_1}\right)\hat\sigma^2} },$$

$$\hat\sigma^2=\frac{1}{k \cdot r - 1} \sum_{i=1}^k \sum_{k=1}^r (x_{ij} - \bar {x_{ij}})^2,$$

$$\bar {x_{ij}}=\frac{1}{k \cdot r} \sum_{i=1}^k \sum_{k=1}^r x_{ij},$$

where $n_1$ is the number of instances used for training, $n_2$ is the number of instances used for testing, $k$ is the number of folds, $r$ is the number of repetitions, $x_{ij}$ is the difference in performance between the two models.

[1]: R. R. Bouckaert and E. Frank, ‘Evaluating the Replicability of Significance Tests for Comparing Learning Algorithms’, in Advances in Knowledge Discovery and Data Mining, vol. 3056, H. Dai, R. Srikant, and C. Zhang, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2004, pp. 3–12. doi: 10.1007/978-3-540-24775-3_3.

Answer 2

Not an entire solution, but some thoughts I have about this situation.

(Repeated) cross validation has (at least) two different sources of variance:

• variance uncertainty due to case (I'll use n for number of statistically independent cases), and
• variance due to model instability.

Looking at more surrogate models (n_repeats * k) will reduce the part of the variance uncertainty on the final estimate that is due to model instability. But the total number of tested cases stays the same (n) after the first complete run. That part of the variance uncertainty can only be reduced by more cases, more surrogate models cannot possibly help.

---

There's a further consideration: cross validation estimates for generalization error are often used as approximation of the generalization error of the model trained on all cases. This is the case when the task at hand is building a model on the data set at hand for application/production use. As opposed to comparing the performance of training algorithms for the given type of data (in that case, there's the problem that only part of the relevant variance components can be assessed by cross validation experiments - see Y. Bengio, Y. Grandvalet, No Unbiased Estimator of the Variance of K-Fold Cross-Validation, J. Mach. Learn. Res. 5 (2004) 1089–1105.)

For the production-use scenario, we say that the variation we observe between our estimate and any single surrogate model could serve as approximation for the variance due to instability in the training between the average performance at n and the single model we then train on the full data set. That means, while we can say that relevant variance uncertainty due to the finite number of tested cases does down with $\frac{1}{n}$, no such reduction can be claimed for the model-instability part of the variance.