Test the hypothesis that the performances of k machine learning models on the same test set is the same

Question

I have $k$ Machine Learning models, trained on the same training set, and I want to test the hypothesis that their performance on a fresh test set is the same. See EDIT 1 below for what I mean with performance. I'm looking for frequentist answers. Also, I'd prefer an answer that refers to a paper or preprint on this topic (testing the null that $k$ ML models have the same performance on a finite-sample test set). The models are Computer Vision models.

1. Is there a distribution-free, exact test that allows me to test the above hypothesis?
2. If there's no such test, is there at least a distribution-free, approximate test?
3. As a worst-case scenario, I could be contented with a parametric, asymptotic test. But it seems to me a terrible option: this is ML, so we know nothing about the models, and we know nothing about the test set data distribution (apart from the fact that it's an iid sample). Yeah, I know "if the size of the test set $n_{test}\to\infty$ we can assume that...". Point is, $n_{test}$ isn't going anywhere. Acquiring and labelling the test set was very expensive, and I'm not getting new labeled data anytime soon. Thus, I'd rather get an answer that doesn't rely on the size of the test set being very large.

EDIT 1: the models are Computer Vision models, meaning that the test set $S=\{(X_i,Y_i)\}_{i=1}^{n_{test}}$ is made of couples (image, structured label). The problem is that, for each image $X_i$, the corresponding label $Y_i$ is not just a categorical variable (we're not doing image classification here), but it's a list of $n_i$ categorical variables (the classes of each of the $n_i$ objects presents in image $X_i$) and $4n_i$ real numbers (the coordinates of the bounding boxes of each object: see here for details). Thus, it's not so simple to say if, for image $X_i$, the prediction $f(X_i)$ of the model is right (1) or wrong (0). I naively thought of solving this using some aggregate metric, since for each of the $k$ models, we can compute a loss over the whole test set, which is a nonnegative number. However, as correctly noted by Jacques Wainer, this approach doesn't work for hypothesis testing. Let's consider then a different setting, where each of the $k$ models (each of the $k$ treatments, in statistical parlance) returns a real number in $[0,1]$ for each of the $n_{test}$ elements of the test set. Thus, for a single metric, our dataset is a matrix $M$ of $k\times n_{test}$ numbers in $[0,1]$. See EDIT 3 for the multiple metrics case.

EDIT 2: there's something that I think was overlooked in comments & answers. I think the $k\times n_{test}$ entries of $M$ are not iid. All models are trained on the same training set, and are tested on the same images, thus I would expect the model results for the a given image $f_1(X_i),\dots,f_k(X_i)$ to be highly correlated.

EDIT 3 (hopefully the final one): let's consider now more than one metric for each image (the multiple hypotheses case I mentioned in earlier versions of this question). For $d$ metrics, I'd get $d$ matrices $M_1,\dots, M_d$ of size $k\times {n_{test}}$. As duly noted by Christian Hennig, different metrics for the same image are likely to be correlated (see also EDIT 2 above), and thing soon get very complicated. So let's forget about multiple hypothesis testing. I'll accept Christian Hennig's suggestion:

decide in advance to use only one test as "major test of interest", interpret that one, and give only a rough "exploratory" interpretation of all the other tests in relation to the major test

Answer 1

Statistical tests, any of them, work on sets of data, not on a single data (for each of the $k$ ML models). Thus, if you compute any measure on the whole test set (let us call them aggregating measures) be it precision, recall, or IoU, you will get a single number for each of the ML algorithms, and there is no statistical test that can receive a single number (for each treatment) and compute a p-value.

So, you cannot compute singe measure for each algorithm on the test set. The only set of data that you have is whether, for each data point in the data set, a particular algorithm got the prediction on that data right or wrong. Thus for each algorithm you have a set of binary data (0 or 1, correct or incorrect, and so on) on each of the data points in the test set.

These measures are paired, or in the parlance of statistical tests, it is blocked - for the same data point in the test set you have the corresponding binary outcome (right or wrong) for each algorithm.

Therefore, you want a test for binary variables (right or wrong), on multiple treatments, and blocked.

The only one I know is the Cochran's Q test. The test is distribution free (but I am not sure if it is exact).

If the p-value is high enough, you can conclude that all algorithm are equally correct, and thus (I believe) any summary measure, such as precision, recall, accuracy, will be "statistically equivalent" (there is no such a thing, but I think given that the Q test tells you that there is no statistical significant difference among the output of all algorithms I believe one can conclude that there is "no difference" for the aggregating measures).

Answering the EDITs:

EDIT1: if the output for each data point is a number (for example between 0 and 1 as you suggested but this works for any number) then you are in luck. What you have is a set of numbers and not of 0/1 and there are many more statistical tests for non-binary number data.

The usual test in machine learning is the one proposed by Demsar (Statistical Comparisons of Classifiers over Multiple Data Sets) in https://www.jmlr.org/papers/volume7/demsar06a/demsar06a.pdf Remember that in your case the multiple datasets of the paper is the multiple datapoints in your test set. Demsar proposed a Friedman test followed by the Nemenyi test to determine which algorithm is significantly different than the others. Since you are hoping that all algorithms are equivalent, if you are lucky the Friedman test will result in the p-value high enough (be mindful of the caveats listed in the comments for my answer). There are implementations of these tests in both Python and R (at least).

Garcia e Herrera ( An Extension on" Statistical Comparisons of Classifiers over Multiple Data Sets" for all Pairwise Comparisons.) https://www.jmlr.org/papers/volume9/garcia08a/garcia08a.pdf proposed other post-hoc tests (beyond the Nemenyi test).

EDIT 2: The data used in the tests are i.i.d. The fact that the algorithms are trained on the same set is not a problem. The conclusion will take that into consideration - your conclusion is that the algorithms when trained on the that same training set, and tested on the that same test set are statistically significantly different or are not statistically significantly different. Your conclusions will be for that particular training and test sets.

EDIT 3: I don't know about multiple aggregating metrics. But first, you are right that they are likely not i.i.d Second, there will be very few data for the statistical test. Say you will use 5 or 10 aggregating metrics, that will leave you with only 5 or 10 data for each algorithm. With so few data, the statistical tests will not likely find that the differences are significative!

Answer 2

Under the assumption that the $k$ ML algorithm performances are "the same", which I will interpret as "independent and identically distributed", on any given test, they are also exchangeable within the test. We can therefore construct a permutation test that preserves potential differences between the tests. Permutation tests are distribution-free, exact, and most powerful conditional upon the data, given the exchangeability of the elements being permuted. They have the additional advantage of allowing any test statistic to be used, as we do not need to be able to calculate the distribution of the statistic under a parametric null hypothesis.

Let us assume an appropriate test statistic $T(X)$, where $X$ is the $k \times n$ matrix of the test results. Within each column (test), the observations are assumed exchangeable (across algorithms) under the null hypothesis; we can therefore in principle permute their row (algorithm) indices without changing the distribution of $T$. We should not permute across the columns, as we are not assuming that the test results are exchangeable given a particular algorithm - the test has an influence on the score. A permutation test would calculate all possible combinations of permutations for each test, and for each combination $k$ would calculate the test statistic $T_k$. We would then compare the test statistic $T_0$ as calculated on the observed data with the $T_k$ test statistics calculated on the permuted data, and find a p-value in an obvious way based upon what fraction of the $T_k$ test statistics exceeded $T_0$.

As a simple example, we can construct a five-"algorithm" ten-"test" problem, where algorithms 1-4 have the same within-test distributional performances and algorithm 5 is mildly better (in this case, has a smaller result on the test.) The tests themselves are slightly different as well. We construct the matrix so that the test results are in $(0,1)$, close to your problem's results range:

X <- matrix(0,5,10)
for (test in 1:10) {
  X[1:4, test] <- rbeta(4, 5-0.2*test, 5+0.2*test)
  X[5, test] <- rbeta(1, 4-0.2*test, 6+0.2*test)
} 
colnames(X) <- paste("Test", 1:10)
rownames(X) <- paste("Algo", c("A","B","C","D","E"))
round(X,4)
       Test 1 Test 2 Test 3 Test 4 Test 5 Test 6 Test 7 Test 8 Test 9 Test 10
Algo A 0.7210 0.3176 0.3773 0.2571 0.3365 0.4364 0.5966 0.2157 0.2142  0.2992
Algo B 0.4153 0.4240 0.4542 0.1841 0.6346 0.4591 0.3750 0.1338 0.5906  0.2124
Algo C 0.7269 0.3499 0.5192 0.6179 0.4444 0.4069 0.2565 0.1840 0.4877  0.2929
Algo D 0.2334 0.8319 0.2045 0.3491 0.2584 0.2575 0.5532 0.3079 0.3876  0.4806
Algo E 0.2510 0.3357 0.3252 0.3784 0.1769 0.1322 0.5008 0.2208 0.3242  0.2482

Now to construct a test statistic. A plausible statistic might involve normalizing each column (as the tests are known to be different,) taking the row means of the normalized observations, and seeing whether some function of the row means, e.g. the minimum, was larger or smaller (whichever is appropriate) than we expect. The construction of a suitable test statistic can be quite important to your overall results, so it is worth spending some time thinking about.

The complexity of generating each possible permutation exactly once is too great (for me, in this case), so we will adopt a simpler approach; we will generate many, many permutations of the data instead:

X_norm <- X
for (test in 1:10) {
  X_norm[,test] <- X_norm[,test] / mean(X_norm[,test])
}

test_statistic <- function(x) min(rowMeans(x))

observed_statistic <- test_statistic(X_norm)

permutation_results <- rep(0, 10000) 
for (i in seq_along(permutation_results)) {
  for (test in 1:10) {
    X_norm[, test] <- sample(X_norm[, test])
  }
  permutation_results[i] <- test_statistic(X_norm)
}

hist(permutation_results)
abline(v = observed_statistic)
mean(observed_statistic >= permutation_results)
[1] 0.0745

We interpret the final number (0.0745) as the (approximate) permutation p-value of the test statistic.

We have taken a few shortcuts in the code that do not affect the results, for example, as we are permuting each column's values, we don't need to re-normalize the data at each permutation step (which we would have to do were we, say, sampling from the column with replacement,) and we can just permute the previously-permuted columns without going back to the original data at each permutation step.

Note that, in this example, the overall test procedure isn't that powerful in an absolute sense, as we only have a few tests and algorithms, and four of the five algorithms had the same expected performance. (As it happens, in 1,000 repeats of the experiment, we rejected the null hypothesis 35.4% of the time at the 95% level of confidence and 16.6% of the time at the 99% level of confidence.) Your mileage will vary! But the ability to choose a custom test statistic combined with the distribution-free, exact, and most powerful (conditional upon the data) properties of permutation tests make them an extremely useful tool in the applied statistician's toolbox.

References for permutation tests aren't as common as one might hope; they are usually discussed as part of more general nonparametric statistics books. Permutation, Parametric, and Bootstrap Tests of Hypotheses by P. Good is my favorite among the more specialized books.

Answer 3

Have a look at the Diebold-Mariano test (DM-Test): here, the forecast accuracy of 2 (arbitrary) models can be compared using their losses, e.g. MSE or RMSE. Let $L_{A,t}$ be the loss of model A at time t. If the difference between the losses $\Delta_{A,B,t}=L_{A,t}-L_{B,t}$ is significantly greater than zero then the forecast from B is better than A. Otherwise not. Let $\hat \varepsilon_A=[\hat \varepsilon_{A,T+1} \dots \hat \varepsilon_{A,T+N}]'$ be the vector of forecasting errors for model A for N-step-ahead forecasts. Define e.g. $L_{A,t}=|\hat \varepsilon_{A,t}|$ (MAE) or $L_{A,t}=|\hat \varepsilon_{A,t}|^2$ (MSE). Then $\bar \Delta_{A,B}=\frac{1}{N}\sum_{i=1}^N\Delta_{A,B,i}$ is the mean difference. The DM test statistic is

$$ t_{DM} = \frac{\bar \Delta_{A,B}}{\sigma(\bar \Delta_{A,B})} \overset{a}{\sim} \mathcal{N}(0,1) $$

under the null that the difference is not significant (i.e. the 2 models are the same), where $\sigma(\bar \Delta_{A,B})$ is the standard deviation of $\bar \Delta_{A,B}$.

Usually, the DM-Test is used to pairwise compare multiple models in a 2D matrix and can easily be extended to forecasting vectors (instead of point-forecasts). It is frequently used in electricity price forecasting. See Figure 3 of this paper (and the corresponding sections) for a very good example of how the test is employed. Also see the dm.test() R function.