Functional Factor Modeling of Brain Connectivity

There exist a number of imaging modalities that allow researchers to measure the physiological changes that accompany neuronal activation. Among the most widely used techniques is functional magnetic resonance imaging (fMRI). FMRI is based on the haemodynamic response wherein blood delivers oxygen to active brain regions at a greater rate than to inactive regions. Blood-oxygen-level-dependent (BOLD) imaging uses magnetic fields to measure relative oxygenation levels across the brain, allowing researchers to use BOLD signal as a proxy for neuronal activation. During an fMRI experiment, a subject is placed in a scanner that collects a temporal sequence of three-dimensional brain images. Each image is partitioned into a grid of three-dimensional volume elements, called voxels, and contains BOLD measurements at each voxel (Lindquist, 2008). The resulting data are both large and complex. The data from a single session typically contains hundreds of brain images, each with BOLD observations at more than 100,000 voxels; large fMRI datasets include hundreds of such scans. Moreover, these data exhibit spatiotemporal dependencies that complicate analysis.

A common goal of fMRI analyses is to characterize functional connectivity, or the statistical dependencies among remote neurophysiological events (Friston, 2011). This is of particular interest in resting-state fMRI, wherein subjects do not perform an explicit task while in the scanner. In this work, we study functional connectivity in resting-state data from the PIOP1 dataset of the Amsterdam Open MRI Collection (AOMIC) (Snoek, van der Miesen, Beemsterboer, van der Leij, Eigenhuis, and Steven Scholte, 2021), which contain six-minute scans collected from 210 subjects at rest. To do so, we develop a new approach that integrates two fields of statistics: factor analysis and functional data analysis.

Techniques used to study resting-state functional connectivity fall in two categories: seed methods and whole-brain methods. To conduct a seed-based connectivity analysis, researchers first select a seed voxel (or region), then correlate the time series of that seed with the time series of all other voxels (or regions), resulting in a functional connectivity map (Biswal, Zerrin Yetkin, Haughton, and Hyde, 1995). Although readily interpretable, these functional connectivity maps only describe dependencies of the selected seed region (Van Den Heuvel and Pol, 2010). Principal component analysis (PCA) is a rudimentary whole-brain approach to functional connectivity that finds the collection of orthogonal spatial maps (a.k.a., principal components) that maximally explains the variation – both between- and within-voxel – in a set of brain images (Friston et al., 1993). The orthogonality constraint, however, means that the low-order subspace is expressed by dense spatial maps containing many high-magnitude weights that complicate interpretation. Independent component analysis (ICA; Beckmann and Smith, 2004), perhaps the most popular whole-brain approach, offers a more parsimonious representation of this subspace. After computing the principal spatial components and associated temporal components from an fMRI scan, ICA rotates these components until the spatial set is maximally independent. Although the resulting spatial maps are typically more parsimonious than those of PCA, ICA may fragment broad areas of activation into multiple maps with highly correlated time courses in its pursuit to maximize spatial independence (Friston, 1998). This means ICA findings are not robust to dimension overspecification.

To address the shortcomings of the above approaches, we develop a method based on exploratory factor analysis. The goal of factor analysis is to describe the covariance relationships between many observed variables using a few unobserved random quantities. Classical factor analyses use the orthogonal factor model (OFM), which posits that $M$ zero-mean observed quantities $\underline{X}\in\mathbb{R}^{M}$ having covariance $\boldsymbol{C}\in\mathbb{R}^{M\times M}$ are the sum of a common component, given by linear combinations of $K<M$ unobserved uncorrelated factors $\underline{F}\in\mathbb{R}^{K}$ with zero-mean and unit variance, and an idiosyncratic component, given by a vector of uncorrelated errors $\underline{\epsilon}\in\mathbb{R}^{M}$ with zero-mean. The loading matrix $\boldsymbol{L}\in\mathbb{R}^{M\times K}$ specifies the coefficients of these linear combinations in the common component:

If $\boldsymbol{D}\in\mathbb{R}^{M\times M}$ denotes the diagonal covariance matrix of $\underline{\epsilon}$, then $\boldsymbol{C}$ decomposes as

The OFM is identifiable only up to an orthogonal rotation since, for $\underline{F}^{*}=\boldsymbol{R}^{T}\underline{F}$ and $\boldsymbol{L}^{*}=\boldsymbol{L}\boldsymbol{R}$, where $\boldsymbol{R}\in\mathbb{R}^{K\times K}$ is a rotation matrix, we have the equivalences

Practitioners address this rotation problem by choosing from a suite of factor rotation algorithms, like quartimax (Neuhaus and Wrigley, 1954) and varimax (Kaiser, 1958) rotation, that bring the loading matrix closer to a so-called simple structure, where in (i) each variable has a high loading on a single factor but near-zero loadings on other factors, and (ii) each factor is described by only a few variables with high loadings on that factor while other variables have near-zero loadings on that factor. This structure provides a parsimonious representation of the common component and one can opt for a rotation procedure that is robust to dimension overspecification.

To grasp the value of factor analysis in functional connectivity, consider a naive application wherein an OFM is fit to a vector of voxel-wise BOLD observations. Used in conjunction with factor rotation, the OFM would provide $K$ simple spatial maps amenable to interpretation. This sparse and low-dimensional representation of the entire brain makes factor analysis an attractive alternative to seed and ICA methods which, respectively, are not whole-brain and lack such interpretable results. However, there is a problem; the OFM crucially assumes that its errors are uncorrelated, a condition violated in this application given that errors for nearby voxels will be correlated. Such autocorrelation may dissipate as $K$ grows, but this has two undesirable effects. First, the parsimony of the low-order model is lost. Second, a high-order common component will capture variation at small scales – possibly even between adjacent voxels – an outcome at odds with the purpose of functional connectivity studies: to study correlations between distant brain regions. In a more appropriate factor model, the common component would be low-rank and capture only large-scale variation while its error term would permit variation at a smaller scale.

To develop this better suited factor model, it helps to view fMRI data in the functional data analysis (FDA) framework. Data are said to be functional if we can naturally assume that they arise from some smooth curve. FDA is then the statistical analysis of sample curves in such data (Ramsay and Silverman, 2005; Wang et al., 2016; Kokoszka and Reimherr, 2017). Since it is reasonable to view BOLD observations of a brain image as noisy realizations from some smooth function over a three-dimensional domain, fMRI data is an example of multidimensional functional data. Though not widely adopted by the neuroimaging community, fMRI has been studied using FDA tools such as functional principal component analysis (FPCA; Viviani et al., 2005), functional regression (Reiss et al., 2015, 2017), and functional graphical models (Li and Solea, 2018). To adapt factor analytic techniques to fMRI data, we can extend classical factor methods to this functional setting.

Early extensions of the OFM come from studies of large-dimensional econometric panel data, in which $M$ time series (e.g., one time series per asset) are observed at the same $n$ points in time. When $M$ is large, traditionally reliable time series tools (e.g., vector autoregressive models) require estimation of too many parameters. Factor models are a practical alternative as they provide more parsimonious parameterizations. These factor models should, however, possess two characteristics to accommodate properties inherent to economic time series. First, they ought to be dynamic as problems involving time series are dynamic in nature. Second, they should permit correlated errors as uncorrelatedness is often unrealistic with large $M$. Sargent and Sims (1977) and Geweke (1977) innovated on the OFM by making latent factors dynamic in their dynamic factor model (DFM). However, these early DFMs were exact in that they assumed uncorrelated errors. Factor models did not permit such error correlation until Chamberlain and Rothschild (1983) introduced their approximate factor model (AFM), which required only that the largest eigenvalue of the error covariance be bounded. Although this AFM was static, much subsequent study focused on factor models that were both approximate and dynamic (Forni et al., 2000; Stock and Watson, 2002; Bai, 2003).

Factor models are well-studied in large dimensions, but the literature on infinite-dimensional factor analysis is scarce. Hays et al. (2012) were the first to merge ideas from FDA with factor analysis in their dynamic functional factor model (DFFM) which they used to forecast yield curves. Liebl (2013) proposed a different DFFM which, unlike the model of Hays et al. (2012), placed no a priori constraints on the latent process of factor scores. Kowal et al. (2017) presented a Bayesian DFFM for modeling multivariate dependent functional data. None of these early functional factor models permit any degree of error dependence, a structure that certainly arises when representing potentially infinite-rank observations with low-rank objects. Otto and Salish (2022) were the first to address this issue using their approximate DFFM wherein the error covariance is left unrestricted. However, as this is a model for completely-observed functional data, one must either observe these data at infinite resolution (a practical impossibility) or obtain functional representations of discretely-observed sample curves via spline smoothing (Ramsay and Silverman, 2005) prior to model estimation.

Though spline smoothing (and other forms of pre-smoothing) may be an effective first step in many analyses of functional data, it is not in factor analysis. Analyses that begin by smoothing functional data assume that a functional observation $X(s)$ is the sum of two uncorrelated components: a smooth signal $Y(s)$ and a white noise term $\epsilon(s)$ whose covariance kernel is zero off the diagonal (or some infinitesimally narrow band). When observed at discrete points,

the smoothed curve $\tilde{X}(s)$ is defined as

When $\epsilon(s)$ is white noise, $\tilde{X}(s)$ is a good approximation to $Y(s)$. If, instead, the signal $Y(s)$ and the error $\epsilon(s)$ are defined as in our theorized factor model – the former capturing only large-scale variation and the latter permitting small-scale variation – then minimizing the “uncorrelated” objective of (1) yields a smoothed curve $\tilde{X}(s)$ contaminated by short-range variation of the error term, making it a poor approximation to $Y(s)$. Using $\tilde{X}(s)$ as a stand-in for $Y(s)$ in subsequent analyses is ill-advised. Thus, a factor model for functional data should not presume completely-observed sample curves, but rather deal directly with discretely-observed data.

As it turns out, an orthogonal factor structure indeed underlies discretely-observed functional data. Hörmann and Jammoul (2020) formalized this observation, noting, however, that the assumption of uncorrelated errors is rarely justified. They suggest that a more realistic factor model would permit error correlations that “taper to zero with increasing lag.” But this poses an identifiability question: is it possible to distinguish such error dependence from that residing in the common component? Descary and Panaretos (2019) show that this is, indeed, possible if one assumes the covariance of the common term is smooth with low rank and that of the error term is banded with potentially infinite rank. These assumptions give rise to a distinctly functional re-characterization of the common and idiosyncratic components of the OFM: the former becomes “global”, capturing dependencies between distant regions of the domain, while the latter becomes “local”, describing dependencies between nearby regions. As this is the precise perspective put forth by our theorized factor model for functional connectivity, we let the ideas of Descary and Panaretos (2019) guide the development of this model.

In this paper, we develop a factor analytic approach for discretely-observed multidimensional functional data that is well-suited to the study of functional connectivity in fMRI. Our methodology is inspired by the framework of Descary and Panaretos (2019), which we build upon in three key respects. First, we extend their conditions for identifying the global and local components of a full covariance to multidimensional functional data, and show that these conditions become less restrictive as the dimensionality grows. Second, we improve upon their global covariance estimator by incorporating a roughness penalty. Third, we develop a post-processing procedure – inspired by methods from classical factor analysis – that facilitates interpretation of global estimates. Using this new approach, we uncover a rich covariance structure in the 210 resting-state scans from the AOMIC-PIOP1 dataset.

The rest of this paper proceeds as follows. In Section 2, we introduce functional factor models for both completely- and discretely-observed functional data, and present identifiability conditions for both. In Section 3, we develop estimation methods for the discretely-observed functional factor model. We then take a brief detour, in Section 4, to review an ICA-based approach used throughout the neuroimaging community. In Section 5, we use simulations to compare our estimator to that of this ICA-based approach and several other comparators. Section 6 contains an application of our approach to the AOMIC resting-state data which we contrast with an application of the ICA-based estimator. We conclude, in Section 7, with some final remarks and future work.

This section considers two data observation paradigms: (i) data observed completely (e.g., a brain image with infinite spatial resolution, making it a function), and (ii) data observed discretely on a grid (e.g., a brain image with finite spatial resolution, making it a tensor). Given that our methodology forgos the usual data smoothing step in many analyses of functional data, the completely-observed paradigm is not practically possible. We instead use this first paradigm as an intermediate step in the development of our theory for the second. Unless otherwise specified, let $a$ or $A$ denote scalars, $\underline{a}$ or $\underline{A}$ denote vectors, $\boldsymbol{A}$ denote matrices, $\mathcal{A}$ denote tensors, and $\mathscr{A}$ denote operators.

In this section, we formulate estimation in the discrete observation paradigm. The methodology is split in two phases: initial loading estimation and post-processing. In the first phase, we define an estimator for $\mathcal{G}$, from which we derive estimates for the number of factors $K$ and the loading tensors $\mathcal{L}_{k}$. To estimate the $\mathcal{L}_{k}$ through $\mathcal{G}$, we must assume that these $\mathcal{L}_{k}$ are orthogonal. This, of course, is merely a mathematical convenience; there is no reason to suspect that a collection of brain maps describing connectivity structures is orthogonal. We depart from this assumption during post-processing, wherein the initial loading estimates are brought into a more interpretable form via rotation and shrinkage. This section emphasizes practical implementation of our method; Section 3.3 of Supplement A develops theory for these estimators and provides derivations of their asymptotic properties.

In Sections 5 and 6, we compare our approach to, among others, one based on ICA, the most popular whole-brain connectivity method. In doing so, we rely on the MELODIC function (Beckmann and Smith, 2004) of the FMRIB Software Library (FSL; Jenkinson et al., 2012), arguably ICA’s most widely used implementation for connectivity analyses. This section details the MELODIC model and its estimation.

FSL’s MELODIC function estimates the probabilistic ICA model which assumes that observations for each voxel are generated from a set of statistically independent non-Gaussian sources via a linear mixing process corrupted by additive Gaussian noise:

Here, $\underline{x}_{m}$ denotes the $n$-dimensional vector of observations at voxel $m$, $\underline{s}_{m}$ denotes the $K$-dimensional vector of non-Gaussian sources, $\boldsymbol{A}$ denotes the ($n\times K$)-dimensional mixing matrix, $\underline{\mu}$ denotes the mean of the observations $\underline{x}_{m}$, and $\underline{\eta}_{m}$ denotes Gaussian noise. The model assumes there are fewer sources than observations (i.e., $K<n$), and that the noise covariance is voxel-dependent (i.e., $\underline{\eta}_{m}\sim N(\underline{0},\sigma^{2}\boldsymbol{\Sigma}_{m})$).

The goal of estimation is to find the $(K\times n)$-dimensional unmixing matrix $\boldsymbol{W}$ such that $\hat{\underline{s}}_{m}=\boldsymbol{W}\underline{x}_{m}$ is a good approximation to $\underline{s}_{m}$ for each voxel. Prior to invoking MELODIC, it is common to spatially smooth the data. We do so using the Gaussian filter of FSL’s FSLMATHS utility (Jenkinson et al., 2012), tuning the kernel width via cross-validation. Before estimation, MELODIC prepares the multi-session data by (i) temporally concatenating the preprocessed slices to form observation vectors $\underline{x}_{m}$ for each voxel, (ii) temporally whitening each $\underline{x}_{m}$ so that the noise covariance is isotropic at each voxel location, then (iii) normalizing each $\underline{x}_{m}$ to zero mean and unit variance. By default, the utility also reduces the dimension of the concatenated data via MELODIC Incremental Group PCA (MIGP; Smith et al., 2014); however, we switch off this feature as Sections 5 and 6 consider datasets of sufficiently reduced sizes. After this preprocessing, one can show that the maximum likelihood (ML) mixing matrix estimate is

where $\boldsymbol{U}_{K}$ and $\boldsymbol{\Lambda}_{K}$ contain the first $K$ left singluar vectors and singular values of $\boldsymbol{X}$, the $n\times M$ matrix of preprocessed data, and $\boldsymbol{Q}$ is some $K\times K$ rotation matrix. From $\hat{\boldsymbol{A}}$, the ML source estimates are obtained via generalized least squares,

where $\tilde{\underline{x}}_{m}=(\boldsymbol{\Lambda}_{K}-\sigma^{2}\boldsymbol{I}_{K})^{-1/2}\boldsymbol{U}_{K}^{T}\underline{x}_{m}$. Thus, to fix $\hat{\boldsymbol{A}}$, one chooses the $\boldsymbol{Q}$ that maximizes the non-Gaussianity of the source estimates $\hat{\underline{s}}_{m}$.

For a random variable $s$, one measure of non-Gaussianity is negentropy,

where $H(y)=-\mathbb{E}_{y}[\log y]$ is the entropy of a random variable $y$ and $s_{\text{gauss}}$ is a Gaussian having the same variance (or covariance if $s$ is a random vector) as $s$. Since a Gaussian random variable has the largest entropy among all random variables of equal variance, negentropy is a non-negative function that equals zero if and only if $s$ is Gaussian. As estimation of negentropy is difficult, it is often approximated using the contrast function

where $\nu$ is a standard normal and $G$ is a general non-quadratic function. To estimate the $k$th source $s_{mk}$ of the $m$th voxel MELODIC uses a fixed-point iteration scheme (Hyvärinen and Oja, 2000) to choose the $k$th row $\underline{q}_{k}^{T}$ of $\boldsymbol{Q}$ so that $J(\hat{s}_{mk})=J(\underline{q}_{k}^{T}\tilde{\underline{x}}_{m})$ is maximized for some domain-specific choice of $G$. By default, MELODIC sets $G$ equal to the pow3 function. To estimate $K$ sources, the procedure is repeated $K$ times under the constraint that the vectors $\underline{q}_{k}$ are mutually orthogonal.

To assess the efficacy of our methodology, we conduct three simulation studies. In the first, we compare the accuracy of the post-processed estimator of (4) (denoted by FFA) to that of four other estimators. In the second, we assess the relative interpretability of these estimators. In the third, which may be found in Section 4 of Supplement A, we explore how the scree plot approach for selecting the number of factors behaves in different settings. We begin, in Section 5.1, by describing the data generating procedure for these studies.

We now consider the AOMIC-PIOP1 dataset which includes six-minute resting-state scans for 210 healthy university students. Prior to each subject’s resting-state session, a structural image was acquired at a resolution of 1 mm isotropic while the functional image was acquired at a spatial resolution of 3 mm isotropic and a temporal resolution of 0.75 seconds (resulting in a single-subject data matrix with dimension $65\times 77\times 60\times 480$). Including all $\sim$ 300,000 voxels in this analysis would necessitate estimation and in-memory storage of a prohibitively large spatial covariance. We avoid this by focusing on one axial slice (at $z=30$). Preprocessing of anatomical and functional MRI were performed using a standard Fmriprep pipeline (Esteban et al., 2019), as detailed in Snoek et al. (2021).

Our strategy for discovering the connectivity patterns present in these resting-state data was to prewhiten the voxel time courses in each subject’s preprocessed scan (which are, in general, nonstationary and autocorrelated), then treat the multi-subject collection of slices as an i.i.d. sample $\mathcal{X}_{1},\dots,\mathcal{X}_{n}$, $n=100800$, which we used to fit the FFM. Although fMRI data is inherently nonstationary and autocorrelated, many connectivity studies that correlate voxel time series skip this essential prewhitening step, risking the discovery of spurious correlations (Christova et al., 2011; Afyouni et al., 2019). Of course, prewhitening only addresses temporal correlation within voxel time courses, not between time courses. There are several ways to prewhiten voxel time courses, including those based on Fourier (Laird et al., 2004) and wavelet (Bullmore et al., 2001) decompositions. Taking another common approach, we performed prewhitening by fitting a non-seasonal AutoRegressive Integrative Moving Average (ARIMA) model to each voxel time course using the R function forecast::auto.arima (Hyndman and Khandakar, 2008), then extracting the residuals. Figure 5 from Section 6.1 of Supplement A demonstrates the effects of prewhitening on voxel time courses. Taking after Christova et al. (2011), we form samples $\mathcal{X}_{1},\dots,\mathcal{X}_{n}$ from the scaled ARIMA residuals.

We begin analysis of these samples by computing the empirical covariance $\hat{\mathcal{C}}$ in equation (3). Based on visual inspection and the fact that the sample size is quite large, we opt not to smooth (i.e., we set $\alpha=0$), then plot the non-increasing function $g(j)=\|\mathcal{A}\circ(\hat{\mathcal{C}}-\hat{\theta}_{j})\|_{F}^{2}$ for $j=1,\dots,100$ (see Figure 6 of Section 6.1 in Supplement A). Based on this plot, we set $\hat{K}=12$ and use the corresponding $\hat{\mathcal{L}}_{1},\dots,\hat{\mathcal{L}}_{12}$ as the initial loading estimates. To choose the shrinkage parameter $\underline{\kappa}=(\kappa_{1},\dots,\kappa_{12})$ used to compute the post-processed loadings $\tilde{\mathcal{L}}_{1},\dots,\tilde{\mathcal{L}}_{12}$, we define a manageable search grid then follow the tuning procedure described in Section 3 of Supplement A. Figure 9a displays the rich covariance structure captured by these post-processed loadings. In addition to anatomical segmentation (e.g., grey matter segmentation in loadings 1, 2, and 6; white matter segmentation in loadings 3 and 9; and ventricle segmentation in loadings 5 and 10), these 12 loadings contain known resting-state networks. The prominent bilateral structure in the first loading captures parts of the auditory resting-state network, including co-activation in Heschl’s gyrus and the posterior insular (see $7_{20}$ of Figure 1 in Smith et al., 2009). Next, frontal activation in the second loading includes areas associated with resting-state executive control, like the anterior cingulate (see $8_{20}$ of Figure 1 in Smith et al., 2014). Lastly, the sixth loading contains activation in the occipital lobe, an area associated with visual processing (see $2_{20}$ of Figure 1 in Smith et al., 2009).

To provide a methodological comparison, we also performed group ICA by temporally concatenating the preprocessed slices then running ICA on the resulting data to produce spatial maps common to all subjects (Calhoun et al., 2001). As in Section 5, we do so using FSL’s MELODIC function. For this analysis, we fit 12-, 25-, and 50-independent component (IC) models. Transformed components (for interpretability) for the 12-IC model are shown in Figure 9, and those for the 25- and 50-IC models are shown in Figures 8, 9, and 10 of Section 5.2 in Supplement A. To facilitate comparison with FFA, we also estimate the FFM with 25 and 50 factors (see Figures 5, 6, and 7 in Section 5.2 of Supplement A) even though the plots of Figure 6 in Section 6.1 of Supplement A suggest a smaller number of factors.

We first compare the 12-component models (see Figure 9). Several prominent structures are captured by both FFA and ICA: dorsal correlations appear in the 6th loading, and in the 4th and 9th ICs; frontal correlations appear in the 2nd loading and the 7th IC; bilateral correlations appear in the 1st loading and in the 8th IC; similar ventricle correlations appear in the 5th loading and in the 1st IC; and brain outlines appear in the 4th and 7th loadings, and in the 11th IC. Despite the many common features, there exist some discrepancies (e.g., the spikes of activation in the 10th and 12th IC do not display prominently in any single loading).

However, as observed in Section 5.3, the most practically meaningful difference is that MELODIC lacks FFA’s robustness to dimension misspecification. Consider, for instance, the similar ventricle correlations in loading 5 and IC 1 of the order-12 plots in Figures 10a and 10b. In the FFMs with 25 and 50 factors, this structure is preserved. However, this structure, captured in a single IC of the 12-IC model, is splintered across 2 ICs in both the 25- and 50-IC models, with the components of the latter containing less of the original structure than those of the former. Similar fragmentation occurs for other structures found in the 12-IC model.

One way to combat ICA’s instability is through multi-scale approaches. For example, Iraji et al. (2023) estimate common spatial maps by first generating 100 half-splits of their multi-subject dataset, then fitting IC models of order 25, 50, 75, 100, 125, 150, 175, and 200 to each data split. This produces 100 sets of 900 components from which they select the 900 with the highest average spatial similarity (calculated by Pearson correlation) across the 100 sets. They then identify those ICs among the 900 that are most distinct from each other (spatial similarity ¡ 0.8). In Section 6.2 of Supplement A, we adapt this procedure to the AOMIC data, showing that the final subset of estimated ICs does, in fact, preserve low-order structure. However, this requires more computing resources than a single-model ICA approaches. Moreover, this output lacks some of the interpretability afforded by that from a single-model approach: estimated ICs from this multi-scale algorithm are cherry-picked from many IC models, consequently collapsing the straightforward independence-based interpretation arising from the assumptions of a single IC model.

This work recognizes the central contribution of Descary and Panaretos (2019) – a framework that decouples the global and local variation of a functional observation – as the proper scaffolding upon which to build a factor analytic approach for functional data. Our work extends their framework to multidimensional functional data, enhances the estimator of the global covariance via a roughness penalty, then appends a post-processing procedure that improves interpretability. The result is a factor analytic approach tailored to the study of functional connectivity in fMRI.

Three characteristics of our methodology make it uniquely suited to the study of functional connectivity. First, it outputs spatial maps possessing simple structure, which facilitates interpretation of the latent factors that drive correlations in the data. It is a hallmark of factor analytic methods and has, in fields like education and psychology, made multivariate factor analysis a popular alternative to less legible techniques, like PCA. We demonstrate the utility of this structure in the neuroimaging context. Second, our approach is not preceded by data smoothing like other functional connectivity methods. Given that connectivity studies aim to characterize correlations between distant brain regions, model estimation should exclude local variation, yet data smoothing forbids this. It smears local variation into the global signal, thereby corrupting subsequent attempts to estimate this signal. Third, our approach is robust to dimension overspecification. This model feature, not found in popular ICA-based methods, is critical in fMRI where it is difficult to correctly specify the dimension of the signal.

The limitations of our methodology present several opportunities for future work. One weakness of the approach is that its application does not consider all spatial information available in an fMRI scan. Although the model remains valid for 3-dimensional brain volumes, inclusion of a third spatial dimension renders the existing estimation method computationally impractical. Additionally, the model ignores the temporal information present in fMRI data. We attempted to mitigate the impact of these temporal dynamics in our analysis through prewhitening, but such dependencies ought to be acknowledged by the model. Future work on functional factor modeling of large-scale functional or spatiotemporal data should address these concerns.

Supplement A: Functional Factor Modeling of Brain Connectivity

The FFM is identifiable up to an orthogonal rotation. To see this, let $\mathfrak{L}\in[0,1]^{D}\times\mathbb{N}_{K}$, defined by $\mathfrak{L}_{\underline{s},k}=l_{k}(\underline{s})$, be a set of loading functions indexed by $k$. For a $K\times K$ rotation matrix $\boldsymbol{R}$, define the rotation operator $\mathcal{M}_{\boldsymbol{R}}\mathrel{\mathop{\ordinarycolon}}[0,1]^{D}\times\mathbb{N}_{K}\to[0,1]^{D}\times\mathbb{N}_{K}$ by $(\mathcal{M}_{\boldsymbol{R}}\mathfrak{L})_{\underline{s},k}=(\boldsymbol{R}\mathfrak{L}_{\underline{s},\cdot})_{k}$. That is, $\mathcal{M}_{\boldsymbol{R}}$ rotates the loading functions by rotating the $K$-dimensional vector $(l_{1}(\underline{s}),\dots,l_{K}(\underline{s}))$ by $\boldsymbol{R}$ at each fixed $\underline{s}\in[0,1]^{D}$. If $\underline{f}=(f_{1},\dots,f_{K})$ is a vector of factors and $\boldsymbol{R}$ is some rotation matrix, then the following equivalences hold:

and

Thus, when the global component $Y$ and its covariance $\mathscr{G}$ are known, it is impossible to distinguish the loading functions of $\mathfrak{L}$ from those of $\mathcal{M}_{\boldsymbol{R}}\mathfrak{L}$.

This section includes statements and proofs of results described in the manuscript, as well as those for auxiliary results that enrich our theory.