Learning Robust Hierarchical Patterns of Human Brain across Many fMRI Studies
Year: 2021

Abstract:
Multi-site fMRI studies face the challenge that the pooling introduces systematic non-biological site-specific variance due to hardware, software, and environment. In this paper, we propose to reduce site-specific variance in the estimation of hierarchical Sparsity Connectivity Patterns (hSCPs) in fMRI data via a simple yet effective matrix factorization while preserving biologically relevant variations. Our method leverages unsupervised adversarial learning to improve the reproducibility of the components. Experiments on simulated datasets display that the proposed method can estimate components with higher accuracy and reproducibility, while preserving age-related variation on a multi-center clinical data set.

heading: 1 Introduction
text: Multi-site fMRI studies have gained a lot of interest over the last decade [1, 2]. One reason for this is the necessity to evaluate a hypothesis in multiple settings/sites and make the hypothesis result generalizable to a diverse population. Also, the pooling of data is essential when studying rare disorders or neurological conditions where the aim is to generalize the results to diverse populations [3, 4]. However, the data pooling often results in the introduction of non-biological systematic variance due to differences in scanner hardware and imaging acquisition parameters [5]. This additional variability can lead to spurious results and a decrease in statistical power. The variability can also result in hindrance in the estimation of true biological changes or in inferring non-biological differences as biological because of the correlation between site effects and biological predictors. Many studies working with multi-site data fMRI have reported considerable variability due site or scanner effects [6, 7, 1].
The non-biological variability introduced due to inter-scanner and inter-protocol variations can affect the estimation of the common features derived from fMRI [8], such as functional connectivity [5] or sparse hierarchical factors (this paper). These features were used for the study of the brain’s function during aging [9], of various neurological disorders [10, 11], and tasks [12]. The non-biological variability can considerably reduce these features’ reproducibility across different datasets and hence their utility as biomarkers for diseases that disrupt functional connectivity. Thus the removal of non-biological variance introduced by pooling of the data is essential for many neuroimaging studies.
Independent Component Analysis [13], Non Negative Matrix Factorization [14], Sparse Dictionary Learning [15] and Sparse factor Modeling [16, 17] are some of the commonly used methods for estimating interpretable components that capture the complexity of the human brain function at rest. However, these methods do not account for the hierarchical organization of the human brain [18, 19]. Recently, deep learning based methods [20, 21] have been developed to estimate hierarchical networks. However, they fail to capture overlapping components and the existence of positively and negatively correlated nodes in a component that are shown to capture the underlying structure better [22, 23]. Hierarchical Sparse Connectivity Patterns (hSCPs) overcome these limitations by
35th Conference on Neural Information Processing Systems (NeurIPS 2021)
capturing interpretable hierarchical, sparse, and overlapping components using a linear deep matrix factorization approach. The method not only captures shared patterns but also the subject-specific weights of these patterns, thus capturing the heterogeneity of brain activity patterns across individuals. In this paper, we focus on robust estimation of hSCPs [23–25] in a multi-site regime.
Many existing methods to reduce site effects are based on an empirical Bayes method ComBat [26], which was developed to remove ‘batch effects’ in genetics and has been applied for harmonizing different measures derived from structural [27, 28] and functional MRI [8]. But these harmonization methods can not be directly applied to the hSCP method because of the loss in the structure of captured heterogeneity (more details in section 2.1). See Appendix A for extended discussion on related work.
Contribution. In this paper, we develop a new model that is robust to site-effects in the estimation of sparse hierarchical connectivity pattern components (rshSCP). For this, the method learns site-specific features and global space, storing the information about the scanner and site, and uses these features to reduced site effects in the components. We also use an existing adversarial learning approach [25] on top of our method to improve the reproducibility and generalizability of the components across components from the same site. We formulate the method as a non-convex optimization problem which is solved using stochastic gradient descent. Experiments on simulated and real datasets show that our method can improve the split-sample and leave one site reproducibility of the components while retaining age-related biological variability in the data, thus capturing informative heterogeneity.

heading: 2 Method
text: In this section, we first present a brief overview of the hierarchical Sparse Connectivity Pattern (hSCP) and its adversarial formulation. We then describe our main method and the joint formulation incorporating adversarial learning. We follow the notation in [24] and [25]. The set of symmetric positive definite matrices of size P × P is denoted by SP×P++ . Matrix A with all the elements greater than or equal to 0 is denoted by A ≥ 0. JP denotes P × P matrix with all elements equal to one. P × P identity matrix is denoted by IP and element-wise product between two matrices A and B is denoted by A ◦B.

heading: 2.1 Introduction to hierarchical Sparse Connectivity Patterns
text: Hierarchical Sparse Connectivity Patterns (hSCP) [24] is a hierarchical extension of Sparse Connectivity Patterns (SCPs), first defined by [23] to estimate sparse functional patterns in the human brain using fMRI data. Let there be N number of subjects or participants, and each subject’s BOLD fMRI time series has T time points and P nodes representing regions of interest. The input to hSCP are correlation matrices Θn ∈ SP×P++ where ith and jth element of the matrix is the correlation between time series of ith and jth node. hSCP then outputs a set of shared hierarchical patterns following the below equations:
Θn ≈W1Λn1W>1 , . . . Θn ≈W1W2 . . .WKΛnKW>KW>K−1 . . .W>1 ,
where Λnk is a diagonal matrix having positive elements storing relative contribution of the components for the nth subject at kth level, K is the depth of hierarchy and P > k1 > . . . > kK i.e. each successive level in the hierarchy has less number of components than the previous one. In the above formulation, W1 ∈ RP×k1 stores k1 components at the bottom most level, and each successive multiplication by W2, W3, . . ., WK linearly transforms to a lower dimensional space of k2, k3, . . ., kK dimension. LetW = {Wr | r = 1, . . . ,K} be the set storing sparse components shared across all subjects and D = {Λnr | r = 1, . . . ,K;n = 1 . . . , N} be set storing subject specific diagonal matrix with Λnr ≥ 0. The hierarchical components are estimated by solving the below optimization problem:
min W,D H(W,D, C) = N∑ n=1 K∑ r=1 ‖Θn − ( r∏ j=1 Wj)Λ n r ( r∏ j=1 Wj) >‖2F
s.t. ‖wrl ‖1 < τr, ‖wrl ‖∞ ≤ 1, trace(Λnr ) = 1, Λnr ≥ 0, ∀n, r, l; Wj ≥ 0, j = 2, . . . ,K,
(1)
where C = {Θn | n = 1, . . . , N}, l = 1, . . . , kr, n = 1, . . . , N and r = 1, . . . ,K. L1, L∞ and trace constraints help the problem to identify a decomposition which can provide reproducible components. More details can be found in the original paper [24]. We will be denoting above constraint set as ΩW = {W | ‖wrl ‖1 < τr, ‖wrl ‖∞ ≤ 1,Wj ≥ 0, j = 2, . . . ,K} and Ψ = {Λ | trace(Λnr ) = 1,Λ n r ≥ 0}. A detailed description of the method can be found in [24].
Can we use standard harmonization approaches? These methods reduce site effects by adjusting for additive and multiplicative effects for each feature in data separately and use emperical Bayes estimates the model parameters. These methods can be used in the case of hSCP in two ways. First, harmonization can be directly applied to each element of the correlation matrices, which is the input of hSCP. This will reduce site effects from each element of the correlation matrix, thus from the complete input, but the final matrix that does not necessarily follow the essential property of a correlation matrix i.e., positive definiteness. For similar reasons, COMBAT can not be directly applied to time series; if applied, it can change the inference derived from the correlation matrix. Second, harmonization can be directly applied to Λ to remove site effects. To understand this, we look at the hSCP formulation at one level:
Θn ≈ k∑ l=1 dnl wlw > l ≈WΛnW>,
where dnl are non-zero elements storing the subject-specific information, which can be affected by the variability introduced by the site. In this model, harmonizing each feature across different sites will change the relative contribution of the components in each subject’s functional structure, which will remove the interpretability of the subject-wise weights, which is not desirable. Instead, a two step optimization procedure can be used to incorporate ComBat with hSCP (ComBat hSCP). We first run hSCP and use ComBat on the extracted Λn to get harmonized subject specific information ∆n ∈ Rk1×k1 for each subject. We then re-fitted W using the below decomposition-
Θn ≈W(∆n + S)W>.
We added a diagonal shift matrix S ∈ Rk1×k1 such that ∆n + S is positive for each subject and performed the optimization to estimate W and S. We show through experiments that this baseline two step optimization procedure is not optimal and performs worse than vanilla hSCP.

heading: 2.2 Adversarial Learning in hSCP
text: Adversarial learning has shown to achieve state-of-the-art performance of various matrix factorization approaches [29, 30]. Recently, Sahoo and Davatzikos [25] demonstrated that incorporating adversarial learning in the estimation problem of hSCP can improve the reproducibility of the hierarchical components. The method is based on perturbation of input data Θn to learn stable components robust to adversaries. First the input data is perturbed to generate new data Γn = Θn + σJP where σ is the standard deviation of the data. We consider this perturbation as rank one perturbation which will transform eigenvalues of each subject’s correlation matrix differently. Since the addition of ones does not change the symmetric property of the matrix, we just scale the matrix such that the diagonal matrix contains 1 and by Weyl’s inequality about perturbation [31] it will be positive definite. We have experimented with randomly positively scaled rank one perturbation, but the results were worse than the vanilla hSCP. Also, using rank one perturbation, it is easier to control eigenvalues (and noise added) of the modified matrix than using rank k perturbation. The perturbed set of components W̃1 are then estimated using the new data and are ensured to be close to W1 by solving the below minimization problem:
A(Ŵ1) = α‖W̃1 −W1‖2F + N∑ n=1 ‖Γn − W̃1ΛnW̃>1 ‖2F . (2)
In the above optimization problem, the first term constraints W̃1 to be close to W1 and second term is used for learning the components using the perturbed data. Aim of the learner is to estimate Λn and W1 by minimizing the below cost function:
D(W,Λ) = N∑ n=1 ‖Θn − W̃1ΛnW̃>1 ‖2F + β N∑ n=1 ‖Θn −W1ΛnW>1 ‖2F , (3)
for a fixed W̃1. Learner first estimates subject specific information Λn using perturbed weight matrix and then use it to learn W1. Below is the complete optimization problem for adversarial learning of hSCP:
min W,D H(W̃,D, C) + βH(W,D, C)
s.t. W̃r = arg min Ŵr α‖Ŵr −Wr‖2F +H(W̃,D,P) r = 1, . . . ,K
W̃r,Wr ∈ Ω, D ∈ Ψ,
(4)
where W̃ = {W̃r | r = 1, . . . ,K} and P = {Γn | n = 1, . . . , N}. We next discuss our main method which aims to reduce the site effects.

heading: 2.3 Robust to site hSCP
text: Estimating hSCPs in multi-site data can introduce non-biological variances in the components and the subject-specific information. One of the typical approaches would be to use harmonization methods mentioned previously, but it would lead to a loss in the structure of these features, which in turn will lose interpretability. Instead of removing site effects after estimating the components, we jointly model the sparse components and the site effects, and estimate robust to site hSCP (rshSCP). We first look at the case when there is only one level in the hierarchy, which then can be extended to multiple levels. Let there be total S sites, Is be the set storing subjects from site s and y ∈ RN×S be the one-hot encoded site labels We hypothesize that there is a space V ∈ RP×P storing site and scanner information for all the possible available data, and for each site s, we have space Us ∈ RP×P storing site-specific information for s = 1, . . . , S. Based on the above hypothesis, we decompose the correlation matrix Θn of n ∈ Is to jointly estimate the hSCPs, Us and V as:
Θn ≈ WΛnW>︸ ︷︷ ︸ decomposition of
subject components + UsV︸ ︷︷ ︸ decomposition of site components , (5)
where Us is constrained to be a diagonal matrix and L1 sparsity constraint is used for V to prevent overfitting. In addition to estimating the site effects, we train the model such that the predictive power of subject-specific information D for predicting site is reduced, which can assist in removing site information. For this, we train a differentiable classification model F (ζ,D) parameterized by ζ with input Λn that return site predictions ŷ ∈ RN×S . These predictions indicate the probabilities that each of N inputs belongs to each of S site labels. The classification model is trained by optimizing for ζ such that the cross-entropy loss L(ζ,D,y) between the predictions ŷ and the true site labels y is minimized:
ζ∗ = arg min ζ − 1 N N∑ n=1 S∑ s=1 yn,s log ŷn,s. (6)
Using this model, we modify Λn such that its predictability power reduces. We achieve this by maximizing the above loss with respect to D. This will result in a minimax game, where the classifier learns to minimize the cross-entropy or the surrogate classification loss, andD is adjusted to maximize the loss. The joint optimization problem can be written as:
max ζ min W,D,U,V
S∑ s=1 ∑ n∈Is ‖Θn −WΛnW> −UsV‖2F − γL(ζ,D,y)
s.t. W ∈ Ω, D ∈ Ψ, ‖vp‖1 < µ, p = 1, . . . , P,
(7)
where U = {Us|s = 1, . . . , S}, vp is the pth column of V.

heading: 2.4 Complete Model
text: We can combine the above formulation (7) at multi level with the adversarial learning (4) to jointly model hSCPS and site effects. Let
G(W,D, C) = S∑ s=1 ∑ n∈Is K∑ r=1 ‖Θn − ( r∏ j=1 Wj)Λ n r ( r∏ n=1 Wn) > −UsrVr‖2F , (8)
then the joint optimization problem can be written as:
max ζ min W,D,U,V J(W̃,W,D, C) = G(W̃,D, C) + βG(W,D, C) + γL(ζ,D,y)
s.t. W̃r = arg min Ŵr α‖Ŵr −Wr‖2F +G(W̃,D,P) r = 1, . . . ,K
W̃r,Wr ∈ Ω D ∈ Ψ, ‖vp‖1 < µ, p = 1, . . . , P,
(9)
The optimization problem defined above is a non-convex problem that we solve using alternating minimization. Complete algorithm and the details about the optimization are described in Appendix B. Note that the random initialization of the variables can result in a very different final solution that might be far from the ground truth. One such solution for U and V would be the identity matrix since all the correlation matrices have one as their diagonal element, which can drastically change the final components. It might also be possible that V might store highly reproducible components since they are present in most individuals, leading to a decrease in reproducibility of hSCPs. We prevent these cases by using svd− initialization [24, Algorithm 2] forW and D, where, in the starting, most of the variability associated with data is stored inW and D. In this way, we can prevent V from storing highly reproducible components during initial iterations. We initialize Us and V using the below equation:
Usr =  1 |Is| ∑ n∈Is Θn − ( r∏ j=1 Wj)Λ n r ( r∏ n=1 Wn) > Jp  ◦ Ip, Vr = 1 P JP . (10)
This complete initialization procedure ensures that the algorithm starts with the majority of variability in the data stored in W and D, and Us start from the residual variance left in site s after the svd− initialization procedure. We show in the next sections that this simple strategy, though suboptimal, can help estimate reproducible components with diminished site effects. All the code is implemented in MATLAB and will be released upon publication.

heading: 3 Experiment

heading: 3.1 Simulated Dataset
text: One level. We first generate simulated dataset at one level to evaluate the performance of our model against the vanilla hSCP. We simulate data with p = 50, k1 = 10, S = 4 with 200, 300, 400 and 500 number of participants in each site. We generated sparse shared components W1 with percentage of non-zeros equal to 60% and each element sampled from N (0, 1). We then generate correlation matrix for nth subject belonging to sth site using:
Θn = (W1 + E n 1 ) Λ n (W1 + E n 2 ) > + UsV + En2 , (11)
where Us is a diagonal matrix with positive elements sampled from N (1, .1), V is a random matrix sampled from wishart distribution, each element of Λn is sampled from N (4, 1) and En1 is the noise matrix added to the components whose each element is sampled from N (0, .1) and En2 is added to ensure that the final matrix is positive definite. However the diagonal elements of Θn are not equal to 1. To make them 1, we extract diagonal elements D of Θn and get the new correlation matrix as D1/2ΘnD1/2. We used a feed-forward neural network for the classification model with two hidden layers. The networks contain the following layers: a fully connected layer with 50 hidden unites, dropout layer with rate 0.2, ReLU, a fully-connected layer with 4 hidden units and a softmax layer. Optimal value of hyperparameters α, β, µ and τ1 are selected from [0.1, 1], [1, 5], [0.1, 0.5, 1] and 10[−2:2]. The criterion for choosing the best hyperparameter is maximum splitsample reproducibility. The split sample reproducibility is the normalized inner product between the components estimated on two random equal splits of the data. Split sample reproducibility tries to answer the question of whether the components are generalizable across subjects from the same sites or not. We compared different methods for estimation of hierarhical components- hSCP, ComBat hSCP, hSCP with adversarial learning (Adv. hSCP), rshSCP, rshSCP with adversarial learning (Adv. rshSCP), rshSCP and Adv. rshSCP with random initialization (rshSCP w/ rand. and Adv. rshSCP w/ rand.). Table 2 shows the reproducibility of the components generated from different methods. It is computed over 15 runs in all the experiments. We used accuracy of the estimated components as a
performance measure. It is defined as the normalized inner product between ground truth components and estimated components. All the experiments were run on a four i7-6700HQ CPU cores single ubuntu machine.
Accuracy. Table 1 displays the accuracy of different methods on the simulated dataset. Here, accuracy is defined as the average correlation between estimated components and the ground truth components. From the results, we can see that the rshSCP with adversarial learning can significantly improve the components’ accuracy and the reproducibility of the components. The baseline (ComBat hSCP) performs worse than vanilla hSCP. One reason for this might be that the harmonized Λ extracted using ComBat might not necessarily result in optimal highly reproducible W. This result bolsters our method that we need a joint optimization procedure to obtain W and Λ with reduced site effects. The results using random initialization instead of using the initialization strategy mentioned in the previous section indicates that random initialization brings significant variability to performance. On average, it performs worse than our strategy, but there might be instances where the random initialization can perform better, which might suggest that there might be some better strategy for initialization. Also, for V, there is an optimal sparsity value, which achieves the best result. If V is dense, then it might remove essential information that might reduce reproducibility, and if it is too sparse, then we might not have desired effects to make the model robust. The results showing the variation in the accuracy and reproducibility with the sparsity of V are in Appendix D.1.
Site prediction. To check if the estimated subject information (Λ) has reduced predictive power to predict the site to which the subject belonged, we performed a 5 fold cross-validation using SVM with RBF kernel. We also ran our experiment using two different feed forward networks with two different architectures: (a) a fully connected layer with hidden units, dropout layer with the rate
, ReLU, a fully-connected layer with hidden units and a softmax layer and (b) a fully connected layer with hidden units, dropout layer with the rate , ReLU, a fully-connected layer with hidden units, dropout layer with rate , ReLU, a fully-connected layer with hidden units and a softmax layer. Our model leads to a decrease in average cross-validation accuracy from 97.6% to 67% for SVM, 98.1% to 67.3% for neural network with architecture (a) and 98.2% to 66.9% for neural network with architecture (b). This suggests that our model can reduce the prediction capability to predict site. More details about the result are in Appendix D.1.
Two level. Under the same settings as defined above, we generate correlation matrix from two level components with k2 = 4 using:
Θn = W̃1W̃2Λ nW̃>2 W̃ > 1 + U sV + En3 ,
W̃1 = W1E n 1 , W̃2 = W2 + E n 2
(12)
where each element of W2 is sampled from N (0, 1), the percentage of non-zeros equal to 40%, En1 and En2 is the noise added to the components whose each element is sampled from N (0, .1) and En3 is added to ensure that the final matrix is positive definite. Table 3 shows the accuracy for different values of k1 and k2. Selection of hyparameter is same as in the previous paragraph. We can see that the proposed method estimates most accurate ground truth components. Reproducbility results and random initialization results are availble in Appendix D.1.

heading: 3.2 Real Dataset
text: Data. We collected functional MRI data from 5 different multi-center imaging studies- 1) Baltimore Longitudinal Study of Aging (BLSA) [32, 33], the Coronary Artery Risk Development in Young Adults study (CARDIA) [34], UK BioBank (UKBB) [35], Open access series of imaging studies (OASIS) [36] and Aging Brain Cohort Study (ABC) from Penn Memory Center [37]. Although UK Biobank has more than 20000 scans, we only used 2023 randomly selected scans to avoid estimating the results that would be heavily influenced by the UK Biobank. We projected the data into a lower-dimensional space such that the number of nodes in each subject’s data was 100. Table 4 summarizes the number of participants in each site and age distribution. More details on the dataset and the preprocessing pipeline are given in Appendix C. CARDIA data is divided into three parts because of the acquisition at three different sites.
Reproducibility. Since we don’t have access to ground truth here, we compare the methods based on the split sample and leave one site reproducibility. Leave one site out reproducibility is defined as the similarity between components derived from the site s and all sites except s. Split sample reproducibility tries to answer the question of whether the components are generalizable to other sites or not. For estimating rshSCP with only one site, we used V estimated from all sites except s since the idea behind V was to store information about the site/scanner from various sites. This would also help analyze the generalization power of V. The optimum value of the hyperparameters is selected from the range defined in section 3.1. τ1 and τ2 are selected from 10[−2:2] based on maximum split-sample reproducibility. The criterion for choosing the best value is the maximum split sample reproducibility. Table 5 shows the split sample reproducibility for varied values of k1 and k2 = 4. Leave one site out reproducibility results are shown in Table 6. Table 13 and Table 14 in D.2 shows split sample reproducibility and leave one site out reproducibility respectively at two-level for k2 = 6. The results demonstrate that the proposed method can significantly improve the split sample reproducibility and leave one site out reproducibility. For the remaining paper, we focus on
the comparison between components learned using adversarial learning from hSCP and Adv. rshSCP.
Site prediction. We performed the same experiment under the same settings as mentioned in the previous section to check Λ has reduced predictive power to predict the site. Using SVM, our model leads to a decrease in average cross-validation accuracy from 51% to 32%. Using the first neural network architecture defined in section 3.1, the cross-validation accuracy for hSCP model is 59.3% and for the rshSCP is 33.6%. Using the second architecture, the cross-validation accuracy for hSCP model is 58.7% and for the rshSCP is 33.4%. This suggests that our model can reduce the prediction capability to predict site.
Age prediction We used subject specific information (Λ) to predict age of each subject as the metric to check if our method is able to preserve age related biological variation. Table 7 shows 10 fold cross validation mean absolute error (MAE) using random forest . From the table, we can see that the proposed method has comparable performance as the hSCP suggesting that it preserves age related biological variance. More details about the experiment are given in Appendix D.2.

heading: 3.3 Analysis of components
text: A robust method should be able to reduce non-biological variability caused by site and scanner while retaining biological variability. In this study, we look at brain aging-related associations and leave analysis with other variables for future work. We also discuss the difference between the components with and reduced site effects. We selected the subjects with age greater than 60 to find an association between brain aging and the components derived from hSCP and rshSCP. The total number of subjects having an age greater than 60 is 2746. We first computed Spearman correlation between subject-specific information (Λ) and their respective age and used 0.05 as the significance level for the hypothesis test of no correlation against the alternative hypothesis of a nonzero correlation. We derived 10 fine-scale (1 − 10) and 4 coarse-scale components (I-IV) because of the high split sample reproducibility and easier interpretation of each component. The correlation and p-values are displayed in Table 15 and 16 in Appendix D.2 for the hSCP and rshSCP.
We first compare components from the two methods. Figure 1 shows the components derived from hSCP and the proposed method. The first row of the figure displays the components with
anti-correlation between Default Mode Network (DMN) and Dorsal Attention Network (DAN). The component derived using hSCP has a part of the visual area positively associated with DMN, but the opposite is true, as shown by the previous sparse connectivity patterns [23]. On the other hand, the component with reduced site effects is cleaner since it does not include that relation. This component has a negative correlation with age which has been previously shown in resting-state fMRI and task-based fMRI [38]. The magnitude of anti-correlation has been connected to individual differences in task performances in healthy young adults [39]. However, in the case of older adults, the behavioral implications of reduced anti-correlation remain unclear. The second row of the Figure 1 displays another set of components for comparison. The components stores information about the anti-correlation between DMN and sensorimotor, which aligns with the previous literature [40]. But the addition of a positive correlation of DMN with visual areas will cause misleading inference since it contradicts the previous SCPs and studies. Hence making an inference without removing the site effects can be misleading. Discussion of the association between other relevant components and aging is in Appendix D.2 showing biological interpretability of the components in brain aging. From the results, we can see that there is an increase (or decrease in anti-correlation) in connectivity between different networks in the aging brain. This suggests that there is a reorganization of the aging brain aligning with the previous findings [41]. This can serve as a base to explore rshSCP as a biomarker of neurodegenerative diseases.
Figure 2 displays one of the hierarchical components with coarse-scale component storing relation between different fine-scale components comprising DMN, sensorimotor, and visual areas, previously studied by [40]. These findings give evidence that even after removal site effects, the components can have a meaningful interpretation. The results indicate that our approach can extract robust informative patterns without using traditional seed-based methods that are dependent on the knowledge of the seed region of interest.

heading: 4 Conclusion
text: In this work, we have presented a method for estimating site effects in hSCP. We formulated the problem as a minimax non-convex optimization problem and solved it using AMSgrad. We also propose a simple initialization procedure to make the optimization procedure deterministic and improve the performance on an average on a simulated and real dataset. Experimentally, using a simulated dataset, we showed that our method accurately estimates the ground truth compared to the vanilla method with better reproducibility. On the real dataset, we show that the proposed method can capture components with a better split sample and leave one site out reproducibility without losing biological interpretability and information. We also show that without removal of site effects, we can have a noisy estimate of sparse components resulting in misleading downstream analysis.
Below we mention some directions for future research. First, it would be interesting to consider the framework for the analysis of task-induced activity to investigate the extent of site effect and corrections on underlying networks activated by the task. Second, one could look at the changes in the associations of hSCPs with various clinical variables such as Mini-mental score, Digit Span Forward score, etc., after removing site effects. Third, we can also look low-dimensional modeling of V along with sparse constraints which has been used several robust matrix factorization problems. Since we have only shown age related biological preservation, future studies will focus on whether the proposed method preserves components associated with other demographic, clinical phenotypes, and pathological biomarkers.
There are few weaknesses of our proposed model, which also adds directions for future work. First, our method only captures linear site effects, it would be interesting to see if explicitly capturing nonlinear site effects can improve the performance of the model. Second, the result of the optimization algorithm depends on the initialization procedure, which has been shown to perform well on the simulated dataset and real dataset but can be sub-optimal.

heading: 5 Broader Impact
text: In this work, we provide a new method to diminish site effects in hSCPs and robustly estimate the components in multi-center studies. The formulation used in the method is not limited to hSCP. It can be easily extended to various matrix factorization approaches such as Independent Component Analysis, Non-negative Matrix Factorization, Dictionary Learning, etc., to improve the reproducibility of functional networks/components. Our work not only has broader applicability in terms of methods used for estimation of components but also to different types of neuroscience data, which includes EEG, MEG, etc.
Meanwhile, it should be noted that any unsupervised machine learning model has risk associated with it in terms of accuracy of the model; our work is no departure from this. It can lead to a negative impact since the reproducibility of the components is not 100%, which can lead to misleading clinical and biological interpretability of some of the components. This can be prevented by analyzing the components having high reproducibility as shown in section 3.2 and 3.3.

Review summaries:
This paper proposes a combined technique based on matrix factorization and adversarial models to reduce the batch effect problem in fMRI data. The new model, called rshSCP, is designed to be robust to the site-effects in estimation of sparse connectivity pattern components. The idea is to decompose correlation matrices into a set of shared hierarchical pattern of subject components as well as site related components. Then an existing adversarial model is applied in order to increase the reproducibility of them. They also applied a simple initialization method for optimization procedure and evaluated their approach on both simulated and real datasets.

This paper introduces the method site-robust hierarchical sparse connectivity patterns (rshscp) and its adversarial extension. It builds on sparse connectivity patterns and its adversarial extension introduced in references 22 and 23. These propose a factorization of correlation matrices extracted for functional connectivity analysis. The factorization consists of a symmetric product of  sparse factors around a central diagonal weighting matrix. The present contribution introduces a site-dependent component into the optimization problem and adds a classifier to predict site, which is then adversarially optimized to make the learned structure lose information about site. In synthetic and real-data experiments it is shown that the method is indeed capable of discarding site information while retaining other relevant information.   

This paper deals with the important problem of how to pool fMRI data from multiple sites for connectivity estimation. The authors propose a simple extension to adversarial hierarchical sparse connectivity pattern (hSCP) estimation to remove site effects. The argument against adopting methods, such as the widely-used COMBAT, is to preserve the manifold structure of connectivity matrices. The method is tested on simulated data with higher accuracy shown. The method is also tested on real data with higher reproducibility demonstrated and age effects retained.

This paper proposes a method called rshSCP (robust-to-site hSCP) to estimate hierarchical sparse connectivity patterns (hSCPs) from multi-site fMRI data while reducing the impact of site-specific variation due to scanner and experiment details rather than biological effects. Their approach is applied on top of the previous hSCP method as well as an extension that uses adversarial learning to improve reproducibility of estimated patterns (Adv. hSCP). The proposed method includes an additive component to the decomposition of the correlation matrix in hSCP, which models site effects as a matrix factorization with diagonal site-specific and full shared (across-site) matrix components. The method then additionally applies an adversarial classifier (a two hidden-layer feed-forward neural network) to predict the site from the subject-specific components of the hSCP decomposition. The entire optimization is formulated as a minimax problem, where the classifier maximizes its ability to predict site information while the remaining parameters are optimized to reconstruct the correlation matrices while minimizing the classifier's accuracy. The authors also propose an SVD-based initialization which improves reproducibility of the estimated components. Experiments on simulated data as well as fMRI data from multi-center imaging studies show that rshSCP improves accuracy (on simulated data) as well as reproducibility (on simulated and real data) of the estimated components, both when applied with hSCP and Adv. hSCP. The proposed method also leads to a reduction in accuracy when predicting site from the estimated subject-specific components using an SVM, as desired. Predicting age from the estimated subject-specific components with a random forest and visualizing the shared components shows that rshSCP preserves relevant biologically-related variation.