A hypothesis-driven method based on machine learning for neuroimaging data analysis

A hypothesis-driven method based on machine learning for neuroimaging data analysis J.M. Gorriz a,b,c, ⇑ , R. Martín-Clemente e , C.G. Puntonet d , A. Ortiz f , J. Ramírez a , SiPBA group a , J. Suckling b a DaSCI Institute, University of Granada, Spain b Dpt. of Psychiatry, University of Cambridge, UK c ibs.Granada, Granada, Spain d Dpt. Computer Architecture and Technology, University of Granada, Spain e Dpt. Signal Theory and Communications, University of Seville, Spain f Dpt. Communication Engineering, University of Málaga, Spain article info Article history: Received 15 June 2022 Revised 18 July 2022 Accepted 3 September 2022 Available online 9 September 2022 Keywords: General Linear Model Linear Regression Model Support Vector Regression permutation tests Magnetic Resonance Imaging Random Field Theory abstract There remains an open question about the usefulness and the interpretation of machine learning (ML) approaches for discrimination of spatial patterns of brain images between samples or activation states. In the last few decades, these approaches have limited their operation to feature extraction and linear classification tasks for between-group inference. In this context, statistical inference is assessed by randomly permuting image labels or by the use of random effect models that consider between-subject variability. These multivariate ML-based statistical pipelines, whilst potentially more effective for detecting activations than hypotheses-driven methods, have lost their mathematical elegance, ease of interpretation, and spatial localization of the ubiquitous General linear Model (GLM). Recently, the estimation of the conventional GLM parameters has been demonstrated to be connected to an univariate classification task when the design matrix in the GLM is expressed as a binary indicator matrix. In this paper we explore the complete connection between the univariate GLM and ML-based regressions. To this purpose we derive a refined statistical test with the GLM based on the parameters obtained by a linear Support Vector Regression (SVR) in the inverse problem (SVR-iGLM). Subsequently, random field theory (RFT) is employed for assessing statistical significance following a conventional GLM benchmark. Experimental results demonstrate how parameter estimations derived from each model (mainly GLM and SVR) result in different experimental design estimates that are significantly related to the predefined functional task. Moreover, using real data from a multisite initiative the proposed ML-based inference demonstrates statistical power and the control of false positives, outperforming the regular GLM. Ó2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). 1. Introduction Whole-brain analyses in neuroimaging (NI), comprising a large number of independent statistical tests, have been traditionally conducted with classical statistics, either hypothesis testing or Bayesian inference, and the univariate General Linear Model (GLM) [8]. These hypothesis-driven methods gained their popularity due to the ease of interpretation and function localization across experimental designs [9]. However, they usually rely on assumptions that are frequently violated; e.g. homogeneity, Gaussianity, etc. [30] and, consequentially, inflated type I error rates have become problematic and a key contributor to the replication crisis [8,23]. Furthermore, technological advances are increasing spatial and temporal resolutions as well as the range of available measurements of anatomy and physiology; a true exemplar of the curse of dimensionality [1]. In this context, analyses of contemporary large image repositories retain the difficulties associated with small sample sizes. One of them is the inflated false-positives observed across experiments as a consequence of the multiple comparison problem. This problem is partly solved by over conservative approaches, such as Bonferroni or Random Field Theory (RFT) corrections [7]. One promising solution for the aforementioned problems in NI is machine learning (ML) [11,12]. ML provides us with highdimensional relationships between datasets that are empirically established based on data-driven methods [36]. Estimating depenhttps://doi.org/10.1016/j.neucom.2022.09.001 0925-2312/Ó2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). ⇑ Corresponding author at: DaSCI Institute, University of Granada, Spain. E-mail address: [email protected] (J.M. Gorriz). Neurocomputing 510 (2022) 159–171 Contents lists available at ScienceDirect Neurocomputing journal homepage: www.elsevier.com/locate/neucom dencies in regression or classification tasks using statistical learning theory (SLT), unlike classical statistics, characterizes the actual relationships or effects with a limited dataset. NI in particular has embraced ML as a technology to deliver diagnostic and prognostic classification of neurological and psychiatric disorders [4,19,42,21,20]. Nevertheless, the mainstay of NI studies are observational and mechanistic, seeking to identify regional betweengroup differences in brain structure and function. 1.1. Related Works Data-driven analysis methods based on ML [27,39,40,33,13] have demonstrated their ability for detecting activations in fMRI data, outperforming conventional hypothesis-driven approaches, i.e. the standard GLM inference based on random effect models (see Fig. 1). The core idea on these agnostic (model-free) approaches is to perform an accurate feature extraction based on a fixed-complexity MLE classifier, e.g. linear support vector machine (SVM), between predefined groups. They all share the same characteristic processing pipeline of a data-driven multivariate approach that enhances detection ability within a classification task. Consequently, the Statistical Parametric Maps (SPM) derived from the GLM parameters are replaced by the spatial discriminance map (SDM) 1 obtained by ML [37,31]; or some other specific feature extracted at the training stage, e.g. distance to the separating hyperplane [39], replaces the explanatory variables in the GLM, previous to the statistical inference stage. These data-driven maps are then deployed in conventional pipelines (hybrid approaches) for statistical inference. However, some approaches depart from the p-value- based frequentist, Bayesian or permutation analyses, and introduce the concept of the probability of the worst case in neuroimaging [13]. Despite the popularity of ML as a solution for a wide range of complex problems, there remains an open question about its usefulness for statistical inference. Mainly, what is the statistical significance of its performance in a classification/regression task? Moreover, several works in NI [37,38,11] have revealed the shortcomings of predictive inference methods based on ML, e.g. crossvalidation failure results in poor estimates. Another example is the inclusion of covariates and nuisance variables 2 in the analysis using ML. This procedure, that is massively used in the conventional NI analysis, is usually overlooked when using ML tools [19,42,21,20,25]. At most, the aforementioned hybrid approaches make use of a fitting process in the ‘‘ML space” in combination with the classical GLM, which includes these variables [17]. Efforts with ML around these issues are increasing with continuous output variables ([5], with remarks in [29]) rather than the more typical categorical classifications. However, the classification task is just a particular case of the regression problem with discrete labels, thus exploring general ML methods for linear regression, such as SV regression (SVR), is currently relevant in addition to its use as a simple extension of SVM [42]. Recently, a connection between both domains, that is, the standard GLM and classification tasks by ML, has been established [14] using binary experimental design matrices. The aim was to formalize the relationship that has been evidenced in several neuroimaging applications in the extant literature. In this paper, we show a complete and novel connection between the classical GLM, including random effect models, and the ML framework in the estimation of optimum regression parameters. 1.2. Aims Overall, the main achievement of the article is the formalization of the use of ML in the context of hypothesis-driven neuroimaging statistical inference by deriving its formal relationship to the GLM; the method almost ubiquitously used for mass univariate testing of imaging datasets. Outputs from both approaches are then compared in a variety of commonly encountered applications. In doing so, we leverage one of the most commonly used packages for mass Fig. 1. NI analysis and methods. An univariate hypothesis-driven method based on the regression parameters obtained by ML (dash line) is proposed. Note the different meaning among the activation maps prior to probability thresholding, e.g. by means of a multiple comparison correction (RFT, Bonferroni, etc.) in classical approaches or selecting accuracies above random performance in SDM. In particular, SPM obtained by hypothesis testing indicates the likelihood of the observation given the null hypothesis H 0 (data likelihood) whilst the Bayesian framework considers the posterior probability map of the observation (effect likelihood). SDM and current multivariate ML approaches for image classification provide the (class-membership) probability of the observation given the alternative hypothesis H 1 (prevalence). MANCOVA: Multivariate analysis of covariance, MANOVA: Multivariate analysis of variance. 1 mainly based on accuracy or prevalence, that is, the proportion of a population who have a specific characteristic or effect versus control subjects 2 they are usually referred to as confounds, i.e. variables that may blur the effect that is being sought. J.M. Gorriz, R. Martín-Clemente, C.G. Puntonet et al. Neurocomputing 510 (2022) 159–171 160 univariate testing: SPM [8], and within it the methods for spatially extended statistics and p-value corrections for multiple comparisons. In summary, in this paper: We formalize the use of ML within the univariate GLM pipeline for statistical inference, We give an explanation on why hybrid approaches, massively used in this context, are actually working, e.g. the main regressor (experimental condition) in the design matrix could be replaced by a data-driven property [39]. The connection is established by selecting the SVR in the inverse domain (SVR-iGLM) that allows us to include all the covariates (age, sex, etc.) in the analysis at once. Subsequent analyses based on frequentist inference and permutation testing are carried out to calculate the level of significance in between-group testing of SPMs with a multiple comparisons correction. 2. Theory 2.1. The General Linear Model and its statistical framework The GLM [9] is defined for a single observation level, e.g. in a between-subject comparison, as: y¼Xhþ  ð1Þ where yis the N1 observation vector with units of time, signal change, volume, etc., is the N1 vector of errors that is assumed to be Gaussian distributed, Xis the NMmatrix containing the explanatory variables or constraints, and his the M1 vector of parameters explaining the observations y. Note that: i) for a hierarchical observation model each level requires prior estimation at the previous levels; and ii) in terms of ML, Xare the multidimensional labels or regressors acting on the observations y. In the classic GLM, his usually estimated by a maximum likelihood criterion based on the Gaussianity assumption and is given by: ^ h¼ðX t C 1  XÞ 1 X t C 1  yð2Þ where C  is the covariance matrix of errors. Inferences on this estimate 3 determine the components of h, and the relationship between classical GLM and ML-based prevalence inferences can be obtained using a linear compound specified by a contrast weight vector c, and writing a T statistic as: T¼c t ^ h ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c t Co v ð^ hÞc qð3Þ where Covð^ hÞ¼ðX t C 1  XÞ 1 . This T statistic gives us the probability of observing the ML estimation under H 0 , and when it is small enough, e.g. p<0:05, the linear compound is considered significantly different from zero. As an example, given a set of two parameters in h¼½h 1 ;h 2  t , if we select c¼½11we are assessing how large is the first parameter with respect to the second; i.e. the difference h 1 h 2 . Thus, if the T statistic suggests a small probability, the contrast is statistically significant with observations generated from different sources. A similar procedure could be established based on a Bayesian estimation and inference to handle complex hierarchical observational models. This framework would be based on the Expectation Maximization (EM) algorithm for parameter estimation, along with known priors and a priori probability models, with the aim of evaluating the posterior probability (ppm). By thresholding the ppm, relationships between this and the frequentist approach can be established for both their similarities (statistical power) and differences (specificity) [9]. 2.2. Machine learning and the inverse GLM The GLM can be interpreted as the inverse problem of regressing the observations onto the conditions (see Fig. 2). Instead of assuming the model in Eq. 1, the inverse GLM is defined as: X¼y x þ^  ð4Þ where xis a set of (1 M) parameters that best explains the design matrix given the observations, and ^ is noise with unknown pdf. We can readily see that: ðX^  Þ x t ð xx t Þ 1 ¼yð5Þ where we assume that the inverse of the norm ðxx t Þ 1 exists. After some manipulations we finally get: y¼X~ hþ~  ð6Þ where we define: ~ h¼ x t ð xx t Þ 1 ;~  ¼ ^ x t ð xx t Þ 1 ð7Þ Therefore, solving the multiple regression problem in equation 4, e.g. using a ML approach, is equivalent to estimating the parameters of the GLM. If we carefully examine Eq. 7, each column cof the design matrix (c¼1;...;Mincluding covariates and nuisance variables) can be described by: x 1;c . . . x N;c 2 6 6 4 3 7 7 5 ¼y x c þ^  c ð8Þ The set of parameters xcould be determined by a classical approach such as Parametric Empirical Bayes (PEB), assuming a Gaussian model for the noise ^ . However, we prefer to use the ML approach based on the regularized risk minimization and linear Fig. 2. Illustration example of the connection between GLM and iGLM. On the left observations in yand the best parameters of the GLM (h 1 ;h 2 ) given the experimental conditions (blue and red) are shown. On the right, given the set of observations, we show the parameters (x 1 ;x 2 ) that better explain the classes (red and blue) and the equivalent GLM parameters (^ h 1 ;^ h 2 ). 3 Here, we refer to voxelwise inference since we use a threshold uto classify voxels ias ‘‘active” if TiPu. Clusterwise inference uses a cluster-forming threshold to define contiguous suprathreshold regions [28]. J.M. Gorriz, R. Martín-Clemente, C.G. Puntonet et al. Neurocomputing 510 (2022) 159–171 161 regressors. For the reader’s convenience, further mathematical discussion is available on 7. Thus, given the set of observations yand the experimental conditions and covariates in X, we regress Mindependent linear equations as: X¼ywþB¼^ y^ wð9Þ where ^ w¼x B  , the NMmatrix B¼½b t ...b t  t and bis the 1 M vector of biases. In view of Eq. 8, we realistically assume under this approach that the bias does not depend on the experimental condition or covariate realization. Once the set of parameters and biases fw;bgare estimated by a suitable procedure (see Appendix), such as least squares (LS) or SVR, we can calculate the observation yby simply inverting Eq. 9 as: y est ¼ðXBÞw t ðww t Þ 1 ¼~ X~ hð10Þ where ~ his the ML estimated vector of parameters that combines the weight of multiple predictors to explain the observation, and ~ X¼ðXBÞis the adjusted design matrix. Instabilities that could arise from the scalar inversion in Eq. 10, that depend on the selected ML algorithm, are easily solved by bounding its value between 1 and 1. It is worth mentioning here that our model does not use any data reduction techniques, e.g. principal component analysis (PCA), preserving the function localization of the univariate GLM. However, a similar description of the method can be given in terms of the difference of signed distances by replacing the main reference function in the GLM by yw jjwjj 2 . In this way we a give an explanation of the multivariate approach proposed in [39] in terms of the solution in the iGLM. 2.3. Support Vector Regression In general, the core idea of SVR [32] is to do (non-) linear regression in a feature space F: fðyÞ¼wUðyÞþb;U:R n !F;w2Fð11Þ where denotes dot product and, in the case of the linear regression, Uis simply the identity function. We determine wfrom the data by minimizing the sum of the empirical risk, e.g. -insensitive loss function [2], and a complexity term proportional to its norm: R reg ðfÞ¼1 2jjwjj 2 þCX N i¼1 jx i fðy i Þj  ð12Þ This minimization can be transformed into a uniquely solvable quadratic programming problem [34] that provides the vector of parameters win terms of the samples or support vectors: w¼X N i¼1 ð a i  a  i ÞUðy i Þð13Þ where a i ;a  i are Lagrange multipliers; that is, the solutions of the quadratic programming problem. Finally, the bias term bcan be computed by determining the prediction error on the margin d i ¼fðy i Þx i ¼signða i a  i Þand taking the average of differences as b¼<fðy i ÞwUðy i Þ>. 2.4. Model Equivalence After performing the regression in the iGLM domain by SVR and assuming that the explanatory matrix Xcontains two experimental conditions, i.e. indicator variables that refer to class membership, and a set of continuous covariates or nuisance variables, then we have X¼½X b X c and the observations can be estimated, following Eq. 10, by: Fig. 3. Simulated data with noisy observations example with h¼½101 t , CNR¼1 and N¼1000. We evaluated the GLM in Eq. 1given the experimental conditions (main regressors) with a covariate and the Gaussian noise assumption. J.M. Gorriz, R. Martín-Clemente, C.G. Puntonet et al. Neurocomputing 510 (2022) 159–171 162 y est ¼X b ~ h b þX c ~ h c B~ hð14Þ where ~ h¼½ ~ h t b ~ h t c  t . Eq. 14 can be used to approximate a set of parameters that best explain the observation vector with maximum and minimum influence of covariates and nuisance variables ranged in the interval ½0;1as (see an example in Fig. 2) ~ h est ¼~ h b þ½0;1~ h c B~ hð15Þ 3. Materials and Methods 3.1. Synthetic data We generated synthetic data with the aim of modeling different scenarios in an fMRI time-series analysis with inserted activations using a block-design (baseline and task) paradigm [39]. For this purpose, we simulated one dimensional observation vectors y (Eq. 1) with different contrast to noise ratios (CNR: 0:25;0:5;0:75;1) and sample sizes (N:100–1000). The design matrix Xwas the canonical HRF convolved boxcar function for fMRI simulated data. This matrix contained an exponentially decaying function fðtÞ¼ð1 t N Þ 0:5 representing habituation during the fMRI task, or a covariate to simulate the effect of age when spatially testing brain activations, and as shown in Fig. 3.AnN-dimensional Gaussian noise vector v was randomly drawn with zero mean and controlled variance ( r 2 X =CNR). Finally, a vector of observations was constructed by adding the noise to the design matrix with ideal parameters h¼½10c v  t , where c v is a constant that modulates the exponentially decreasing covariate in the time-series. The use of synthetic data allowed us to estimated the noise covariance matrix by averaging a set of 100 noise realizations as an ideal comparison with the GLM benchmark. This is useful to simulate different performances of the ReML estimation that depends upon the observation vector yand some parametrization of the covariance components. The latter is the procedure in SPM12 that obtains the noise covariance matrix C  in Eq. 3. 3.2. A structural MRI dataset including covariates: ADNI Data were obtained from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). The ADNI database contains T1-weighted structural MRI scans acquired at 1.5 T and 3.0 T from patients with Alzheimer’s disease (AD), Mild Cognitive Impairment (MCI), and cognitively normal controls (NC) at multiple time points. Here we only included structural MRI collected at 1.5T. The original database contains more than 1000 T1- weighted MRI images in total, although for this study only the first MRI examination of each participant was included, resulting in 417 structural images in the sample. Demographic data is summarized in Table 1. The dataset was processed using the standard neuroimaging methods and protocols implemented in the SPM software (www.- fil.ion.ucl.ac.uk/spm/), including registration to MNI space by spatial normalization and segmentation, to generate maps of grey matter (GM) volume [8]. Following the recommendation of the National Institute on Aging and the Alzheimer’s Association for the use of imaging biomarkers [22], we considered the group comparison NC vs. AD for establishing a clear framework for comparing statistical paradigms. Age, sex and intracraneal volume (ICV) were included as covariates in between-group modelling as representing the most common set in the extant literature [17]. All the covariates were standardized with zero mean and standard deviation equal to one. We selected two regions of interest using the 116-area automated anatomical labeling template [35]: one relevant area in AD, the bilateral hippocampus (denoted Hippocampus L and Hippocampus R for left and right hemispheres, respectively) with 2559 voxels; and the cerebellum (Cerebelum9 L and Cerebelum9 R of the atlas), which is considered not relevant to AD pathology, with 3027 voxels. 3.3. Statistical Analyses To assess the performance of the methods presented in this paper, we used the benchmark proposed in SPM12 for statistical inference; that is, the GLM and RFT with FWE correction and p¼0:05 for second-level statistical inference. First, we estimated the best set of parameters by regressing the design matrix, or the observations, using several configurations: i) the ideal ML method with synthetic data, where C  is estimated from noise realizations, ii) Restricted (Re) ML, iii) LS, and iv) SVR. Then, we connected the ML-based estimates (LS and SVR) with the corresponding set of parameters in the GLM space, as shown in Eq. 15. To compare the inferences of each regression method, we assumed the same noise model and evaluated the T statistic in Eq. 3on the set of parameters. Finally, we thresholded the resulting T-maps, e.g. derived from the SVR-iGLM method, by a detection threshold based on RFT. A permutation analysis was also adopted to provide an alternative statistical inference based on a non-parametric approach [27]. By randomly permuting the experimental conditions 1000 times, we calculated the non-parametric T statistic based on the contrast of the estimates, thus avoiding the estimation of the denominator of Eq. 3).The probability of observed contrast was then calculated relative to the distribution of permuted contrasts representing the null distribution. If this value was less than a selected threshold, e.g. p¼0:05, then we rejected the null hypothesis. 4. Experimental results 4.1. Synthetic data: Estimating functional tasks With the synthetic data we estimated the parameters by all the methods followed by permutation inference. We additionally plotted the non-normalized statistic (contrast) and evaluated the set of parameters, as a classifier, in the label domain, wy> <0. The observations were artificially drawn from the same Gaussian sample distribution (Fig. 3) with varying sample size (i.e. time points) and CNR. We regressed both explanatory variables and observations to obtain the experimental parameters for each model h;wgiven the simulated covariate. All these estimations were employed to calculate the regressed observed variables using Eqs. 1 and 10, given the explanatory matrix and the estimated parameters, and illustrated in Figs. 4 and 5. Only for extremely Table 1 Demographic details of the MRI ADNI dataset, with group means and their standard deviations Status N Age Sex (M/F) ICV(10 5 ÞMMSE NC 229 75.95.0 119/110 15.31.6 29.11.0 AD 188 75.37.5 99/89 15.51.8 23.22.2 J.M. Gorriz, R. Martín-Clemente, C.G. Puntonet et al. Neurocomputing 510 (2022) 159–171 163 noisy observations (CNR<0:5) and considerably high sample sizes did the ReML (and the ideal ML) outperform ML approaches in MSE, even under the Gaussian assumption (Fig. 4). However, in the label domain, the set of parameters derived from the ML approaches provided the highest classification accuracy in the associated classification task, recapitulating prior findings [13].In summary, the ML parameter fitting in the GLM does not imply an accurate regression in the iGLM. The p-value of the estimated contrast (difference between parameters h 1 h 2 ) was less than p¼0:05 for all methods. The probability of observation was p¼1=1001 in all cases when testing the null hypothesis. Thus, no false positive were detected during the simulated task, although the SVR estimation based inference provided results closest to the nominal false positive rate, whilst the remainder were over-conservative (Fig. 6). 4.2. Empirical data: a case-control design with the ADNI Dataset In this section we show the inference derived from the two methodologies in each domain. We regressed on the observations and on the design matrix including covariates for age, sex and ICV (Table 1). Then, we constructed the spatially extended statistical processes, generating maps of significance, using GM estimated from the MRI ADNI dataset[13]. We compared the SPM (a twosample T-statistic similar to Eq. 3), where significance is first individually assessed at each voxel, and then on clusters with p¼0:05 FWE corrected based on RFT. The univariate test based on the inverse GLM in Section 2.2 was also conducted. In short, we will assess the detection ability (true positive detection) and the control of type I error (false positives) of the methods in a semi-controlled environment. To this purpose we selected two baseline regions in late AD, the hippocampus, that is primarily and generally affected in early AD; and the cerebellum, that remains unaffected. Finally, we will explore the detection ability within an Omnibus test (whole-brain analysis) using the multiple comparison correction provided by RFT and the optimum threshold derived from the baseline regions. Fig. 4. Top: Example estimation with N¼1000, CNR¼1. Bottom: Averaged mean squared error (MSE) with different activation CNR and sample sizes (CNR: 0:25;0:5;0:75;1 and N:100,200,...1000). Fig. 5. Distribution of observations (in grey) separated by condition and estimations h(symbols in blue and red) for all the analysed methods including the covariate effect as shown in Eq. 15. Fig. 6. Shadeplots (average and standard deviation) showing contrast estimations (h 1 h 2 ) for each method under the null hypothesis. They were calculated from 1000 random permutations and CNR¼1. J.M. Gorriz, R. Martín-Clemente, C.G. Puntonet et al. Neurocomputing 510 (2022) 159–171 164 Fig. 7. Contrasts and T-statistics in the hippocampus with RFT FWE-corrected threshold at p¼0:05 and varying sample sizes. Note that the ML approaches are upper bounded for visualization purposes and the threshold plotted T thres ¼5:2758 is that obtained for N¼350. J.M. Gorriz, R. Martín-Clemente, C.G. Puntonet et al. Neurocomputing 510 (2022) 159–171 165 4.2.1. Hippocampal analysis to assess statistical power First, we performed the univariate statistical analysis on GM in the left and right hippocamppal ROIs, which are typically the earliest sites of atrophy associated with AD [18]. Thus, the group comparison (advanced stage of the disease) AD vs. NC can be considered as a true positive (TP) region. We compared the GLM contrast with that obtained by the ML-based approach using the same parametrization of the noise covariance in Eq. 3, the standard approach in SPM12. SPMs were collected voxelwise and the inference method, based on a two-sample T-test with the linear compound c¼½1;1;0;0;0, was thresholded by means of the RFT FWE rate correction at p¼0:05, as previously. Figs. 7 and 8 show the contrasts (shadeplots with increasing sample size) derived after model estimation, and the T-statistics obtained by the use of the same normalization term in Eq. 3. In this simulation, we employed the T-scores obtained by the standard SPM approach to bound the values obtained by our method. We readily see the higher statistical power (T>T thres ¼5:2758 for N¼350) of the ML methods and the strong dependence of the standard SPM on the sample size. Both methods converged as sample size increased. The LS approach was more unstable with fewer significant voxels than the SVR-based procedure (N¼ 100;2314 v s2341;N¼150;2300 v s:2324;N¼200;2300 v s:2329;N¼ 250;2334 v s:2355;N¼300;2355 v s:2372); N¼350;2344 v s:2360). 4.2.2. Type I error control in a putatively-null cerebellar region We repeated the experiments in cerebellar ROIs (Cerebellum 9 left and right [35]) to evaluate the ability of the ML-based inference methods to control false positive rates using the same inference strategy used for thresholding the T-maps; i.e. FWE correction based on RFT in the standard GLM. Fig. 9 illustrates the overconservative voxelwise inference with FWE correction based on RFT, although the proposed methods based on ML provided a significant number of tests closer to the expected value, i.e. p¼0:05. In other words, by chance the number of FPs should be around 5%of the total number of voxels. Parametric voxelwise inference is known to be valid but conservative, often falling below the nominal rate [6]. 4.2.3. Global analysis in contrast images with limited sample size Fig. 10 shows the vector of parameters in the selected axial slice Z = 47 derived from ReML and SVR in the analysis of the whole brain. The ML approach implemented in SPM finds very few significant relationships between the covariate effects (age, sex and ICV) and the observations, compared with the experimental conditions. On the contrary, the SVR yields stronger connections, mainly for the sex covariate. From these parameter images, contrast images are then derived and, subsequently, inference on their size relative to the estimate of their standard error is made in large target regions, as depicted in Figs. 7 and 9. Selecting the aforementioned target regions and plotting the tscore histograms we can approximately compute the optimum threshold at a given level of significance a using the Neyman- Pearson lemma, e.g. at a ¼0:05, SVR: t thresh ¼10:03 and SPM: T thresh ¼2:03. By comparing the latter value with the one used in previous sections (RFT correction) we can readily see the overconservative nature of the standard voxelwise inference (Fig. 11). Extending these thresholds to the whole volume we obtain the activation maps for the SVR and standard SPM approaches as shown in Fig. 12. (See Fig. 13). 5. Discussion In the context of neuroimaging statistical inference, there is an increasing trend to incorporate exploratory methods into wellestablished GLM-based data analysis. Not only data preprocessing techniques, such as independent or principal component analysis (PCA) [26], but also multivariate ML approaches have been widely used in classification tasks to replace the predefined design matrix in the regular GLM pipeline [39] or to provide novel statistical maps of prevalence [27,13]. Existing multivariate approaches based on ML, including stages for smoothing or orthogonal decompositions, such as PCA [27,39,40,13], have provided promising results in (f) MRI-data analysis where there is a trade-off between sensitivity and computational cost. In this paper we present a novel univariate methodology for (f) MRI image analysis based on the optimum performance in limited Fig. 8. Group-level statistical analysis results in the hippocampus of the ADNI MRI data. The univariate analysis was conducted using the standard GLM (top) and SVRRFT approach (bottom). Note T-maps scales for comparison. The T-maps were thresholded at u>5:5704 (equivalent to an uncorrected p<1:1802e 7 ) with sample size N¼100. J.M. Gorriz, R. Martín-Clemente, C.G. Puntonet et al. Neurocomputing 510 (2022) 159–171 166 Fig. 9. Contrasts and T-statistics for the cerebellum with FWE-corrected threshold t thres ¼5:2758 at p¼0:05. J.M. Gorriz, R. Martín-Clemente, C.G. Puntonet et al. Neurocomputing 510 (2022) 159–171 167