Connectivity Regression

Summary

Assessing how brain functional connectivity networks vary across individuals promises to uncover important scientific questions such as patterns of healthy brain aging through the lifespan or dysconnectivity associated with disease. In this article, we introduce a general regression framework, Connectivity Regression (ConnReg), for regressing subject-specific functional connectivity networks on covariates while accounting for within-network inter-edge dependence. ConnReg utilizes a multivariate generalization of Fisher’s transformation to project network objects into an alternative space where Gaussian assumptions are justified and positive semidefinite constraints are automatically satisfied. Penalized multivariate regression is fit in the transformed space to simultaneously induce sparsity in regression coefficients and in covariance elements, which capture within network inter-edge dependence. We use permutation tests to perform multiplicity-adjusted inference to identify covariates associated with connectivity, and stability selection scores to identify network edges that vary with selected covariates. Simulation studies validate the inferential properties of our proposed method and demonstrate how estimating and accounting for within-network inter-edge dependence leads to more efficient estimation, more powerful inference, and more accurate selection of covariate-dependent network edges. We apply ConnReg to the Human Connectome Project Young Adult study, revealing insights into how connectivity varies with language processing covariates and structural brain features.

1. INTRODUCTION

Functional connectivity, the study of temporally coincident brain activity, seeks to link geographically distinct brain regions via measures of statistical dependence (Toga 2015). Since the brain is known to sub-specialize regions for distinct responsibilities, functional connectivity sheds insight into which brain regions coordinate together by assessing mutual brain activation measured by blood oxygenation scores (BOLD) from functional MRI (fMRI) images across timepoints (Eickhoff et al. 2015). fMRI is typically classified as task-based or resting state; while task-based fMRI seeks to identify brain regions involved in a particular activity, resting state fMRI explores the intrinsic sub-specialization of brain regions at rest and how they organize together for core cognitive functions (Zhang et al. 2016b; Lv et al. 2018). In practice, functional connectivity is often represented as a network/graph based on a correlation, covariance, or precision matrix, with an entry denoting dependencies between brain regions assessed from BOLD scores (Stulz et al. 2018; Mohanty et al. 2020). For functional connectivity derived from resting state fMRI, active research areas involve either the accurate estimation of these networks or an assessment of how functional connectivity varies with covariates in order to determine the biological link between phenotype and the brain’s natural functional organization (Zhu et al. 2019). The core scientific aim of this article is the latter; we seek to develop a rigorous statistical framework to link phenotype to functional connectivity networks derived from resting state fMRI.

The Human Connectome Project (HCP) aimed to construct a comprehensive network mapping the human brain to provide insight into anatomic and functional connectivity and aid research into disorders such as autism, Alzheimer’s, and schizophrenia. Using a variety of noninvasive imaging methods, the project collects extensive connectome information for healthy, diseased, and longitudinally sampled patients. Among its large scale goals to comprehensively map brain circuits, a major application of the HCP is to assess the link between patient connectivity structure and phenotypic traits (Van Essen et al. 2013). The link between brain connectivity and phenotype has the potential to provide deeper nuanced insight into brain function. For example, contrasting overall connectivity of diseased and healthy populations highlights circuitry changes potentially driving diseases (Zhao et al. 2020). Analogously, a comparison of connectivity structure across ages helps scientists understand the biological underpinnings of the maturation process (Lawrence et al. 2019).

While these examples compare global connectivity differences, a more ambitious goal is to link phenotypic changes explicitly with circuits between specific brain regions. Given rapid advances in neuroscientific knowledge, finding precise region-specific links with phenotype can better integrate existing biological knowledge and determine the relationship between structure and function. To this end, in this article, we develop a general regression approach, Connectivity Regression (ConnReg), that relates subject-specific network connectivity objects to a set of subject-specific covariates and then characterizes which network circuits are driving these differences. We begin with a simple example in Fig. 1 to illustrate this goal and to describe one of the key conceptual underpinnings of our approach: inter-edge association that engenders a type of second order dependence.

Motivated by functional data analysis (FDA) approaches, we undertake a joint modeling approach that accounts for the internal structure of the network object. In the setting of FDA, the internal structure is the within-function correlation structure of the functional data object |$ Y(t) $| in that observations at nearby domain points |$ t $| are similar to one another; models incorporating smoothness in |$ t $| and accounting for correlation in the between-function residuals lead to more efficient estimation and more powerful joint inference (Morris 2015). In a similar spirit, for these connectivity network objects the internal structure involves inter-edge dependence structures, as illustrated in Fig. 1. In the toy example, note that individuals with strong connectivity in their vision-memory circuits also have strong connectivity in their memory-emotion and vision-emotion circuits (and vice-versa). This suggests that these circuits operate in concert with each other, inducing an association among the corresponding edges, a type of “correlation of correlations.” We wish to estimate and account for this type of second order association in our modeling that, as we will demonstrate, tends to produce more efficient estimation and inference and also reveals additional insights into the connectivity dynamics.

1.1. Network-object regression: regression approaches treating entire subject-specific networks as units of observation

In contrast to our approach, most preceding network regression approaches assume subject-specific vectors as their units of observation. Many of these methods focus on differences in network structure between groups, or other discrete covariates (Xia et al. 2015; Peterson et al. 2015; Liu et al. 2017; Durante and Dunson 2018). Ni et al. (2022) have published a review paper on Bayesian graphical models for biomedical applications which covers some methods on group graphs and their applications. The graphical model literature has also produced methods for the general regression of network edges for observed binary data (Cheng et al. 2014), directed acyclic graphs (Ni et al. 2019), and undirected graphs (Wang et al. 2022). For regression of covariance matrices on covariates, much work assumes a natural ordering of the multivariate dimensions, including autoregressive conditionally heteroscedastic models that take sequential fits of covariance across time series (Engle and Kroner 1995) and Cholesky-decomposition based models (Pourahmadi et al. 2007). Among methods for unordered variates, Chiu et al. (1996) performed joint regression of the mean and covariance on a multivariate vector, utilizing the matrix log-transformation to model covariance elements. Hoff and Niu (2012) took a random effects approach and parameterized heteroscedastic covariance by the inner product of random effect regression coefficients.

In our setting, we have sufficient data from each subject such that we are able to obtain an entire |$ K\times K $| subject-specific connectivity network across |$ K $| brain regions as opposed to a |$ K $|-dimensional vector. There exists recent literature modeling subject-specific networks as units of observation, making various tradeoffs among interpretability, scalability, and modeling higher order structure like inter-edge dependence that can be leveraged in this setting. Seiler and Holmes (2017) focused on covariance matrices and introduced 2 models. Both models implicitly assumed a Gaussian distribution on observed vectors and restricted their attention to discrete covariates representing group identities; these models also do not seek to directly model replicated network edges. Zhao et al. (2021a) proposed a 2-step procedure where group-ICA was applied to entire functional connectivity networks followed by a covariate-assisted principal component procedure. The authors’ model identifies global changes in functional connectivity based on covariates and mitigates multiplicity issues via the use of common linear projections; their approach does not focus on producing interpretable results for assessing network edge differences in the original space of networks. Xia et al. (2020) describe a modeling framework that takes community structure amongst edges into account to potentially improve inference, produces an interpretable result and is solved efficiently via a convex objective function. However, the method does not seek to directly examine the effects of covariates on edges and focuses on the specific case where sub-networks are known a priori; inter-edge dependence is treated as fixed and not estimated as an unknown parameter. Tomlinson et al. (2022) proposed a regression framework that regresses dissimilarity summary statistics between subject-specific networks on continuous and discrete traits and develop a variety of statistical hypothesis tests to accompany their regression approach. While the method provides a framework to determine a covariate’s impact on differences between subject-specific networks, the authors do not focus on assessing inter-edge dependence or characterizing individual edge effects. The authors’ work was highlighted in a review article that contrasted their distance-based framework, which relates covariates to measurements that summarize spatial or topological similarities between subject-specific networks, to complementary approaches which related covariates to the probability and strength of individual network edges (Simpson et al. 2024). To provide a comprehensive summary of approaches treating networks as the units of observation, we compare how our approach differs from recent network-response regression approaches in Fig. 2 below.

In this article, we propose a new framework for regressing functional connectivity networks, represented by correlation-matrix objects, on covariates. Our approach first projects subject-specific correlation matrices to Fisher correlations using a matrix logarithm transformation which theoretically justifies Gaussianity while ensuring positive semi-definiteness (Archakov and Hansen 2021). Using a multivariate regression approach with specification of sparsity-inducing penalties on both regression coefficients and off-diagonal precision matrix elements, we then regress connectivity on predictors while estimating and accounting for inter-edge dependence across network edges. We show this strategy leads to greater efficiency and statistical power based on the principles of seemingly unrelated regression. Furthermore, we determine which covariates affect connectivity using multiplicity-adjusted permutation tests that control experiment-wise error rates and then for selected covariates rank the edges in importance of driving these differences using a false discovery rate guided stability selection. Through extensive simulation studies we validate the performance of our permutation test and provide empirical evidence that desired edge-selection operating characteristics can be controlled via appropriately calibrated stability-selection thresholds. We apply our framework to functional connectivity network data derived from HCP (Van Essen et al. 2013). Our results demonstrate that the key findings in the transformed Fisher space comprise the dominant results in the original correlation space, and thus are interpretable in practice. We find, along with the substantial presence of second order dependence, language processing ability and 2 key anatomical areas to be significantly associated with changes in functional connectivity. The rest of the article is organized as follows. Section 2 describes ConnReg statistical framework and inferential procedures. Section 3 summarizes the simulation results aimed to assess the performance of ConnReg against competing methods. Section 4 shows the results of our framework applied to data from HCP and Section 5 concludes the article with a discussion of the properties and scope of our proposed framework.

2. CONNECTIVITY REGRESSION

2.1. Overview and notations

For each subject |$ i\,=\,1,\ldots, n $|⁠, we have a set of |$ p $| covariates |$ \boldsymbol{X_{i}}=\{X_{ij},j\,=\,1,\ldots, p\} $| and |$ K\times K $| network matrix |$ \boldsymbol{R_{i}} $| with off-diagonal elements |$ R_{ikk^{\prime}} $| providing some measure summarizing functional connectivity between brain region |$ k $| and |$ k^{\prime} $| for a prespecified set of |$ K $| brain regions. Our goal is to build a general regression framework for relating the network object |$ \boldsymbol{R_{i}} $| to covariates |$ \boldsymbol{X_{i}} $|⁠, identifying which covariates are associated with the connectivity network, and then characterizing which network circuits vary by each covariate. Although there are various options for choice of connectivity measure, here we focus on correlation, so that |$ \boldsymbol{R_{i}} $| is a |$ K\times K $| subject-specific correlation matrix. A common approach is to perform separate regressions for each pairwise association |$ R_{ikk^{{}^{\prime}}} $| (Scheinost et al. 2015; Zhang et al. 2016a). To effectively achieve this goal, Fisher’s transformation |$ \frac{1}{2}\log\left({1\,+\,R_{ikk^{\prime}}\over 1\,-\,R_{ikk^{\prime}}}\right) $| is often used to transform the pairwise correlations to the real line and engender Gaussian assumptions in the regression analyses (Fisher 1915). This approach is insufficient for our setting since we are interested in performing a joint regression on the entire network object that accounts for its internal structure and constraints. For correlation matrix network objects, the constraints require |$ \boldsymbol{R_{i}} $| to be a symmetric positive semi-definite matrix with 1’s down the diagonal and pairwise correlations |$ R_{ikk^{\prime}}\in[-1,1] $| in the off-diagonals. The internal structure includes inter-edge dependencies, higher order associations among the off-diagonal elements of |$ \boldsymbol{R_{i}} $| indicating which network circuits tend to be associated with each other. Subject-specific network objects enable this type of dependence to be estimated and taken into account.

For our setting, we propose utilizing a multivariate generalization of Fisher’s transformation, the matrix logarithm transformation applied to correlation matrices. We transform the entire  |$ M\,=\,K\times(K-1)/2 $| dimensional vector of correlations |$ \boldsymbol{\rho_{i}}=\mbox{vecl}\left(R_{ikk^{\prime}}\right) $|⁠, with vecl|$ (\cdot) $| the lower diagonal region of a symmetric matrix, into an M-dimensional vector of Fisher correlations |$ \boldsymbol{z_{i}}=[z_{i1},\ldots, z_{iM}]^{\prime}\in\mathcal{R}^{M} $|⁠, with each real vector mapping back to a legitimate correlation matrix. We will specify our regression model |$ E(\boldsymbol{z}|\boldsymbol{X}) $| in the Fisher space, which we detail next, with a fast inverse mapping producing a unique and valid predicted correlation network object |$ \boldsymbol{R} $| for any given |$ \boldsymbol{X} $|⁠.

2.2. Fisher space transformation

The matrix exponential of a matrix |$ \boldsymbol{R_{i}} $| is given by |$ e^{\boldsymbol{R_{i}}}=\sum_{k\,=\,0}^{\infty}\frac{\boldsymbol{R_{i}^{k}}}{k!} $| (Higham 2008). When |$ \boldsymbol{R_{i}} $| is a |$ K\times K $| positive semi-definite matrix, this can be written in terms of spectral decomposition |$ e^{\boldsymbol{R_{i}}}=\boldsymbol{V_{i}}e^{\boldsymbol{D_{i}}}\boldsymbol{V_{i}^{{}^{\prime}}} $|⁠, with |$ \boldsymbol{D_{i}} $| a diagonal matrix of eigenvalues |$ \lambda_{i1},\ldots,\lambda_{iK} $| for |$ \boldsymbol{R_{i}} $| corresponding to set of eigenvectors |$ \{\boldsymbol{V_{ik}},k\,=\,1,\ldots, K\} $|⁠.

In the positive semi-definite case, for a given matrix |$ \boldsymbol{R_{i}} $|⁠, there is corresponding matrix logarithm function logm given by:

with log|$ (\boldsymbol{D_{i}})=\{\log(\lambda_{1}),\ldots,\log(\lambda_{K})\} $|⁠. Archakov and Hansen (2021) introduced an algorithm for computing a matrix logarithm transformation for correlation matrices, with |$ \boldsymbol{Z}_{i}=\hbox{logm}(\boldsymbol{R_{i}}) $| being a |$ K\times K $| matrix with lower triangular elements |$ \boldsymbol{z}_{i}=\mbox{vecl}(\boldsymbol{Z}_{i})\in\mathcal{R}^{M},M\,=\,K(K-1)/2 $|⁠. When |$ K\,=\,2 $|⁠, |$ \boldsymbol{z}_{i} $| is a scalar and identical to Fisher’s transform. As a result, the matrix logarithm transform of a correlation matrix can be viewed as a multivariate generalization of Fisher’s transformation, and so we refer to the transformed space as the Fisher space and the elements of |$ \boldsymbol{z}_{i} $| as Fisher correlations or Fisher edges.

Fisher correlations have several key beneficial modeling properties. First, there is a |$ \frac{K(K-1)}{2}\leftrightarrow\frac{K(K-1)}{2} $| mapping between |$ \boldsymbol{R} $| and |$ \boldsymbol{z} $|⁠, with every correlation matrix |$ \boldsymbol{R} $| generating a unique vector |$ \boldsymbol{z} $| of Fisher correlations, and every real vector |$ \boldsymbol{z}\in\mathcal{R}^{M} $| inducing a unique symmetric positive semidefinite |$ K\times K $| correlation matrix |$ \boldsymbol{R} $| via the efficient exponentially fast inverse mapping algorithm |$ \boldsymbol{R}=g^{-1}(\boldsymbol{z}) $| (Archakov and Hansen 2021). This allows us to fit unconstrained regression models in the Fisher space that induce valid predicted correlation matrices. Second, unlike the Cholesky transformation, this transform is invariant to the ordering of the elements, so can be used with unordered categories. Third, the Fisher correlations are asymptotically Gaussian, so the use of normal-likelihood based models is justified in this space; justifiable Gaussian assumptions enable within network inter-edge associations to be fully captured parameterically by a covariance or precision matrix. Based on these properties, our approach is to specify a multivariate Gaussian regression model |$ E(\boldsymbol{z}|\boldsymbol{X}) $| in the Fisher space that induces for each |$ \boldsymbol{X} $| a predicted positive semidefinite correlation matrix |$ \boldsymbol{R}=g^{-1}(\boldsymbol{z}) $|⁠.

2.3. Sparse Fisher space multivariate regression model

Given an |$ M $|-dimensional vector of Fisher correlations |$ \boldsymbol{z}_{i} $| for each subject |$ i\,=\,1,\ldots, n $|⁠, we assume the following multivariate regression model:

Here we assume linear regression and define the regression as |$ \boldsymbol{\mu}_{i}(\boldsymbol{X}_{i})=\mbox{E}(\boldsymbol{z}_{i}|\boldsymbol{X}_{i})=\boldsymbol{\beta}_{0}+\boldsymbol{\beta}\boldsymbol{X}_{i} $|⁠, with |$ \boldsymbol{\beta}_{0} $| an |$ M $|-vector of intercepts, |$ \boldsymbol{\beta}=\{\boldsymbol{\beta}_{1},\ldots,\boldsymbol{\beta}_{p}\} $| an |$ M\times p $| matrix of regression coefficients corresponding to |$ p $|-dimensional covariate vector |$ \boldsymbol{X}_{i} $|⁠. The regression coefficient |$ \beta_{mj} $| indicates the partial linear effect of covariate |$ X_{j} $| on Fisher edge |$ m $|⁠. Predicted connectivity networks (correlation matrices) for any given |$ \boldsymbol{X}=\boldsymbol{x} $| can be obtained by |$ \boldsymbol{R}=g^{-1}\{\boldsymbol{\mu}(\boldsymbol{X}_{i}=\boldsymbol{x})\} $|⁠. While it is valid to constrain |$ \boldsymbol{\Omega} $| diagonal and fit |$ M $| parallel single response models, we are primarily interested in a model allowing general |$ \boldsymbol{\Omega} $| that can estimate and account for inter-edge dependence within connectivity networks, which as we will show, can lead to substantially more efficient estimation and inference as well as provide additional biological insights into the dynamics of brain connectivity networks. The off-diagonal elements of |$ \boldsymbol{\Omega} $| contain partial correlations among the transformed Fisher correlation edges, in some sense an association of associations. This allows our model to capture overdispersion in the subject-specific connectivity network circuits, introducing inter-edge dependence to account for network circuits that move in concert with one another. We refer to this concept as second-order dependence, since by capturing statistical dependence among network edges we are capturing variability not accounted for in the original functional connectivity networks, which we leverage for potential gains in selection accuracy and statistical efficiency. Thus, multivariate regression with general |$ \boldsymbol{\Omega} $| enables us to estimate and account for this second-order dependence in the networks.

2.4. Motivating sparsity—connection to seemingly unrelated regressions

Zellner (1963) introduced the concept of seemingly unrelated regressions (SURs), where M separate regression equations |$ \boldsymbol{z}_{mi}=\boldsymbol{x}_{mi}^{T}\boldsymbol{\beta}_{m}+\boldsymbol{\epsilon_{mi}} $| for |$ m\,=\,1,\ldots, M $| and observations |$ i\,=\,1,\ldots, n $| have correlated errors |$ \boldsymbol{\epsilon_{mi}} $|⁠.

While each regression could validly be fit separately, accounting for residual dependence in estimation led to increased efficiency in coefficient estimation in the case where each regression had a non-identical set of predictors. In our setting, we start with the same set of predictors |$ \boldsymbol{x} $| across M and utilize a sparsity prior on |$ \boldsymbol{\beta} $| to induce different predictor sets across |$ M $| to satisfy the conditions necessary for potential gains in estimation efficiency. For |$ \boldsymbol{\Omega} $|⁠, we similarly utilize a sparsity prior to both enable meaningful potential biological interpretations and to ease computational burdens associated with estimating a high dimensional object |$ \{ $|dim|$ (\boldsymbol{\Omega})=M(M-1)/2\,=\,O(K^{4}))\} $|⁠.

Recent literature has focused on the joint estimation of sparse coefficient matrix |$ \boldsymbol{\beta} $| and precision matrix |$ \boldsymbol{\Omega} $| in the context of multivariate regression. These methods have been in both the Bayesian paradigm via hierarchical models with a hyper-inverse Wishart prior on the covariance matrix (Bhadra and Mallick 2013; Consonni et al. 2017) and in a penalized optimization framework (Rothman et al. 2010; Lee and Liu 2012; Yin and Li 2013; Cai et al. 2013). In this article, we utilize the multivariate spike-slab lasso (mSSL), a penalized optimization approach that reduces to solving a series of penalized likelihood problems with a penalty that adapts until convergence (Deshpande et al. 2019). Analogous to the well-known Lasso’s connection to a posteriori optimization with a Laplace prior, the mSSL has a connection to a posteriori optimization with separate priors that each consist of a mixture of Laplace distributions on coefficients and precision elements. The method is fit with an Expectation Conditional Maximization (ECM) algorithm, producing sparse point estimates of both |$ \boldsymbol{\beta} $| and |$ \boldsymbol{\Omega} $| that have empirically been shown to produce results with excellent estimation accuracy (Deshpande et al. 2019). We next outline ConnReg|$ {}_{\rm{mSSL}} $|⁠, an implementation of our ConnReg framework in which we use the mSSL to induce sparsity in both |$ \boldsymbol{\beta} $| and |$ \boldsymbol{\Omega} $| while accounting for second order dependence. Given that we use an optimization-based approach that produces sparse point estimates, we describe how we build an inferential procedure that obtains multiplicity adjusted P-values for which covariates affect functional connectivity networks and produces bootstrap-based stability selection scores to characterize which network edges are driving those effects.

2.5. Inducing sparsity via the multivariate spike-slab lasso

Conceptually, the mSSL uses sparsity priors to jointly produce estimates of coefficients |$ \boldsymbol{\beta} $| and inter-edge precision |$ \boldsymbol{\Omega} $|⁠. The mSSL’s penalized multivariate optimization framework equates to placing separate spike-slab priors on the elements of |$ \boldsymbol{\beta} $| and |$ \boldsymbol{\Omega} $| (Deshpande et al. 2019). The spike-slab lasso prior is a mixture of 2 Laplace distributions with a large and small scale parameter:

|$ \delta_{mj} $| is a 0-1 indicator of whether the regression coefficient for covariate |$ j $| and Fisher edge |$ m $| is substantially large enough to be included in the model. For |$ \boldsymbol{\beta}_{mj} $|⁠, large scale parameter |$ \lambda_{1}^{-1} $| is designed to capture substantially large coefficients and small scale parameter |$ \lambda_{0}^{-1} $| should capture negligible ones to be forced to 0. Analogously, |$ \kappa_{mm^{\prime}} $| is a 0-1 indicator of whether precision element |$ (m, m^{\prime}) $| is substantially large enough to be included in the model. For the elements of |$ \boldsymbol{\Omega} $|⁠, large scale parameter |$ \xi_{1}^{-1} $| is designed to capture substantially large precision elements and small scale parameter |$ \xi_{0}^{-1} $| should capture negligible elements to force to 0. Through the use of separate hyperparameters for the active and null set, the mSSL avoids overshrinkage issues for elements in the active set commonly encountered in optimization-based strategies. As is typical for spike-slab frameworks, Beta-Bernoulli priors are placed on elements of |$ \eta $| and |$ \delta $|⁠:

2.6. Model fitting

To produce a final sparse estimate for |$ \boldsymbol{\beta} $| and |$ \boldsymbol{\Omega} $|⁠, the mSSL uses an ECM algorithm to create a sequential series of optimized point estimates for both coefficients |$ \boldsymbol{\beta} $| and inter-edge precision matrix |$ \boldsymbol{\Omega} $|⁠, where each solution was “warm-started” at the previous solution. The sequence is defined via ladders for |$ \lambda_{0} $| and |$ \eta_{0} $|⁠: |$ L_{\lambda}=\{\lambda_{0_{1}} \lt\ldots\lt\lambda_{0_{L}}\} $| and |$ L_{\eta}=\{\eta_{0_{1}} \lt\ldots\lt\eta_{0_{L}}\} $| and a final solution is obtained when estimates of |$ \boldsymbol{\beta} $| and |$ \boldsymbol{\Omega} $| are determined to have empirically converged (Deshpande et al. 2019). For each step along the ladder, the optimization strategy is equivalent to maximizing a posterior mode for a given specification of hyperparameters (⁠|$ a_{\theta} $|⁠, |$ b_{\theta} $|⁠, |$ a_{\eta} $|⁠, |$ b_{\eta} $|⁠, |$ \lambda_{1} $|⁠, |$ \lambda_{0} $|⁠, |$ \xi_{1} $|⁠, |$ \xi_{0} $|⁠). The algorithm first holds (⁠|$ \boldsymbol{\Omega},\boldsymbol{\eta} $|⁠) constant and updates covariate coefficient matrix (⁠|$ \boldsymbol{\beta},\boldsymbol{\eta} $|⁠) via refined coordinate ascent and a Newton algorithm. Then, given these estimates, the algorithm updates |$ \boldsymbol{\eta} $| (closed form) and inter-edge dependence matrix |$ \boldsymbol{\Omega} $|⁠, which reduces to solving a graphical lasso problem (Friedman et al. 2008).

Given a procedure for finding sparse estimates of |$ \boldsymbol{\beta} $| describing covariate-edge associations and |$ M\times M $| matrix |$ \boldsymbol{\Omega} $| describing edge-edge associations, we now describe ConnReg’s procedure for performing inference.

2.7. Inference

ConnReg’s inferential aims are to rigorously select which covariates are associated with changes in functional connectivity networks and to characterize which particular network edges drive significant covariate effects on functional connectivity. To select covariates, we propose a permutation test that produces multiplicity adjusted P-values for which covariates are significantly associated with functional connectivity networks. Our method permutes rows of our predictor matrix |$ X $| to generate a null distribution |$ N_{max}^{G} $|⁠, for a given regression method and then computes a test statistic |$ N_{j} $| on unperturbed data. Null |$ N_{max}^{g} $|⁠, for |$ g\,=\,1,\ldots,G $| permutations, is computed by first summing indicators for non-zero covariate-network edge combinations for each covariate. Given scores for |$ J $| covariates (⁠|$ N_{j}^{g} $|⁠, |$ j\,=\,1,\ldots,J $|⁠), |$ N_{max}^{g} $|⁠, the maximum of these scores, is then used for the null to control for multiplicity. To determine which network edges drive these covariate effects, we produce stability selection scores with a computationally tractable bootstrap-based procedure that gives network edges a natural ordering in terms of importance when characterizing a covariate effect. As a secondary benefit, our stability selection procedure can be applied to precision matrix |$ \boldsymbol{\Omega} $| to create a measure of confidence for which edges display inter-edge within network dependence. Our permutation test and stability selection procedure is fully outlined in Algorithms 1 and 2 below.

2.8. Assessment of results in space of correlation matrices

Given the lack of an exact one-to-one transformation between the original space and Fisher space, we compare differences in predicted connectivity networks at distinct values of a covariate to measure its effect in the original space of correlation matrices, holding all other covariates constant: |$ \boldsymbol{R}_{\Delta}=g^{-1}\{\boldsymbol{\mu}(\boldsymbol{X}_{i}=\boldsymbol{x}_{2})\} $|  |$ -\boldsymbol{g}^{-1}\{\boldsymbol{\mu}(\boldsymbol{X}_{i}=\boldsymbol{x}_{1})\} $|⁠. This assessment of covariate-driven graph differences between networks is an active research area with many applications in neuroscience, as reviewed in recent work by Tomlinson et al. (2022). For example, one could compare differences in predicted graphs for a dichotomous case–control covariate while holding other covariates at their mean values. For our approach, we will see that the graph differences observed in the Fisher space comprise the predominate effects observed in the original space.

3. SIMULATION STUDIES

We conduct simulation studies to investigate the performance of our proposed methods. To reflect uncertainty in the estimation of subject-specific networks, we estimate each of |$ n $| “observed” subject-specific correlation matrices from 1,000 draws of a multivariate normal distribution |$ \hbox{MVN}(0,\boldsymbol{\Gamma}_{sim}^{*}) $|⁠, where |$ \boldsymbol{\Gamma}_{sim}^{*} $| is a known |$ q\times q $| subject-specific correlation matrix constructed to vary by covariates. |$ \boldsymbol{\Gamma}_{sim} $| is generated in the Fisher space for all subjects, where Gaussian assumptions are justified, correlation matrix constraints are automatically satisfied, and inter-edge dependence among |$ m=\frac{q(q-1)}{2} $| network edges can be fully controlled. The Fisher-space generative model used to produce |$ n\times m $| dataset |$ \boldsymbol{\Gamma}_{sim} $| is given by |$ \boldsymbol{\Gamma}_{sim}=\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{E} $|⁠. Each row of |$ \boldsymbol{\Gamma}_{sim} $| induces a unique |$ q\times q $| true subject-specific correlation matrix |$ \boldsymbol{\Gamma}_{sim}^{*} $|⁠. Subject-specific residuals |$ \boldsymbol{E} $| are generated by |$ MVN(0_{M},\boldsymbol{\Lambda}_{\rho}) $|⁠, where |$ \boldsymbol{\Lambda}_{\rho} $| directly controls the level of inter-edge dependence in the dataset. We assume |$ \boldsymbol{\Lambda}_{\rho} $| follows a Toeplitz structure: |$ \boldsymbol{\Lambda}_{\rho}=(\boldsymbol{\rho}^{|k-k^{{}^{\prime}}|})^{M}_{k, k ^{{}^{\prime}}=1} $|⁠. Covariates X are produced from |$ n $| draws of |$ MVN(0_{p},\boldsymbol{I}_{p}) $|⁠.

In this simulation, we fix the number of edges to |$ m\,=\,45 $| (corresponding to a |$ 10\times 10 $| correlation network) and the number of covariates to |$ p\,=\,30 $| while varying sample size (⁠|$ n\,=\,150,250 $|⁠) and level of inter-edge dependence (⁠|$ \rho\,=\,0.9, 0.7, 0.5, 0.2 $|⁠). For each setting of |$ n $| and |$ \rho $| considered, we simulate 100 datasets (900 datasets in total). All results for |$ n\,=\,150 $| are included in Appendix A.1 as well as a setting with inter-edge dependence structure resembling our real data (ie non-Toeplitz structure including positive and negative second-order correlations). For every simulated dataset, a sparse |$ \boldsymbol{\beta} $| is generated as follows. First, 30 |$ \% $| of edges are randomly selected to have any signal with covariates. Next, for each selected edge 30 |$ \% $| of covariates are randomly selected to have non-zero signal with that edge. Five covariates are then randomly selected to have 0 effect on any edges. For selected elements of |$ \boldsymbol{\beta} $|⁠, |$ \frac{1}{3} $| are randomly drawn from |$ \hbox{runif}[0.1,0.5] $|⁠, |$ \frac{1}{3} $| are randomly drawn from |$ \hbox{runif}[0.5,1] $|⁠, and |$ \frac{1}{3} $| are randomly drawn from |$ \hbox{runif}[1,1.5] $| to assess our method’s robustness to capturing signals of varying strength.

We next assess our method’s performance on our inferential framework. To determine the impact of accounting for second-order dependence, we compare our recommended approach against 2 modified versions of our framework that run single response sparse regressions on each edge (using the well-known LASSO (Tibshirani 1996) and Sorted |$ L_{1} $| Penalized Estimation (SLOPE) (Bogdan et al. 2015)) that do not take the presence of inter-edge dependence into account. We also assessed the performance of the method of Zhao et al. (2021b) on our simulated data but did not include results for this method due to the high Type I error rate. The authors’ method focuses on identifying principal linear projections associated with both subject-specific covariance matrices and covariates, so may be susceptible to false positives in a multiple covariate setting seeking to identify significant covariate-network edge associations because of the way it borrows strength across networks. More information on the SLOPE’s optimization procedure is available in Appendix A.1.

3.1. Assessing inferential framework

Table 1 illustrates the performance of our permutation-based hypothesis test and our bootstrap-based edge selection procedure. We compare our recommended approach (⁠|$ \hbox{ConnReg}_{\rm{mSSL}} $|⁠) against 2 modified versions of our framework (⁠|$ \hbox{ConnReg}_{\rm{SLOPE}} $|⁠, |$ \hbox{ConnReg}_{\rm{LASSO}} $|⁠) that do not take second order dependence into account. For our permutation-based hypothesis test, we assess power and experiment-wise type 1 (⁠|$ T_{1} $|⁠) error. To assess edge-level selection, 200 bootstrap samples were run for each simulated dataset considered. With the false positive rate fixed at |$ \hbox{FP}=0.01 $|⁠, we compare the true positive rate (⁠|$ \hbox{TP}_{01} $|⁠) across all methods. Likewise, we compare a partial area under the receiver operating characteristic curve (⁠|$ \hbox{\rm{pAUC}}_{01} $|⁠) up to |$ \hbox{FP}=0.01 $| across all methods. Finally, we compare estimation accuracy for regression coefficient |$ \boldsymbol{\beta} $| across all 3 methods. A full explanation of how |$ \hbox{\rm{pAUC}}_{01} $| is computed, additional simulations testing the sensitivity of our inferential procedure to the number of bootstraps and permutations, and additional simulations testing our framework’s performance in settings with both varying levels of sparsity in coefficient matrix |$ \beta $| and varying levels of signal strength in coefficient matrix |$ \beta $| can be seen in Appendix A.1.

|$ \hbox{ConnReg}_{\rm{mSSL}} $| demonstrates higher power and lower experiment-wise |$ T_{1} $| error at |$ \rho\,=\,0.2 $|⁠, |$ \rho\,=\,0.5 $|⁠, |$ \rho\,=\,0.7 $| and |$ \rho\,=\,0.9 $| (with differences becoming more stark as |$ \rho $| increases). |$ \hbox{ConnReg}_{\rm{SLOPE}} $| also produces competitive results at |$ \rho\,=\,0.2 $| and |$ \rho\,=\,0.5 $|⁠, suggesting its utility in settings where our framework is applied to data with low inter-edge dependence or with high dimension where computational speed is important. For edge selection, |$ \hbox{ConnReg}_{\rm{mSSL}} $| displayed superior true positive rates at the edge level while maintaining a low false positive rate of 0.01 at |$ \rho\,=\,0.5,0.7,\hbox{and}0.9 $|⁠. |$ \hbox{ConnReg}_{\rm{mSSL}} $| also demonstrated superior pAUC for all settings of |$ \rho $| considered. Finally, each method also demonstrated different stability selection thresholds needed to obtain a false positive rate of 0.01; |$ \hbox{ConnReg}_{\rm{mSSL}} $| only required stability thresholds under 0.1 for all values of |$ \rho $| while comparatively the |$ \hbox{ConnReg}_{\rm{LASSO}} $| required thresholds over 0.9, suggesting that when run without stability selection |$ \hbox{ConnReg}_{\rm{LASSO}} $| tends to be liberal and |$ \hbox{ConnReg}_{\rm{mSSL}} $| conservative in edge selection. The conservative selection properties of the mSSL are also reflected in the low type I error observed by |$ \hbox{ConnReg}_{\rm{mSSL}} $| across all settings. Additional simulations run in Appendix A.1 reveal similar trends to those seen in the main body of the manuscript, with particularly stark differences in power between |$ \hbox{ConnReg}_{\rm{mSSL}} $| and other marginal approaches in settings where |$ \beta $| is extremely sparse. This improvement in power is also particularly noticeable in simulations where non-sparse entries of |$ \beta $| have low magnitude. This phenomena demonstrates the benefits of a joint approach, where by the principles of SURs leveraging shared structure across responses can improve the ability to capture low-to-moderate signal in the presence of substantial noise (Zellner 1963; Deshpande et al. 2019).

In terms of estimation accuracy, |$ \hbox{ConnReg}_{\rm{mSSL}} $| showed superior average |$ \beta_{MSE} $| at all levels of |$ \rho $| considered. Estimation accuracy of |$ \hbox{ConnReg}_{\rm{SLOPE}} $| and |$ \hbox{ConnReg}_{\rm{LASSO}} $| can be improved by using a 2-step procedure where selection is first performed followed by ordinary least squares (OLS) to reduce bias—|$ \hbox{ConnReg}_{\rm{mSSL}} $| does not have an ad-hoc 2-step alternative. Estimation accuracy results using these 2-step procedures are shown in Appendix A.1 and find that |$ \hbox{ConnReg}_{\rm{LASSO}} $| and |$ \hbox{ConnReg}_{\rm{SLOPE}} $| perform competitively at low levels of |$ \rho $|⁠. For |$ \hbox{ConnReg}_{\rm{LASSO}} $|  |$ \lambda $| was selected by 10-fold cross-validation minimizing MSE.

In summary, simulation results provide evidence for the validity of our inferential procedure. For our designed permutation hypothesis test, in the presence of inter-edge dependence |$ \hbox{ConnReg}_{\rm{mSSL}} $| dominates in terms of power and experiment-wise |$ T_{1} $| error compared to approaches running independent edge-specific regressions. Similarly, for our bootstrap-based edge selection procedure in the presence of inter-edge dependence our recommended approach shows superior true positive rates and higher pAUC at fixed false positive rates. Finally, the superior estimation accuracy of |$ \boldsymbol{\beta} $| by |$ \hbox{ConnReg}_{\rm{mSSL}} $| in the presence of inter-edge dependence demonstrates the tangible benefits of accounting for second order associations when regressing subject-specific networks on covariates. Our framework also performed well with faster parallel-edge specific regression approaches when inter-edge dependence was low, suggesting it can effectively be used to perform rigorous inference for applications with low inter-edge dependence or for applications where the practitioner wants to scale up to networks with very high dimensions. We note that results from all simulations conducted with sample size |$ n\,=\,150 $| and with dependence structure matching our Human Connecctome Project data set (see Appendix A.1) demonstrate the same beneficial properties for our approach as those observed in the main body of this article for |$ n\,=\,250 $|⁠.

4. ANALYZING THE HCP YOUNG ADULT STUDY

We apply ConnReg to characterize heterogeneity in connectivity patterns of 1,003 subjects taken from the HCP Young Adult Study (Van Essen et al. 2013). Subjects are healthy individuals between the ages of 22–35 years old with fMRI BOLD score time series summarized for 15 aggregated functional brain regions (Glasser et al. 2016; Robinson et al. 2018). Detailed information on data retrieval, pre-processing steps, and a comprehensive anatomical breakdown of each aggregated functional brain region is discussed in Appendices A.2.1 and A.2.2. We considered 7 covariates: 2 cognitive tests related to language processing ability, and the structural area of 5 brain regions related to cortical hubs connecting intrinsic functional connectivity sub-networks. Cortical hubs are hypothesized to be the means through which distinct sub-networks communicate with one another (Buckner et al. 2009; Liu et al. 2018). From the pre-processed time series, we form |$ 15\times 15 $| subject-specific networks and transform to the Fisher space to produce a |$ 1,003\times 105 $| stacked matrix of network edges. For hypothesis tests, 250 permutations were used to generate null distributions. An initial empirical assessment of our data (see Appendix A.2.3) indicated a high presence of inter-edge dependence within subject-specific networks. Similarly, an empirical comparison of within network inter-edge dependence in the Fisher and original space of correlation networks show highly similar patterns, in terms of both magnitude and direction. This high degree of similarity suggests that it may be reasonable to directly interpret second-order associations found in the Fisher space. Our permutation test took 12.76 min and the bootstrap-based stability selection took 14.00 min run on a personal computer with a 3.3 GHz Dual-Core Intel Core i5 processor and 16 GB of memory.

4.1. Linking phenotype with functional connectivity networks

The 7 covariates considered are grouped into language processing assessments (“ReadEng,” “PicVocab”) and anatomical areas of intrinsic functional connectivity hubs (“Precuneus,” “Postcentral Gyrus,” “Cuneus,” “Posteriorcingulate,” “Temporalpole”). We find ReadEng, Precuneus area, and Postcentral Gyrus area (P-value: <0.004) to be significantly associated with changes in subject-specific functional connectivity networks and all other covariates to not be significantly associated (P-values: >0.996). Given covariates found to significantly vary with functional connectivity, we ran our bootstrap-based stability selection procedure to identify which edges were driving these effects. Figure 3 shows heatmaps demonstrating stability selection results over 200 bootstrap samples for covariate-network edge associations and second order edge-edge associations in the Fisher space. Stability selection results for second order associations validate the empirical examination performed at the beginning of this section, as many within-network edge-edge associations are found to be associated over 85% of the time (ie stability selection scores >0.85). A comprehensive overview of first and second order stability selection results is shown in Appendix A.2.3.

Holding all other covariates at mean levels, to characterize a covariate’s impact on functional connectivity we assess graph differences for predicted networks at the |$ X_{0.1} $| percentile and the |$ X_{0.9} $| percentile in both the Fisher space and the space of correlation networks. Figure 4 displays graph differences for our 3 flagged covariates. Results are overlaid with stability selection proportions conducted for each covariate-edge combination in Fisher space. We observe that salient features are preserved in both Fisher space and correlation space with respect to direction and magnitude of edge differences. Predicted network differences are sparse in the Fisher space by construction, and approximately sparse in the original correlation space but with some near zero magnitude edges appearing when transformed from the Fisher space. The results provide strong evidence that the primary edges driving functional connectivity changes can be identified in Fisher space both in terms of magnitude and direction (increasing or decreasing)—the dominant edge differences identified in Fisher space in Fig. 4 are the same sign (color), similar magnitude (thickness), and the same direction (solid vs. dashed) when transformed to the original space of correlation networks for each significant covariate.

For ReadEng, prominent edges involve nodes that correspond to the cerebellum, amygdala, brain stem, thalamus, and putamen. Combined with known functions of these subcortical structures, recent research implicating the cerebellum’s role in language recognition indicates that predicted graph differences may link language recognition both with emotional response and with episodic memory and emotional expression (Koziol et al. 2014; Basinger and Hogg 2019; Torrico and Munakomi 2019). For Precuneus Area, the most prominently flagged edges involve nodes that correspond to the hippocampus, cerebellum, amygdala, thalamus, diencephalon, and brain stem. Given the function of the precuneus, edges seem to link episodic memory and motor movement, addiction, and emotion (Fauchon et al. 2019; Tomasi and Volkow 2020; Berkovich-Ohana et al. 2020). The most prominent edge flagged for Postcentral Gyrus Area involves nodes that correspond to the brain stem, cerebellum, and hippocampus. Given the function of the postcentral gyrus, the majority of edges seem to link memory with motor movement, with the association between node 5 and node 6 possibly linking emotional response with episodic memory (Fu et al. 2019; DiGuiseppi and Tadi 2020). A more comprehensive biological interpretation of prominent flagged edges for all significant covariates is given in Appendix A.2.3.

4.2. Interpretation of second-order dependence

Given the presence of substantial inter-edge dependence, Fig. 5 illustrates the mostly highly correlated edges (those with stability selection = 1 and partial correlation |$ \geq $| 0.7). The cerebellum, heavily represented in node 8, has a central role in regulating associations between anatomical areas and appears to act as a “hub” of strong second order associations. A comprehensive overview of all potential second order associations, in terms of stability selection and partial correlation, is provided in Appendix A.2.3 along with a more comprehensive interpretation of the second order associations highlighted in Fig. 5.

5. DISCUSSION

In this article, we introduce a principled framework for regressing entire correlation networks on covariates, connectivity regression (ConnReg), motivated by applications in functional connectivity. We utilize a sparse multivariate linear model in a transformed Fisher space that enables unconstrained modeling on the real line, justifiable Gaussian assumptions, and automatically satisfied correlation matrix constraints. Via the precision matrix in the transformed Fisher space, our framework is able to fully estimate and account for within network inter-edge dependence, which is shown via simulations to lead to gains in estimation efficiency and selection accuracy. We introduced a permutation hypothesis test that produces multiplicity-adjusted P-values for which covariates vary with functional connectivity networks and a bootstrap-based stability selection procedure that characterizes which network edges drive significant covariate effects. Results in the Fisher space are shown to be highly concordant with those transformed back to the space of correlation networks, aiding interpretability.

The fact that connectivity regression of a |$ p\times p $| network requires assessment of an |$ \mathcal{O}(p^{4}) $| second-order association matrix raises computational challenges to apply |$ \hbox{ConnReg}_{\rm{mSSL}} $| to very large networks. However, assessment of functional connectivity networks represented by lower dimensional embeddings is an active area of contemporary research (Casanova et al. 2021; Gallos et al. 2021; Gao et al. 2021), so networks of small to moderate size are clearly of interest. Also, |$ \hbox{ConnReg}_{\rm{SLOPE}} $| runs much faster, not requiring calculation of the |$ \mathcal{O}(p^{4}) $| second-order association matrix, and does not sacrifice too much performance relative to |$ \hbox{ConnReg}_{\rm{mSSL}} $| when the magnitude of second order dependence is not too high, so could be an alternative to consider for very high-dimensional networks. In future work, we will seek to improve the computational efficiency of our joint modeling algorithm to enable the benefits of accounting for second order associations while scaling to larger networks. Other future extensions can include the relaxation of linearity assumptions to allow connectivity to vary smoothly with continuous covariates while maintaining computational feasibility or the incorporation of random effects to adjust for inter-individual correlation structure when multi-level designs are encountered.

In conclusion, we view this manuscript’s framework and the idea of fully leveraging within network inter-edge dependence as fundamental building blocks for new methods to be developed for the improved assessment of the link between functional connectivity and phenotypes. In the context of existing literature, our method seeks to determine how whole networks, with their dependence structure preserved, vary across covariates. As an alternative to existing distance-based approaches that relate covariates to measures summarizing similarities between subject-specific brain networks (Simpson et al. 2024), our method assesses this question at both the whole-network and edge-specific (ie specific network connection) level, finding which covariates are associated with changes in subject-specific functional connectivity networks (via multiplicity adjusted statistical inference) and determining which network edges are most associated with changes in covariates (bootstrap-based stability selection). In addition, our method also directly assesses the association among network edges themselves (ie second order dependence) and leverages these associations for gains in statistical efficiency for relating covariates to functional connectivity networks. Using shared structure to improve efficiency can be seen in other fields such as machine learning, where research in multi-task learning seeks to perform multiple tasks in parallel while using a shared representation across tasks (Caruana 1997).

SOFTWARE

The Supplementary Materials for Connectivity Regression contain both additional simulation results and additional analyses conducted on the data from the Human Connectome Project (HCP) Young Adult Study analyzed in this article. Software used to conduct the analyses in this manuscript can be found at https://github.com/nmd1994/Connectivity-Regression.

SUPPLEMENTARY MATERIAL

Supplementary material is available at Biostatistics Journal online.

FUNDING

This work was partially supported by R01-CA244845, R01-CA178744, UL1-TR001878, P30-CA046592, R01-MH112847, and R01-MH123550.

CONFLICT OF INTEREST

None declared.

References