Open Access
14 June 2023 Random matrix theory tools for the predictive analysis of functional magnetic resonance imaging examinations
Derek Berger, Gurpreet S. Matharoo, Jacob Levman
Author Affiliations +
Abstract

Purpose

Random matrix theory (RMT) is an increasingly useful tool for understanding large, complex systems. Prior studies have examined functional magnetic resonance imaging (fMRI) scans using tools from RMT, with some success. However, RMT computations are highly sensitive to a number of analytic choices, and the robustness of findings involving RMT remains in question. We systematically investigate the usefulness of RMT on a wide variety of fMRI datasets using a rigorous predictive framework.

Approach

We develop open-source software to efficiently compute RMT features from fMRI images and examine the cross-validated predictive potential of eigenvalue and RMT-based features (“eigenfeatures”) with classic machine-learning classifiers. We systematically vary pre-processing extent, normalization procedures, RMT unfolding procedures, and feature selection and compare the impact of these analytic choices on the distributions of cross-validated prediction performance for each combination of dataset binary classification task, classifier, and feature. To deal with class imbalance, we use the area under the receiver operating characteristic curve (AUROC) as the main performance metric.

Results

Across all classification tasks and analytic choices, we find RMT- and eigenvalue-based “eigenfeatures” to have predictive utility more often than not (82.4% of median AUROCs > 0.5; median AUROC range across classification tasks 0.47 to 0.64). Simple baseline reductions on source timeseries, by contrast, were less useful (58.8% of median AUROCs > 0.5, median AUROC range across classification tasks 0.42 to 0.62). Additionally, eigenfeature AUROC distributions were overall more right-tailed than baseline features, suggesting greater predictive potential. However, performance distributions were wide and often significantly affected by analytic choices.

Conclusions

Eigenfeatures clearly have potential for understanding fMRI functional connectivity in a wide variety of scenarios. The utility of these features is strongly dependent on analytic decisions, suggesting caution when interpreting past and future studies applying RMT to fMRI. However, our study demonstrates that the inclusion of RMT statistics in fMRI investigations could improve prediction performances across a wide variety of phenomena.

1.

Introduction

Though first developed to describe the fluctuation of nuclei energy levels in quantum physics,1,2 Random matrix theory (RMT) has been shown to have extremely broad potential. In small-scale physical systems, RMT universalities have been observed in quantum chaotic systems, complex nuclei, atoms, molecules and disordered mesoscopic systems;17 at larger scales, RMT has been applied to atmospheric physics,8 stock cross-correlations,9 social networks,10 random networks,11 network-formation in liquids,12,13 and amorphous clusters.1416 Within biological systems, RMT has also been used to successfully model aspects of amino acid functional relationships,17 synchrony in epileptic seizures,18 and in protein–protein interactions both in different species19 and breast cancer.20 RMT has also been used to guide statistical decisions in principal components analyses2123 and, more recently, has provided insights into the behaviors of deep neural networks.24,25

At its most rudimentary, RMT describes the expected behavior of the eigenvalues—also often called the spectra or levels—of a number of classes of random matrices.1,26,27 A random matrix is a matrix where each entry is a random variable. Typically, each random variable is independent and identically distributed (iid), however, modern extensions allow for some dependence.28 A class or ensemble of random matrices can be defined by the type of iid distribution involved. For example, the Gaussian orthogonal ensemble (GOE1,2) is a “class” or “ensemble” of random matrices that comprises orthogonal matrices with each entry being sampled from a standard Gaussian, and the Marchenko–Pastur distribution describes the limiting behavior of the singular values of iid rectangular random matrices where the iid distribution is arbitrary.1 RMT often finds that in the infinite limit, the expected eigenvalue distribution of a number of classes of random matrices can be described quite precisely.1,26,27

Of course, given a particular non-random matrix or instantiation of a random matrix, one can, in general, only estimate the likelihood that it is a member of a certain class or ensemble since certain RMT ensembles (such as the GOE) can yield arbitrary matrices with extremely low probability. Thus RMT also includes statistical measures (“spectral observables”) and tools to help compare empirical observations to theory.1,2

When or if RMT has real-world explanatory potential, this will most likely be when dealing with a complex system of many (hundreds or more1) interacting components, rather than a system that is too small for statistical regularities to reliably surface. If such a system has a matrix representation, and the eigenvalues and spectral observables of this representation have distributions similar to those predicted by RMT, it suggests that the system is either highly random or chaotic. By contrast, if the observed spectra deviate significantly from RMT-predicted spectra, this suggests otherwise. A number of studies have used RMT to make such interpretations and comparisons between systems.8,10,11,19,20,2932

1.1.

RMT and Neurobiological Signals

In the human brain, each neuron, collection of neurons, or region of interest (ROI) is a potentially interacting component in a complex system. RMT may have potential in describing the totality of these interactions, provided that measurements of functioning can be obtained with sufficient spatial and temporal resolution to speak to some neurobiological or neuropsychological phenomenon of interest.

For an imaging modality like functional magnetic resonance imaging (fMRI), where changes in the blood-oxygenation-level-dependent (BOLD) signals are related to neural activity, RMT may be an ideal starting point, as each voxel time course (or collection of such time sources, i.e., ROIs) can be considered an interacting component of the system. Likewise, in functional connectivity analyses, statistical relationships of the BOLD ROI time courses are investigated in the hope of gaining insights into brain function (see Refs. 33 and 34 for reviews, but also Refs. 35 and 36 for challenges facing functional connectivity analyses). In this framework, each connection or correlation can be considered a system component, and the eigenvalues of such a correlation matrix can be examined from the perspective of RMT.

The earliest study taking this approach demonstrated that spectra of the correlations between electroencephalographic signals closely resemble those of the GOE.29 In fMRI, RMT has been used to evaluate the quality of whole brain features extracted from fMRI data,37,38 and in diffusion MRI to aid in the selection of the number of components to employ in principal-component reduction analysis and denoising.22,23,38

RMT has also been used in ROI-based fMRI functional connectivity studies to investigate differences between rest and task states,30 between subjects with and without attention-deficit hyperactive disorder (ADHD),31 between pain and non-pain states,32 and between left-sided versus right-sided motor imagery.39 Although no differences were found in the latter study, across the first three studies, the spectra of resting or low-attention states exhibited properties close to the GOE. These findings suggest that certain aspects of psychological processes might be characterized, in part, by features computed from the eigenvalues of fMRI correlation matrices, and that these features might vary in an interpretable manner across psychological processes. If this is the case, RMT could aid in understanding the functioning of the human brain.

1.2.

Eigenvalue Features

The basic insight of RMT is thus that eigenvalues alone may provide interesting information about highly complex systems. However, real systems are usually noisy and involve a mixture of random and non-random components and interactions. RMT is statistical in nature, describing only the expected behavior of the spectra of iid random matrices. To this end, a number of summary statistics or “spectral observables”1,2 can be computed from empirically observed spectra, with these spectral observables sometimes being better suited for further analysis, or for the comparison of empirical observations to theory. Two such summary statistics that have been popular8,10,11,19,20,2932 are the “spectral rigidity” and “level number variance” (“rigidity” and “level variance” for short; details in Sec. 3).

However, from a predictive standpoint, these summary statistics may mask predictively useful information. If the basic RMT insight is that the eigenvalues alone can provide understanding of a system, then those eigenvalues (or other simple transformations of them) ought also to be predictively useful. We examine a number of such features in this study and refer to both RMT-derived features, such as the rigidity and level variance, and non-RMT-derived features collectively as “eigenfeatures” in this study.

1.3.

Functional Connectivity Reduction

The functional connectivity—often, the matrix of correlations of various collections of voxels of an fMRI scan—is a priori a useful representation of the fMRI data. However, the full correlation matrix between all voxels is itself often computationally infeasible to work with, and thus is reduced in various ways prior to being used in analyses.

Typical reductions include the use of a “seed” voxel or collection of voxels:40,41 the mean signal over that ROI is correlated with all other N ROIs of interest to generate a reduced functional connectivity matrix (or image) with one correlation value at each non-seed ROI. This sort of reduction reduces the functional connectivity to N values but necessarily biases the analysis to the seed ROI.

Likewise, one can work with ROI mean signal reductions, or independent component analysis reductions,40,41 and work with the N×N matrix of the correlations of these reductions. This analysis is, a priori, sensitive to the choice of ROIs, and, in the case of anatomical atlases/parcellations for ROIs, may also lack theoretical justification.

RMT suggests a potentially useful reduction in the form of the eigenvalues of the voxelwise functional connectivity matrix. An N×t matrix of N time series of length t, and where tN will have a symmetric correlation (or covariance) matrix with t1 non-zero, positive real-valued eigenvalues. These eigenvalues can be computed highly efficiently via transposition of the voxelwise correlation or covariance matrix.

1.4.

Limitations of Previous Work

Previous studies applying RMT to functional connectivity data31,32,39 took largely descriptive/explanatory approaches, noting only differences in RMT across subgroups and/or conditions. For example, Wang et al.31 and Gu et al.39 used a Kolmogorov–Smirnov test, and Wang et al. used a t-test to provide evidence of differences in various RMT metrics. However, statistically significant differences in a metric do not necessarily imply practical significance or predictive utility,42 nor do they imply generalization or replicability.43 Cross validation, by contrast, attempts to more directly assess these properties.44

In addition, the previous papers do not provide publicly accessible code to reproduce results. The extraction of the spectral rigidity and level variance is computationally demanding and mathematically and algorithmically non-trivial and additionally requires an “unfolding” procedure.1,2 Unfolding is an exponential fitting procedure that involves a number of highly subjective decisions regarding outliers and the flexibility of the fitting function. These decisions are, unfortunately, often poorly documented or even entirely missing from method descriptions, despite being known to often dramatically impact RMT conclusions.4549

In addition to the flexibility in the implementation and application of RMT to empirical data, there is also flexibility introduced by analytic choices made in the complex preprocessing pipelines of fMRI.50 These “researcher degrees of freedom,”51 in combination with the absence of reproducible code and data availability, and predominance of descriptive and explanatory approaches, may raise doubts about the basic robustness and practical utility of past RMT-based functional connectivity analyses.

What is missing is a rigorous, systematic, reproducible investigation of the predictive value of RMT metrics across a wide variety of data, analytic choices, RMT features, and preprocessing decisions, with RMT features being compared to simpler alternative baseline predictors. The current study aims to remedy this.

2.

Datasets

2.1.

Overview

We selected fMRI data publicly available on the OpenNeuro platform.52 Selection criteria were somewhat subjective, but we required each dataset to (1) comprise 10 or more human subjects, (2) have all fMRI runs have the same spatial and temporal resolutions, (3) have fMRI images that can be split into various classes for a number of binary classification tasks, and (4) allow reasonable classification of each run using all of the run 4D voxel data only. That is, for this last point, classification of a run should be reasonable in the absence of run-specific task or event timings or details and should not require using only some subset of the total set of 3D volumes (e.g., those associated only with some task or event timings). This yielded seven datasets total.

In order not to inundate readers with dataset details and because this is an exploratory investigation of RMT, which is not committed to any specific theory (In fact, datasets and, later, the classification tasks were chosen without investigation into any findings or publications associated with the data, both in order to reduce bias in our decisions, and because the original authors analytic/statistical approaches and intentions have very limited relevance to our whole-brain, voxelwise, multiverse predictive approach.), we refer the reader to the original publications and/or data releases when greater detail is needed and highlight only the most basic aspects of each dataset here. Dataset scan parameters and acquisition details are summarized in Table 1, and sample and subgroup sizes are summarized in Table 2.

Table 1

Included fMRI dataset details. Name, identifier for paper. Dimensions listed as M×N×P, indicate P axial slices each with dimensions M×N. TR, time of repetition (s). Volumes, number of 3D volumes per run. n_scans, total number of 4D images in dataset (number of subjects times number of runs per subject).

NameDimensionsVoxel size (mm)TRVolumesn_scans
Aging74 × 74 × 323.0 × 3.0 × 4.02.030062
Bilingual100 × 100 × 721.8 × 1.8 × 1.80.8882390
Depress112 × 112 × 252.0 × 2.0 × 5.02.510072
Learn64 × 64 × 363.0 × 3.0 × 3.02.0195432
Osteo64 × 64 × 363.4 × 3.4 × 3.02.530074
Park80 × 80 × 433.0 × 3.0 × 3.02.4149552
Attention128 × 128 × 701.5 × 1.5 × 1.53.030090

Table 2

Sizes and other details for classification task subgroups. ID, Identifier for paper; subgroup, name of subgroup used in classification task; subjects, number of subjects; task, fMRI task; ANT, ANT.53

IDSubgroupSubjectsTaskScans per subject
AgingOlder28Resting-state1
Younger34Resting-state1
BilingualBilingual59Resting-state1
Monolingual33Resting-state1
DepressDepression51Resting-state1
Control21Resting-state1
LearnTask24Learn image sketches16
Rest24Resting-state2
OsteoDuloxetine19Resting-state1
Pain37Resting-state1
Nopain20Resting-state1
ParkParkinson’s25ANT12
Control21ANT12
AttentionVigilant11Resting-state2
Non-vigilant11Resting-state2
Trait-attentive11Resting-state2
Trait-non-attentive11Resting-state2
Task-attentive11Resting-state2
Task-non-attentive11Resting-state2

2.2.

Datasets and Classification Tasks

The “aging” dataset54 included rs-fMRI data for 34 subjects with a mean age of 22 years (range: 18 to 32 years) and 28 subjects with a mean age of 70 years (range: 61 to 80 years). For our analysis, we make the binary classification task the prediction of subject age-group membership, i.e., younger v older.

The “bilinguality” dataset55,56 examined English and Spanish-speaking monolinguals and multilinguals during a prolonged resting state. The original study grouped participants into three subgroups: early versus late bilinguals versus monolingual controls.56 For our analysis, we make the binary classification task monolingual v bilingual.

The “depression” dataset57,58 includes rs-fMRI scans from non-depressed controls and mildly or moderately depressed subjects. In our analysis, the binary classification task is depress v control.

The “learning” dataset59,60 has both rs-fMRI and task fMRI scans available for all subjects, and task v rest is chosen as the binary classification task for our analysis. The task is complex and difficult to summarize here adequately, but it involved multiple phases where subjects were presented with various related abstract images, and then later tested for their memory of certain aspects of the learned images. Thus the classification task might also be considered learning v rest.

The “osteo” dataset61 includes whole-brain rs-fMRI scans of healthy, pain-free controls (nopain condition), and individuals with knee osteoarthritis. Osteoarthritic patients were treated for two weeks with either placebo (pain condition) or duloxetine (duloxetine condition). We use the three binary classification tasks nopain v duloxetine, nopain v pain, and pain v duloxetine for our analysis.

The “park” dataset62 includes subjects with non-demented Parkinson’s disease and healthy controls performed a number of repetitions of the attention network task (ANT53) during scans. The classification task for this dataset is Parkinson’s versus controls, i.e., park v ctrl.

The “attention” dataset63 is a high-resolution rs-fMRI dataset including a battery of psychological measures for each subject. Since previous studies employing RMT sometimes interpreted their findings with respect to attentional processes,31,32 we divided subjects into various high versus low attention binary classification tasks based on median splits of subsets of the metrics available in this study.

The vigilant v nonvigilant task was formed based on self-report questionnaire items involving “vigilance.”64 The trait_attend v trait_nonattend task was formed using PANAS-X65 items related to self-reported wakefulness and attention over the past weeks. Finally, the task_attend v task_nonattend classification task was constructed using the scores of the Conjunctive Continuous Performance Task.66 This is a behavioral performance task that requires the subject to quickly and selectively respond to only certain visual stimuli. Details allowing exact reproduction of these splits are available in this paper’s source code.

3.

Methods

3.1.

fMRI Preprocessing

We limited preprocessing to a simple pipeline of, in order: (1) brain extraction, (2) slice timing correction, (3) motion correction, and (4) registration to the MNI152.67,68 2mm isotropic, asymmetric template version C made available through TemplateFlow.69

Brain extraction was performed first with BET70 and then with ANTs71 to clean up residue left behind by BET. BET results were visually inspected for each image in each dataset, and BET fractional intensity threshold (-f argument) values were modified to ensure brain extraction was acceptable. Motion-correction was performed with FSL’s MCFLIRT72 and slice time correction with FSL’s slicetimer73 tool. Finally, functional images were registered directly to the MNI152 template via ANTs.71

We limited preprocessing methods to these steps because data required for other typical preprocessing steps (e.g., field intensity maps and physiological measurements) was missing, but also because all pipeline intermediates were saved to later compare the effect of increasing degrees of preprocessing on feature predictive utility (see Sec. 3.6), and it was important to limit the number of preprocessing steps to prevent excessive computational costs and comparisons.

3.2.

Feature Extraction

When dealing with large inputs, any predictive analysis must reduce that input into predictively useful features. Any method that reduces a larger set of data is potentially valid as a feature extraction method, and RMT, which summarizes matrix data by summarizing the eigenvalue distributions, may be useful for feature extraction. However, different feature extraction methods may have different advantages and disadvantages in terms of the computational complexity, amount of summary or selection involved, and in interpretability. RMT is both mathematically and computationally complex, and thus RMT-based features should be compared to alternate features that may be easier to understand or compute.

Thus for each preprocessing degree and fMRI image, features were extracted from all N non-constant brain voxels, where N may differ for each image. After this process, each scan yields an N×t matrix M, where t is the number of volumes acquired. All features extracted summarize this matrix M.

3.2.1.

Raw eigenvalues

The computing time and memory requirements for calculating the full N×N matrix M of Pearson correlation coefficients (and subsequent eigenvalues) are too large to be tractable. However, since transposition does not change the eigenvalues, we can use the transpose to efficiently compute these eigenvalues (see Appendix A). These raw eigenvalues (feature eigs in Table 3), most directly test the idea that eigenvalues alone are useful functional connectivity features.

Table 3

Feature groupings for summarization. eigs, non-RMT eigenfeatures and their combinations with other non-RMT eigenfeatures. rmt, RMT eigenfeatures and their combinations. tseries, baseline timeseries reductions.

Coarse groupingFine groupingFeature id
eigseigseigs
eigs maxeigsminmax10
eigsminmax20
eigsminmax5
eigs middleeigsmiddle10
eigsmiddle20
eigsmiddle40
eigs smootheigs + eigs_smooth
eigs + savgol
eigs_savgol
eigs_smooth
rmtrmt + eigseigs + levelvar
eigs + rigidity
eigs + rigidity + levelvar
eigs + unfolded
eigs + unfolded + levelvar
eigs + unfolded + rigidity
rmt onlylevelvar
rigidity
rigidity + levelvar
unfolded
unfolded + levelvar
unfolded + rigidity
unfolded + rigidity + levelvar
tseriesLocationT-max
T-mean
T-med
T-min
T-p05
T-p95
ScaleT-iqr
T-rng
T-rrng
T-std

However, not all eigenvalues may have predictive utility. We also examine as features the central eigenvalues (defined as the middle 10%, 20%, or 40% of the sorted eigenvalues, eigsmiddle in Table 3) and tail eigenvalues (defined as the first and last 5%, 10%, and 20% of each tail of the spectrum, eigsminmax in Table 3).

3.2.2.

RMT features

To compare the observed eigenvalue distribution to some RMT theoretical predictions, the eigenvalues must be unfolded.1,2 The details and motivation behind the unfolding procedure are documented well elsewhere.2 Practically, the unfolding process can be viewed as a smoothing and rescaling procedure where the originally observed spectrum λ1λn is smoothly mapped to an unfolded spectrum e1,,en, such that di/n1, if di=ei+1ei. For most empirical data, the true smoothing function is not known, and so must be approximated.1,2 Typically, a polynomial is chosen.45

Given empirically observed eigenvalues Λ, the spectral rigidity Δ3(L) is calculated for any positive real value L<max(Λ) as

Eq. (1)

Δ3(L)=minA,B1Lcc+L(η(λ)AλB)2c,
where η(λ) is the number of unfolded eigenvalues less than or equal to λ, ·c denotes the average with respect to all starting points c, and where A and B denote the slope and intercept, respectively, of the least squares fit of a straight line to η(λ) on [c,c+L].2 Viewing the unfolded eigenvalues as a timeseries, the spectral rigidity for a value L is the “average non-linearity” of all intervals of length L over the series.

The level number variance Σ2(L) or level variance, for short, is closely related to the spectral rigidity1 and is calculated as

Eq. (2)

Σ2(L)=η2(L,c)cη(L,c)c2,
where η(L,c) is the number of unfolded eigenvalues in [c,c+L], and where c, L, and ·c are as above.2 Viewing the unfolded eigenvalues as an irregular timeseries, the level number variance is the variation of the number of samples in all intervals of length L over the series.

As eigenfeatures, we compute the unfolded eigenvalues, and the spectral rigidity and level number variance for all L{1,2,,20}. These are unfolded, rigidity, and levelvar, respectively, in Table 3.

Trimming variants. As the unfolding procedure operates on the sorted eigenvalues and involves fitting a smooth polynomial or exponential to these values, extreme values are often omitted from fitting.48,49 The number of eigenvalues to trim is subjective, and unfortunately, with a large number of different source matrices, it is too destructive to naively use the same hard criterion (e.g., largest 10%, largest three) for all unfoldings.

We develop and test four trimming variants: no trimming, precision-based, largest, and middle, with less subjective criteria for determining trimming thresholds. The details of these trimming procedures are straightforward and can be found in Appendix B, or the source code.

Unfolding variants. Following any trimming, the eigenvalues are fit with a polynomial. Because the choice of degree has been known to dramatically impact certain analyses involving the spectral rigidity or level variance,4549 we compute the rigidity and level variance with all possible combinations of our trimming procedures and unfolding polynomial degrees of 3, 5, 7, and 9.

3.2.3.

Smoothed eigenfeatures

When computing RMT eigenfeatures, the polynomial fitting during unfolding has a smoothing effect, as do the additional averaging c operations. However, it is possible other, simpler transformations of the eigenvalues (e.g., uniform or Savitsky–Golay smoothing) might have equal or greater predictive utility. We thus also test smoothed variants of the eigenvalues as predictive features: the sorted eigenvalues are smoothed with either a uniform (moving average) filter or Savitsky–Golay, using window sizes of 3, 5, 7, and 9 to yield 8 total features (eigs_smooth and eigs_savgol in Table 3).

3.2.4.

Feature combinations

A combination of RMT and non-RMT eigenfeatures could be more predictive than either alone. We test this with a variety of eigenfeature combinations, with a focus on combinations that involve the simplest features (e.g., raw eigenvalues and unfolded eigenvalues) and then features involving additional processing (e.g., smoothed eigenvalues, level variance, and rigidity). Combined features are formed by concatenation so that if we have features f1,,fn with dimensions p1,,pn, then the combined feature is [f1;;fn] with pi dimensions. The final combinations chosen can be found in Table 3, where the “+” symbol indicates concatenation.

3.2.5.

Slicing variants

The largest, middle, or smallest eigenvalues could be most useful in characterizing any given system. For example, if the smallest eigenvalues correspond to random/noise aspects of a system but differences in the nature of the system noise most differentiate between systems, then the smallest eigenvalues may have the most predictive utility. Likewise, the L value in each RMT eigenfeature defines the degree of locality in which we summarize the spectrum: with small L-values, the spectral rigidity summarizes the non-linearity of eigenvalues that are relatively close to each other in magnitude. At large values of L, the rigidity summarizes the long-range non-linearity of the spectrum.

Since predictive utility may vary with summary locality or with the region of the original spectrum, we investigate various front, middle, and end contiguous slices of each eigenfeature in all analyses, where the size of each slice is either the first or last 5%, 10%, or 20% of the full eigenfeatures for non-middle slices, or the middle 10%, 20%, or 40% of the full eigenfeature for the middle eigenvalues.

For combined features, slicing variants are also computed, with slicing performed first on each feature to be combined. For example, if we have features f1,,fn with dimensions p1,,pn, then the sliced max 25% feature for each component feature, f^i is the last pi/4 elements of fi, and the combined sliced feature is [f^1;;f^n].

3.2.6.

Baseline features

If RMT or eigenvalue-based features fail to predictively outperform features that are simple summary statistics of the fMRI data, it is difficult to justify the greater computational and interpretational complexities of the former. We compute 10 simple statistical reductions of the voxel time series (“tseries” features in Table 3) to help assess the relative value of RMT features. That is, each statistic reduces the fMRI data along the voxel dimension, yielding a feature of t dimensions (e.g., the “mean” feature, T-mean, is the usual global mean signal).

The baseline reductions used as statistical summaries were: robust measures of location (mean, max, and min); non-robust measures of location (median, 95th percentile, and 5th percentile); non-robust measures of scale (standard deviation and the range—T-rng in Table 3); and robust measures of scale (interquartile range and difference between 95th percentile and 5th percentile—T-rrng in Table 3).

We do not take slice variants of these baseline features since the baseline features are still time series. However, we do also evaluate a number of smoothing degrees of these features, to account for noise and to be somewhat similar to the various smoothing variants for the eigenfeatures. Each baseline feature is tested with a degree of uniform smoothing, where the size of the smoothing window is either 1 (no smoothing), 2, 4, 8, or 16.

3.3.

Classifiers

We use a variety of standard machine-learning classifiers available in the Scikit-learn74 Python library to solve each classification task. We use a gradient-boosted decision tree, random forest classifier, support vector classifier (SVC) with radial basis function, and k-nearest neighbors classifiers with k equal to 3, 5, and 9 (KNN3, KNN5, and KNN9, respectively), in all cases with the default hyperparameter values. We originally also attempted to test a simple logistic regression classifier but found that this model frequently failed to converge for a number of features, even after increasing iterations significantly.

3.4.

Preprocessing Levels

Each preprocessing intermediate in the (1) brain extraction, (2) slice time correction, (3) motion correction, and (4) template registration pipeline was saved and used for entirely separate feature analyses. For example, all features were extracted from an fMRI image prepared with one of four levels or degrees of preprocessing, where preprocessing level k includes preprocessing steps 1 to k, inclusive, and starting at k=1.

3.5.

Normalization

Because of the exponential distribution of the eigenfeatures, normalization is somewhat non-trivial, and a simple standardization or min–max normalization is ineffective. We instead first apply a logarithm to all eigenfeatures and then test each classifier with un-normalized or min–max normalized versions of the log-features. Baseline features are also tested with un-normalized and min–max normalized versions, but without any log transform.

3.6.

Multiverse Analysis

We perform a “multiverse analysis75” to assess the overall predictive potential of the various eigenfeatures across all previously mentioned analytic choices. That is, for each combination of comparison task, classifier, and analytic choices, we use fivefold cross validation and use the mean area under the receiver operating characteristic curve (mAUROC) across folds as our metric to evaluate feature predictive utility.

The area under the receiver operating characteristic curve (AUROC) metric was chosen because it handles class imbalances present in the various datasets and is naturally normalized and interpretable such that mAUROC<0.5 indicates predictive performance worse than guessing, and mAUROC>0.5 indicates positive predictive utility.76 However, we also collected as performance metrics the accuracy, adjusted accuracy (proportion of samples in the entire dataset largest class minus accuracy), and the F1 score. A complete table with these additional metrics is available with the data for this study.

Each non-baseline feature is thus evaluated, for each dataset classification task and each classifier, exactly 1280 times: there are 4 preprocessing levels × 2 normalization methods × 4 trimming choices × 4 unfolding degrees × 10 slice choices. Each baseline feature is evaluated exactly 40 times for each comparison task and classifier: 4 preprocessing levels × 2 normalization methods × 5 smoothing window sizes. In total, with all the datasets and features examined in this study, this yields 2,053,920 mAUROC values to summarize, with five main analytic factors—preprocessing, normalization, trimming, unfolding/smoothing degree, and slicing—to consider.

3.7.

Software Release

Because the computation of these statistics is non-trivial and to make the unfolding procedure more transparent and reproducible, we developed and release a separate, open-source Python library empyricalRMT.77 The library allows for efficient, parallel computation of the spectral rigidity and level variance via Monte Carlo methods, and automatically ensures convergence of the statistics to user-specified tolerances. The library also makes available various other functions useful for empirical RMT analyses, such as unfolding and trimming functions, and plotting facilities for classic RMT ensembles.

4.

Results

4.1.

Overview

Summarizing this quantity of data requires some caution. Measures of location or scale, even when robust, are, for the most part, misleading and uninformative when distributions are skewed. In most cases, we find skewed mAUROC distributions, and so have chosen primarily to present our findings visually, with kernel density estimates. Because aggregating across datasets with different class imbalances can be extremely misleading, even with a metric like the mAUROC,78 we restrict such summaries to Table 4 only, and note that Table S1 in the Supplementary Material is intended as a supplement only to the distribution plots: no conclusions should be drawn about feature performance from the table values alone.

Table 4

Numerical summaries of feature mAUROCs across predictable comparisons, and all combinations of analytic choices, sorted by 95% percentile (robust max) value. Bold values indicate column “best” values, when reasonable.

FeatureMeanMin5%50%95%maxstd
eigs + eigs_smooth0.5830.1880.4080.5670.7880.9060.111
eigs + savgol0.5800.1880.4080.5630.7850.9130.110
Unfolded0.5680.1500.4000.5550.7750.9310.108
Unfolded + levelvar0.5690.1500.4000.5550.7750.9250.108
Unfolded + rigidity0.5680.1500.4010.5550.7750.9250.108
Unfolded + rigidity + levelvar0.5680.1500.4010.5550.7750.9310.108
eigs + rigidity + levelvar0.5640.2020.4010.5480.7630.9370.105
eigs + unfolded0.5640.2020.4010.5490.7630.9190.105
eigs + unfolded + levelvar0.5640.2020.4000.5480.7630.9330.105
eigs + unfolded + rigidity0.5640.2020.4010.5480.7630.9200.105
eigs + levelvar0.5630.2020.4000.5480.7580.9140.104
eigs0.5630.2020.4010.5480.7560.8750.104
eigs + rigidity0.5630.2020.4000.5480.7560.8820.104
eigs_smooth0.5690.1880.4170.5560.7550.9070.102
T-p050.5260.2490.3400.5000.7540.8560.110
eigs_savgol0.5660.1880.4100.5530.7480.9090.102
eigsminmax200.5560.2450.4110.5440.7380.9060.098
eigsminmax100.5500.2310.4090.5380.7370.9060.095
eigsminmax50.5490.2130.4000.5390.7290.8880.096
T-mean0.5570.2080.4270.5580.7170.8140.085
levelvar0.5340.0750.3670.5270.7120.8790.101
eigsmiddle400.5470.3240.4190.5330.6970.7850.084
Rigidity + levelvar0.5290.1540.3770.5230.6960.8880.094
Rigidity0.5290.1540.3780.5230.6940.8880.094
T-rrng0.5450.2810.3840.5420.6880.7590.090
eigsmiddle200.5450.3210.4190.5330.6860.7980.080
T-iqr0.4940.1750.3120.4900.6820.7590.099
eigsmiddle100.5410.3230.4110.5300.6810.7550.081
T-std0.5130.2520.3920.5000.6790.7290.086
T-rng0.5380.2790.3920.5390.6780.7790.092
T-max0.5400.2750.3850.5460.6780.7540.095
T-med0.5390.2210.3990.5410.6650.7020.077
T-p950.5240.2940.4020.5140.6530.6930.077
T-min0.5040.3210.4260.5000.5820.6980.046

In addition, a large number of interactions are possible between each analytic factor in this study. For example, it could be that trimming impacts the overall mAUROC distribution only at a certain preprocessing level and for certain feature slicing. However, as there are too many potential interactions to present, we limit our presentation to main effects.

To aid in summarizing the performances of the large number of different features, we also summarize patterns of results across two abstract groupings of similar features (“coarse” and “fine”) shown in Table 3. Thus, for example, figures depicting the coarse feature grouping “eigs” in fact depict all observed mAUROC values for the fine feature groupings of raw, max, middle, and smoothed eigenvalues, and the fine feature grouping “eigs smooth,” for example, depicts all observed mAUROC values for both the eigs_savgol and eigs_smooth features.

4.2.

Classification Tasks

4.2.1.

Non-predictable comparisons

It is important, to reduce figure clutter and complexity, to exclude a classification task or an analytic factor from summaries if there is no evidence that the classification task is solvable in general, or if, when restricting the analytic factor to a particular instance, there is no evidence of solvability. For example, if a particular trimming procedure were to render all classification tasks unsolvable, it would be better to note this and exclude the associated mAUROCs from further visualizations, rather than to have the mAUROC distributions diluted by the bad procedure.

We take as lack of evidence of solvability an mAUROC distribution that is either (1) roughly symmetric and has mean and median close to 0.5, i.e., performance appears random or (2) with median and mode <0.5. Based on the mAUROC distributions for either coarsely grouped (Fig. 1) or finely grouped features (Fig. S1 in the Supplementary Material), neither bilinguality nor depression was predictable by either baseline or eigenfeatures. Additionally, the pain v duloxetine comparison in the osteo dataset, and the trait_attend v trait_nonattend conditions (“WeeklyAttend” in Fig. 1 and Fig. S1 in the Supplementary Material) in the “attention” dataset were also not meaningfully predictable by any feature. As such, we exclude these classification tasks from further figures and discussion. We did not, however, find that any analytic factor resulted in unsolvability.

Fig. 1

AUROC distributions across gross feature groupings and comparison tasks.

JMI_10_3_036003_f001.png

4.2.2.

Predictable comparisons

As visible in Fig. 1, when coarsely summarizing features, eigenfeatures were more likely to be predictively useful than not, and except for in the task_attend v task_nonattend comparison, were also more predictively useful than the baseline features. Eigenfeatures most strongly and consistently demonstrated predictive utility in the aging dataset older v younger classification task, the osteo dataset nopain v duloxetine task, and in the attention vigilant v nonvigilant comparison.

4.3.

Largest mAUROCs

Considering the various analytic choices as tunable parameters, it makes sense to examine the largest portion of mAUROCs as an indication of the maximum predictive potential of the eigenfeatures. In this case, it is clear that eigenfeatures using RMT features almost always had the highest potential predictive utility (Fig. 2). Figure S2 in the Supplementary Material shows that this was primarily due to either the “rmt only” or “rmt + eigs” features (see Table 3). However, RMT eigenfeatures were also most likely to cause poor performance and overfitting (indicated by an mAUROC<0.5; Fig. S3 in the Supplementary Material).

Fig. 2

Distributions of largest 500 mAUROCs by coarse feature grouping.

JMI_10_3_036003_f002.png

Examining these RMT features more closely, it is clear that these features’ performance distributions differ mostly due to the unfolding procedure. That is, combined features that used the unfolded eigenvalues plus some other RMT eigenfeature tended to have visually indistinguishable mAUROC distributions to those using the unfolded eigenvalues alone (Figs. S7 and S8 in the Supplementary Material). Instead, the mAUROC distributions of these features differed mostly in the tails (Figs. S2 and S3 in the Supplementary Material).

4.4.

Effect of Preprocessing

At a coarse level of feature grouping, slice time correction followed by motion correction tended to slightly increase the predictive utility of the eigenfeatures (Fig. 3) relative to brain extraction only. Subsequent registration following these steps did not generally impact mAUROC distributions further, except in the task_attend v task_nonattend comparison, where registration reduced the predictive utility of the eigenfeatures (Fig. 3, second-last column). When examining features more finely, it is clear that preprocessing most impacts the mAUROC distribution of the largest and central eigenvalues (“eigs middle” and “eigs max” in Fig. S4 in the Supplementary Material).

Fig. 3

mAUROC distributions by preprocessing degree.

JMI_10_3_036003_f003.png

4.5.

Effect of Classifier

Within any feature grouping (coarse or fine) and within a classification task, mAUROC distributions were generally similar across classifiers (Figs. S5 and S6 in the Supplementary Material). Additionally, these figures also show that within each classification task, choice of classifier does not result in dramatic changes to the rough rank ordering of features. For example, if features are ranked on predictive utility, using the median, modal, or bulk of the mAUROC values, this rank ordering appears to remain similar across classifiers.

In the attention task_attend v task_nonattend condition, eigenfeatures were modally predictive only with the RF classifier, whereas in the aging data, an SVC was least likely to have mAUROC<0.5. Overall, however, differences in the mAUROC distributions due to the classifier were small and inconsistent.

4.6.

Normalization

There was no visible effect of feature normalization on mAUROC distributions.

4.7.

Effect of Trimming

Trimming based on numerical precision (see Sec. 3.2.2 and Appendix B) did not result in meaningfully different mAUROC distributions in any case (Fig. S7 in the Supplementary Material). However, trimming away the largest, or both largest and smallest eigenvalues, generally had a significant positive effect on the predictive quality of RMT features, most especially for RMT features involving the unfolded eigenvalues. When employing these trimming methods, these features were consistently more predictive than not (Fig. S7 in the Supplementary Material).

4.8.

Effect of Unfolding Degree

The choice of polynomial unfolding degree significantly impacted the mAUROC distributions for most classification tasks and most RMT features, and most significantly for the level variance features (Fig. S8 in the Supplementary Material). Overall, Fig. S8 in the Supplementary Material weakly suggests that either smaller (degree = 3) or larger (degree = 9) unfolding degrees tend to yield the most predictively useful RMT features. However, when restricting to the most predictive RMT features (those including the unfolded eigenvalues), it seems clear from Fig. S8 in the Supplementary Material that the largest unfolding degree of 9 produces the most favourable mAUROC distributions.

4.9.

Effect of Slicing

The best slice size and location depended complexly on the classification task and feature, and few general summaries can be made of these interactions. However, features including the full spectrum (e.g., raw eigenvalues, smoothed eigenvalues, and their combinations with RMT eigenfeatures) were slightly more predictive when using the largest portions (“max-XX,” rows in Fig. S9 in the Supplementary Material), and usually least predictive when features primarily involved the smallest or central portions.

4.10.

Choice of Summary Metric

We note briefly that most of the above findings regarding the impacts of analytic factors, and rank ordering of feature predictive utilities, are similar when using the adjusted accuracy (see Sec. 3.6), instead of the mAUROC (see Figs. S10–S18 in the Supplementary Material). However, if comparison task predictability is defined as requiring adjusted accuracies to be more positive than negative, then the learning and Parkinson’s datasets do not appear to be predictable with any feature (Fig. S10 in the Supplementary Material).

5.

Discussion

In this study, eigenfeatures inspired by RMT and derived from the eigenvalues of the full, whole-brain voxelwise fMRI correlation matrix were found to have predictive utility across a wide variety of phenomena and analytic choices. Compared to simple baseline reductions of the fMRI data, these eigenfeatures had more consistent predictive utility and a higher maximum predictive potential (Figs. 1 and 2). In addition to evidence from previous studies,29,31,32 this suggests RMT may be a useful analytic and theoretical tool for understanding functional connectivity.

However, eigenfeature mAUROC values observed in this study were highly sensitive to the overall analytic procedure, and there was no single analytic choice (e.g., choice of trimming procedure, unfolding polynomial degree, number of preprocessing steps, or feature slicing) that ensured, for any combination of feature and classification task, that all other analytic choices resulted in mAUROC values >0.5. In addition, the mean, median, and modal mAUROCs were generally close to 0.5 and adjusted mean and median accuracies also tended to be close to zero. Thus we find limited evidence that functional-connectivity-based eigenfeatures have general, “out of the box” predictive utility, with general utility likely requiring either careful tuning, or different preprocessing decisions and analytic choices than those examined here.

Nevertheless, in all datasets, there were combinations of analytic choices that resulted in cross-validated mean prediction performances well beyond mere guessing (Figs. S2 and S11 in the Supplementary Material and Table 5). Whether or not these should be considered to have practical relevance depends on one’s goals, however, we note that with small datasets of rs- or task-fMRI data, binary, subject-level classification using whole-brain features is generally challenging.

Table 5

Top three robust maximum (95th percentile) mAUROC and adjusted accuracy (acc+) values for each predictable classification task and fine feature grouping, sorted by mAUROC. A dash indicates that the fine feature grouping for that row was not in the top three, i.e., that the top three performing features differed depending on the performance metric.

Classification taskSource featuremAUROCacc+
Aging: younger versus oldereigs smooth0.8340.226
rmt + eigs0.8280.212
eigs0.8230.211
Learning: rest versus taskrmt + eigs0.6380.007
eigs0.6360.007
eigs max0.631
eigs middle0.012
Osteo: nopain versus duloxetinermt only0.8190.259
eigs max0.8120.226
rmt + eigs0.806
eigs smooth0.212
Osteo: nopain versus paintseries loc0.7620.104
tseries scale0.697
eigs smooth0.684
eigs middle0.069
eigs0.034
Parkinsons: ctrl versus parktseries scale0.6870.107
tseries loc0.6570.063
rmt only0.5900.027
TaskAttention: task_attend versus task_non-attendtseries loc0.6660.135
rmt only0.6600.121
tseries scale0.6560.124
Vigilance: vigilant versus non-vigilanteigs middle0.7330.184
eigs smooth0.726
rmt + eigs0.7190.162
eigs0.162

For example, deep learning methods improve upon guessing by 17% for autism79 or 3% to 30% for ADHD,80 16% for severe depression,81 and 23% for obsessive compulsive disorder.82 Manual feature engineering with more separable conditions (e.g., schizophrenia) can result in classification accuracies well above 90%,83 and with larger data, sophisticated custom feature extraction methods can achieve near perfect accuracies at classifying task versus rest.84 However, for functional connectivity data and with classical machine learning algorithms (such as SVC), we in general only expect large prediction accuracies (e.g., >80%) when the group functional connectivities are already strongly separated (e.g., Cohen’s d>1.0).85 In this study, the (robust) maximum improvements upon guessing are shown in Table 5 and vary from 3% to 26%.

It is somewhat surprising that the eigenfeatures examined here ever have net cross-validated predictive utility. The reduction of the functional connectivity matrix to the sorted t1 eigenvalues uses all brain voxels (including gray matter voxels) and destroys a large amount of information (radically different matrices can have identical sorted eigenvalues). The subsequent small-degree polynomial fit used in the unfolding procedure further reduces variance in the raw data, and all eigenfeatures, due to the eigenvalue sorting, are monotonically increasing curves (or approximately monotonic). All such curves could likely be fit near-perfectly with 3 to 5 parameters, i.e., the inherent dimensionality of these features is quite small. Interpreting the eigenvalues as the magnitudes of the principal components of the standardized data, this suggests that a rough summary of the magnitudes of the principal components can often be surprisingly predictive in fMRI.

We speculate that the unfolded eigenvalues may have predictive utility in part because of their smoothing and rescaling effect (see also Sec. 4.3). Figure 4, which depicts the raw eigenvalues and unfolded eigenvalues with different trimming procedures, shows how the raw eigenvalues have strongly exponential distributions, even with logarithmic axes. This is due to the magnitude of the largest eigenvalues, and the unfolded and trimmed feature distribution in Fig. 4 are far less pathologically distributed.

Fig. 4

Raw eigenvalues (first subplot) and unfolded eigenvalues (polynomial degree 9; last 3 subplots) for osteo dataset duloxetine v nopain classification task, brain extraction, and slice time correction only. Median observed mAUROC = 0.643 (range = 0.346 to 0.925).

JMI_10_3_036003_f004.png

5.1.

Limitations

As it is unfortunately typical of fMRI research,86 the number of subjects in each dataset was quite small (Table 2). With such small numbers of subjects, randomization cannot be expected to effectively control group differences, and it is possible the predictive utility of the eigenfeatures was due to capitalization of such differences. For example, eigenfeatures were generally predictive in the aging dataset, and it is quite possible for randomization failures to introduce age imbalances across classes.

Likewise, one of the other more predictable classification tasks was the duloxetine v nopain task in the osteoarthritis dataset. An active drug could introduce any number of physiological confounds relevant to fMRI,87 but we could not control for such effects due to the absence of physiological recording in most datasets.

In general, much larger fMRI datasets are needed to adequately control for and test what exactly the connectivity eigenvalues actually predict and to test if these predictions generalize to larger, different populations.

6.

Conclusion

Eigenvalue-based features inspired by RMT and extracted from fMRI functional connectivity matrices were found to have predictive utility across a wide variety of datasets and classification tasks. However, the predictive utilities were modest and highly dependent on preprocessing steps and other fitting and feature selection procedures.

Given the sensitivity to these decisions observed in this study and considering the lack of consensus regarding the preprocessing of fMRI or other complex biological signal data, RMT should probably not be currently considered a tool with “out-of-the-box” utility in these domains. Although RMT likely still has potential to yield insights in a variety of contexts, our results suggest that, absent strong evidence otherwise, it is not necessarily safe to assume that these insights will generalize broadly beyond the specific preprocessing and analytic pipelines involved.

Further research might establish specific sets of analytic choices that allow RMT to consistently extract useful information in a wide variety of contexts. However, it is also possible that each unique context might require a specific combination of choices. In this latter case, there should be strong theoretical justification for the use of RMT, and for each of the various analytic choices involved in its use: the variability of results seen in this study suggests that there may often be a set of analytic choices that provide favourable results, and which can be weakly justified by plausible, but ultimately post hoc justifications. Provided that these precautions are followed and that future studies employing RMT also carefully investigate the sensitivity of any findings to such analytic decisions, then there is likely considerable untapped potential for RMT in the analysis of fMRI.

7.

Appendix A: Correlation Eigenvalues via Transposition

Let X be a real n×p matrix with n>p. Let

Z=norm(X)=(X1X¯1||XpX¯p),
where Xi denotes column i of X. Denote the (unordered) set of eigenvalues of X as eigs(X), and let r=(p1)1. Denote the covariance matrix of X as cov(X). Then:
eigs(cov(X))=eigs(r·ZZ)=r·eigs(ZZ)=r·eigs((ZZ))=r·eigs(ZZ).

Denote the correlation matrix of X as corr(X), and let

Y=standardize(X)=((X1X¯1)/σ1||(XpX¯p)/σp).

Then

eigs(corr(X))=eigs(cov(Y))=r·eigs(YY).

8.

Appendix B: Trimming Procedures

We implement three trimming procedures: precision-based, largest, and middle trimming. The source code is the definitive reference for the procedures, but we describe the motivations belows.

In precision-based trimming, we trim away any eigenvalues that are close enough to zero to be considered a result of numerical error due to floating point representation. There are two thresholds we consider, including the one used by NumPy88 for the determination of matrix rank, and those recommended by LAPACK89 in their user guide on the error bounds for symmetric eigenproblems and related additional details. We trim each matrix eigenvalues to whichever threshold is largest for the matrix in question.

For largest trimming, we must determine a threshold in which to separate “large” from “small” eigenvalues. The eigenvalues for our data tended to grow exponentially, so we instead looked at thresholding on the logarithms. We then used k-means with k=2 on the precision-trimmed eigenvalues and took the largest cluster (which also always had the smaller mean) as the “largest” eigenvalues to trim away. “middle” trimming simply reflects the threshold found by the largest trim method, e.g., if the largest trim method removes the last n precision-trimmed eigenvalues, then we also trim the first n smallest eigenvalues remaining after precision-trimming.

We chose one-dimensional k-means partly due to efficiency and simplicity and because of the general relation to classical thresholding methods such as the Otsu method.90

Disclosures

J.L. is the founder of Time Will Tell Technologies, Inc. The authors have no relevant financial interests in the manuscript and no other potential conflicts of interest to disclose.

Data, Materials, and Code Availability

All raw fMRI data used in this study are publicly available on the OpenNeuro platform.52,91 All codes used specifically to perform preprocessing and analyses for this paper is publicly available online in the paper GitHub repository.92 This repository also includes all produced prediction data in a single table available in the file all_produced_prediction_data.json.xz available in this paper repository. Code used to compute the spectral observables and perform unfolding is publicly available in the empyricalRMT library.77

Acknowledgments

This work was supported by the Natural Science and Engineering Research Council of Canada’s Canada Research Chair Grant (Grant No. 231266) to J.L., Natural Science and Engineering Research Council of Canada Discovery Grant to J.L., a Canada Foundation for Innovation and Nova Scotia Research and Innovation Trust Infrastructure Grant (Grant No. R0176004) to J.L., a St. Francis Xavier University Research Startup Grant to J.L. (Grant No. R0168020), and a St. Francis Xavier University UCR Grant to J.L.

References

1. 

M. L. Mehta, “Random matrices,” Pure and Applied Mathematics, 142 3rd ed.Academic Press, Amsterdam; San Diego, California (2004). Google Scholar

2. 

T. Guhr, A. Müller-Groeling and H. A. Weidenmüller, “Random-matrix theories in quantum physics: common concepts,” Phys. Rep., 299 189 –425 https://doi.org/10.1016/S0370-1573(97)00088-4 PRPLCM 0370-1573 (1998). Google Scholar

3. 

T. A. Brody et al., “Random-matrix physics: spectrum and strength fluctuations,” Rev. Mod. Phys., 53 385 –479 https://doi.org/10.1103/RevModPhys.53.385 RMPHAT 0034-6861 (1981). Google Scholar

4. 

C. W. J. Beenakker, “Random-matrix theory of quantum transport,” Rev. Mod. Phys., 69 731 –808 https://doi.org/10.1103/RevModPhys.69.731 RMPHAT 0034-6861 (1997). Google Scholar

5. 

O. Bohigas, R. U. Haq and A. Pandey, “Higher-order correlations in spectra of complex systems,” Phys. Rev. Lett., 54 1645 –1648 https://doi.org/10.1103/PhysRevLett.54.1645 PRLTAO 0031-9007 (1985). Google Scholar

6. 

D. Wintgen and H. Marxer, “Level statistics of a quantized cantori system,” Phys. Rev. Lett., 60 971 –974 https://doi.org/10.1103/PhysRevLett.60.971 PRLTAO 0031-9007 (1988). Google Scholar

7. 

A. Pandey and S. Ghosh, “Skew-orthogonal polynomials and universality of energy-level correlations,” Phys. Rev. Lett., 87 024102 https://doi.org/10.1103/PhysRevLett.87.024102 PRLTAO 0031-9007 (2001). Google Scholar

8. 

M. S. Santhanam and P. K. Patra, “Statistics of atmospheric correlations,” Phys. Rev. E, 64 016102 https://doi.org/10.1103/PhysRevE.64.016102 PLEEE8 1539-3755 (2001). Google Scholar

9. 

V. Plerou et al., “Random matrix approach to cross correlations in financial data,” Phys. Rev. E, 65 066126 https://doi.org/10.1103/PhysRevE.65.066126 PLEEE8 1539-3755 (2002). Google Scholar

10. 

S. Jalan et al., “Uncovering randomness and success in society,” PLoS One, 9 e88249 https://doi.org/10.1371/journal.pone.0088249 POLNCL 1932-6203 (2014). Google Scholar

11. 

J. N. Bandyopadhyay and S. Jalan, “Universality in complex networks: random matrix analysis,” Phys. Rev. E, 76 026109 https://doi.org/10.1103/PhysRevE.76.026109 PLEEE8 1539-3755 (2007). Google Scholar

12. 

S. Sastry, N. Deo and S. Franz, “Spectral statistics of instantaneous normal modes in liquids and random matrices,” Phys. Rev. E, 64 016305 https://doi.org/10.1103/PhysRevE.64.016305 PLEEE8 1539-3755 (2001). Google Scholar

13. 

G. S. Matharoo, M. S. G. Razul and P. H. Poole, “Spectral statistics of the quenched normal modes of a network-forming molecular liquid,” J. Chem. Phys., 130 124512 https://doi.org/10.1063/1.3099605 (2009). Google Scholar

14. 

S. K. Sarkar, G. S. Matharoo and A. Pandey, “Universality in the vibrational spectra of single-component amorphous clusters,” Phys. Rev. Lett., 92 215503 https://doi.org/10.1103/PhysRevLett.92.215503 PRLTAO 0031-9007 (2004). Google Scholar

15. 

G. S. Matharoo, S. K. Sarkar and A. Pandey, “Vibrational spectra of amorphous clusters: universal aspects,” Phys. Rev. B, 72 075401 https://doi.org/10.1103/PhysRevB.72.075401 PRBMDO 1098-0121 (2005). Google Scholar

16. 

G. S. Matharoo, “Universality in the vibrational spectra of amorphous systems,” (2008). Google Scholar

17. 

P. Bhadola and N. Deo, “Targeting functional motifs of a protein family,” Phys. Rev. E, 94 042409 https://doi.org/10.1103/PhysRevE.94.042409 PLEEE8 1539-3755 (2016). Google Scholar

18. 

I. Osorio and Y.-C. Lai, “A phase-synchronization and random-matrix based approach to multichannel time-series analysis with application to epilepsy,” Chaos, 21 033108 https://doi.org/10.1063/1.3615642 (2011). Google Scholar

19. 

A. Agrawal et al., “Quantifying randomness in protein–protein interaction networks of different species: a random matrix approach,” Physica A, 404 359 –367 https://doi.org/10.1016/j.physa.2013.12.005 PHYADX 0378-4371 (2014). Google Scholar

20. 

A. Rai, A. V. Menon and S. Jalan, “Randomness and preserved patterns in cancer network,” Sci. Rep., 4 6368 https://doi.org/10.1038/srep06368 (2015). Google Scholar

21. 

S. B. Franklin et al., “Parallel analysis: a method for determining significant principal components,” J. Veg. Sci., 6 (1), 99 –106 https://doi.org/10.2307/3236261 JVESEK 1100-9233 (1995). Google Scholar

22. 

J. Veraart et al., “Denoising of diffusion MRI using random matrix theory,” NeuroImage, 142 394 –406 https://doi.org/10.1016/j.neuroimage.2016.08.016 NEIMEF 1053-8119 (2016). Google Scholar

23. 

M. Ulfarsson and V. Solo, “Dimension estimation in noisy PCA with SURE and random matrix theory,” IEEE Trans. Signal Process., 56 5804 –5816 https://doi.org/10.1109/TSP.2008.2005865 ITPRED 1053-587X (2008). Google Scholar

24. 

C. H. Martin and M. W. Mahoney, “Implicit self-regularization in deep neural networks: evidence from random matrix theory and implications for learning,” J. Mach. Learn. Res., 22 (1), 165:7479 –165:7551 https://doi.org/10.48550/arXiv.1810.01075 (2021). Google Scholar

25. 

C. H. Martin, T. Peng and M. W. Mahoney, “Predicting trends in the quality of state-of-the-art neural networks without access to training or testing data,” Nat. Commun., 12 4122 https://doi.org/10.1038/s41467-021-24025-8 NCAOBW 2041-1723 (2021). Google Scholar

26. 

G. W. Anderson, A. Guionnet, O. Zeitouni, “An introduction to random matrices,” Cambridge Studies in Advanced Mathematics, 118 Cambridge University Press, New York (2010). Google Scholar

27. 

The Oxford Handbook of Random Matrix Theory, Oxford University Press, Oxford; New York (2011). Google Scholar

28. 

J. P. Bouchaud, M. Potters, “Financial applications of random matrix theory: a short review,” The Oxford Handbook of Random Matrix Theory, Oxford University Press( (2015). Google Scholar

29. 

P. Šeba, “Random matrix analysis of human EEG data,” Phys. Rev. Lett., 91 198104 https://doi.org/10.1103/PhysRevLett.91.198104 PRLTAO 0031-9007 (2003). Google Scholar

30. 

R. Wang et al., “Spectral properties of the temporal evolution of brain network structure,” Chaos, 25 123112 https://doi.org/10.1063/1.4937451 (2015). Google Scholar

31. 

R. Wang et al., “Random matrix theory for analyzing the brain functional network in attention deficit hyperactivity disorder,” Phys. Rev. E, 94 052411 https://doi.org/10.1103/PhysRevE.94.052411 PLEEE8 1539-3755 (2016). Google Scholar

32. 

G. S. Matharoo and J. A. Hashmi, “Spontaneous back-pain alters randomness in functional connections in large scale brain networks: a random matrix perspective,” Physica A, 541 123321 https://doi.org/10.1016/j.physa.2019.123321 PHYADX 0378-4371 (2020). Google Scholar

33. 

A. M. Bastos and J.-M. Schoffelen, “A tutorial review of functional connectivity analysis methods and their interpretational pitfalls,” Front Syst Neurosci, 9 175 https://doi.org/10.3389/fnsys.2015.00175 (2016). Google Scholar

34. 

M. P. van den Heuvel and H. E. Hulshoff Pol, “Exploring the brain network: a review on resting-state fMRI functional connectivity,” Eur. Neuropsychopharmacol., 20 519 –534 https://doi.org/10.1016/j.euroneuro.2010.03.008 EURNE8 0924-977X (2010). Google Scholar

35. 

S. Noble, D. Scheinost and R. T. Constable, “A decade of test-retest reliability of functional connectivity: a systematic review and meta-analysis,” NeuroImage, 203 116157 https://doi.org/10.1016/j.neuroimage.2019.116157 NEIMEF 1053-8119 (2019). Google Scholar

36. 

D. J. Lurie et al., “Questions and controversies in the study of time-varying functional connectivity in resting fMRI,” Network Neurosci., 4 30 –69 https://doi.org/10.1162/netn_a_00116 (2020). Google Scholar

37. 

M. Voultsidou, S. Dodel and J. M. Herrmann, “Feature evaluation in fMRI data using random matrix theory,” Comput. Vis. Sci., 10 99 –105 https://doi.org/10.1007/s00791-006-0037-6 (2007). Google Scholar

38. 

A. A. Vergani, S. Martinelli and E. Binaghi, “Resting state fMRI analysis using unsupervised learning algorithms,” Comput. Methods Biomech. Biomed. Eng.: Imaging Vis., 8 252 –265 https://doi.org/10.1080/21681163.2019.1636413 (2019). Google Scholar

39. 

L. Gu et al., “Random matrix theory for analysing the brain functional network in lower limb motor imagery,” in 42nd Annu. Int. Conf. of the IEEE Eng. in Med. & Biol. Soc. (EMBC), 506 –509 (2020). https://doi.org/10.1109/EMBC44109.2020.9176442 Google Scholar

40. 

S. E. Joel et al., “On the relationship between seed-based and ICA-based measures of functional connectivity,” Magn. Reson. Med., 66 (3), 644 –657 https://doi.org/10.1002/mrm.22818 MRMEEN 0740-3194 (2011). Google Scholar

41. 

D. V. Smith et al., “Characterizing individual differences in functional connectivity using dual-regression and seed-based approaches,” NeuroImage, 95 1 –12 https://doi.org/10.1016/j.neuroimage.2014.03.042 NEIMEF 1053-8119 (2014). Google Scholar

42. 

A. Lo et al., “Why significant variables aren’t automatically good predictors,” Proc. Natl. Acad. Sci. U. S. A., 112 13892 –13897 https://doi.org/10.1073/pnas.1518285112 (2015). Google Scholar

43. 

V. Amrhein, F. Korner-Nievergelt and T. Roth, “The earth is flat (p>0.05): significance thresholds and the crisis of unreplicable research,” PeerJ, 5 e3544 https://doi.org/10.7717/peerj.3544 (2017). Google Scholar

44. 

T. Yarkoni and J. Westfall, “Choosing prediction over explanation in psychology: lessons from machine learning,” Perspect. Psychol. Sci., 12 1100 –1122 https://doi.org/10.1177/1745691617693393 (2017). Google Scholar

45. 

A. A. Abul-Magd and A. Y. Abul-Magd, “Unfolding of the spectrum for chaotic and mixed systems,” Physica A, 396 185 –194 https://doi.org/10.1016/j.physa.2013.11.012 PHYADX 0378-4371 (2014). Google Scholar

46. 

S. M. Abuelenin, “On the spectral unfolding of chaotic and mixed systems,” Physica A, 492 564 –570 https://doi.org/10.1016/j.physa.2017.08.158 PHYADX 0378-4371 (2018). Google Scholar

47. 

R. Fossion, G. T. Vargas and J. C. L. Vieyra, “Random-matrix spectra as a time series,” Phys. Rev. E, 88 060902 https://doi.org/10.1103/PhysRevE.88.060902 PLEEE8 1539-3755 (2013). Google Scholar

48. 

S. M. Abuelenin and A. Y. Abul-Magd, “Effect of unfolding on the spectral statistics of adjacency matrices of complex networks,” Procedia Comput. Sci., 12 69 –74 https://doi.org/10.1016/j.procs.2012.09.031 (2012). Google Scholar

49. 

I. O. Morales et al., “Improved unfolding by detrending of statistical fluctuations in quantum spectra,” Phys. Rev. E, 84 016203 https://doi.org/10.1103/PhysRevE.84.016203 PLEEE8 1539-3755 (2011). Google Scholar

50. 

D. B. Parker and Q. R. Razlighi, “The benefit of slice timing correction in common fMRI preprocessing pipelines,” Front. Neurosci., 13 821 https://doi.org/10.3389/fnins.2019.00821 1662-453X (2019). Google Scholar

51. 

J. P. Simmons, L. D. Nelson and U. Simonsohn, “False-positive psychology: undisclosed flexibility in data collection and analysis allows presenting anything as significant,” Psychol. Sci., 22 1359 –1366 https://doi.org/10.1177/0956797611417632 1467-9280 (2011). Google Scholar

52. 

C. J. Markiewicz et al., “The OpenNeuro resource for sharing of neuroscience data,” eLife, 10 e71774 https://doi.org/10.7554/eLife.71774 (2021). Google Scholar

53. 

J. Fan et al., “The activation of attentional networks,” NeuroImage, 26 471 –479 https://doi.org/10.1016/j.neuroimage.2005.02.004 NEIMEF 1053-8119 (2005). Google Scholar

54. 

C. N. Wahlheim et al., “Intrinsic functional connectivity in the default mode network predicts mnemonic discrimination: a connectome-based modeling approach,” Hippocampus, 32 (1), 21 –37 https://doi.org/10.1002/hipo.23393 (2022). Google Scholar

55. 

C. E. Gold, “Exploring the resting state neural activity of monolinguals and late and early bilinguals,” (2018). Google Scholar

56. 

C. E. Gold, “Exploring the resting state neural activity of monolinguals and late and early bilinguals,” Brigham Young University, (2018). Google Scholar

57. 

D. D. Bezmaternykh et al., “Resting state with closed eyes for patients with depression and healthy participants,” OpenNeuro, https://doi.org/10.18112/openneuro.ds002748.v1.0.5 (2021). Google Scholar

58. 

D. D. Bezmaternykh et al., “Brain networks connectivity in mild to moderate depression: resting state fMRI study with implications to nonpharmacological treatment,” Neural Plast., 2021 1 –15 https://doi.org/10.1155/2021/8846097 (2021). Google Scholar

59. 

A. C. Schapiro et al., “Human hippocampal replay during rest prioritizes weakly learned information and predicts memory performance,” Nat. Commun., 9 3920 https://doi.org/10.1038/s41467-018-06213-1 NCAOBW 2041-1723 (2018). Google Scholar

60. 

A. C. Schapiro et al., “Human hippocampal replay during rest prioritizes weakly learned information and predicts memory performance,” OpenNeuro, (2020). https://doi.org/10.18112/openneuro.ds001454.v1.3.1 Google Scholar

61. 

P. Tétreault et al., “Brain connectivity predicts placebo response across chronic pain clinical trials,” PLoS Biol., 14 e1002570 https://doi.org/10.1371/journal.pbio.1002570 (2016). Google Scholar

62. 

T. M. Madhyastha et al., “Dynamic connectivity at rest predicts attention task performance,” Brain Connect, 5 45 –59 https://doi.org/10.1089/brain.2014.0248 (2015). Google Scholar

63. 

K. J. Gorgolewski et al., “A high resolution 7-Tesla resting-state fMRI test-retest dataset with cognitive and physiological measures,” Sci. Data, 2 140054 https://doi.org/10.1038/sdata.2014.54 (2015). Google Scholar

64. 

K. J. Gorgolewski et al., “A correspondence between Individual differences in the brain’s intrinsic functional architecture and the content and form of self-generated thoughts,” PLoS One, 9 e97176 https://doi.org/10.1371/journal.pone.0097176 POLNCL 1932-6203 (2014). Google Scholar

65. 

L. A. Clark and D. Watson, “The PANAS-X: manual for the positive and negative affect schedule - expanded form,” https://iro.uiowa.edu/esploro/outputs/other/9983557488402771 (1994). Google Scholar

66. 

L. Shalev et al., “Conjunctive Continuous Performance Task (CCPT)—a pure measure of sustained attention,” Neuropsychologia, 49 2584 –2591 https://doi.org/10.1016/j.neuropsychologia.2011.05.006 NUPSA6 0028-3932 (2011). Google Scholar

67. 

V. Fonov et al., “Unbiased average age-appropriate atlases for pediatric studies,” NeuroImage, 54 313 –327 https://doi.org/10.1016/j.neuroimage.2010.07.033 NEIMEF 1053-8119 (2011). Google Scholar

68. 

V. Fonov et al., “Unbiased nonlinear average age-appropriate brain templates from birth to adulthood,” NeuroImage, 47 S102 https://doi.org/10.1016/S1053-8119(09)70884-5 NEIMEF 1053-8119 (2009). Google Scholar

69. 

R. Ciric et al., “TemplateFlow: FAIR-sharing of multi-scale, multi-species brain models,” Neuroscience, 19 (12), 12 https://doi.org/10.1038/s41592-022-01681-2 (2022). Google Scholar

70. 

S. M. Smith, “Fast robust automated brain extraction,” Hum. Brain Mapp., 17 143 –155 https://doi.org/10.1002/hbm.10062 (2002). Google Scholar

71. 

B. B. Avants et al., “A reproducible evaluation of ANTs similarity metric performance in brain image registration,” NeuroImage, 54 2033 –2044 https://doi.org/10.1016/j.neuroimage.2010.09.025 NEIMEF 1053-8119 (2011). Google Scholar

72. 

M. Jenkinson et al., “Improved optimization for the robust and accurate linear registration and motion correction of brain images,” NeuroImage, 17 825 –841 https://doi.org/10.1006/nimg.2002.1132 NEIMEF 1053-8119 (2002). Google Scholar

73. 

M. Jenkinson et al., “FSL,” NeuroImage, 62 782 –790 https://doi.org/10.1016/j.neuroimage.2011.09.015 NEIMEF 1053-8119 (2012). Google Scholar

74. 

F. Pedregosa et al., “Scikit-learn: machine learning in Python,” J. Mach. Learn. Res., 12 (85), 2825 –2830 (2011). Google Scholar

75. 

S. Steegen et al., “Increasing transparency through a multiverse analysis,” Perspect. Psychol. Sci., 11 702 –712 https://doi.org/10.1177/1745691616658637 (2016). Google Scholar

76. 

J. N. Mandrekar, “Receiver operating characteristic curve in diagnostic test assessment,” J. Thorac. Oncol., 5 1315 –1316 https://doi.org/10.1097/JTO.0b013e3181ec173d (2010). Google Scholar

77. 

DM-Berger, “Stfxecutables/empyricalRMT: V1.1.1,” (2022). Google Scholar

78. 

J. Brabec et al., “On model evaluation under non-constant class imbalance,” in Computational Science – ICCS 2020, 74 –87 (2020). Google Scholar

79. 

M. Bengs, N. Gessert and A. Schlaefer, “4D spatio-temporal deep learning with 4D fMRI data for autism spectrum disorder classification,” (2020). Google Scholar

80. 

A. Riaz et al., “DeepFMRI: end-to-end deep learning for functional connectivity and classification of ADHD using fMRI,” J. Neurosci. Methods, 335 108506 https://doi.org/10.1016/j.jneumeth.2019.108506 JNMEDT 0165-0270 (2020). Google Scholar

81. 

R. Ramasubbu et al., “Accuracy of automated classification of major depressive disorder as a function of symptom severity,” Neuroimage Clin., 12 320 –331 https://doi.org/10.1016/j.nicl.2016.07.012 (2016). Google Scholar

82. 

Y. Takagi et al., “A neural marker of obsessive-compulsive disorder from whole-brain functional connectivity,” Sci. Rep., 7 7538 https://doi.org/10.1038/s41598-017-07792-7 (2017). Google Scholar

83. 

W. Du et al., “High classification accuracy for schizophrenia with rest and task fMRI data,” Front. Hum. Neurosci., 6 145 https://doi.org/10.3389/fnhum.2012.00145 (2012). Google Scholar

84. 

S. Zhang et al., “Characterizing and differentiating task-based and resting state fMRI signals via two-stage sparse representations,” Brain Imaging Behav., 10 21 –32 https://doi.org/10.1007/s11682-015-9359-7 (2016). Google Scholar

85. 

C. Dansereau, “Statistical power and prediction accuracy in multisite resting-state fMRI connectivity,” NeuroImage, 149 220 –232 https://doi.org/10.1016/j.neuroimage.2017.01.072 NEIMEF 1053-8119 (2017). Google Scholar

86. 

B. O. Turner et al., “Small sample sizes reduce the replicability of task-based fMRI studies,” Commun. Biol., 1 1 –10 https://doi.org/10.1038/s42003-018-0073-z (2018). Google Scholar

87. 

K. Murphy, R. M. Birn and P. A. Bandettini, “Resting-state FMRI confounds and cleanup,” NeuroImage, 80 349 –359 https://doi.org/10.1016/j.neuroimage.2013.04.001 NEIMEF 1053-8119 (2013). Google Scholar

88. 

C. R. Harris et al., “Array programming with NumPy,” Nature, 585 357 –362 https://doi.org/10.1038/s41586-020-2649-2 (2020). Google Scholar

89. 

E. Anderson et al., LAPACK Users’ Guide, 3rd ed.Society for Industrial and Applied Mathematics, Philadelphia, PA (1999). Google Scholar

90. 

D. Liu and J. Yu, “Otsu method and K-means,” in Ninth Int. Conf. Hybrid Intell. Syst., 344 –349 (2009). Google Scholar

91. 

R. A. Poldrack and K. J. Gorgolewski, “OpenfMRI: open sharing of task fMRI data,” NeuroImage, 144 259 –261 https://doi.org/10.1016/j.neuroimage.2015.05.073 NEIMEF 1053-8119 (2017). Google Scholar

92. 

DM-Berger, “DM-Berger/random-matrix-fMRI: paper release,” (2022). Google Scholar

Biography

Derek Berger received his BSc degree with a specialization in psychology from the University of Toronto, Toronto, Ontario, Canada, in 2015. He is currently employed as a research assistant at St. Francis Xavier University, Antigonish, Nova Scotia, Canada. His research interests include reproducibility in machine- and deep-learning research for biomedical imaging and signal data, and the development of software for the automation of machine learning pipelines for predictive modeling.

Gurpreet S. Matharoo received his BSc (Hons.) and MSc degrees in physics from the University of Delhi, Delhi, India, in 1998 and 2000, respectively, and his PhD in computational and statistical physics from Jawaharlal Nehru University, New Delhi, India, in 2006. He is a research consultant at ACENET. His current research focuses on applications of advanced data analytics and high-performance computing tools in health, biomedical, and biophysical systems.

Jacob Levman received his BASc degree in computer engineering from the University of Waterloo, his MASc degree in electrical and computer engineering from Ryerson University, and his PhD in medical biophysics from the University of Toronto. He completed postdoctoral fellowships in imaging research at Sunnybrook Hospital, biomedical engineering at the University of Oxford, and neuroscience at Boston Children’s Hospital. As an associate professor of computer science at St. Francis Xavier University and a visiting scientist at Massachusetts General Hospital of Harvard Medical School, his research includes computational technologies, medical imaging, patient monitoring, and neuroscience.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 International License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Derek Berger, Gurpreet S. Matharoo, and Jacob Levman "Random matrix theory tools for the predictive analysis of functional magnetic resonance imaging examinations," Journal of Medical Imaging 10(3), 036003 (14 June 2023). https://doi.org/10.1117/1.JMI.10.3.036003
Received: 15 November 2022; Accepted: 22 May 2023; Published: 14 June 2023
Advertisement
Advertisement
KEYWORDS
Matrices

Functional magnetic resonance imaging

Analytics

Feature extraction

Binary data

Brain

Voxels

Back to Top