{'fraction_non_alphanumeric': 0.03411595927079324, 'fraction_numerical': 0.01578081808320795, 'mean_word_length': 4.680779169151262, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian Modeling of Effective and Functional Brain Connectivity using Hierarchical Vector Autoregressions

Analysis of brain connectivity is important for understanding how information is processed by the brain. We propose a novel Bayesian vector autoregression (VAR) hierarchical model for analyzing brain connectivity in a resting-state fMRI data set with autism spectrum disorder (ASD) patients and healthy controls. Our approach models functional and effective connectivity simultaneously, which is new in the VAR literature for brain connectivity, and allows for both group- and single-subject inference as well as group comparisons. We combine analytical marginalization with Hamiltonian Monte Carlo (HMC) to obtain highly efficient posterior sampling. The results from more simplified covariance settings are, in general, overly optimistic about functional connectivity between regions compared to our results. In addition, our modeling of heterogeneous subject-specific covariance matrices is shown to give smaller differences in effective connectivity compared to models with a common covariance matrix to all subjects.


Introduction
The use of functional magnetic resonance imaging (fMRI) to investigate brain connectivity dates back to seminal papers from the mid-1990:s (Friston 1994, Biswal et al. 1995, and during the last decade the interest has increased dramatically ). Brain connectivity is an important tool in understanding how information is processed by the brain, with applications in both non-clinical and clinical settings. On the clinical side, connectivity is frequently investigated in association with different kinds of neuropsychological conditions, e.g. schizophrenia (Lynall et al. 2010), ADHD (Konrad and Eickhoff 2010), epilepsy (Morgan, Abou-Khalil and Rogers 2015) and autism spectrum disorder (ASD) (Easson, Fatima and McIntosh 2019). The most common approach for studying brain connectivity using fMRI is so-called resting state fMRI (rs-fMRI), and the analyses are mainly divided into functional and effective connectivity. Functional connectivity pertains to the investigation of undirected associations between brain regions, whereas effective connectivity refers to directed associations (Friston 1994(Friston , 2011. Our work targets both functional and effective connectivity, which is new in the literature, and the possibility to perform both group-level inference, which is most common and important in practice, as well as subject-specific inference. In an rs-fMRI scan, the subject is typically scanned for around ten minutes or less, without specific instructions or tasks to complete. The brain is divided into three-dimensional pixels, called voxels, with a side of about 3-4 mm. Within each voxel, the blood oxygenation level depend (BOLD) signal is measured, which is a proxy measure for neural activity. The BOLD signal in every voxel is usually sampled at about 0.5-2 seconds intervals; for a rigorous and detailed introduction to fMRI, see Buxton (2009). The data obtained are therefore voxelwise time series, where the time series length is usually 200-1000 observations. The number of voxels tend to be in the hundreds of thousands, and some kind of dimension reduction is therefore typically used to divide the brain into at most a few hundred functional regions using a brain atlas, see e.g. Power et al. (2011) or Glasser et al. (2016). A single time series for each region is then constructed by e.g. averaging the time series for all voxels within the region. Our work assumes the rs-fMRI data has been preprocessed to yield such regionwise time series.
The overarching statistical problem in functional and effective connectivity analysis is to estimate suitable measures of association between brain regions. The perhaps simplest solution, which is still widely used in practice, is to calculate the pairwise Pearson correlation, or partial correlation, between every pair of regions. These approaches have the obvious drawback of completely disregarding autoregressive dependencies, both within and between regions. Afyouni, Smith and Nichols (2019) point out that the standard error of the sample correlation coefficient becomes biased from autocorrelation, which implies that the commonly used Fisher transformation can not stabilise the variance.
Attention has therefore turned to time series modeling, e.g. wavelet expansions (Zhang et al. 2014) and vector autoregressive (VAR) models either in the time domain (Chiang et al. 2017, Goebel et al. 2003, or in the frequency domain (Cassidy, Rae and Solo 2015;Cassidy et al. 2018). Dynamic Causal Models (DCM, Friston, Harrison and Penny 2003), a state-space model with ambitious neurophysiological modeling have been extended from task fMRI data to analyze connectivity from resting-state data (Friston et al. 2014).
We propose a Bayesian hierarchical VAR model for both effective and functional connectivity that accounts for autoregressive dependencies within and between brain regions. The hierarchical setting allows for comparison of grouplevel inference in a straightforward manner. Our results differ substantially from other comparable VAR models in the brain connectivity literature. In general, the results from more simplified covariance settings overestimate functional connectivity between regions compared to our results. We also observe smaller differences in EC, which we suspect are mainly due to heterogeneous subject-specific covariance matrices, which other approaches with a common covariance matrix for all subjects can not account for. We fit these more complex and computationally very demanding models by integrating out the subject-specific parameters to obtain posterior inference on the group-level parameters using highly efficient Hamiltonian Monte Carlo (HMC) sampling. The article is organized as follows. In Section 2, we define our proposed model with subject-specific covariance matrices, and models with more simplified covariance settings (Chiang et al. (2017) and Gorrostieta et al. (2012Gorrostieta et al. ( , 2013). Our Bayesian setting is described and posterior inference is derived in Section 3. Group-level inference on a real rs-fMRI data set with controls and individuals diagnosed with ASD is presented in Section 4. Concluding remarks are given in Section 5.

A Bayesian VAR hierarchical modeling for Brain Connectivity
This section describes a Bayesian vector autoregressive (VAR) hierarchical model for brain connectivity applied to a group of subjects. The model allows for subject-specific VAR parameters centered around a group-level VAR. We also discuss two special cases of our model which have been used for effective connectivity by Chiang et al. (2017) and Gorrostieta et al. (2012Gorrostieta et al. ( , 2013. Section 3 proposes an efficient posterior sampling algorithm for the model that combines analytical marginalization with HMC.

The hierarchical VAR model
Let Y rst be the fMRI BOLD signal for subject s in region r at time t, where r = 1, . . . , R, s = 1, . . . , S, and t = 1, . . . , T . The Bayesian VAR (BVAR) model of order L can be defined for where N k () and IW () denote the k-dimensional multivariate normal and inverse Wishart distributions, respectively. Thus, the user needs to specify prior precision matrices P 0 , Ψ 0 , P s for s = 1, . . . , S , and the degrees of freedom ν 0 > R + 1. We elaborate on our chosen prior specification in Section 3.3. We label the general hierarchical VAR in (1) as Model 1. The global parameters B and Σ in (1) are of main interest, in particular comparing the connectivity implied by B and Σ between groups of subjects; see the application in Section 4 where a group of subjects diagnosed with autism spectrum disorder (ASD) are compared to healthy controls.
We also consider the following two submodels of Model 1. First, Model 2 assumes a common covariance matrix for all subjects, i.e. Σ s = Σ for all s. This model is clearly nested in Model 1, since Model 3 makes the additional simplification also made by Chiang et al. (2017) and Gorrostieta et al. (2012Gorrostieta et al. ( , 2013; that the common covariance matrix Σ is diagonal with prior independent diagonal elements following a conjugate inverse gamma (IG) distribution (Press 2005), Hence, Model 3 does not allow estimation of functional connectivity from nonzero off-diagonal elements of Σ.

Bayesian Inference
The Bayesian approach updates a prior distribution for all model parameters with observed data through the likelihood function to a posterior distribution where B 1:S = {B 1 , . . . , B S }, and Σ 1:S is defined analogously. The object of main interest is the marginal posterior of the group-level parameters which is obtained by integrating out the subject-specific parameters from (2). The marginal posterior distribution in Model 2, i.e. the model with common Σ for all subjects, can be derived in closed form. However, the posterior for the general model where all parameters are subject-specific is not tractable. We propose an efficient HMC algorithm to sample from the marginal posterior in (3). The analytical result for Model 2 is also exploited for determining the prior hyperparameters in Model 1.

Posterior inference when Σ is common to all subjects
Let X s = (Y s,t−1 , . . . , Y s,t−L ) be the set of covariates in the BVAR model. The likelihood function of (B s , Σ) for each subject s is given by (see Appendix A for details) The marginal likelihood function of (B, Σ) for all subjects is given by Multiplying this marginal likelihood with the prior distribution of (B, Σ), the posterior distribution of (B, Σ) becomes (see Appendix A for details) where c 0 does not depend on B and Σ, Conditional on Σ the posterior distribution of B is given by i.e. a matrix-Normal distribution with posterior meanB as a weighted aver- where ν n = ν 0 + Sn − q. Hence, the marginal posterior distribution of Σ is an Inverse-Wishart distribution with ν n degrees of freedom and scale matrix Ψ n in Model 2. This implies for Model 3 that the marginal posterior distribution of each element rr in the diagonal matrix Σ follows an inverse gamma distribution as (Press 2005) Σ rr |Y ∼ IG ν n − 2R 2 , ν n Ψ nrr 2 .

Posterior inference for the hierarchical VAR with subject-specific Σ s
Replacing Σ with Σ s in Equation (4) gives the likelihood function of (B s , Σ s ) for each subject s. Then, the marginal likelihood function of (B, Σ, ν) becomes (see Appendix B for details) where c 1 does not depend on B, Σ, and ν. The posterior distribution of (B, Σ, ν) is intractable and high-dimensional, so we use the HMC algorithm with hyperparameters tuned adaptively using the No-U-Turn Sampler (NUTS) (Hoffman and Gelman, 2014) to sample from the posterior. We implement the algorithm in the probabilistic programming language Stan, see Appendix C for the Stan model specification. To monitor convergence to the posterior, we run three parallel MCMC chains until the diagnostic convergence measureR in Gelman and Rubin (1992) is close to 1.

Prior specification
Let max s 2 r be the maximum sample variance in region r for all subjects. We choose a non-informative prior for Σ by letting Ψ 0 in (1) be a diagonal matrix with elements Ψ 0rr = max s 2 r , r = 1, . . . , R, and a low degree of freedom ν 0 = R + 2. Following Litterman (1986), it is common practice in the BVAR literature to impose heavier shrinkage on higher lag orders. To implement this effect we let P −1 0 = λD and P −1 s = κ s D for each subject s, where the diagonal elements of D for lag l are given by l 2s2 r −1 , wheres 2 r is the mean of the subjects' sample variances in region r.
The values of λ and κ s are obtained from an empirical Bayes approach by maximizing the analytical, tractable, marginal likelihood of Y in Model 2. This is expected to be a good approximation to the optimal hyperparameters for Model 1 since λ and κ s are not related to Σ or Σ s , which is the aspect that differs between Models 1 and 2. Integrating out (B, Σ) from the posterior distribution in Equation (5), the marginal likelihood of the data Y can be written as a function of λ and κ = (κ 1 , . . . , κ S ) as where c 2 does not depend on κ and λ. Optimizing this function with respect to λ and κ, gives the estimated values of the hyperparameters in the prior precision matrices P 0 and P s of all models, respectively.

Brain Connectivity in resting-state fmri data
In Section 4.1, we describe the data used for group comparisons between healthy controls and individuals diagnosed with ASD. Effective and functional connectivity results are presented and compared between the models in Section 4.2. In Section 4.3, we present a brief overview of computational time requirements for analyzing the data with different number of time lags and number of regions considered.

Description of data and ROI selection
We use data from ABIDE 1 (Di Martino et al., 2014) preprocessed 2 (Craddock et al., 2013) consisting of resting state fMRI data from 539 individuals diagnosed with ASD and 573 healthy controls; we use randomly selected subsets of 20 controls and 20 ASD patients from the data collected at New York University. The fMRI data were collected using a 3 T Siemens Allegra scanner using a TR of 2 seconds. Each fMRI dataset contains 180 time points. No motion scrubbing has been performed, but the first four volumes were dropped in the processing to obtain 176 time points. The ABIDE Preprocessed data have been processed with four different pipelines, and we use the data from the CCS (connectome computation system) pipeline here. We use the data preprocessed without global signal regression and without bandpass filtering, as bandpass filtering will substantially change the autoregressive structure and we prefer to model it. Interested readers are referred to ABIDE preprocessed for preprocessing details. As all the preprocessed data are freely available, other researchers can reproduce our findings. We select the ROI:s guided by Easson, Fatima and McIntosh (2019), as their rs-fMRI dataset also included ASD patients and healthy controls. We include regions belonging to networks which are active during resting-state scans for both groups, as well as there being some indication of between-group differences in network configuration. The present analyses use ten regions (five in each hemisphere) belonging to the Default-Mode Network (DMN) and ten regions (also here five in each hemisphere) belonging to the Sensory-Motor Network (SMN). More details on the selected regions are given in Appendix D. To make graphs more readable, we refer to the 20 brain regions by numbering them as R1-R20 instead of naming them in the graphs (see Appendix D for a full list of region locations). Regions R1-R10 belong to the DMN and regions R11-R20 to the SMN.

Results on effective and functional connectivity
We present results on both effective (EC) and functional (FC) connectivity for the three models in Section 2, using the data described in Section 4.1. The EC and FC results are presented for two time lags (L = 2) for each BVAR model, which is the optimal number of lags for Model 1 by the widely applicable information criterion (WAIC, Vehtari, Gelman and Gabry (2017)), see Table 1. Models with two time lags for rs-fMRI data have also been suggested previously in the literature (Chiang et al. 2017, Gorrostieta et al. 2012. In addition, note that Model 1 is superior to Model 2 and 3 for each lag and group with substantially lower values of WAIC. Hence, the results clearly suggest that the heterogeneous subject-specific covariance matrices in Model 1 are indeed needed for modeling this data. . Results on EC corresponds to posterior inference on the AR-coefficients in our Bayesian VAR models. For illustration purposes, we apply thresholds on the posterior distribution of the AR-coefficients, where the threshold is applied on both size of the coefficient and whether a 90 or 95 % credible interval includes zero or not. We show results on EC for each of the ASD and control groups separately, as well as the group differences in EC at different time lags. Results for the separate groups are seen in Figures 1 (control group) and 2 (ASD group). Each subgraph gives a visual network description of the thresholded directed connections between regions. The general pattern across groups and models is that there are considerably more connections at one time lag than at two time lags, despite having a stricter threshold at one time lag. There are more connections within a network than between networks, as expected for the definition of the network, and most of the coefficients are positive. For a given time lag, there are some differences between the models. Model 2 and 3 yield more connections than Model 1 for lag 1, while Model 1 yields some additional, mostly negative, connections compared to Model 2 and 3. Figure 3 illustrates differences in EC between the groups as the difference in corresponding AR-coefficients. Results from Models 2 and 3 indicate substantially more group differences than Model 1 for both time lags, while for a given model the number of differences is greater for lag 1. In the figure, the differences for the two lags look comparable, but note that we use a more leniant threshold for lag 2 for illustrative purposes (otherwise there would have been only one connection for the differences of lag 2). In general, most of the differences are within-network, as for the group-specific connections. There are also some discrepancies in EC between Models 2 and 3, but much less than the corresponding discrepancies between any of these models and Model 1.
Results for FC are shown in Figure 4. We only compare Models 1 and 2, since Model 3 has a diagonal covariance matrix and therefore does not allow estima- tion of FC. We follow the same procedure as for EC by considering FC for each group separately, as well as the difference in FC between groups. The figures were constructed in a similar manner as for the figures for EC, but FC is undirected such that connection lines between regions are undirected. It is clear that the number of functional connections is considerably different between the models in the figure, where Model 2 yields many more functional connections than Model 1. Hence, it is important to account for heterogeneous subject-specific covariance matrices in Model 1 compared to a common covariance matrix in Model 2 in order to obtain accurate inference on functional connectivity. The overestimation of functional connections in Model 2 also implies an overestimation of group differences, resulting in many spurious functional connections between the groups from Model 2. We elaborate further on comparing results from the different models by pairwise comparisons of the posterior means of the parameters. In Figure 5, ARcoefficients and the elements of the FC correlation matrix are compared between the models. The comparison between Model 1 and 3 was omitted as it is virtually undistinguishable from the Model 1-Model 2 comparison. We only present this comparison for the control group, as the results for the ASD group were very similar. The posterior means of the AR coefficients from Models 1 and 2 are quite close (upper two panels in the left column of Figure 5). The posterior standard deviations of the respective coefficients are slightly higher for Model 2 compared to Model 1 (upper two panels in the middle column), but the t-ratios are quite similar (upper two panels in the right column). Moving on to the bottom row, the posterior means of FC correlations are fairly similar, although Model 1 yields slightly lower posterior mean correlations compared to Model 2. This means that the differences in FC between the Models are mainly due to the standard deviation of the respective posterior distributions being much smaller for Model 2, which is evident in the t-ratios as well. Table 2 presents computation times for applying the HMC algorithm to the posterior of Model 1 for different numbers of brain regions; the regions are randomly selected for this purpose, since interest here is only on the computing times. Our approach works fine for around 20-30 brain regions compared to previous VAR modeling of 5-6 regions, e.g. Chiang et al. (2017); Gorrostieta et al. (2012Gorrostieta et al. ( , 2013, but becomes very computational demanding for more than p = 30 regions and at least two time lags.

Conclusions
We propose a novel Bayesian VAR hierarchical model for brain connectivity that accounts for autoregressive dependencies within and between brain regions, and apply it to an existing, openly available rs-fMRI data set. Compared to existing Bayesian VAR hierarchical models for this purpose, we incorporate more  Table 2 Computation times (hours) for Model 1 with different lag lengths L and number of regions p. The data consists of S = 20 randomly selected healthy controls and the length of each time series within each region is T = 176. The computer analyses were run using two parallel MCMC chains with 200 warm-up and 500 sampling draws for each chain on two CPUs with a 2.7GHz processor. N/A (not applicable) means that the HMC algorithm took too long time to converge and was therefore interrupted.
flexible, subject-specific, covariance modeling that estimates both effective and functional connectivity simultaneously. By using the information criteria WAIC, we show that our proposed model is superior to special cases of our model with a common covariance matrix for all subjects. Similar simplified, diagonal, covariance matrices has been used previously for effective connectivity by Chiang et al. (2017) and Gorrostieta et al. (2012Gorrostieta et al. ( , 2013. We are also able to handle 20-30 brain regions compared to previous VAR modeling in the literature of typically 5-6 regions.
Overall, our flexible model displayed the most conservative results with respect to the number of effective and functional connections deemed to be nonzero from common thresholds. This is especially true for functional connectivity where the special case with a common covariance matrix for all subjects substantially overestimates the number of non-zero connections. We find that the standard deviations of the corresponding posteriors to functional connectivity are much lower for the special model case, which implies that between-subject variation is underestimated in such models.
We suggest some future extensions of our work. Flexible VAR modeling for effective and functional connectivity is not yet ready for large-scale brain connectivity, which typically involves hundreds of brain regions. Our derived, analytical result for the posterior inference of the model with a common covariance matrix for all subjects can be directly applied, but with the obvious drawback of biased inference for effective and, especially, functional connectivity. Another possibility can be to extend our modeling to handle longitudinal data, e.g. repeated rs-fMRI scans over time.
Appendix A: Posterior distribution of (B, Σ) for Model 2 with common covariance matrix Σ Let X s,t = (Y s,t−1 , . . . , Y s,t−L ) be the set of covariates in the BVAR model. The likelihood function of (B s , Σ) for each subject s is given by Completing the squares of B s , the likelihood function of (B s , Σ) can be written as

The marginal likelihood function of (B, Σ) for all subjects becomes
and c 0 is a constant that does not depend on B and Σ. Rewriting this marginal likelihood on a quadratic form of B and then multiplying the marginal likelihood with the prior distribution of (B, Σ), the posterior distribution of (B, Σ) becomes Rewriting on a quadratic form of B, the posterior distribution is finally given by Appendix B: Marginal likelihood of (B, Σ, ν) for Model 1 with subject-specific covariance matrix Σ s Replacing Σ with Σ s in Equation (4) gives the likelihood function of (B s , Σ s ) for each subject s as The marginal likelihood function of (B, Σ) for all subjects becomes where c 1 is a constant that does not depend on B and Σ. Rewriting on a quadratic form of B in the determinant, the marginal likelihood of (B, Σ) can be expressed as Sum logdet = Sum logdet + log determinant(Part s); } target += S*( lmgamma(p,0.5*(T+nu)) -lmgamma(p,0.5*nu) ); target += 0.5*S*nu*log determinant(nu*Sigma) -0.5*(nu+T)*Sum logdet; }

Appendix D: ROI information
Information on our selected ROI:s in the Default-Mode Network (DMN) and Sensory-Motor Network (SMN) is given below in the following order: abbreviation in the manuscript, type of network the ROI is classified to, volume of the ROI, (x, y, z)−coordinates for the ROI center of mass, and AAL atlas annotation (Tzourio-Mazoyer et al. 2002).  Control group EC: Effective connections, measured by AR coefficients, for the first (left column) and second (right column) time lag for each of our three models (rows one, two and three, respectively). The regions in the DMN are R1-R10 (brown), regions R11-R20 (orange) are the SMN regions. Thresholds for the first time lag were set to the following: connections for which all posterior draws lie above or below zero were kept. For the second time lag we retained connections whose 95% credible intervals do not include zero. The arc length for each region corresponds to the sum of the absolute value of coefficients for the model with the highest sum after thresholding. Direction of connections are indicated with an arrowhead towards the "receiving" region, thickness of the arrow indicate connection strength (coefficient size), and the color of the arrow indicate a positive (red) or negative (blue) coefficient.  ASD group EC: Effective connections, measured by AR coefficients, the first (left column) and second (right column) time lag for each of our three models (rows one, two and three, respectively). The regions in the DMN are R1-R10 (brown), regions R11-R20 (orange) are the SMN regions. Thresholds for the first time lag were set to the following: connections for which all posterior draws lie above or below zero were kept. For the second time lag we retained connections whose 95% credible intervals do not include zero. The arc length for each region corresponds to the sum of the absolute value of coefficients for the model with the highest sum after thresholding. Direction of connections are indicated with an arrowhead towards the "receiving" region, thickness of the arrow indicate connection strength (coefficient size), and the color of the arrow indicate a positive (red) or negative (blue) coefficient. Difference in EC, Control group EC -ASD group EC: Effective connection difference, measured by AR coefficient difference, for one (left column) and two (right column) time lags, for each of our three models (rows one, two and three, respectively). The regions in the DMN are R1-R10 (brown), regions R11-R20 (orange) are the SMN regions. Thresholds for the first time lag was set to the following: connections whose 95% credible intervals do not include zero, were kept. For the second time lag we retained connections whose 95% credible intervals do not include zero. The arc length for each region corresponds to the sum of the absolute value of coefficient difference for the model with the highest sum after thresholding. Direction of connections are indicated with an arrowhead towards the "receiving" region, thickness of the arrow indicate size of coefficient difference, and the color of the arrow indicate a positive (red) or negative (blue) coefficient difference. FC for both groups and the group difference in FC, measured by the estimated FC correlation matrix between regions and the FC correlation matrix difference between regions (columns) for Models 1 and 2 (rows). The regions in the DMN are R1-R10 (brown), regions R11-R20 (orange) are the SMN regions. The following thresholding was used in the plot: for the separate groups, the absolute posterior mean correlation was thresholded at 0.35, and all posterior draws lie above or below zero. When considering the group difference, the absolute posterior mean difference had to exceed 0.05 and the 95% credible intervals exclude zero. The arc length for each region corresponds to the sum of the absolute value of the correlations -or the correlation differences -for the model with the highest sum after thresholding. Thickness of the line indicate size of correlation/correlation difference, and the color of the arrow indicate a positive (red) or negative (blue) coefficient/coefficient difference.