Statistical test for comparing performance metrics of two regression models on a *single* validation dataset?

Question

What is an appropriate way to apply a statistical hypothesis test for evaluating model performance metrics (e.g. MAE, MAPE, MSE) between two regression models based on a single holdout (validation) dataset?

Possible approaches

What I'd like to do is just apply a paired t-test on the observations directly (at least for measures like absolute error (AE) or absolute percent error (APE) where there are meaningful corresponding observation level measures). The problem with this though seems like we would be at risk of violating some of the assumptions of a paired t-test, e.g.

• dependent variable might not be normally distributed (as they would be when comparing the average of the metrics of cross-validaon splits, see [Context], which we can trust in many metrics due to central limit theorem)
• will not be safe from outliers

I guess though how much do I need to worry about these (isn't the t-test pretty robust)? Or what adjustments can we make to the test to still make it valid? Would using a Wilcoxan signed-rank test be a better way to go... (tweet about this)?

Preference is for approaches that are simple, reasonably appropriate / generalizable and use performance on holdout data that do not depend on building the model multiple times...

Context

A common model comparison technique based on predictions on hold-out data is to build multiple models each using k-fold cross validation, calculate the average performance for each on the k holdout sets for each model, and then compare performance between the model's using a paired t-test (e.g. from the {tidymodels} team at Rstudio see Feature Engineering..., 3.7 or Tidy Models with R, 11.3).

I use this procedure frequently and love how simple it is, however a few concerns I have with it are:

1. cannot be used when doing a simple training-validation split for evaluation (as only have one hold-out set)
2. the power of the test seems like it could be dependent on the number of splits, k
3. Post on Statistical Significance Tests for Comparing Machine Learning Algorithms critiques the methodology for a lack of independence (due to the training sets being largely similar between k models) though still notes the approach is fine and may just have a little high type 1 errors -- this point doesn't bother me too much as I primarily am using the test as an indicator.

This question is in regards to point 1. where we are doing evaluation using a single validation set. I found a lot of cross-validated posts on statistical tests in model comparisons, but most are in the context of cross-validation or other resampling procedures.

Related / follow-up questions

• Similarly, I recently tweeted about adding confidence intervals on model performance estimates in a validation set using simple methods (based on calculating the standard error), is this approach fine or problematic?
• In the case when doing k-fold cross-validation would it be OK to do a paired t-test on the errors of the individual hold-out predictions from the training set (i.e. run the test on the sample of n observations, rather than the sample of k average errors of each split). However the point in 3 and the concerns from the bullet points give me even greater pause here... (Also, similar to above note, some metrics cannot be calculated at the observation level, e.g. in the classification context the AUC under ROC curve).

Answer 1

A problem with evaluation on a single validation set is that there is sample-to-sample noise that is ignored by this. Based on resampling in general, this can be the largest source of variation.
-- Max Kuhn (paraphrased)

That said, you could bootstrap the predictions to get confidence intervals for the performance metrics on the validation set. Example below bootstraps the difference in absolute errors between two models (per original question).

library(tidyverse)
library(tidymodels)

mpg <- mutate(mpg, cyl = as.factor(cyl))

set.seed(213)
split <- initial_split(mpg)
train <- training(split)
test <- testing(split)

mod_tree <- decision_tree(mode = "regression") %>% 
  set_engine("rpart") %>% 
  fit(cty ~ hwy + displ + drv + cyl, data = train)

mod_lm <- linear_reg(mode = "regression") %>% 
  set_engine("lm") %>% 
  fit(cty ~ hwy + displ + drv + cyl, data = train)

preds <- bind_cols(
  predict(mod_tree, test) %>% 
    rename(.pred_tree = .pred),
  predict(mod_lm, test) %>% 
    rename(.pred_lm = .pred),
  select(test, cty)
) %>% 
  mutate(across(contains(".pred"), list(resid = ~.x - cty))) %>% 
  mutate(diff_abs_resids = abs(.pred_lm_resid) - abs(.pred_tree_resid))

### Bootstrap confidence interval on diff of absolute errors
n_preds <- nrow(preds)

set.seed(1234)
resamples <- tibble(diff = map_dbl(1:10000, ~mean(sample(preds$diff_abs_resids, size = n_preds, replace = TRUE))))
# see gist to MUCH slower alternative approach to producing `resamples`:
# https://gist.github.com/brshallo/1349cc19e616e5d1df03cc7efe77e67a

conf_95 <- summarise(resamples, band = quantile(diff, c(0.025, 0.95)))

resamples %>% 
  ggplot(aes(x = diff))+
  geom_histogram()+
  geom_vline(xintercept = mean(resamples$diff))+
  geom_vline(aes(xintercept = band, colour = "95% confidence"), data = conf_95)+
  labs(title = "Difference Absolute Errors Linear Model vs. Decision Tree",
       subtitle = "Linear Model Has Significantly Smaller Absolute Errors (10000 Resamples)",
       caption = "Black line represents mean difference.")
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Similar to results we see from running a paired t-test on the absolute errors
t.test(abs(preds$.pred_lm_resid), abs(preds$.pred_tree_resid), paired = TRUE)
#> 
#>  Paired t-test
#> 
#> data:  abs(preds$.pred_lm_resid) and abs(preds$.pred_tree_resid)
#> t = -3.1764, df = 57, p-value = 0.002407
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#>  -0.7652232 -0.1734593
#> sample estimates:
#> mean of the differences 
#>              -0.4693413

^{Created on 2021-04-27 by the reprex package (v1.0.0)}

This result seems to suggest significant differences, however (per quote from Max) it is sensitive to sample-to-sample variation. As an example, if you change the seed before splitting the data from set.seed(213) to set.seed(214) the results will no longer appear significant (p-value on paired t-test goes from 0.002 to 0.185).