# Statistical Agnostic Regression: a machine learning method to validate
regression models

J.M. Gorriz, J. Ramirez, F. Segovia, F. J. Martinez-Murcia, C. Jiménez-Mesa
and J. Suckling.
Data Science and Computational Intelligence Institute
University of Granada
Granada, Spain<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Regression analysis is a central topic in statistical modeling, aiming to
estimate the relationships between a dependent variable, commonly referred to
as the response variable, and one or more independent variables; i.e.,
explanatory variables. Linear regression is by far the most popular method for
performing this task in several fields of research such as prediction,
forecasting, or causal inference. Classical methods for solving linear
regression problems such as Ordinary Least Squares, Ridge, or Lasso
regressions, are often the foundation for more advanced machine learning (ML)
techniques that have been successfully applied although without a formal
definition of statistical significance. At most, permutation or analyses based
on empirical measures (e.g., residuals or accuracy) have been conducted that
leverage the greater sensitivity of ML estimations for detection. In this
paper, we introduce Statistical Agnostic Regression (SAR) for evaluating the
statistical significance of an ML-based linear regression based on
concentration inequalities of the actual risk (expected loss) using the
analysis of the worst case. To achieve this goal, we define a threshold that
establishes that there is sufficient evidence with a probability of at least
$1-\eta$ to conclude that there is a linear relationship in the population
between the explanatory (feature) and the response (label) variables.
Simulations in only two dimensions demonstrate the ability of the proposed
agnostic test to provide an analysis of variance similar to that from the
classical $F$ test for the slope parameter.

_K_ eywords Ordinary Least Squares, $\cdot$ K-fold cross-validation $\cdot$
linear support vector machines $\cdot$ statistical learning theory $\cdot$
permutation tests $\cdot$ Upper Bounding..

##  1 Introduction

Ordinary Least Squares (OLS) is the most popular method to perform linear
regression analysis due to its optimal statistical properties assuming the
linear model:

$y=\upbeta_{1}^{T}\mathbf{x}+\beta_{0}+\epsilon$ (1)

where $y$ is the response variable or observation, $\mathbf{x}$ is the
$P\times 1$ explanatory variable or predictor, $\upbeta_{1}$ and $\beta_{0}$
are unknown parameters, slope and intercept, respectively, that define the
linear regressor (hyperplane) and $\epsilon$ is a random variable with zero
mean and variance $\sigma^{2}$. The model above can be rewritten if we define
$\upbeta=[\upbeta_{1},^{T}\beta_{0}]^{T}$ and
$\hat{\mathbf{x}}=[\mathbf{x}^{T},1]^{T}$ as:

$y=\upbeta^{T}\hat{\mathbf{x}}+\epsilon$ (2)

Using a set of observations $\mathbf{y}=[y_{1},\ldots,y_{N}]^{T}$ and the
predictor matrix
$\mathbf{X}=[\hat{\mathbf{x}}_{1},\ldots,\hat{\mathbf{x}}_{N}]^{T}$, the OLS
estimator is obtained by minimizing the sum of squares
$||\mathbf{y}-\mathbf{X}\upbeta||^{2}_{2}$. Following the Gauss-Markoff
theorem, the best linear unbiased estimator of any linear function is obtained
by the OLS [Zelen62].

In pattern recognition using machine learning (ML), the regression algorithm
is formulated by incorporating a regularizer $||\upbeta||^{2}$ into the
minimization of the risk functional. In this context, the objective is to find
a linear function $f(\hat{\mathbf{x}})=\upbeta^{T}\hat{\mathbf{x}}$ that
minimizes the expected loss or actual risk
$\mathcal{R}(f)=E[\mathcal{L}(f,\mathbf{x},y)]$, considering two terms. One
term is akin to the OLS approximation and is referred to as the empirical
risk, while the other is associated with model complexity. This insight is
fundamental to the statistical learning theory (SLT), which serves as a strong
foundation for all current artificial intelligence (AI) approaches [Grohs22].

Various approaches have been proposed to validate regression models. [Snee77]
recommended comparing model predictions and coefficients with theoretical
expectations, collecting new data, and employing techniques such as data
splitting or cross-validation [Kohavi95]. In [Kleijnen96] the use of
bootstrapping for validating metamodels was introduced, identified as a potent
methodology whilst [Miller91], concentrated on logistic regression models and
put forth a comprehensive approach that incorporated measures of goodness-of-
fit, calibration, and refinement. Finally, in [Oredein11] a comparative
analysis of the predictive index accuracy between data splitting and residual
resampling bootstraps was conducted, concluding that bootstrapping yields a
more precise estimate of the coefficient of determination.

The main drawback of current AI approaches for classification and regression
[Lecun15, Gorriz2020] is their lack of rigorous analysis of significance in
comparison to classical approaches such as the analysis of residuals in OLS
linear regression using hypothesis testing. These approaches often limit their
analyses to the use of permutation testing on empirical measures derived
during the training stages, such as p-value analysis using cross-validation
(CV). For instance, numerous commentaries in the neuroimaging literature
[Varoquaux18, Eklund16, Jollans2019, Gorriz18, Braga04] highlight the high
variability of performance across CV folds in various analytic designs, with
clear implications for predictive inference.

Another explored possibility is to consider classical analysis on the
residuals produced by these approaches, which heavily relies on the
assumptions of classical statistics, such as Gaussianity, as indicated in
Equation 1. Nonetheless, deep feature extraction analyses are typically
conducted in low-dimensional spaces using linear classifiers [vandermaaten08,
Jimenez22, Gorriz2021]. More critically, in pattern recognition, and
particularly in regression problems, there is a significant concern about
formulating machine learning analyses exclusively based on learning curves
derived from selected loss functions [Viering23] that merely demonstrate that
the learning algorithm is converging to a potentially unreliable solution.

In this paper, we propose the use of SLT to formulate a statistical test for
assessing the significance of ML regression models. Initially, we establish an
upper bound on the expected loss of a (linear) support vector regressor under
the worst-case scenario. Subsequently, we compare the actual risk with that
obtained under the null hypothesis, similar to the 50% rate in a
classification problem, meaning there is no linear relationship between the
predictors and the observed variables. Whenever the corrected risk is less
than this threshold, we reject the null hypothesis and conclude that, with at
least a probability $\eta$, there is a linear relationship between the
predictor and observation.

##  2 Background in Support Vector regression: theory and practice

The Support Vector Regression (SVR) algorithm can be generalized to the case
of regression estimation [Vapnik95], where the sparseness property is
preserved by the definition of the $\epsilon$-insensitive loss function:

$|y-f(\mathbf{x})|_{\epsilon}=\max\\{0,|y-f(\mathbf{x})|-\epsilon\\},$ (3)

where the parameter $\epsilon$ is automatically computed [Scholkopf95]. The
expected loss to be minimized is the regularized risk functional:

$\frac{1}{2}||\upbeta||^{2}+\frac{c}{N}\sum_{i=1}^{N}|y_{i}-f(\mathbf{x_{i}})|_{\epsilon}$
(4)

Observe that $c$ is the trade-off constant between complexity and training
error. The minimization of equation 4 is equivalent to a constrained
optimization problem detailed elsewhere [Burges98]. It is important to note
that the choice of loss function depends on the specific regression problem at
hand [Huber64]. In this paper, we will employ the two most commonly used
losses: the squared loss (also known as $\mathcal{L}_{2}$) and the absolute
value loss or $\mathcal{L}_{1}$ (the insensitive loss for $\epsilon$=0).
However, previous studies have indicated that when additive noise deviates
from Gaussianity, better approximations to the regression problem are achieved
by using estimators based on alternative loss functions, beyond the quadratic
loss function employed in OLS [Huber64, Vapnik95]. To make the two losses
comparable, we will employ the rounded quadratic loss in the case of the OLS
algorithm whilst ML models are using the least modulo method, e.g.,
$\epsilon$-insensitive loss. Finally, we point out that the quadratic problem
associated with equation 4 is connected to support vector classification
[Scholkopf95]. Furthermore, the hyperplane that is constructed by this
minimization, and lies close to as many of the data points as possible, is a
self-supervised training strategy to classify the samples outside the
$\epsilon$-tube.

###  2.1 Flatness and the empirical loss

As an example, and connected to the novel test proposed in the next section,
let’s assume we are provided with a dataset of observations $y$ and predictors
$x$ in 2 dimensions. A learning algorithm fits a linear function
$f(x)=\beta_{1}x+\beta_{0}$ using a loss function that penalizes the
difference between the observation and prediction. Starting from an i.i.d. 2D
Gaussian distribution, we can define scaling and rotation transforms
$\mathbf{T}=\mathbf{S}\mathbf{R}$, where $\mathbf{S}=\text{diag}(1,1-\tau)$,
where $\tau\in(0,1)$ is the correlation level, and
$\mathbf{R}=\begin{pmatrix}\cos(\theta)&-\sin(\theta)\\\
sin(\theta)&\cos(\theta)\end{pmatrix}$. These transforms are applied to the
data to create a non-diagonal covariance distribution, i.e. a linear
relationship (see figures 1(a) and 1(b)).Then, we assess the individual losses
and the expected loss (average) of the problem in the same figures.

From these examples, we readily see that with decreasing correlation level,
the flatness property of the OLS and SVR solutions dominate over the
minimization of the empirical risk that, with uncorrelated data, provides an
extreme value. For instance, within the framework of the OLS method, the
expression for the squared loss function, denoted as $\mathcal{L}_{2}$, is
given by $(f(\mathbf{x})-y)^{2}$. Recall that $\boldsymbol{\beta}$ is
determined by the OLS solution
$\boldsymbol{\beta}=(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{y}$
when the rank of the design matrix $\mathbf{X}$ is equal to the number of
observations $N$ [Zelen62]. The solution becomes null
($\boldsymbol{\beta}=\mathbf{0}$) when $\mathbf{X}^{T}\mathbf{y}=\mathbf{0}$,
indicating a perfect orthogonality issue in the dataset, which is the opposite
of multicollinearity. Therefore, the empirical loss simplifies to
$\mathcal{R}_{N}=(\mathbf{y}-\mathbf{X}\upbeta^{T})(\mathbf{y}-\mathbf{X}\upbeta)=\mathbf{y}^{T}\mathbf{y}$.

In this hypothetical case, the effect of sampling and sample size is depicted
in figures 1(c)-1(e), where we additionally illustrate individual and averaged
losses converging to the theoretical value. This theoretical value could serve
as a threshold to establish the significance level of a regression problem:
regressions with losses lower than this threshold, with at least a certain
probability, allow us to conclude that there is enough evidence of a linear
relationship in the population. Under the null hypothesis, the expected loss
can be computed on a mesh-grid of uniformly distributed points ($\mathbf{x}$,
$y$) within a standardized hypercube of $P+1$ dimensions centered at the
origin. After some calculus with $P=1$, we find (see appendix)
$\mathcal{R}=\frac{b}{2}+\frac{a^{2}}{6b}$ for $\mathcal{L}_{1}$ and
$\mathcal{R}=\frac{b^{2}+a^{2}}{3}$ for $\mathcal{L}_{2}$, where $a$ and $b$
are the maximum values for $\hat{y}$ and $y$, respectively. These values can
be easily approximated under the null hypothesis, as illustrated in figures
1(c) and 1(e). It is worth mentioning that when using the $\mathcal{L}_{1}$
loss function, widely employed in predictive models for linear regression,
such as K-fold and leave-one-out CV, the results are overly optimistic. This
is because they provide convergent values lower than the theoretical expected
loss under $H_{0}$.

(a) Regression fitting for $\tau=0.9$.

(b) Regression fitting for $\tau=0.1$.

(c) Averaged (empirical) loss under $H_{0}$ ($\mathcal{L}_{1}$).

(d) Individual losses in the volume $\mathcal{V}$ of pairs $\\{\hat{y},y\\}$
for $\mathcal{L}_{1,2}$

(e) Averaged (empirical) loss under $H_{0}$ ($\mathcal{L}_{2}$)..

Figure 1: In this 2D example of regression fitting, we explore effects ranging
from large to very small. Observe the flatness of the linear functions in
subfigure 1(b). In the regression examples of subfigures 1(a) and 1(b), we
present several representations: y vs. x; the result of regression; y vs.
$\hat{y}$; and $\mathcal{L}$ vs. (y, $\hat{y}$).

##  3 A formal test for the error estimation of ML algorithms

In this section, we present a method for testing the null hypothesis
$H_{0}:\upbeta=0$ against the alternative hypothesis $H_{A}:\upbeta\neq 0$,
utilizing a common measure employed by ML researchers, namely the expected
loss $\mathcal{R}$. Under $H_{0}$, with the linear regressor
$f_{N}\in\mathcal{F}$, we anticipate the loss value to be comparable to that
obtained with uncorrelated variables. Conversely, this value is expected to be
lower when the regression coefficients are significant. The test is formulated
as:

$\begin{array}[]{c}H_{0}:\mathcal{R}(f_{N})=E[\mathcal{L}_{1,2}(f_{N},\mathbf{x},\mathbf{y})|\mathbf{X}\perp\mathbf{y}]\\\
H_{A}:\mathcal{R}(f_{N})<E[\mathcal{L}_{1,2}(f_{N},\mathbf{x},\mathbf{y})|\mathbf{X}\perp\mathbf{y}]\end{array}$
(5)

Here, the actual risk $\mathcal{R}(f_{N})$ is obtained by analyzing the worst
case with at least a probability $1-\eta$; i.e., by adding an upper bound
$\Delta(N,\mathcal{F})$ to the empirical risk $\mathcal{R}_{N}$:

$\mathcal{R}(f_{N})=\mathcal{R}_{N}(f_{N})+\Delta(N,\mathcal{F})$ (6)

Note that sub-indices stand for the selected loss, i.e. $\mathcal{L}_{1}$ or
$\mathcal{L}_{2}$. The proposed test, named Statistical Agnostic Regression
(SAR), assesses the significance of the regression by comparing the actual
risk with that obtained under the assumption of orthogonality between
predictors and observations.

The former test is similar to the classical F-test (also known as analysis of
variance), defined as the ratio $F^{*}=MSR/MSE$, with $P$ and $N-P-1$ degrees
of freedom, where $MSE$ is the mean square error and $MSR$ is the mean square
due to regression [Chatterjee06]. The difference lies in the linear model
assumptions, i.e., $E[MSE]=\sigma^{2}$. If these assumptions are fulfilled,
the $F^{*}$ transforms into an $F$ distribution, and the regression problem
becomes a test of hypothesis or the analysis of p-values. Even if the linear
model assumptions are not fulfilled, we can test our ML measures in the same
manner [Reiss15] and compare them with the test provided in equation 9 (see
the experimental section). Moreover, previous approaches to achieve
statistical significance using ML and statistical agnostic theory in
classification tasks [Gorriz23] are naturally extended in this regression
analysis.

##  4 Materials and Methods

###  4.1 Upper Bounding the expected loss under the worst-case scenario

The primary objective of SLT is to establish a structured approach for
tackling the challenge of statistical inference, as highlighted in works such
as [Vapnik82, Haussler92]. A notable accomplishment of SLT is the development
of straightforward and robust confidence intervals that effectively delimit
the actual risk $\mathcal{R}(f)$ [Boucheron13]. Our specific focus lies in the
estimation of risk derived from an empirical quantity, ensuring a probability
of at least $1-\eta$. This estimation is articulated through the concentration
inequalities (CI) presented as:

$\mathcal{R}(f_{N})\leq\mathcal{R}_{N}(f_{N})+\Delta(N,\mathcal{F})$ (7)

Here, $f_{N}$ is carefully determined to prevent overfitting, wherein we
restrict the class of functions by setting
$f_{N}=\underset{f\in\mathcal{F}}{\arg\min}\mathcal{R}_{N}(f)$. The term
$\Delta(N,\mathcal{F})$ serves as an upper bound for the true risk. In the
worst-case scenario, the inequality transforms into an equality. This
deviation or inequality can be understood through various perspectives in
classical probability theory [Vapnik82] offering insights into how closely the
sum of independent random variables (empirical risk) are to their expectations
(actual risk).

###  4.2 A Probably Approximately Correct Bayesian bound

In this work, we leverage a significant advancement in the field rooted in
Probably Approximately Correct (PAC)-Bayesian theory [MacAllester2013].
Specifically, we assess a dropout bound inspired by the recent success of
dropout training in deep neural networks. The bound, as represented in
equation 7, is articulated with respect to the underlying distribution $Q$,
which samples the function $f$ from the set of "rules," denoted as
$\mathcal{F}$.

For any constant $\lambda>1/2$ and a class of linear regressors,
$\upbeta\in\mathbb{R}^{P+1}$, selected according to the distribution $Q$, we
establish that with a probability of at least $1-\eta$ over the sample draw,
the following CI are valid for all distributions $Q$ with dropout rate
$\delta$:

$\begin{array}[]{l}\mathcal{R}(\upbeta)\leq\mathcal{R}_{N}(\upbeta)+\\\
\min_{1\leq i\leq
k}\frac{1}{2\lambda_{i}-1}\left(\mathcal{R}_{N}(\upbeta)+\frac{2\lambda_{i}^{2}\mathcal{L}_{max}}{N}\left(\mathcal{D}(Q,Q_{u}))+\ln\frac{k}{\eta}\right)\right)\end{array}$
(8)

Here, $\mathcal{D}(Q,Q_{u})=\frac{1-\delta}{2}||\boldsymbol{\upbeta}||^{2}$
represents the Kullback-Leibler divergence from $Q$ to the uniform
distribution $Q_{u}=\mathcal{N}(0,1)^{n}$, $\lambda\in(1/2,10)$ can assume $k$
different values and $\mathcal{L}_{max}$ is an outlier threshold
[MacAllester2013]. Other approaches to formulate upper bounds can be
considered based on more general assumptions [Gorriz19, Vapnik82].

###  4.3 ANOVA and the analysis of residuals in ML

In classical regression [Chatterjee06] results are tested for significance by
assessing the analysis of variance of the residuals $r=y_{i}-\hat{y}_{i}$. As
aforementioned, the $F^{*}$-statistic is compared to an F-distribution with
$P=1$ numerator and $N-2$ denominator degrees of freedom and the p-value, or
the probability that we get this statistic as large as we did under the null-
hypothesis, is determined. This is also known as the formal test for the slope
parameter $\upbeta_{1}$ and, if the linear model assumptions are fulfilled,
can be extended to evaluate the residuals obtained by ML methods.

Another possibility for drawing conclusions about the population (and not only
the particular observed sample) is by the use of confidence intervals.
Confidence intervals and hypothesis tests are two different ways of learning
about the values of population parameters ($\beta_{0}$ and $\upbeta_{1}$).
Both approximations can be conducted on ML residuals as well to demonstrate
the reliability of the results, however this is rarely found in the ML
literature and there is a trend to exclusively show different performance
scores based on empirical measures; e.g. accuracy in learning curves.

Similar to the permutation analysis proposed for group comparisons
[Bullmore99], we can test the power of the proposed SAR method by comparing a
set of Monte Carlo simulations with the expected loss of the null
distribution, following the ideas presented in previous sections. We estimate
the $1-\beta$ by counting the number of times the corrected empirical loss was
less than the averaged power of the observation.

###  4.4 Heteroscedascity and the Breusch-Pagan test

When the residuals are heterogeneously distributed along the explanatory
variable (predictor), we encounter the issue of heteroscedasticity in the
data. The Breusch-Pagan (BP) test [Breusch79] checks the null hypothesis
(homoscedasticity):

$\begin{array}[]{c}H_{0}:\mathcal{R}(f_{N})\text{ is constant}\\\ \end{array}$
(9)

by re-fitting the explanatory (independent) variable to the squared residuals
as observation variables. A simple way to assess this condition is by using
the test statistic $T=N\cdot R^{2}$, where $R$ is the coefficient of
determination of the auxiliary fitting, and follows a $\chi^{2}$ distribution
with $P$ degrees of freedom under $H_{0}$ [Koenker81]. If the p-value
associated with it is less than the level of significance, we reject $H_{0}$.
The homoscedasticity condition is assumed in the linear regression model, and
when it fails, further analyses based on the F-test are not only imprecise
(biased) but also invalidated from a statistical point of view. ML residuals
could be tested in the BP test as well, to assess the ability of data-driven
models to detect heteroscedasticity in regression problems with limited sample
sizes.

###  4.5 Synthetic, realistic & real datasets

We evaluated linear SVR models using the validation techniques: K-fold, LOO,
Resubstitution, and SAR, along with the OLS method, on various synthetic and
realistic datasets. First, we assessed the 2D problem with Gaussian-
distributed variables ranging from perfectly uncorrelated to strongly
correlated signals. To control the degree of correlation, we drew samples
$\mathbf{Z}$ from a 2D normal distribution, i.e., $\mathcal{N}(0,1)$.
Subsequently, we computed a random linear transformation matrix $\mathbf{T}$
and modified its singular value decomposition (SVD) diagonal matrix using a
scaling transform similar to the earlier simulated dataset (1). Finally, we
applied this transformation matrix to the original data, obtaining
$\hat{\mathbf{Z}}=\mathbf{Z}\mathbf{T}$ that generated observations and
predictors (see figure 2). We fitted our models on this data with two losses
$\mathcal{L}_{1}$ and $\mathcal{L}_{2}$, several samples sizes
$N=10,20,\ldots,300$ and tested for significance using the F-test and the SAR
test. We drew up to $R=100$ realizations (sampling) from the ideal
distribution and then averaged the regression results to study the effect of
sample size [Friston13]. Whenever CV methods split the dataset into folds, we
analyzed the variability of the performance measures accordingly. This
represents a realistic situation with a fixed realization or dataset instead
of having an ensemble of $R$ realizations.

Figure 2: Data transformed by rotation and scaling with a non-diagonal
covariance matrix, assuming a Gaussian distribution.

To enhance the realism of our data, we employed the procedure used in
[Gorriz19] to segment the data into $N_{c}$ clusters. By selectively removing
some of these clusters, we introduced non-Gaussian characteristics. Under this
experimental setting, the F-test and, in particular, the OLS method are no
longer optimal. This allows us to truly test the robustness of the SAR test
developed in Section 3. We conducted tests for up to four experimental
settings by randomly removing a set of $N_{c}-N_{g}$ clusters, with $N_{g}=3$.
In figure 3 we plotted the remaining clusters with different colors to
illustrate the complete data generation.

Figure 3: Gaussian data transformed by rotation and scaling, along with
cluster pruning. Colors simply indicate the applied transformations to the
data and identify the clusters that were removed.

Finally, we assessed how the BP test worked on the residuals obtained from the
set of methods previously described. Specifically, we evaluated the ability of
the resulting models to detect heteroscedasticity versus sample size. For this
purpose, we employed the dummy dataset provided in the XLSTAT software
[Soft19] designed to compare a homoscedastic model to another with strong
heteroscedasticity. In particular, we generated a dependent variable (“Size”)
based on an independent variable (“Age”), where the residuals are defined as
the product of the independent variable by the random normal error (figure
4(a)). In this case, the residuals were obviously correlated with age, and the
problem (figure 4) was apparently linear, but heteroscedastic.

(a) Datasets.

(b) Tests for linearity.

(c) Qqplot revealing non gaussianity.

(d) Explanatory variable vs. residuals.

Figure 4: Dataset with heteroscedascity and increasing sample size in figure
4(a) and SAR and F tests on linearity in figure 4(b). Assumptions needed to
perform the F-test are not fulfilled as shown in the qqplots in figure 4(c)
and residuals vs. explanatory variable plot in figure 4(c).

####  4.5.1 Cancer dataset in multiple dimensions

We tested our multivariate methods on real data set with data downloaded from
US National Cancer Institute and the US Census American Community Survey to
explore the linear relationship between between cancer mortality rate and
several predictors (socio-economic status) in US counties. See
https://data.world/nrippner/cancer-linear-regression-model-
tutorial\\}forafulldescriptionofthedata(seefigure\ref{fig:6a}).\par%****␣template.tex␣Line␣250␣****FollowingtheOLSanalysis(figure\ref{fig:cancer})ofthisdatasetandeliminatingmulticollinearitybyapplyingvarianceinflationfactors,weemployeduptosixindependentpredictorsincluding:Bothmaleandfemalereportedbelowpovertylinepercapita(All\\_Poverty\\_PC),medianincomeofallethnicities(Med\\_Income),malesandfemaleswithandwithouthealthinsurancepercapita(All\\_With$/$out\\_PC),lungcancerincidencerate(Incidence\\_Rate),andthepopulationestimatein2015(Pop\\_estimate).Atotalof$N=2809$samplesizedatasetwascreatedandrandomlydownsampledforrandomizationanalysis.\par\begin{figure*}\centering\begin{subfigure}{\textwidth}\includegraphics[width=\textwidth]{datareal}\caption{Cancerdataset.}\label{fig:6a}\end{subfigure}\hfill%\begin{subfigure}{0.49\textwidth}%\includegraphics[width=\textwidth]{BP2}%\caption{Testsforlinearity.}%\label{fig:5b}%\end{subfigure}%\hfill\begin{subfigure}{0.49\textwidth}\includegraphics[width=\textwidth]{qqplot}\caption{Qqplotrevealingnongaussianity.}\label{fig:6b}\end{subfigure}\hfill\begin{subfigure}{0.49\textwidth}\includegraphics[width=\textwidth]{residuals}\caption{Explanatoryvariablevs.residuals.}%****␣template.tex␣Line␣275␣****\label{fig:6c}\end{subfigure}\par\caption{CancerDatasetwithmultiplepredictors(6)infigure\ref{fig:6a}.Thefirstcolumnrepresentsthemortalityrateastheobservablevariable.Infigure\ref{fig:6b}werepresenttheqqplotfollowingtheanalysisgiveninthereferenceprovidedinthetext.AssumptionsneededtoperformtheF-
testarenotfulfilledasshownintheqqplotsinfigure\ref{fig:6b}andresidualsvs.explanatoryvariableplotinfigure\ref{fig:6c}.}\label{fig:cancer}\end{figure*}\par%\section{Results}\label{sec:experiments}\par\subsection{Gaussiandata:OLSisthegoldstandard}\parInthissection,linearlytransformedGaussiandataisconsidered,wheretheOLSmethodisthegoldstandard.Weaveraged$R=100$realizationsandusedthesameloss($\mathcal{L}_2$)withstandardMLalgorithms,whichperformedquitewell(SARwasconservativebutconvergedtoOLSasthesamplesizeincreased).However,therewashighvariabilityintheexpectedlossfromfoldswhenthenumberofsampleswaslimited($N<50$).Itisworthnotingthatwhenusingthe$\mathcal{L}_1$loss,thetypicallossinMLapplications,thecorrectionappliedtotheresubstitutionerrormadeSARconvergetoOLS,unliketheotherMLvalidationmethods.\begin{figure*}\centering\includegraphics[width=0.49\textwidth]{EL_rmse}\includegraphics[width=0.49\textwidth]{EL_eps}\caption{Expectedlosses($\mathcal{L}_2$ontheleftand$\mathcal{L}_1$ontheright)andanalysisofvariancefromfolds.Notethatthevariancefromfoldsisthatexperiencedinrealscenarioswithasinglesamplerealization.}\label{fig:datanG}%\end{figure*}ThemethodsarenowtestedwiththeformalF-
testfortheslopetoassesstheirabilitytodetectaslightlinearrelationshipwithacorrelationlevelequalto0.1.TheresultsrevealoptimisticbehavioroftheMLmethods,exceptfortheSARmethodwhichstruckatrade-
offbetweenOLSandMLmethods.ItisessentialtorecallthatOLSisthegoldstandardinthiscase,andtheconservativebehavioroftheSARallowsustoobtainoutputssimilartoML.InFigure\ref{fig:Ftest},wedisplaythep-
valueinthischallengingcasewherethenullhypothesisisnotrejected.Atthebottomofthefigure,weprovidetheSARtestbasedonthedetectionthresholdthatmimicsclassicaltestsfortheslopebasedonasignificancelevel.\begin{figure}\centering\includegraphics[width=0.49\textwidth]{Fstat}\includegraphics[width=0.49\textwidth]{tests}%****␣template.tex␣Line␣300␣****\caption{TheF-
testfortheslope$\upbeta_1$inthecaseofacorrelationlevelequalto0.1.WeutilizedtheresidualsderivedfromeachmethodandaveragedtheF-
statisticanditsvariabilityover$R=100$repetitions.ItisnoteworthythatSARprovidesatrade-
offbetweenMLmethodsandOLS,exhibitinglessvariabilityintherepetitions.TheSARtestisobtainedwithaprobabilityofatleast$\eta=0.5$.}\label{fig:Ftest}%\end{figure}TheoptimisticresultsobtainedbyCVmethodswithalimitedsamplesize,whichcouldindicatetheirbetterabilityforlinearitydetection,areindeedaconsequenceofpoorcontroloffalsepositives(FPs).Todemonstratethisissue,wedesignedaputativetaskconsistingofaregressionproblemwithnocorrelationatall,i.e.,acorrelationlevelequalto0.Inthiscase,weanalyzetherealcasewithalmostnorepetitions($R=2$),wherethevariabilityoferrorsinfoldsaffectstheevaluationoftheF-
statistic.Inthisscenario,weobserveanon-smoothbehavioroftheF-
statisticcurves,butonlyMLCVlinescrossbelowthesignificancelevel(e.g.,$N=100$),unliketheSARtestortheF-
testfortheOLSestimation,whicharealwaysabovethedetectionthreshold,controllingtherateofFPs(seeFigure\ref{fig:Ftest2}).Moreover,toillustratethisissue,weperformedapoweranalysisinFigure\ref{fig:power1}withawidevarietyofcorrelationlevels,showingthatCVmethodsareinflatingtherateofFPs(detectionratesabovetheOLSmethodforlowcorrelationlevels).\begin{figure}\centering\includegraphics[width=0.49\textwidth]{Fstat2}\includegraphics[width=0.49\textwidth]{tests2}\caption{TheF-
testfortheslope$\upbeta_1$inthecaseofacorrelationlevelequalto$0.1$.WeutilizedtheresidualsderivedfromeachmethodandaveragedtheF-
statisticanditsvariabilityover$R=2$repetitions.}\label{fig:Ftest2}%\end{figure}\begin{figure}\centering\includegraphics[width=0.49\textwidth]{power1}\caption{PoweranalysisderivedfromtheF-
testinGaussiandata.TheselectedlossforMLmethodsis$\mathcal{L}_1$.ObservehowtheresultsobtainedbytheSARfeaturesarebyfarthemostsimilartotheonesobtainedbytheOLS.TheSARtestbasedonexpectedlossismoreconservativethantheothermethods(correlationlevelsbelow$0.37$arenotconsideredsignificantbasedontheworst-
caseanalysiswithaprobabilityofatleast$\eta=0.5$).}\label{fig:power1}%\end{figure}\par\subsection{Non-
Gaussiandata:largeCross-ValidationresidualsandFalsePositives}Fornon-
Gaussiandistributeddata,theOLSmethodisnotthebestestimatorbutbiased,andananalysisderivedfromtheGaussianityassumption,i.e.theF-
test,canleadtounreliableresults.HerewefindthetrueutilityoftheSARtest:tovalidatethesignificanceofapriorregressionapproach.Weanalyzedtwocases,asshowninFigure\ref{fig:datanG},withcorrelationlevels$0.37$and$0.1$.Bothexhibitedthesimilarp-
valuesforthestandardMLapproachesandOLS.However,thisisnotalwaysgoodnews,especiallyfor$R=100$asshowninFigure\ref{fig:Ftest3}.\parFromapproximately$N=100$,allmethodsrejectedthenullhypothesisandonemightbetemptedtoconcludethatthereissufficientevidenceforalinearrelationshipinthedata.UnliketheSARtest,thevaluesoftheexpectedlossconvergetothedetectionthresholdbutalwaysremainaboveit.Inrealisticscenarioswith$R=2$,allmethodsprovidedap-
valueabove$0.05$.So,forexample,resultsfromasinglelaboratorywithonlytwosamplescannotbeextrapolatedtotherestoftherealizations.Inotherwords,dependingontherepetition,wemayhaveasignificantresult,orthenullhypothesiscannotberejected.Asshowninthesamefigure,theSARtest,beingconservative,isrobustinthisscenariocomparedtotheF-
test.\begin{figure}\centering\includegraphics[width=0.49\textwidth]{tests3}%****␣template.tex␣Line␣325␣****\includegraphics[width=0.49\textwidth]{tests4}\caption{TheF-
testfortheslope$\upbeta_1$inthecaseofcorrelationlevelsequalto$0.37$and$0.1$withnon-
Gaussiandata.WeemployedtheresidualsderivedfromeachmethodandaveragedtheF-
statisticanditsvariabilityin$R=100$and$2$repetitionsforcomparisonpurposes.}\label{fig:Ftest3}%\end{figure}\parAnotherexampleisdepictedinFigure\ref{fig:Ftest4},wheretheSARtestvalidatesthesignificanceofthetestsfromapproximately$N=150$.NotethattheSARtestisrobustagainstthenumberofrepetitionsusedfromitscomputation,andfor$R=100$withsamplesizesabove$N=50$thesignificantregionsonaveragearenotfulfilledwithonly$R=2$repetitionsleadingtopoorreplication.Again,weperformedapoweranalysis,includingaputativetaskwithacorrelationlevelequaltozerotoassesshowthemethodscontroltherateoffalsepositives(FPs)\cite{Rosenblatt16}.Inthiscase,therateofFPsisapproximatelyequaltothelevelofsignificance$\alpha$.Figure\ref{fig:power2}showsthattheMLmethodsareinflatingFPs(greensolidline)abovethelevelofsignificance,i.e.,$1-\beta>\alpha$.RecallthattheresultsofOLSunderconditionsofnon-
Gaussianityshouldbesubjecttoscrutiny(exceptforthecase$\tau=0$).\par\begin{figure}\centering\includegraphics[width=0.49\textwidth]{tests5}\includegraphics[width=0.49\textwidth]{tests6}\caption{TheF-
testfortheslope$\upbeta_1$inthecaseofcorrelationlevelsequalto$0.37$and$0.1$withnon-
Gaussiandata.WeemployedtheresidualsderivedfromeachmethodandaveragedtheF-
statisticanditsvariabilityin$R=100$and$2$repetitionsforcomparisonpurposes.}\label{fig:Ftest4}%\end{figure}\par\begin{figure}\centering\includegraphics[width=0.49\textwidth]{power2}\caption{PoweranalysisderivedfromtheF-
testinnon-
Gaussiandata(exceptforthegreensolidlinethatrepresentsuncorrelateddata).TheselectedlossforMLmethodsis$\mathcal{L}_2$.ObservehowtheresultsobtainedbytheSARfeaturesarebyfarthemostsimilartothoseobtainedbytheOLSforlowcorrelationlevels.However,forhighlevels,theperformanceofSARisatrade-
offbetweenMLmethodsandOLS.TheSARtestbasedonexpectedlossismoreconservativethantheothermethods(correlationlevelsbelow$0.37$arenotconsideredsignificantbasedontheworst-
caseanalysiswithaprobabilityofatleast$\eta=0.5$).}\label{fig:power2}%\end{figure}\par\subsection{Heteroscedascitydata:superiordetectionusingMLtechniques}\parInthissection,wetestedtheabilityofthemethodstodetectheteroscedasticityusingtheBPtestdescribedinSection\ref{sec:matmet}.Theproblemanalyzedhereisclearlylinear,butthedistributionofresidualsisnothomogeneous(seeFigure\ref{fig:BP1}).Asimple’qqplot’analysisofthedatarevealsthestrongheteroscedasticityofthedata,\parWeevaluatedtheBPtestontheresidualsobtainedwithallmethodsandcomparedtheaveragedandstandarddeviation(std)P-valuewithincreasingsamplesizein$R=100$repetitions.TheBPtestusingtheSARresidualsprovidesafasterdetectionofheteroscedasticityonaverageandalowerstandarddeviationoftheP-
value,asshowninFigure\ref{fig:BP2}.\begin{figure}\centering\includegraphics[width=0.49\textwidth]{BP5}\caption{TheBPtestontheresidualsobtainedbyMLmethodsandtheOLS.}\label{fig:BP2}%\end{figure}\par\subsection{MultivariateLinearRegression:SARfeatureseffectivelydetectlinearity.}\parWerepeatedthetestforlinearityandthepoweranalysiswithincreasingdimensions,fromonetosixpredictors,usingtheCancerdataset.Weobservedtwoprimaryeffects.First,thevarianceofCVmethodsisextremelylargeforsmallsamplesizes($N=10-50$),especiallywiththe$\mathcal{L}_2$loss,althoughitissignificantfortherestofthesimulations.Thiseffectdependsonthepredictoraddedintheregressionanalysis.Notably,thelastpredictorthatwasadded,’Pop\\_estimate’(predictor6),increasesthevarianceoftheestimatorsresultinginnon-
reliableestimateswithintrainingfolds,mainlywithsmallsamplesizes(Figure\ref{fig:datarealloss}).ThereasonforthisanomalousoperationwaspreviouslyshowninFigure\ref{fig:6c}inlightblue.Here,wereadilyseeintheSARtest(bottomofthefigure)thatevidenceforalinearrelationshipinthedata,undernon-
idealconditions,isachievedwithonly10samples(theriskisalwayslessthanthethreshold).\par\begin{figure*}\centering\includegraphics[width=0.49\textwidth]{reallossrmse}\includegraphics[width=0.49\textwidth]{reallossepsins}\caption{Expectedlosses($\mathcal{L}_2$ontheleftand$\mathcal{L}_1$ontheright)andanalysisofvariancefromfolds.Againthevarianceusingthe$\mathcal{L}_2$lossisgreaterthanthe$\mathcal{L}_1$lossforMLmethodsusingCV.However,withrealdatathepresenceoflargeoutliersinfoldsresultsininadmissiblelossestimations.}\label{fig:datarealloss}%\end{figure*}\parThesecondeffectworthcommentingonisthedetectionabilityoftheSARmethodwhenthenumberofpredictorsincreases,asshowninFigure\ref{fig:7b}.Moreover,theSARmethodconvergestotheOLSmethodwiththetwotestedlosses,unliketheCV-
basedMLmethods.ThissuspiciousfeaturefoundfortheSARtest(andforanyotherMLmethod)canbeexplainedbyoverfittingsinceitgoesagainstitscommonbehavior.Todiscardthispossibility,weincreasedthedimensionsorthenumberofpredictorsusedinthesimulation,asshowninFigure\ref{fig:power2},representingtheputativetaskwherethereisnocorrelation.Inthistask,weadditionallyemployedageneralupperboundbasedontheassumptionofsamplesingeneralposition(i.g.p.),asproposedin\cite{Gorriz19}.TheanalysisforthetwoboundsisgiveninFigures\ref{fig:7c}and\ref{fig:7d}.Asshowninthisanalysis,whenweincreasethenumberofdimensions(onlyinthissection),theboundproposedin\cite{MacAllester2013}isoverlyoptimisticcomparedtothatproposedin\cite{Gorriz19}when$N<100$.Inthiscase,PACapproachesshouldbereformulated,andanalyticalboundssuchasthoseproposedin\cite{Vapnik82,Gorriz19}arevalidsolutionstoformulatetheSARtestastheyeffectivelycontroltheFPrate.\par\begin{figure*}\centering%****␣template.tex␣Line␣375␣****\begin{subfigure}{0.49\textwidth}\includegraphics[width=\textwidth]{F_SAR_test_real}\caption{TestingthelinearrelationshipinCancerdatasetwith6predictorsusing$\mathcal{L}_1$.}\label{fig:7a}\end{subfigure}\hfill\begin{subfigure}{0.49\textwidth}\includegraphics[width=\textwidth]{powerreal}\caption{Poweranalysisversusthe$\\#$predictorsandsamplesize.NotethatresultsareobtainedwiththemostcompetitivelossforMLmethods($\mathcal{L}_1$).}\label{fig:7b}\end{subfigure}\hfill\begin{subfigure}{0.49\textwidth}\includegraphics[width=\textwidth]{SAR}\caption{TestforlinearityintheputativetaskofuncorrelatedsamplesusingPACbound.}\label{fig:7c}\end{subfigure}\hfill\begin{subfigure}{0.49\textwidth}\includegraphics[width=\textwidth]{SARG}\caption{Testforlinearityintheputativetaskofuncorrelatedsamplesusingthei.g.p.bound.}\label{fig:7d}\end{subfigure}\par\caption{PerformancewithincreasingnumberofpredictorsinCancerDataset.}%****␣template.tex␣Line␣400␣****\label{fig:cancer2}\end{figure*}\par%\section{Discussion}%\parInitially,werecapitulatedtheconventionalapplicationofOLSinlinearregression,emphasizingitswidespreadadoptionduetoitsoptimalstatisticalproperties.ThetransitiontoMLintroducesregularizationhighlightingtheshiftinfocustowardsminimizingbothempiricalriskandmodelcomplexity.Critically,weunderscoredthecommondeficiencyinMLapproaches,particularlyinregression,whererigoroussignificanceanalysisisoftenneglected,relyinginsteadonpermutationtestingandCV.TheproposedSARmethodseekstobridgethisgapbyformulatingastatisticaltestgroundedinstatisticallearningtheory.Thisinvolvesestablishinganupperboundontheexpectedlosscomparingittothenullhypothesisandrejectingitifthecorrectedriskislower,signifyingevidenceofalinearrelationship.\parThebackgroundofSupportVectorregressionwasthenexploredpresentingthe$\epsilon$-insensitivelossfunctionandtheregularizedriskfunctional.Thediscussionextendstoexploringtheflatnesspropertyandempiricalloss,highlightingthedifferencesinlossfunctionsemployedinMLmodels.Finally,thenovelSARmethodwasintroducedasaformaltestforassessingtheerrorestimationofMLalgorithms,drawingparallelswithclassicalstatisticaltestsliketheF-
test.ThesectionemphasizedthesignificanceofSARinovercominglimitationsinMLapproachesandextendedtheapplicationtoregressionanalysis.\parTheresultssectionpresentedacomprehensiveexaminationofexperimentaloutcomesfocusingontheperformanceoftheproposedSARmethodincomparisontostandardMLtechniquesacrossdiversescenarios.ForGaussiandistributeddata,whereOLSservesasthegoldstandard,theSARmethodexhibitedconservatismconvergingtoOLSwithincreasingsamplesizes\cite{Noble20}.TestingwiththeF-
testfortheslopeatacorrelationlevelof$0.1$revealedtheoptimisticbehaviorofMLmethods,exceptforSAR,whichstruckatrade-
offwithOLS\cite{Phipson2010}.Analysisofnon-
GaussiandistributeddataemphasizedtheutilityofSARinvalidatingsignificance(reproducibility),especiallyinscenarioswithlimitedexperimentalrepetitions($R=2$)\cite{NAS2019}.ThedetectionofheteroscedasticityusingtheBPtestdemonstratedSAR’sfasterdetectionandlowerstandarddeviationthep-
values.Thisthoroughanalysis,illustratedthroughvariousfigures,underscoresSAR’sversatilityandrobustnessinassessingerrorvariability,validatingsignificance,anddetectingheteroscedasticityacrossdifferentdatascenarios(samplesizes).Tosumup,fromthesetofexperimentsundertakenhere,wedemonstratedthattheproposedSARiseffectiveforassessingthevariabilityoferrorvaluesobtainedbyCVinMLmethods.Moreover,theevidencederivedfromtheseexperimentssuggeststhatSARisarobustapproachtofunctionalregression.\parAsmentionedintheintroduction,thereisincreasingconcernabouttheuseofK-
foldcross-
validationasabaselinemethodformodelselection.AIisproposingmorecomplexalgorithmicarchitecturesthatarevalidatedwiththeeasytoimplementCVmethod,potentiallyleadingtosignificantvariabilityifsamplesizesarelimitedorthedataheterogeneityislarge.Thisraisesthequestionofwhetherthegianthasfeetofclay.Whenthedataallowsustoderiveastableinducerbysplitting(makingaguessfromsomepartsofit),itcanbeaveryrobustmethodtoderivemodelsbeyondthecurrentrealization.However,thisassumptionisnotalwaysmet.Byperformingthisprocedure,wecanexpectthatwhatwelearninthetrainingfoldsisnotnecessarilyusefulforpredictionwiththeremainingsamples.ThisisanalogoustothejudgmentofSolomon;weaskthetwowomentodeterminewhothebabybelongsto,andtheinformationprovidedbyeachpartyshouldbeuseful.Incasethetwowomenreplythesame,Solomonwouldhavetosadlysplitthebabyintwo.Inthebiblicalstory,thereactionofeachofthetwowomentothesuggestionofsplittingthebabywasnotthesame,andSolomongavethebabytotherealmother.ThisisakintowhatwehavedonewithdatausingSLT,althoughwecorrectthisresubstitution(splitintwo)byconsideringtheinformationcontainedinthedata(reactionfromdata).Inourcase,wewouldgivethetwopartstothemoreaffectedwoman.\parThetheoreticalfindingsarealsoapplicabletorealdatainvariousfieldsofresearch.Indeed,themodelsdevisedinthispaper,togetherwiththesimulationofrealisticdatasets,createsuitableexemplarsforcharacterizingperformanceinneuroimagingapplications\cite{Gorgen18,Zhang14,Gorriz2021b,Gorriz2022,Wang07,Wang09}.Asimplecomparisoninclassificationtasksbetweenthesimulateddatasetsandthoseofrealsituationsrevealsthesimilarityintheresultsobtained\cite{Gorriz19}.Nevertheless,thescatterplotsanddatadistributionsprojectedonthedimensionswereclearexamplesdemonstratingthattheconditionstoprovidestableinducers\cite{Kohavi95}werenotmet.Therefore,explorationofalternativevalidationmethodsisapriority.\parInsummary,echoingpracticesfromourpreviouswork,weunderscoretheimportanceofemphasizingnegativeresultstoenhancescientificunderstanding\footnote{Forfurtherinsights,refertothecolumninNatureaddressingthisissue:https://www.nature.com/articles/d41586-019-02960-3}.Whileresearchpapersoftenprioritizepositiveoutcomes,weinfrequentlyassessouralgorithmsinhypotheticaltaskscenarioswherenodiscernibleeffectisexpected.Thisisevidentinourexperimentswithuncorrelateddata.Suchanalysesplayacrucialroleservingtoapproximatethenull-
distributionofthetest-
statisticinpermutationstudies,suchasevaluatingperformanceoraccuracyinaclassificationtaskusingMLtechniques.Forinstance,inthepermutationanalysisofclassificationtasks,theperformancederivedfrompaireddataandlabelsisjuxtaposedwiththatobtainedbyrandomlypermutinggrouplabelsnumeroustimeswiththeanticipateddistributioncenteredaround$50\%$.Iftheperformancedistributionexhibitsnon-
symmetryaroundrandomchanceandbias,conclusionsdrawnfromthetestdatamaybecompromised.Thissuggestsadisparityindatadistributionbetweengroupsunderthenullhypothesis,violatingthei.i.d.assumption.Consequently,theestimationofp-
valuesmightleadtoinaccurateconclusionsatthefamily-
wiselevel\cite{Phipson2010}.\par\par%\section{Conclusions}\parInthispaper,wepresentamethodforvalidatingregressionmodelsinthefieldofmachinelearninganditsapplications.ThemethodisrelatedtotheF-
testofclassicalhypothesistestingforestablishingsignificanceinlinearmodels.WedemonstratethatstandardMLmethodsformodelvalidationtendtooverinflateFPs,thusrequiringtheseapproachestoensuregoodreplicationandextrapolationofresultsinlimitedsamplesizes.WhenSARfeaturesareincorporatedintoclassicalstatisticalframeworks,theyprovideatrade-
offbetweenOLSandMLparadigmswithexcellentcontrolofFPs(aroundthelevelofsignificance).Althoughthesepipelinescouldfacecriticismwhendealingwithnon-
Gaussiandata,wealsoconductedtheformaltestpresentedbasedonCIs.Thistestdiscardedregressionproblemswithlowcorrelationlevelsundertheworst-
casescenarioandprovidedstatisticalsignificancefortherestofthecases.Theuseofthisformaltestisintendedtobecombinedwithclassicalhypothesistesting,allowingustoconfirmtheanalysiswithtechniquesthatarenotformallyvalid,ortoinformtheresearcherthatthereisinsufficientevidencefromtheperspectiveofSLTtoestablishsuchalinearrelationship.Thisconstitutestheformaldefinitionofhypothesistesting:conclusionsaboutthedatacanonlybedrawnwhenwerejectthenullhypothesis;otherwise,cautionmustbeexercisedwiththefindings.\par\par%%StandardCVmethodsincombinationwithMLwereevaluatedtoascertainwhetherstatisticalinferencesmadebydata-
drivenapproachesaresufficientlyconsistent.Theyarefrequentlyclaimedtooutperformconventionalstatisticalapproachessuchashypothesistesting.However,thisimprovementisbasedonmeasures(e.g.accuracy)derivedfromMLclassificationtasksthatdonothaveaparametricdescriptionanddependontheexperimentalsetup;i.e.derivedfromclassification/predictionfoldsofthedatasettoestablishconfidenceintervals.%%Asshowninthispaper,smallsample-
sizedatasetsandlearningfromheterogeneousdatasourcesstronglyinfluencestheirperformanceandresultsinpoorreplication.AnovelstatisticaltestbasedonK-
foldCVandtheUpperBoundoftheactualerror(K-foldCUBV)wasproposedtotackletheuncertainpredictionsofMLandCVmethods.Theanalysisofthe\emph{worstcase}obtainedbya(PAC)-BayesianupperboundforlinearclassifiersincombinationwiththeK-
foldCVestimationisarobustcriteriontodetecteffectsvalidatingaccuracyvaluesobtainedfromMLmodelsandavoidingfalsepositives,complementingtheregularK-
foldmethodforCV.\par\section*{Acknowledgments}ThisresearchispartofthePID2022-137451OB-I00project,fundedbytheCIN/AEI/10.13039/501100011033andbyFSE+.\par\begin{thebibliography}{00}\par\par\bibitem{Zelen62}MarvinZelen.LinearEstimationandRelatedTopics.inSurveyofNumericalAnalysiseditedbyJohnTodd,McGraw-
HillBookCo.Inc.,NewYork,pp.558-577.1962)\par\bibitem{Grohs22}PGrohs,etal.MathematicalAspectsofDeepLearning.CambridgeUniversityPress.ISBN9781009025096.https://doi.org/10.1017/9781009025096.\par\bibitem{Snee77}RonaldD.Snee.ValidationofRegressionModels:MethodsandExamples.Technometrics.Vol.19,No.4(Nov.,1977),pp.415-428\par\bibitem{Kohavi95}R.Kohavi.Astudyofcross-
validationandbootstrapforaccuracyestimationandmodelselection.InternationalJointConferenceonArtificialIntelligence(IJCAI),pp1
--
7,1995.\par\bibitem{Kleijnen96}%****␣template.tex␣Line␣450␣****J.P.C.Kleijnenetal.Validationoftrace-
drivensimulationmodels:regressionanalysisrevisited.ProceedingsWinterSimulationConference..352-359.0-7803-3383-7.Dec1996.\par\bibitem{Miller91}MichaelE.Milleretal..Validationtechniquesforlogisticregressionmodels.StatisticsinMedicine.Volume10,Issue81213-1226.August1991.\par\bibitem{Oredein11}A.IOredeinet.al.OnValidatingRegressionModelswithBootstrapsandDataSplittingTechniques.GlobalJournalofScienceFrontierResearchVolume11Issue6Version1.0September2011\par\bibitem{Lecun15}Y.LeCunetal.Deeplearning.Nature521,436–444(2015).\par\bibitem{Gorriz2020}J.M.Gorriz,etal.Artificialintelligencewithintheinterplaybetweennaturalandartificialcomputation:Advancesindatascience,trendsandapplications.NeurocomputingVolume410,14October237-2702020.\par\bibitem{Varoquaux18}G.Varoquaux.Cross-
validationfailure:Smallsamplesizesleadtolargeerrorbars.NeuroImage180(2018)68-77.\par\bibitem{Eklund16}A.Eklund,etal.Clusterfailure:InflatedfalsepositivesforfMRI.ProceedingsoftheNationalAcademyofSciencesJul2016,113(28)7900-7905.\par\bibitem{Jollans2019}JollansL,etal.Quantifyingperformanceofmachinelearningmethodsforneuroimagingdata.Neuroimage.2019Oct1;199:351-365.\par\bibitem{Gorriz18}J.M.Górriz,etal.AMachineLearningApproachtoRevealtheNeuroPhenotypesofAutisms.Internationaljournalofneuralsystems,1850058.2019.\par\bibitem{Braga04}Braga-
NetoUM,DoughertyER.Iscross-validationvalidforsmall-
samplemicroarrayclassification?Bioinformatics.2004Feb12;20(3):374-80.doi:10.1093/bioinformatics/btg419.PMID:14960464\par\bibitem{vandermaaten08}L.vanderMaatenetal.VisualizingDatausingt-
SNE.JournalofMachineLearningResearch2008vol9,num86,2579--
2605.\par\bibitem{Jimenez22}C.Jimenez-Mesaetal.Anon-
parametricstatisticalinferenceframeworkforDeepLearningincurrentneuroimaging.InformationFusionVolume91,March2023,Pages598-611.\par\bibitem{Gorriz2021}J.M.Gorriz,etal.StatisticalAgnosticMapping:Aframeworkinneuroimagingbasedonconcentrationinequalities.InformationFusionVolume66,February2021,Pages198-212\par\bibitem{Viering23}TomVieringetal.TheShapeofLearningCurves:AReview.IEEETransactionsonPatternAnalysisandMachineIntelligenceJune2023,pp.7799-7819,vol.45DOIBookmark:10.1109/TPAMI.2022.3220744\par\bibitem{Vapnik95}V.Vapnik.ThenatureofstatisticallearningtheorySpringer-
VerlagNewYork,Inc.,(1995)\par\bibitem{Scholkopf95}BernhardScholkopfet.al.LearningwithKernels:SupportVectorMachines,Regularization,Optimization,andBeyond.MITPressISBN:978-0-262-19475-4.December2001\par\bibitem{Burges98}%****␣template.tex␣Line␣500␣****C.J.CBurges.AtutorialonsupportvectormachinesforpatternrecognitionDataMiningandKnowledgeDiscovery,2(2)(1998),pp.121-167\par\bibitem{Huber64}Huber,P.J.(1964).Robustestimationofalocationparameter.Ann.Math.Statist.,35,73–101\par\bibitem{Chatterjee06}Chatterjee,S.andHadi,A.S.RegressionAnalysisbyExample.4thEdition,JohnWiley\&Sons,Hoboken.https://doi.org/10.1002/0470055464.(2006)\par\bibitem{Bullmore99}ETBullmoreetal.Global,voxel,andclustertests,bytheoryandpermutation,foradifferencebetweentwogroupsofstructuralMRimagesofthebrainIEEETransMedImaging(1999)Jan;18(1):32-42.\par\par\bibitem{Reiss15}P.T.Reiss,etal.Cross-
validationandhypothesistestinginneuroimaging:anireniccommentontheexchangebetweenFristonandLindquistetal.Neuroimage.2015August1;116:248-254\par\bibitem{Gorriz23}JuanMGorrizetal.IsK-
foldcrossvalidationthebestmodelselectionmethodforMachineLearning?arXiv:2401.16407\par\bibitem{Vapnik82}V.Vapnik.EstimationdependenciesbasedonEmpiricalData.Springer-
Verlach.1982ISBN0-387-90733-5\par\bibitem{Gorriz19}JMGórriz,etal..Onthecomputationofdistribution-
freeperformancebounds:Applicationtosmallsamplesizesinneuroimaging%****␣template.tex␣Line␣525␣****PatternRecognition93,1-13(2019)\par\bibitem{Haussler92}D.Haussler.DecisiontheoreticgeneralizationsofthePACmodelforneuralnetandotherlearningapplications.InformationandComputationVolume100,Issue1,September1992,Pages78-150\par\bibitem{Boucheron13}S.Boucheronetal.ConcentrationInequalities:ANonasymptoticTheoryofIndependenceISBN:9780199535255OxfordUniversityPress\par\bibitem{MacAllester2013}D.McAllester,APAC-
Bayesiantutorialwithadropoutbound.arXiv10.48550/ARXIV.1307.2118,2013\par\bibitem{Breusch79}T.S.Breuschetal.ASimpleTestforheteroscedasticityandrandomcoefficientvariation.Econometrica,Sep.,1979,Vol.47,No.5(Sep.,1979),pp.1287-1294\par\bibitem{Koenker81}R.Koenker.Anoteonstudentizingatestforheteroscedascity.JournalofEconometrics17.107-l12.North-
HollandPublishingCompany(1981).\par\bibitem{Friston13}K.J.Friston.Samplesizeandthefallaciesofclassicalinference.NeuroImage81(2013)503–504.\par\bibitem{Soft19}Addinsoft,2019.XLSTATstatisticalanddataanalysissolution,LongIsland,NY,USA.<https://www.xlstat.com>\par\bibitem{Rosenblatt16}J.D.Rosenblatt,etal.Better-
than-
chanceclassificationforsignaldetection.Biostatistics(2016).\par\bibitem{NAS2019}NationalAcademiesofSciences,Engineering,andMedicine.(2019).ReproducibilityandReplicabilityinScience.Washington,DC:TheNationalAcademiesPress.https://doi.org/10.17226/25303.\par\bibitem{Noble20}S.Noble,etal.Clusterfailureorpowerfailure?Evaluatingsensitivityincluster-
levelinference.NeuroImage,209,116468,2020.\par\bibitem{Phipson2010}B.Phipsonetal.PermutationP-
valuesShouldNeverBeZero:CalculatingExactP-
valuesWhenPermutationsAreRandomlyDrawn.StatisticalApplicationsinGeneticsandMolecularBiology:Vol.9:Iss.1,Article39.(2010)\par\bibitem{Gorgen18}Gorgen,K.,etal.Thesameanalysisapproach:Practicalprotectionagainstthepitfallsofnovelneuroimaginganalysismethods.NeuroImage,180,19-30.2018.\par\bibitem{Zhang14}Y.Zhangetal.Multivariatelesion-
symptommappingusingsupportvectorregression.HumBrainMapp.2014Dec;35(12):5861-76.\par\bibitem{Gorriz2021b}JMGorriz,etal.AconnectionbetweenpatternclassificationbymachinelearningandstatisticalinferencewiththeGeneralLinearModel.IEEEJournalofBiomedicalandHealthInformatics2021.\par\bibitem{Gorriz2022}JMGorriz,etal.Ahypothesis-
drivenmethodbasedonmachinelearningforneuroimagingdataanalysis.NeurocomputingVolume510,21October2022,Pages159-171\par\bibitem{Wang07}ZWang,etal.Supportvectormachinelearning-
basedfMRIdatagroupanalysis.NeuroImage36(4),1139-1151.2007\par%****␣template.tex␣Line␣575␣****\bibitem{Wang09}ZWang.AhybridSVM–GLMapproachforfMRIdataanalysis.Neuroimage46(3),608-615.2009.\par\par\par\end{thebibliography}\par%\clearpage%\appendix{}\section*{Anoteonthetheoreticallossespresentedincaptionoffigure\ref{fig:example1}}InFigure\ref{fig:example1},weclaimedthatthetheoreticalexpectedlossonauniformdistributionofpairs$\\{y,\hat{y}\\}$was$\mathcal{R}=\frac{b}{2}+\frac{a^2}{6b}$for$\mathcal{L}_1$and$\mathcal{R}=\frac{b^2+a^2}{3}$for$\mathcal{L}_2$,where$a$and$b$arethemaximumvaluesfor$\hat{y}$and$y$,respectively,whenthenumberofpredictorsis$P=1$.Ingeneral,theexpectedlossis:\par\begin{equation}\mathcal{R}=E[\mathcal{L}(f,x,y)]=\int\mathcal{L}(f,x,y)dP(x,y)\end{equation}\parInthissimplecaseforonly2dimensions,assumethat$\\{x,y\\}$areasetofuncorrelatedsamplescenteredattheorigin.Then,thedifference$(y-\beta\cdotx)$isuniformlydistributedaroundzeroaswell,andequallyprobablewithinavolume$\mathcal{V}$.Define$a$and$b$astheboundsfortheuniformlydistributedpairs$\\{\hat{y}=\beta\cdotx,y\\}$containedinthisvolume.Then$\mathcal{V}=2a\times2b$,andtheexpectedlosscanbecomputedas:\par\begin{equation}\label{eq:Ru}\mathcal{R}_u=\frac{1}{\mathcal{V}}\int\mathcal{L}(\hat{y},y)d\hat{y}dy\end{equation}\parwhere$\mathcal{V}=4ab$.If$\mathcal{L}(\hat{y},y)=|y-\hat{y}|$,thisexpressionisequivalentto:\par%****␣template.tex␣Line␣600␣****\begin{equation}\label{eq:Ru1}\mathcal{R}_u=\frac{2}{\mathcal{V}}\int_{-a}^{a}\int_{-b}^{\hat{y}}\mathcal{L}(\hat{y},y)d\hat{y}dy=\frac{2}{4ab}(\frac{1}{3}a^3+b^2a)=\frac{1}{2}b+\frac{1}{6}\frac{a^2}{b}\end{equation}\parIfweselect$\mathcal{L}=(y-\hat{y})^2$,thelossequation\ref{eq:Ru}isequivalentto:\par\begin{equation}\label{eq:Ru2}\mathcal{R}_u=\frac{1}{\mathcal{V}}\int_{-a}^a\int_{-b}^{b}\mathcal{L}(\hat{y},y)d\hat{y}dy=\frac{1}{4ab}(\frac{4}{3}b^3a+\frac{4}{3}ba^3)=\frac{1}{3}(b^2+a^2)\end{equation}\parNotethatweanticipatethesolutiontobeflat($\beta\sim0$)withuncorrelateddata.Subsequently,as$a$approacheszero,thesolutionsconvergetothemeanvalueofthelossevaluatedontheobservedvariable(halfoftheintervalandone-
thirdofthesquare,respectively).Anyalgorithmicdeviationfromthisidealsolutionprovidesthetheoreticalexpectedlossdescribedinequations\ref{eq:Ru1}and\ref{eq:Ru2}.ByexaminingFigure\ref{fig:example1},wecanassesshowsamplingandnon-
idealflatnessaffecttheconvergencetothesetheoreticalvalues.\par\par\end{document}