Test the hypothesis that the performances of k machine learning models on the same test set is the same 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.

*

*Is there a distribution-free, exact test that allows me to test the above hypothesis?

*If there's no such test, is there at least a distribution-free, approximate test?

*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


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.
