Learning Robust Hierarchical Patterns of Human Brain across Many fMRI Studies
Authors: Dushyant Sahoo
Christos Davatzikos
Venue: NeurIPS 2021 Poster (2021)
Keywords: Hierarchical Latent Factor Modeling
Matrix Factorization
Domain Adaptation
fMRI analysis

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. 

Learning Robust Hierarchical Patterns of Human Brain across Many fMRI Studies
Dushyant Sahoo Department of Electrical Engineering University of Pennsylvania saduÊseas. upenn. edu
Christos Davatzikos
Department of Radiology
University of Pennsylvania
christos davatzi upenn edu
Abstract
Multi-site fMRI studies face the challenge that the pooling introduces systematic non-biological site-specific variance due 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 multi-center clinical data set.
Introduction
Multi-site fMRI studies have gained a lot of interest over the last decade [1 One reason for this the necessity to evaluate hypothesis multiple settings/sites and make the hypothesis result generalizable to diverse population. Also, the pooling of data is essential when studying rare disorders or neurological conditions where the aim generalize the results 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 decrease in statistical power. The variability can also result hindrance in the estimation of true biological changes or in inferring non- non- biological differences 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 1.
The nonobiologiaal 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 [0 [] and tasks [12). The non-biological variability can considerably reduce these features' reproducibility across different datasets and hence their utility 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 [6 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 220 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 component that are shown to capture the underlying structure better Hierarchical Sparse Connectivity Patterns (hSCPs) overcome these limitations by
35th Conference on Neural Information Processing Systems (NcurIPS 2021)

---

capturing interpretable hierarchical, sparse, and overlapping components using linear deep matrix factorization approach. The method not only captures shared patterns but also the subject-specine weights of these patterns, thus capturing the heterogeneity of brain activity patterns across individuals. multi-site regime.
In this paper, we foeus on robust estimation of hSCPs [3:25)
Many existing methods to reduce site effects are based on an empirical Bayes method ComBat [26). which was developed to remove 'batch effects' in geneties and has been applied for harmonizing different measures derived from structural [27] (8 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 robust to site-effects in the estimation of sparse hierarchical connectivity pattern components (rshSCP). For this, the method learns site-specife 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 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 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.
2 Method
In this section, we first present a brief overview of the hierarchical Sparse Connectivity Pattern (hSCP) and 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 x P is denoted by P Matrix A with all the elements greater than or equal to 0 denoted by A 2 0. Jp denotes P x matrix with all elements equal to one. P x P identity matrix is denoted by p and element-wise product between two matrices A and B is denoted by A o B.
2.1 Introduction to hierarchical Sparse Connectivity Patterns
Hierarchical Sparse Connectivity Patterns (hSCP) 24 is a hierarchical extension of Sparse Connec- tivity 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 correlation matrices @ SPx 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:
e W'
D W W W W K-1 WI,
W
where Ag is a diagonal matrix having positive elements storing relative contribution of the components for the nith subject at kth level, K is the depth of hierarchy and P >> ki > >... > kk i.c. each successive level in the hierarchy has less number of components than the previous one. In the above formulation, W1 E IR Pxki stores k components at the bottom most level, and each successive multiplication by W2. W....., Wx linearly transforms to a lower dimensional space of le, ka..... kR dimension. Let W
= [W 1 r= 4,...K' be the set storing sparse components shared across
all subjects and D
|r ,,..,,KK K 1 N) be set storing subject specific diagonal matrix with A 2 0. The hierarchical components are estimated by solving the below optimization problem:
N
K
min W.D
H(W,D,C)
= Σ O
W,) II T 1=1
n=
1=1
(I)
s.t.
|w] 11 Tr, |l | 80 < 1, trace (A,) = 1,
A. 2 0, Vn,r,4
W, 2 0, j = 2.... K,
2

---

where C = {@" In=I.....N,111 l...,,k. kr == n .N and == I.....K K. L1. Loc and trace constraints help the problem to identify decomposition which can provide reproducible components. More details can be found in the original paper (2]] We will be denoting above constraint set as flw = (W ||w/|la < Tr, ||w/ll>> 1, W 2 0, j 2...... K) and V = {A
trace( (A)) 1, A 2 03. A detailed description of the method can be found in 22
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 A to remove site effects. To understand this, we look at the hSCP formulation at one level:
-" 2 dw & WA"WT, 1=1
where df 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 A to get harmonized subject specific information A n k xkı for each subject. We then re-fitted W using the below decomposition-
-t ~ W(A" + S.
We added a diagonal shift matrix S R such that A" + 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.
2.2 Adversarial Learning in hSCP
Adversarial learning has shown to achieve state-of-the-art performance of various matrix factorization approaches (0) 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 based on perturbation of input data O learn stable components robust to adversaries. First the input data is perturbed to generate new data r = O" + +JJp where o is the standard deviation of the data. We consider this perturbation as rank one perturbation which will transform cigenvalues of each subject's '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 and by Weyl's inequality about perturbation [1 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, are then estimated using the new data and are ensured be close to W, by solving the below minimization problem:
A A(W) = al Wi W
Σ
T T
W W,A"WIIR 12
In the above optimization problem, the first term constraints W1 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 A and W, by minimizing the below cost function:
N
D(W. Σ O -o-Wai W.A W
1
W A'
(3)
3

---

for fixed W1. Learner first estimates subject specific information A using perturbed weight matrix and then use it to learn WT. Below the complete optimization problem for adversarial learning of hSCP:
min W.D
H(W, + 3HH,,,,W
s.t.
W, arg aarggmin ar W, w. H(W,D,)) r=1...., K
(4)
W,
DEV, V. where W = {W, Ir- ,,..,,,KKK K} and P [T"h n I...,,N, method which aims to reduce the site effects.
Wr. W. EQ,
We next discuss main
2.3 Robust site hSCP
Estimating hSCPs in multi-site data can introduce non non- cooboiooglaaal 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 the case when there only one level in the hierarchy, which then can be extended to multiple levels. Let there be total S sites, I. be the set storing subjects from site and y E R. x be the one-hot encoded site labels We hypothesize that there a space V PxP storing site and scanner information for all the possible available data, and for each site we have space U E PPP storing site-specific information for s = I,...,,, S. Based on the above hypothesis, we decompose the correlation matrix @ of n E Ig jointly estimate the hSCPs, U and V as:
e ~ WA"W
U decomposition of site components
(5)
decomposition of subject components
where U is constrained to be 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-specifie information D for predicting site reduced, which can assist in removing site information. For this, we train differentiable classification model F(C, D) parameterized by 5 with input " that return site predictions y E xRNx 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(S, D,y) y) between the predictions ý and the true site labels y is minimized:
N
5" arg min
1
ΣΣ n=1 s=1
yn,s log in...
(6)
N
Using this model, we modify A x such that its predictability power reduces. We achieve this by maximizing the above loss with respect to D. This will result in minimax game, where the classifier learns to minimize the cross-entropy or the surrogate classification loss, and D is adjusted to maximize the loss. The joint optimization problem can be written as:
S
max
min W.D.U.V
Σ
Σ
1 - WA"WT - UVV - TR((,, D.y)
8=1 NEZ.
s.t.
W ER,
DEV,
|| < 14 p I..,,,
where 21 = {UA = 1, L...,,S,, S}. vp is the pth column of V.
2.4
Complete Model
We can combine the above formulation model hSCPS and site effects. Let
at multi level with the adversarial learning 4 to jointly
K
D(
1
II W W WT)T U.V.
(8)
#=1 NCI, r=1
j=1
n=1

---

then the joint optimization problem can be written as:
max
min
J(W, W, W.D, G(W, (G)))) + +2(())))))
W.D.V
s.t.
W, arg min o| W. C r=
(9)
W.
W W. DEV, 1vpll 14, p
The optimization problem defined above is non-convex problem that we solve using alternating minimization. Complete algorithm and the details about the optimization are described in Appendix Note that the random initialization of the variables can result in 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 most individuals, leading to a decrease reproducibility of hSCPs. We prevent these cases by using svd initialization (24] Algorithm 2] for W and D, where, in the starting, most of the variability associated with data is stored in W and D. In this way, we can prevent V from storing highly reproducible components during initial iterations. We initialize U and V using the below equation:
U
1
Yt II W W WT)
J I
V, -J=
PJP.
(10)
This complete initialization procedure ensures that the algorithm starts with the majority variability in the data stored in W and D, and U start from the residual variance left in site after the svd initialization procedure. We show in the next sections that this simple strategy, though sub- optimal, help estimate reproducible components with diminished site effects. All the code is implemented in MATLAB and will be released upon publication.
Experiment
3.1 Simulated Dataset
One level. We first generate simulated dataset at one level to evaluate the performance our model against the vanilla hSCP. We simulate data with p 50, k 10, S 4 with 200, 300, 400 and 500 number of participants in each site. We generated sparse shared components W, with percentage of non-zeros equal to 60% and each element sampled from ('0, We then generate correlation matrix for with subject belonging to sth site using:
θ = (W, + ET) A" (W1 + W )
UVV + E,
(11)
where U is diagonal matrix with positive elements sampled from (11, .1), V a random matrix sampled from wishart distribution, each element of A" is sampled from W 1) and ET is the noise matrix added to the components whose each element is sampled from AV(0. 1) and E added to ensure that the final matrix is positive definite. However the diagonal elements of @ are not equal to To make them 1, we extract diagonal elements of @" and get the new correlation matrix as D1/2OrD1/2 We used 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 8, and тi are selected from [0.1, 1], [1, 5), [0.1, 0.5, 1] and 10 22:2, The criterion for choosing the best hyperparameter is maximum split- sample reproducibility. The split sample reproducibility 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 shows the reproducibility of the components generated from different methods. computed over 15 runs in all the experiments. We used accuracy of the estimated components
5

---

Table 1: Accuracy of the components on simulated dataset one level.
Method
h1
10
h 12
k1 14
hSCP
0.789
0.787
0.745
0.736
ComBat hSCP
0.763
0.759
0.731
0.718
Adv. hSCP
0.869
0.875
0.862
0.854
rshSCP
0.873
0.865
0.843
0.867
Adv. rshCP
0.903
0.910
0.902
0.908
rshSCP w/ rand.
0.856 + 0.039
0.834 + 0.055
0.824 + 0.031
0.818 + 0.036
Adv. rshSCP. w/ rand.
0.897 + 0.030
0.895 0.036
0.892 + 0.023
0.886 + 0.034
Table 2: Split sample reproducbility simulated dataset one level.
Method
1 =
h = 10
h 12
h1 = 14
hSCP
0.769 + 0.052
0.798 + 0.047
0.739 + 0.053
0.734 + 0.047
ComBat hSCP
0.749 = 0.040
0.750 = 0.052
0.724 0.049
0.719 = 0.052
Adv. hSCP
0.781 + 0.037
0.818
0.780 + 0.034
0.750 0.033
rshSCP
0.825 + 0.039
0.845 0.030
0.826 0.039
0.779 + 0.036
Adv. rshSCP
0.840 + 0.044
0.869 + 0.034
0.815 + 0.035
0.802 + 0.030
rshSCP rand.
0.804 0.085
0.818 0.086
0.780 0.068
0.758 + 0.071
Adv. rshSCP w/ rand.
0.826 = 0.069
0.833 + 0.081
0.801 0.074
0.782 = 0.077
Table 3: Accuracy of the components on hierarchical simulated dataset.
hg
4
hz 6
Method &
8
10
12
14
8
10
12
14
hSCP
0.806
0.801
0.783
0.777
0.797
0.790
0.773
0.766
ComBat hSCP
0.788
0.776
0.743
0.729
0.779
0.766
0.747
0.734
Adv. hSCP
0.875
0.872
0.870
0.864
0.863
0.859
0.849
0.851
rshSCP
0.881
0.876
0.860
0.862
0.874
0.871
0.852
0.859
Adv. rshSCP
0.904
0.909
0.904
0.907
0.902
0.903
0.904
0.902
performance measure. It is defined 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 displays the accuracy of different methods on the simulated dataset. Here, accuracy is defined 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 A 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 A 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, 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 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.
Site prediction. To check if estimated subject information (A) has reduced predictive power to predict the site to which the subject belonged, we performed 5 fold cross- valid-valiotioon 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

---

Table 4: Summary characteristics of the real dataset.
Data Sites
Participants
% of Females
Age Range (Median)
Scanner
BLSA-3T
784
56.5
22, 95 (8)
3T Philips
CARDIAI
199
55.7
42, (52)
3T Siemens Tim Trio
CARDIA2
321
51.4
[43, 61
3T Philips Achieva
CARDIA3
278
55.3
[43, (52)
3T Philips Achieva
UKBB
2023
55.2
[45, (9)
3T Siemens Skyra
OASIS
847
56.0
[42.6, (70
1.5T Siemens Vision
ABC
279
59.1
[23, (5,95((7))
3T Siemens Tim Trio
ReLU, fully-connected layer with hidden units and softmax layer and (b) fully connected layer with hidden units, dropout layer with the rate ReLU, a fully-connected layer with hidden units, dropout layer with rate ReLU, fully-connected layer with hidden units and softmax layer. Our model leads to a decrease in average cross- cross-validation accuracy from 97. 77.%% to 67% for SVM, 98.1% to 67.3% for neural network with architecture (a) and 98.2% 2% 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.I
Two level. Under the same settings as defined above, we generate correlation matrix from two level components with 22 using:
-t W, W.WA W W + UVV + E
(12)
W, WTE,. W₂ W
where each element of W2 is sampled from (0,,)), 1), the percentage of non-zeros equal to 40%, E and E& is the noise added to the components whose each element sampled from (O, .1) and Eg is added to ensure that the final matrix is positive definite. Table shows the accuracy for different values of k and 122. 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.I
3.2 Real Dataset
Data. We collected functional MRT data from 5 different multi-center imaging studies- 1) Baltimore
Longitudinal Study of Aging (BLSA) 32 33 the Coronary Artery Risk Development 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 's data was 100. Table summarizes the number of participants and the preprocessing pipeline are given in Appendix C CARDIA data is divided into three parts because of the acquisition at three different sites.
each site and age distribution. More details on the dataset
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 and all sites except Split sample reproducibility tries to answer the question of whether the components are generalizable other sites or not. For estimating rshSCP with only one site, we used V estimated from all sites except 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 Ti and 72 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. 5 showss the split sample reproducibility for varied values of k and l = Leave one site out reprodueibility results are shown in Table Table [3 and Table [4 in D.2 shows split sample reproducibility and leave one site out reproducibility respectively at two-level for h2 = 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
7

---

Table reipuubbyyl on real dataset (h: =).
Method
h = 10
k 15
h = 20
h1 25
hSCP
0.713 + 0.039
0.707 + 0.038
0.697 = 0.035
0.683 + 0.036
ComBat hSCP
0.673 + 0.049
0.641 + 0.051
0.639 = 0.031
0.611 + 0.038
Adv. hSCP
0.737 = 0.041
0.719 0.033
0.715 = 0.037
0.710 #
rshSCP
0.806 + 0.036
0.768 = 0.032
0.742 = 0.033
0.743 0.044
Adv. rshSCP
0.808 + 0.030
0.772 0.036
0.747 + 0.034
0.746 0.036
Table 6: Leave one site out reproducbility on real dataset(k2 4).
Method
h = 10
k1 15
h1 20
k1 = 25
hSCP
0.652 + 0.038
0.618 + 0.041
0.592 = 0.033
0.571 + 0.035
ComBat hSCP
0.614 + 0.042
0.594 + 0.035
0.542 = 0.041
0.528 + 0.039
Adv. hSCP
0.656 = 0.035
0.629 = 0.039
0.601 = 0.035
0.584 + 0.034
rshSCP
0.712 + 0.034
0.701 + 0.036
0.676 = 0.038
0.665 + 0.034
Adv. rshSCP
0.716 + 0.032
0.709 0.031
0.688 0.034
0.671 + 0.033
Table 7: Mean absolute error (kg 4)
Method
k = 10
k = 15
hn = 20
kj = 25
hSCP
6.490 + 1.485
6.468 + 1.442 6.467
6.425 + 1.412
6.414 + 1.417
Adv. rshSCP
6.494 + 1.501
1.475
6.432 + 1.483
6.409 1.470
the comparison between components learned using adversarial learning from hSCP and Adv. rshSCP.
Site prediction. We performed the same experiment under the same settings mentioned in the previous section to check A has reduced predictive power to predict the site. Using SVM, our model leads to decrease in average cross- vols-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 voss-vvaldation 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 (A) predict age of each subject as the metric to check if our method is able to preserve age related biological variation. Table shows 10 fold cross validation mean absolute error (MAE) using random forest. From the table, 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
see
3.3 Analysis of components
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 (A) 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 I shows the components derived from hSCP and the proposed method. The first row of the figure displays the components with
8

---

(a) Component 5: p - 0.02. In(p) 1.57
(b) Component 5 -0.13, In(p) 28.6
(c) Component 6: p 0.07, hn(p) 9.32
(d) Component 6
0005 - In(p) 5.27
Left column (aa) & (c)) displays the components estimated using hSCP and right column ((b) & (d)) displays the components estimated using rshSCP. Red and blue regions anti-correlated with each other but are correlated among themselves. The colors are not associated with negative or positive correlation since they can be swapped without affecting the final inference.
Figure I:
2
Figure 2: One of the hierarchical components derived from rshSCP comprising of component 2 and 7 fine scale and component II at coarse scale.
anti-correlation between Default Mode Network (DMN) and Dorsal Attention Network (DAN). The component derived using hSCP has 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 IMRI and task-based IMRI (38]. The magnitude of anti-correlation has been connected to individual differences in task performances healthy young adults [39]. However, in the case of older adults, the behavioral implications of reduced amti-correlation remain unclear. The second row of the Figure I 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 positive correlation of DMN with visual areas will cause misleading inference since it contradicts the previous SCPs and studies. Hence making inference without removing the site effects can misleading. Discussion of the association between other relevant components and aging is in Appendix [.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 reorganization of the aging brain aligning with the previous findings This can serve as a base to explore rshSCP as a biomarker of neurodegenerative discases.
F 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 reed-based methods that are dependent on the knowledge of the seed region of interest.

---

4 Conclusion
In this work, we have presented method for estimating site effects in hSCP. We formulated the problem as 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 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 better split sample and leave one site out reproducibility without losing biological interpretability and information. We also show that without removal of site effects, 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 whether the proposed method preserves components associated with other demographie, 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 non- linear 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.
Broader Impact
In this work, we provide a new method 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. can be casily extended to various matrix factorization approaches such as Independent Component Analysis, Non- 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 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
References
[1] Stephanie Noble, Dustin Scheinost, Emily S Finn, Xilin Shen, Xenophon Papademetris, Sarah C McEwen, Carrie E Bearden, Jean Addington, Bradley Goodyear, Kristin S Cadenhead, et al. Multisite reliability of mr mr-based functional connectivity. Neuroimage, 146:959-99 970, 2017.
[2] Adriana Di Martino, Chao-Gan Yan, Qingyang Li, Erin Denio, Francisco X Castellanos, Kaat Alaerts. Jeffrey S Anderson, Michal Assaf, Susan Y Bookheimer, Mirella Dapretto, al. The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture autism. Molecular psychiatry, 19(6);599 667, 2014.
[3] Christian Dansereau, Yassine Benhajali, Celine Risterucci, Emilio Merlo Pich, Pierre Orban, Douglas Amold, and Pierre Bellec. Statistical power and prediction accuracy in multisite resting-state fmri connectivity. Neuroimage, 449:220-232, 2017.
[4] Anisha Keshavan, Friedemann Paul, Mona K Beyer, Alyssa H Zhu, Nico Papinutto. Russell T Shinohara, William Stern, Michael Amann, Rohit Bakshi, Antje Bischof, et al. Power estimation for standaaaddd multisite studies. Neurolmage, 134:281-20 2016.
10

---

[5] Russell T Shinohara, Jiwon Oh, Govind Nair, Peter A Calabresi, Christos Davatzikos, Jimit Doshi, Roland G Henny, Gloria Kim, Kristin A Linn, Nico Papinutto, al. Volumetric analysis from harmonized multisite brain mri study of a single subject with multiple sclerosis. American Journal of Neuroradiology. 38(8):1501-1500 2017.
[6] Alexandre Abraham, Michael P Milham, Adriana Di Martino, R Cameron Craddock, Dimitris Samaras, Bertrand Thirion, and Gael Varoquaux. Deriving reproducible biomarkers from multi-site resting-state data: An autism-based example. Neurolmage, 147:736-745, 2017
[7] Jorge Jovicich, Ludovico Minati, Moira Marizzoni, Rocco Marchitelli, Roser Sala-Llonch, David Bartrés- Faz, Jennifer Amold, Jens Benninghoff, Ute Fiedler, Luca Roccatagliata, al. Longitudinal reproducibility of default-mode network connectivity in healthy elderly participants: multicentric resting-state fiuri study. Neuroimage, 124:442-4545 4544, 2016.
[8] Meichen Yu, Kristin Linn, Philip A Cook, Mary L Phillips, Melvin McInnis, Maurizio Fava, Madhukar Trivedi, Myma M Weissman, Russell T Shinohara, and Yvette I Sheline. Statistical harmonization corrects site effects in functional connectivity measurements from multi-site site fmri data. Human brain mapping, 39 (11):4213 4227,222722222 2018.
[9] Marcus E Raichle. The brain's default mode network. Annual review of neuroscience, 38:433 447, 2015.
[10] Alex Foriio,, Andrew Zalesky, and Edward Bullmore. Fundamentals of brain network analyxis. Academic Press, 2016.
[11] Cornelis J Stam. Modern network science of neurological disorders. Nature Reviews Neuroscience, 15(10): 683-695, 2014.
[12] Ian A Cook, Susan Y Bonkheimer, Laura Mickes, Andrew I. Leuchter, and Anand Kumar. Aging and brain activation with working memory tasks: an fmri study of connectivity. International Journal of Geriatric Psychiatry: A journal of the psychiatry of late life and allied sciences. 22(():3322 342, 2007.
[13] Stephen M Smith, Peter T Fox, Karla L Miller, David C Glahn, P Mickle Fox, Clare E Mackay, Nicola Filippini, Kate Watkins, Roberto Toro, Angela R Laird, Correspondence of brain's functional architecture during activation and rest. Proceedings of the national academy of sciences, (06(31):13040 13045, 2009.
[14] Vamsi K Potluru and Vince D Calhoun. Group learning using contrast nmf: Application to functional and structural mri of schizophrenia. In 2008 IEEE International Symposium Circuits and Systems, pages 1336-1339, IEFF, 2008
[15] Kangjoo Lcc, Sungho Tak, and Jong Chul Ye. A data-driven sparse glm for fini analysis using sparse dictionary learning with mdl criterion. IEEE Transactions Medical Imaging, 107 2010.
[16] Dushyant Sahoo, Nicolas Honnorat, and Christos Davatzikos. Sparse low-dimensional causal modeling for the analysis of brain function. In Medical Imaging 2019: Image Processing, volume 10949, page 109492R. International Society for Optics and Photonics. 2019.
[17] Dushyant Sahoo. Nicolas Honnorat, and Christos Davatzikos. Gpu accelerated extraction of sparse granger causality patterns. In 2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018), pages 604-607. IEEE, 2018.
[18] Gaëlle Doucct, Mikaël Naveau, Laurent Petit, Nicolas Delcroix, Laure Zago, Fabrice Crivello, Gael Jobard, Nathalic Tzourio-Man Mazoyer, Bernard Mazoyer, Emmanuel Mellet, et al. Brain activity at rest: a multiscale hierarchical functional organization. Journal of neurophysiology. 2753))2753--263, 2011.
[19] Murat Demirtas, Joshua B Burt, Markus Helmer, Jie Lisa Ji, Brendan D Adkinson, Matthew F Glasser, David C Van Essen, Stamatios N Sotiropoulos, Alan Anticevic, and John D Murray. Hierarchical heterogeneity across human cortex shapes large-scale neural dynamics. Neuron, 101(6):188-11999) 2019.
[20] Wei Zhang, Shijie Zhao, Xintao Hu, Qinglin Dong, Heng Huang, Shu Zhang, Yu Zhao, Haixing Dai, Fangfei Ge, Lei Guo, al. Hierarchical organization of functional brain networks revealed by hybrid spatiotemporal deep leamingg Brain connectiviry, 10(2):2 82, 2020.
[21] Heng Huang, Xintao Hu, Yu Zhao, Milad Makkie, Qinglin Dong, Shijie Zhao, Lei Guo, and Tianming Liu. Modeling task finri data via deep convolutional autoencoder. IEEE transactions on medical imaging, 37 17)11511-5561, 2017.

---

[22] Michael Andersen, Ole Winther, Lars Kai Hansen, Russell Poldrack, and Oluwasanmi Koyejo. Bayesian structure learning for dynamic brain connectivity. In International Conference on Artificial Intelligence and Statisticss pages 1436-1446. PMLR, 2018.
[23] Harini Eavani. Theodore D Satterthwaite, Roman Filipovych, Raquel E Gur. Ruben C Gur, and Christos Davatzikos. Identifying sparse connectivity patterns in the brain using resting-state fmri. Neuroimage. 105: 286-299, 2015.
[24] Dushyant Sahoo. Theodore D Satterthwaite, and Christos Davatzikos. Hierarchical extraction of functional connectivity components in human brain using resting-state fmri. IEEE Transactions Medical Imaging, 2020.
[25] Dushyant Sahoo and Christos Davatzikos. Extraction of hierarchical functional connectivity components in human brain using adversarial learning. arXiv preprint XXXv:21044102553 2021.
[26] W Evan Johnson, Cheng Li, and Ariel Rabinovic. Adjusting batch effects in microarray expression data using empirical bayes methods, Biostatistics, 8(1):118 1-7, 2007.
[27] Raymond Pomponio, Guray Erus, Mohamad Habes, Jimit Doshi, Dhivya Srinivasan, Elizabeth Mamourian, Vishnu Bashyam, Ilya M Nasrallah, Theodore D Satterthwaite. Yong Fan, et al. Harmonization of large mri datasets for the analysis of brain imaging patterns throughout the lifespan. Neurolmage, 208: 16450, 2020.
[28] Dean-Philipppe Fortin, Drew Parker, Birkan Tunç, Takanori Watanabe, Mark A Elliott, Kosha Ruparel. David R Roalf. Theodore D Satterthwaite, Ruben C Gur, Raquel E Gur, et al. Harmonization of multi-site diffusion tensor imaging data. Neuroinage, 161:149-170, 2017.
[29] Xiangnan He, Zhankui He, Xiaoyu Du and Tat-Seng Chua. Adversarial personalized ranking for recom- mendation. In The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval, pages 355-364, 2018
[30] Lei Luo, Yanfu Zhang, and Heng Huang. Adversarial nonnegative matrix factorization. In International Conference on Machine Learning, pages 6479 6488. PMLR, 2020.
[31] Gilbert W Stewart. Perturbation theory for the singular value decomposition. Technical report, 1998.
[32] Nicole M Armstrong, Yang An, Lori Beason-Held, Jimit Doshi, Guray Erus, Luigi Ferrucci, Christos Davatzikos, and Susan M Resnick. Predictors of neurodegeneration differ between cognitively normal and subsequently impaired older adults. Neurobiology of aging, 75:178-186, 2019.
[33] Susan M Resnick, Dzung I. Pham, Michael A Kraut, Alan B Zonderman, and Christos Davatzikos. Longi- tudinal magnetic resonance imaging studies of older adults: a shrinking brain. Journal of Neuroscience, 23 (8):3295-3301, 2003.
[34] Gary D Friedman, Gary R Cutter, Richard P Donahue, Glenn H Hughes, Stephen B Hulley, David R Jacobs Jr, Kiang Liu, and Peter Savage. Cardia: study design, recruitment, and some characteristics of the examined subjects. Journal of clinical epidemiology. 41(11)11005-11 1116, 1988.
[35] Cathie Sudlow, John Gallacher, Naomi Allen, Valerie Beral, Paul Burton, John Danesh, Paul Downey, Paul Elliott, Jane Green, Martin Landray, et al. Uk biobank: an open access resource for identifying the causes of wide range of complex diseases of middle and old age. Plos med, 12(3)xe100177) 2015.
[36] Daniel S Marcus, Tracy H Wang, Jamie Parker, John G Csernansky, John C Morris, and Randy L Buckner. Open access series of imaging studies (oasis): cross- sectional mri data in young, middle aged, nondemented, and demented older adults. Journal of cognitive neuroscience, 19(9):1498 1498 1507, 2007.
[37] John Pluta, Paul Yushkevich, Sandhitsu Das, and David Wolk. In vivo analysis of hippocampal subfield atrophy in mild cognitive impairment via semi-autonatic segmentation of 12-weighted mri. Journal of Alzheimer's 's disease, 31()):85-4 99, 2012.
[38] R Nathan Spreng, W Dale Stevens, Joseph D Viviano, and Daniel L Schacter. Attenuated anticorrelation between the default and dorsal attention networks with aging: evidence from task and rest. Neurobiology of aging, 45:149-166 160, 2016.
[39] Joseph B Keller, Trey Hedden, Todd W Thompson. Sheeba A Anteraper, John DE Gabrieli, and Susan Chittied--Gabiell Resting-state anticorrelations between medial and lateral prefrontal cortex: association with working memory, aging, and individual differences. Cortex, 64:271-280, 2015.
[40] Fikret Işik Karahanoglu and Dimitri Van De Ville. Transient brain activity disentangles finri resting-state dynamics in terms of spatially and temporally overlapping networks. Nature communications, 6(1):1-10, 2015.
12

---

[41] Jessica S Damoiseaux. Effects of aging functional and structural brain connectivity. Neuroinage, 160: 32-40, 2017.
[42] Emanuele Olivetti, Susanne Greiner, and Paolo Avesani. Adhd diagnosis from multiple data sources with batch effects. Frontiers in systems neuroscience, 6:70, 2012.
[43] Bhim M Adhikari, Neda Jahanshad, Dinesh Shukla, Jessica Turner, Dominik Grotegerd, Udo Dannlowski, Harald Kugel, Jennifer Engelen, Bruno Dietsche, Axel Krug, et al. resting state fmri analysis pipeline for pooling inference across diverse cohorts: an enigma rs- r--fmr protocol. Brain imaging and behavior, 13 (5):1453- 1467, 2019.
[44] Seyed Mostafa Kia. Hester Huijsdens, Richard Dinga, Thomas Wolfers, Maarten Mennes, Ole A An- dreassen, Lars T Westlye, Christian F Beckmann. and Andre F Marquand. Hierarchical bayesian regression for multi-site normative modcling of neuroimaging data. In International Conference Medical Image Computing and Computer-Assisted Intervention, pages 699-709. Springer, 2020.
[45] Johanna MM Bayer, Richard Dinga, Seyed Mostafa Kia, Akhil R Kottaram, Thomas Wolfers, Jinglei Lv, Andrew Zalesky, Lianne Schmaal, and Andre Marquand. Accommodating site variation neuroimaging data using hierarchical and bayesian models. bioRxiv, 2021.
[46] Roberto Vega and Russ Greiner. Finding effective ways to (machine) Iearn finri-based classifiers from multi-site data. In Understanding and Interpreting Machine Learning in Medical Image Computing Applications, pages 32-39. Springer, 2018.
[47] Behnam Gholami, Pritish Sahu, Ognjen Rudovic, Konstantinos Bousmalis, and Vladimir Pavlovic. Unsu- pervised multi-target domain adaptation: An information theoretic approach. IEEE Transactions on Image Processing, 29:3993-4002. 2020.
[48] S Zhao, B Li, X Yue, Gu, P Xu, R Hu, H. Chai, and K Keutzer. Multi domain adaptation semantic segmentation", advances neural information processing systems. Advances neural information processing systems, 2019.
[49] Han Zhao, Shanghang Zhang, Guanhang Wu, José MF Moura, Joao Costeira, and Geoffrey Gordon. Adversarial multiple source domain adaptation. Advances in neural information processing systems, 31: 3559-8570, 2018.
[50] Yitong L., Michael Murias, Samantha Major, Geraldine Dawson, and David E. Carlson. Extracting relationships by multi-domain matching. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 6799-6810, 2018.
[51] Xingchao Peng, Qinxun Bai, Xide Xia, Zijun Huang, Kate Sacnko, and Bo Wang. Moment matching for multi-source domain adaptation. Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 1406-1415, 2019.
[52] Sicheng Zhao, Bo Li, Pengfei Xu, and Kurt Keutzer. Multi-source domain adaptation in the deep learning era: A systematic survey. arXiv preprint arXiv:2002.12166 2020.
[53] Mingliang Wang, Daoqiang Zhang, Jiashuang Huang, Pew-Thian Yap, Dinggang Shen, and Mingxia Liu. Identifying autism spectrum disorder with multi-site fmri via low-rank domain adaptation. IEEE transactions medical imaging. 39(3):644 655, 2019.
[54] Junyi Zhang, Peng Wan, and Daoqiang Zhang. Transport- Trassp based joint distribution alignment for multi-site autism spectrum disorder diagnosis using resting-state fmri. In International Conference Medical Image Computing and Computer-Assisted Intervention, pages 444-453. Springer, 2020.
[55] Xiaoxiao Li, Yufeng Gu, Nicha Dvonek, Lawrence H Staib, Pamela Ventola, and James S Duncan. Multi-site fmni analysis using privacy-preserving federated learning and domain adaptation: Abide results. Medical Image Analysis, 65:101765, 2020.
[56] Haldun Ozaktas. Paradigms of connectivity for computer circuits and networks. Optical Engineering. (17))1563-1561) 1992.
[57] Ian F. Akyildiz, Xudong Wang, and Weilin Wang. Wireless mesh networks: a survey, Computer Networks, 4744)445-487 2005. ISSN 1389-1286. doi: https://dbi..rrg/0.000........ 12.001. http://ttttttteeeeeecceeeeeeeeie httpa/awwwcccccaeeo.oaaaaaaaaaaa
[58] Maximillian Nickel and Douwe Kiela. Poincaré embeddings for learning hierarchical representations.
Advances in Neural Information Processing Systems, 30:6338-633 6348,66372222 2017.
13

---

[59] Osame Kinouchi and Mauro Copelli. Optimal dynamical range of excitable networks at criticality, Nature physics, 2(5):348-351, 2006.
[60] Danielle S Bassett, Daniel L Greenfield, Andreas Meyer-Lindenbers Daniel R Weinberger. Simon W Moore, and Edward T Bullmore. Efficient physical embedding of topologically complex information processing networks in brains and computer circuits. PLoS comput biol, 104))e 2010.
[61] Rishidev Chaudhuri, Alberto Bernacchia, and Xiao-Jing Wang. A diversity of localized timescales in network activity. elife, 3:c01239, 2014.
[62] Luca Ferrarini, Nlya M Veer, Evelinda Bacrends, Marie-José van Tol, Remco J Renken, Nic JA van der Wee, Dirk Veltman, André Aleman, Frans G Zitman, Brenda WJH Pemminx, al. Hierarchical functional modularity in the resting-state human brain. Human bruin mapping, 30(7):2220-223) 2009.
[63] Esraa Al-Sharoa, Mahmood Al-Khassaweneh, and Selin Aviyente. Tensor based temporal and multilayer community detection for studying brain dynamics during resting state fmri. IEEE Transactions on Biomedical Engineering, 66(3):695-700) 2018.
[64] Teddy Akiki and Chadi Abdallah. Determining hierarchical architecture of human brain using subject-level clustering of functional networks. Scientific reports, 9(1):1 15, 2019.
[65] Arian Ashourvan, Qawi Telesford, Timothy Verstynen, Jean M Vettel, and Danielle Bassett. Multi-scalc detection of hicrarchical community architecture in structural and functional brain networks, PLoS One, 14(5))e0215520, 2019.
[66] Annette L Fitzpatrick, Catherine K Buchanan, Richard Nahin, Steven DeKosky, Hal H Atkinson. Michelle C Carlson, and Jeff D Williamson. Associations of gait speed and other measures of physical function with cognition in healthy cohort of elderly persons. The Journals of Gerontology Series A: Biological Sciences and Medical Sciences, 2211))1224-1255 2007.
[67] Liang Zhan, Lisanne M Jenkins, Ouri E Wolfson, Johnson Jonaris GadElkarim, Kevin Nocito, Paul M Thompson, Olusola A Ajilore, Moe K Chung, and Alex D Leow. The significance of negative correlations in brain connectivity. Journal of Comparative Neurology, 525(15):3251
2017.
[68] Michael D Fox, Abraham Z Snyder, Justin I. Vincent, Maurizio Corbetta, David C Van Essen, and Marcus E Raichle. The human brain is intrinsically organized into dynamic, anticorrelated functional networks. Proceedings of the National Academy of Sciences, 102(27):973--967,, 2005.
[69] Xintao Hu, Heng Huang, Bo Peng, Junwei Han, Nian Liu, Jinglei Lv. Lei Guo, Christine Guo, and Tiamming Liu. Lalent source mining firi via restricted boltzmann machine. Human brain mapping. 39 (6):2368-233) 2018.
[70] Qinglin Dong, Fangfei Ge. Qiang Ning, Yu Zhao, Jinglei Lv. Heng Huang, Jing Yuan, Jiang, Dinggang Shen, and Tianming Liu. Modeling hicrarchical brain networks via volumetric sparse deep belicf network. IEEE Transactions on Biomedical Engineering, 67(6):1739-174 2019.
[71] Anastasia Podosinnikova, Matthias Hein, and Rainer Gemulla. Robust Principal Component Analysis as a Nonlinear Figenproblem. PhD thesis, Saarland University, 2013.
[72] Sashank Reddi, Satyen Kale, and Sanjiv Kumnrr On the convergence of adam and beyond. arXiv preprint arXiv:1904.0 0000000..99,.222 2019.
[73] Mark Jenkinson, Christian Beckmann, Timothy EJ Behrens, Mark Woolrich, and Stephen M Smith. Fsl. Neuroimage, 0290).782-790, 2012.
[74] Gholamreza Salimi- Khorshidi, Gwenaëlle Douaud, Christian F Beckmann, Matthew F Glasser, Ludovica Griffanti, and Stephen M Smith. Automatic denoising of functional mri data: combining independent component analysis and hierarchical fusion of classifiers. Neuroimage, 90:449-468, 2014.
[75] Ludovica Griffanti, Gholamreza Khinii-Koosshii,, Christian Beckmann, Edward Auerbach, Gwenaölle Douaud, Claire E Sexton, Enikö Zsoldos, Klaus P Ebmeier, Nicola Filippini, Clare E Mackay, et al. Ica-based artefact removal and accelerated fmri acquisition for improved resting state network imaging. Neuroimage, 95:232-247, 2014.
[76] Stephen M Smith, Aapo Hyvärinen, Gael Varoquaux. Karla L Miller, and Christian F Beckmann. Group-pea for very large fmri datasets. Neuroimage, 101:738-749, 2014.
[77] Yuhui Du and Yong Fan. Group information guided ica for fmri data analysis. Neuroimage, 09:157-197 2013.
14

---

[78] Mark Jenkinson. Measuring transformation error by rms deviation. Studholme, C., Hill, DLG, Hawkes, DJ, 1999.
[79] Christina Triantafyllou, Richard D Hoge, Gunnar Krueger, Christopher Wiggins, Andreas Potthast, Graham C Wiggins, and Lawrence L Wald. Comparison of physiological noise at 1.5 1, 3 and 7 and optimization fmri acquisition parameters. Neuroimage. 26(1):243 250- 2005.
[80] Jonathan D Power, Kelly A Barnes, Abraham Snyder, Bradley L Schlaggar, and Steven Petersen. Spurious but systematic correlations functional connectivity mri networks arise from subject motion. Neuroimage, 59(3):1142 2154, 2012.
[81] Mark Jenkinson. Peter Bannister. Michael Brady, and Stephen Smith. Improved optimizattonn for the robust and accurate linear registration and motion correction of brain images. Neuroimage, 17(2):225-3) 841, 2002.
[82] Chao-Gan Yan, Cameron Craddock, Yong He, and Michael Milham. Addressing head motion dependencies for small-world topologies functional connectomics. Frontiers in human neuroscience, 7: 910, 2013.
[83] Irene E Nagel, Claudia Preuschhof, Shu-Chen Li, Lars Nyberg, Lars Bäckman, Ulman Lindenberger, and Hauke R Heekeren. Load modulation of bold response and connectivity predicts working memory performance in younger and older adults. Journal Cognitive Neuroscience, 23(8):2030-25) 2045, 2011.
[84] Linda Geerligs. Natasha M Maurits, Remco J Renken, and Monicque M Lorist. Reduced specificity of functional connectivity in the aging brain during task performance. Human brain mapping, 35(1):319-3300 2014.
[85] Luiz Kobuti Ferreira, Ana Carolina Brocanello Regina, Natasa Kovacevic, Maria da Graça Morais Martin, Pedro Paim Santos, Camila de Godoi Carneiro, Daniel Shikanai Kerr, Edson Amaro Jr. Authony Randal McIntosh, and Geraldo Busatto, Aging effects on whole- functional connectivity in adults free of cognitive and psychiatric disorders. Cerebral cortex, 2619),385, 2016.
[86] Linda Geerligs, Remco Renken, Emi Saliasi, Natasha M Maurits, and Monicque M Lorist. A brain-wide study of age-related changes in functional connectivity. Cerebral cortex, 1977) 2015.
[87] K Warner Schaic and Sherry Willis. Handbook of the Psychology of Aging. Academic Press, USA, 2021.
[88] Cheryl Grady, Saman Sarraf, Cristina Saverino, Karen Campbell. Age differences in the functional interactions among the default, frontoparietal control, and dorsal attention networks. Neurobiology of aging, 41:159-172, 2016.
15


Reviews:
[{"review_id": "boJ9pBPUQfk", "reviewer": "Reviewer_jmAd", "summary": "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.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Originality\nThe proposed methods, rshSCP and Adv. rshSCP, build upon previous methods, hSCP and its adversarial extension (Adv. hSCP), partly by using the usual trick of adding a loss to discourage prediction of confounds, in this case site effects. The idea at a high level is not particularly original, but the implementation details are useful. A bit surprised domain adaptation and optimal transport are not cited and discussed.\n\nQuality\nThe proposed methods are technically sound, but there are a few concerns regarding assumptions and results.\n1. Is retaining the positive semidefinite property essential for downstream analysis? Comparisons with DL methods that do not retain positive semidefinite property are needed.\n2. Following the theme of modeling the manifold property of correlation matrices, why use l2 loss in (1) as opposed to Gaussian likelihood as in standard Gaussian graphical models or intrinsic distance between covariance?\n3. Will there be an identifiability problem since W1*W2 = W1*A* A’*W2 for A orthogonal?\n4. For estimating rshSCP with only one site s in the reproducibility experiments, would using the same V estimated from all sites except s for both reproducibility assessment splits result in inflated reproducibility values?\n5. How much did site predictability drop using adv. rshSCP compared to adv. hSCP?\n6. Age-related biological variability seems to be retained but just on par with hSCP. Does that imply there is no need to account for site effects? What happens if the two studies with quite different age range were removed? Also, how is the performance of adv. hSCP? Overall, to claim biological variability is retained, need to show results for more than just age.\n7. Not sure what to make of the qualitative analysis in Figures 1 and 2. Are there no exceptions where adv rshSCP pools “wrong” regions into a network at any of the hierarchical levels? \n\nClarity\nThe paper is clearly written but there are some debatable statements in the Introduction.\n1. Is hierarchical representation always better than single level parcellation or is it application dependent?\n2. Is “binarized” components always better than their continuous counterparts, e.g. derived from PCA and ICA? The argument was higher interpretability with binarization, yet in an earlier statement, it was mentioned that ICA, NMF, and dictionary learning provide “interpretable” components. Along this line, allowing for components to overlap matches biology but this is provided by PCA and ICA.\n3. Which biological predictors are likely to be correlated with site effects? If we remove site effects, wouldn’t biological effects be removed as well.\n\nSignificance\nThe problem being tackled is important and the approach in the proposed methods is applicable to other methods though implementation details would likely need to be altered to a certain extent. The results, however, are not particularly impressive. Whether the proposed methods are suitable for EEG and MEG is slightly questionable since temporal dynamics is more important for those modalities.\n\nSummary after Author Responses\nWhile some comments were well addressed, a number of issues remained, such as reproducibility using the same V, how much site predictability dropped using adv. rshSCP compared to adv. hSCP not reported for real data, need for prediction beyond age (which is only correlated with site at 0.2), and using limited subjective prior neuroscience findings to evaluate the extracted components (in a subjective manner). Hence, I will keep my original score.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Limitations are briefly mentioned in the conclusions, and for the current prediction task, there is no negative societal impact.", "review_text": "", "extra_scores": {"ethical_concerns": "There are no ethical issues with this paper.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "J03O9c50tBV", "reviewer": "Reviewer_u9U8", "summary": "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.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I found the paper is well written, the proposed idea interesting, and the motivation is clear. I enjoyed that the proposed approach uses matrix factorization to tackle the site-related batch effect problem in fMRI data. However, the base of the formulations for either the matrix factorization or adversarial learning of hSCP is already introduced in other works. The original novelty of the paper includes decomposing the correlation matrix into both hSCP and site components. \n\nOne primary concern regarding this method could be the reproducibility problem and the lack of theoretical guarantees on sub-optimal solutions. Given that this sparse hierarchical pattern is supposed to estimate a proper human brain functionality using fMRI data, it is valuable to design convex objective functions that lead to a unique optimal point rather than multiple sub-optimal solutions. In the current paradigm, the optimal solution for $W$, $D$, $U$, and $V$ is dependent on the initialization procedure as is also mentioned in the paper. Despite the careful initialization of matrices, comparing the performance of the results in either rshSCP and Adv. rshSCP, with and without random initialization, does not show significant improvement (Table 1 and Table 2). \nMoreover, based on the results demonstrated in the experiment section, more justification is required for higher computational cost in the proposed method, given the slight performance improvement comparing to baselines. \n\nA couple of further comments: \n1) Try to define the dimensionality of matrices explicitly as soon as you introduce them. I am not sure what are the dimensionalities of $U^s$ and $V$? How are they defined? What do their dimensionalities represent, given that we want $U^s V$ to fit the dimensionality of $\\theta^n$? \n2) Designing a [group] regularizer that can capture the biological/functional feasibility of subject-specific components or site-related properties might help with non-convex optimization problems and reproducibility issues. \n\nOverall, I feel like the contribution in this paper is not of great enough significance for this venue right now. Yet, I would be open to increasing my overall score based on the authors' responses since the idea behind the paper is interesting.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors adequately addressed the limitations of the work under either future directions or the weaknesses of the model in the conclusion section. They also provided discussions on different challenges along with the paper and tried to alleviate the corresponding limitations by suggesting and implementing new propositions. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "14"}, "extra_text": {}}, {"review_id": "DoI5xwjugO4", "reviewer": "Reviewer_qpQL", "summary": "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.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Overall, the paper is clear and well-written. The primary contribution is the inclusion of an additional additive term to the hSCP decomposition, as well as an adversarial loss and classifier to remove site-specific variation from the learned components. While an adversarial approach to removing confounding effects is not novel, this seems to be its first application to hSCP and the authors carefully demonstrate the ability of this technique to learn more accurate and reproducible components than versions of hSCP that do not use this technique.\n\nOne interesting comparison that could have been made is applying a harmonization method (e.g., ComBat) to ${\\bf \\Lambda}$ directly (as the authors discuss), then addressing the issue of altered relative contributions of the global components by fixing the harmonized ${\\bf \\Lambda}$ values and re-fitting the ${\\bf W}$ matrices. Determining whether this approach achieves comparable reproducibility and/or interpretability of the learned components would provide a simple baseline to compare to the adversarial approach. However, the authors demonstrate the utility of their proposed method rigorously enough even when compared to other versions of hSCP that I don't feel the lack of this comparison to be a major issue.\n\nAnother choice that could be altered is using an SVM to demonstrate the reduced ability to predict site information from the learned subject-specific components. The authors do not justify this other than to say \"we already used feed forward network\". It is unclear whether a feed forward neural network (such as the architecture used to train the adversarial classifier) would be better at extracting site information even from the learned rshSCP components, and showing the results with the same or similar choice of classifier would strengthen that result.\n\nOther minor comments:\n- The subscript for the second ${\\bf W}$ in the top row of Equation 1 should be $j$, not $n$.\n- If my understanding is correct, the cross-entropy loss term $\\mathcal{L}$ in Equation 7 should have a negative sign (so that the cross-entropy is maximized with respect to the hSCP parameters and minimized with respect to $\\zeta$ in the joint optimization)", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes, the authors discuss the potential negative impact and limitations of their work when considering the validity/interpretability of findings in neuroscience.", "review_text": "", "extra_scores": {"ethical_concerns": "No ethical concerns come to mind.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "CPS5n-WX98S", "reviewer": "Reviewer_8hgE", "summary": "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.\n\n\n", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "While this contribution shows that site information can be discounted, I see several issues with it as it stands, which transcend through all evaluation criteria of interest.\n\nTechnical questions and comments:\n- What is the trace(Lambda) = 1 constraint for? Does this scale correctly with n? It seems like the W need to grow bigger and bigger to accommodate trace(correlation matrix) = n. Have you considered other constraints?\n- Why is Gamma_n created by adding 1 to Theta_n? Why is the adversarial perturbation a matrix of 1s and not e.g. a noise vector? Also, 1_P is not defined in the notation, but J_P is. Is 1_P == J_P or does it mean something else? Furthermore, adding anything to Theta without care or renormalization takes it away from the set of feasible correlation matrices.\n- Eq (5) necessarily creates non-symmetric matrices if U^s is not a multiple of the identity matrix. The diagonal should be applied on both sides if anything. How is this not a problem for the algorithm?\n- How are all the lambda_r set? (also: using lower-case and upper-case lambda for different things is a bit confusing. At the beginning I thought that lower-case lambda were the diagonal entries of upper-case lambda)\n- Is the W_n in formula (1) a typo and should be W_j?\n\nExperiment questions and comments\nThe fact that the proposed method works better under the evaluated metrics than the same method without site information is the absolute minimum bar: If it doesn’t do that, you need to look for another method. The more interesting comparisons are ones that compare to existing literature. While there is a paragraph stating that this is somehow impossible, I do not believe that to be the case.\nTwo things should be done:  1) known preprocessing methods (the ones cited) should be used to remove confounds, followed by the extraction of the hierarchical decomposition 2) Other methods for correlation-matrix extraction should be used in combination with deconfounding methods, e.g. group sparse inverse covariance methods.\nIf evaluation ever becomes a problem, here is a suggested common evaluation metric: Generate BOLD time series. For example, given the BOLD time series in some ROIs or maybe a full hemisphere, use the extracted correlation matrix (+ variances if necessary) to generate BOLD time series for the missing part of the brain, for every voxel. This evaluation method also takes into account the parcellation chosen.\n\nCan this method, being hierarchical, support going down to correlations at the voxel-level? This would be an interesting step forward.\n\nIt would be very interesting to see this method elevated to a general confound elimination method. Sometimes subject age is a confound. It seems like this framework can accommodate this confound in the same way as site is a confound. It would be interesting to see this method predicting a pathological variable with e.g. age, sex and site as confounds. One could evaluate its robustness against choices of test data at very different propensity scores (such as evaluated in the recent https://www.biorxiv.org/content/10.1101/2021.06.16.448764v1).\n\n\nConceptual questions and comments\nThere seems to be some semantic fudging and absence of logical connection around the use of the term “hierarchical”. First of all, if at all applicable at the level of neurons, this does not mean that this term is applicable at the level of voxels and vice versa. More importantly, there is no requirement or indication that incorporating a hierarchical matrix factorization idea into a connectivity estimate captures any of the aforementioned structure. One way of proving that there might be a link would be to show that the particular factorization method introduced here is better than flat methods such as TV+L1-extracted regions -> correlation matrix, followed by conditional BOLD generation.\n\nThe title is phrased way too general. This paper is not about fMRI data analysis, but very specifically about connectivity estimation. At minimum this should be rephrased. \n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "There is unlikely to be any societal impact. The possible societal impacts would be positive for medical diagnosis, but fMRI has not established itself as a tool of choice in any field of medical diagnosis.\n\nThe obvious limitations of any approach related to fMRI \"connectivity\" (first and foremost the neglect of first-order driving effects, second the acknowledgement of all possible other confounding effects) should be required to be mentioned in every paper that addresses it. Sadly, this is at odds with the incentives for publishing, so it would need to become a requirement.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}]

Author rebuttal:

Meta-review: Dear authors,

despite some initial concerns reviewers have now, after discussion and reading
your rebuttal, endorsed the paper (although not strongly). I will therefore suggest
a rather positive outcome for this work.

Best regards,
The AC
Decision: Accept (Poster)