{'fraction_non_alphanumeric': 0.03660012099213551, 'fraction_numerical': 0.020685001628740287, 'mean_word_length': 4.77752386073397, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Statistical power and measurement bias in multisite resting-state fMRI connectivity

Connectivity studies using resting-state functional magnetic resonance imaging are increasingly pooling data acquired at multiple sites. While this may allow investigators to speed up recruitment or increase sample size, multisite studies also potentially introduce systematic biases in connectivity measures across sites. In this work, we measure the inter-site bias in connectivity and its impact on our ability to detect individual and group differences. Our study was based on real multisite fMRI datasets collected in N=345 young, healthy subjects across 8 scanning sites with 3T scanners and heterogeneous scanning protocols. We first empirically show that typical functional networks were reliably found at the group level in all sites, and that the amplitude of the inter-site bias was small to moderate, with a Cohen's effect size below 0.5 on average across brain connections. We then implemented a series of Monte-Carlo simulations, based on real data, to evaluate the impact of the multisite bias on detection power in statistical tests comparing two groups (with and without the effect) using a general linear model, as well as on the prediction of group labels with a support-vector machine. As a reference, we also implemented the same simulations with fMRI data collected at a single site using an identical sample size. Simulations revealed that using data from heterogeneous sites only slightly decreased our ability to detect changes compared to a monosite study with the GLM, and had a more serious impact on prediction accuracy. However, the deleterious effect of multisite data pooling tended to decrease as the total sample size increased, to a point where differences between monosite and multisite simulations were small with N=120 subjects. Taken together, our results support the feasibility of multisite studies in rs-fMRI provided sample size is large enough.


Abstract
Connectivity studies using resting-state functional magnetic resonance imaging are increasingly pooling data acquired at multiple sites. While this may allow investigators to speed up recruitment or increase sample size, multisite studies also potentially introduce systematic biases in connectivity measures across sites. In this work, we measure the inter-site bias in connectivity and its impact on our ability to detect individual and group differences. Our study was based on real multisite fMRI datasets collected in N = 345 young, healthy subjects across 8 scanning sites with 3T scanners and heterogeneous scanning protocols, drawn from the 1000 functional connectome project. We first empirically show that typical functional networks were reliably found at the group level in all sites, and that the amplitude of the inter-site bias was small to moderate, with a Cohen's effect size below 0.5 on average across brain connections. We then implemented a series of Monte-Carlo simulations, based on real data, to evaluate the impact of the multisite bias on detection power in statistical tests comparing two groups (with and without the effect) using a general linear model, as well as on the prediction of group labels with a support-vector machine. As a reference, we also implemented the same simulations with fMRI data collected at a single site using an identical sample size. Simulations revealed that using data from heterogeneous sites only slightly decreased our ability to detect changes compared to a monosite study with the GLM, and had a more serious impact on prediction accuracy. However, the deleterious effect of multisite data pooling tended to decrease as the total sample size increased, to a point where differences between monosite and multisite simulations were small with N = 120 subjects. Taken together, our results support the feasibility of multisite studies

Introduction
Main objective. Multisite studies are becoming increasingly common in restingstate functional magnetic resonance imaging (rs-fMRI). In particular, some consortia have retrospectively pooled rs-fMRI data from multiple independent studies comparing clinical cohorts with control groups, e.g. normal controls in the 1000 functional connectome project (FCP) (Biswal et al., 2010), children and adolescents suffering from attention deficit hyperactivity disorder from the ADHD200 Fair et al., 2012), individual diagnosed with autism spectrum disorder in ABIDE (Nielsen et al., 2013), individuals suffering from schizophrenia (Cheng et al., 2015), or elderly subjects suffering from mild cognitive impairment . The rationale behind such initiatives is to dramatically increase the sample size at the cost of decreased sample homogeneity. The systematic variations of connectivity measures derived using different scanners, called site effects, may decrease the statistical power of group comparisons, and somewhat mitigate the benefits of having a large sample size (Brown et al., 2011;Jovicich et al., 2016). In this work, our main objective was to quantitatively assess the impact of site effects on group comparisons in rs-fMRI connectivity.
Group comparison in rs-fMRI connectivity. We focused in this work on the most common measure of individual functional connectivity, which is the Pearson's correlation coefficient between the average rs-fMRI time series of two brain regions. To compare two groups, a general linear model (GLM) is typically used to establish the statistical significance of the difference in average connectivity between the groups. Finally a p-value is generated for each connection to quantify the probability that the difference in average connectivity is significantly different from zero (Worsley and Friston, 1995). If the estimated p-value is smaller than a prescribed tolerable level of false-positive findings, generally adjusted for the number of tests performed across connections, say α = 0.001, then the difference in connectivity is deemed significant.
Statistical power in group comparisons at multiple sites. The statistical power of a group comparison study is the probability of finding a significant difference, when there is indeed a true difference. A careful study design involves to select a sample-size large enough to reach a given level of statistical power, e.g. 80%.
In the GLM, the statistical power actually depends on a series of parameters (Desmond and Glover, 2002): (1) the sample size (the larger the better); (2) the absolute size of the group difference (the larger the better), and, (3) the intrinsic variability of measurements (the smaller the better). In a multisite (or multi-protocol) setting, differences in imaging or study parameters may add variance to rs-fMRI measures, e.g. the scanner make and model , repetition time, flip angle, voxel resolution or acquisition volume , experimental design such as eyes-open/eyes-closed (Yan et al., 2009), experiment duration (Van Dijk et al., 2010, and scanning environment such as sound attenuation measures (Elliott et al., 1999), or head-motion restraint techniques (Edward et al., 2000;Van Dijk et al., 2012), amongst others. These parameters can be harmonized to some extent, but differences are unavoidable in large multisite studies. The recent work of Yan et al. (2013a) has indeed demonstrated the presence of a significant bias in rs-fMRI measures between sites in the 1000 FCP. Site effects will increase the variability of measures, and thus decrease statistical power. To the best of our knowledge, it is not yet known how important this decrease in statistical power may be.
Sources of variance in rs-fMRI. The relative importance of site effects in rs-fMRI connectivity depends on the amplitude of the many other sources of variance. First, rs-fMRI connectivity only has moderate-to-good test-retest reliability using standard 10-minute imaging protocols (Shehzad et al., 2009), even when using a single scanner and imaging session. Differences in functional connectivity across subjects are also known to correlate with a myriad of behavioural and demographic subject characteristics (Anand et al., 2007;Sheline et al., 2010;Kilpatrick et al., 2006). Taken together, these sources of variance reflect a fundamental volatility of human physiological signals. In addition to physiology, some imaging artefacts will vary systematically from session to session, even at a single site. For example, intensity non-uniformities across the brain depend on the positioning of subjects (Caramanos et al., 2010). Room temperature has also been shown to impact MRI measures (Vanhoutte et al., 2006). Given the good consistency of key findings in resting-state connectivity across sites, such as the organization of distributed brain networks (Biswal et al., 2010), it is reasonable to hypothesize that site effects will be small compared to the combination of physiological and within-site imaging variance.
Multivariate analysis. Another important consideration regarding the impact of site effects on group comparison in rs-fMRI connectivity is the type of method used to identify differences. The concept of statistical power is very well established in the GLM framework, which tests one brain connection at a time (mass univariate testing). However, multivariate methods that combine several or all connectivity values in a single prediction are also widely used and likely affected by the site effects. A popular multivariate technique in rs-fMRI is supportvector machine (SVM) (Cortes and Vapnik, 1995). In this approach, the group sample is split into a training set and a test set. The SVM is trained to predict group labels on the training set, and the accuracy of the prediction is evaluated independently on the test set. Because SVM has the ability to combine measures across connections, unlike univariate GLM tests, we hypothesized that the GLM and SVM will be impacted differently by site effects.
Specific objectives. Our first objective was to characterize, using real data, the amplitude of systematic biases in rs-fMRI connectivity measures across sites, as a function of within-site variance. We based our evaluation on images generated from independent groups at 8 sites equipped with 3T scanners, in a harmonized subset (N = 345) of the 1000 FCP. Our second objective was to evaluate the impact of site effects on the detection power of group differences in rs-fMRI connectivity, as a function of the amplitude of the group difference, sample size, as well as the balancing of groups across sites. We implemented for this purpose a series of Monte Carlo simulations, mixing synthetic data with real data in the 1000 FCP sample. One of the particularity of the 1000 FCP is the presence of one large site of ∼ 200 subjects and 7 small sites of ∼ 20 subjects per site. We were therefore able to implement realistic scenarios following either a monosite or a multisite design (with 7 sites), with the same total sample size. Finally, we repeated the Monte-Carlo using a SVM instead of a GLM, and assessed the impact of site effects on prediction accuracy rather than statistical power.

Imaging sample characteristics
The full 1000 FCP sample includes 1082 subjects, with images acquired over 33 sites spread across North America, Europe, Australia and China. As the 1000 FCP is a retrospective study, no effort was made to harmonize population characteristics or imaging acquisition parameters. Many sites thus featured some outliers characteristics within the sample, such as images acquired at 1.5T or 4T field strengths (5 sites), a population composed mainly of older (4 sites) or Asian (6 sites) participants, samples composed almost exclusively of male or females (8 sites), or partial brain coverage in rs-fMRI. To avoid possible biases in rs-fMRI measures related to such outliers characteristics, a subset of sites was selected based on the following harmonization criteria: (1) 3T scanner field strength, (2) full brain coverage for the rs-fMRI scan, and, (3) a minimum of 15 young or middle aged adult participants, with a mixture of males and females (4) samples drawn from a population with a predominant Caucasian ethnicity. In addition, only young and middle aged participants (18-46 years old) were included in the study, and we further excluded subjects with excessive motion (see next Section). The final sample for our study thus included 345 cognitively normal young adults (150 males, age range: 18-46 years, mean±std: 23.8 ±5.14) with images acquired across 8 sites located in Germany, the United Kingdom, Australia and the United States of America. The total time of available rs-fMRI data for these subjects ranged between 6 and 7.5 min and only one run was available. See Table 1 for more details on the demographics and imaging parameters at each site selected in the study. The experimental protocols for all datasets as well as data sharing in the 1000 FCP were approved by the respective ethics committee of each site. This secondary analysis of the 1000 FCP sample was approved by the local ethics committee at CRIUGM, University of Montreal, QC, Canada.

Computational environment
All experiments were performed using the NeuroImaging Analysis Kit, NIAK 1 (Bellec et al., 2011) version 0.12.18, under CentOS version 6.3 with Octave 2 version 3.8.1 and the Minc toolkit 3 version 0.3.18. Analyses were executed in parallel on the Mammouth supercomputer 4 , using the pipeline system for Octave and Matlab, PSOM (Bellec et al., 2012) version 1.0.2. The scripts used for processing can be found on Github 5 . Prediction was performed using the LibSVM library (Chang and Lin, 2011). Visualization was implemented using Python 2.7.9 from the Anaconda 2.2.0 6 distribution, along with Matplotlib 7 (Hunter, 2007), Seaborn 8 and Nilearn 9 for brain map visualizations.

Preprocessing
Each fMRI dataset was corrected for slice timing; a rigid-body motion was then estimated for each time frame, both within and between runs, as well as between one fMRI run and the T1 scan for each subject (Collins et al., 1994). The T1 scan was itself non-linearly co-registered to the Montreal Neurological Institute (MNI) ICBM152 stereotaxic symmetric template , using the CIVET pipeline (Ad- Dab'bagh et al., 2006). The rigid-body, fMRIto-T1 and T1-to-stereotaxic transformations were all combined to re-sample the fMRI in MNI space at a 3 mm isotropic resolution. To minimize artifacts due to excessive motion, all time frames showing a frame displacement, as defined in Power et al. (2012), greater than 0.5 mm were removed and a residual motion estimated after scrubbing. A minimum of 50 unscrubbed volumes per run was required for further analysis (13 subjects were rejected). The following nuisance covariates were regressed out from fMRI time series: slow time drifts (basis of discrete cosines with a 0.01 Hz highpass cut-off), average signals in conservative masks of the white matter and the lateral ventricles as well as the first principal components (accounting for 95% variance) of the six rigid-body motion parameters and their squares (Giove et al., 2009;Lund et al., 2006). The fMRI volumes were finally spatially smoothed with a 6 mm isotropic Gaussian blurring kernel. A more detailed description of the pipeline can be found on the NIAK website 10 and Github 11 .

Inter-site bias in resting-state connectivity
Functional connectomes. We compared the functional connectivity measures derived from different sites of the 1000 FCP. A functional brain parcellation with 100 regions was first generated using a bootstrap analysis of stable clusters (Bellec et al., 2010b), on the Cambridge cohort of the 1000 FCP (N = 195), as described in Orban et al. (2015). For a given pair of regions, the connectivity measure was defined by the Fisher transformation of the Pearson's correlation coefficient between the average temporal rs-fMRI fluctuations of the two regions. For each subject, a 100 × 100 functional connectome matrix was thus generated, featuring the connections for every possible pair of brain regions.
Inter-site bias. The inter-site bias at a particular connection was defined as the absolute difference in average connectivity between two sites. In order to formally test the significance of the inter-site bias, we used a GLM including age, sex and residual motion as covariates (corrected to have a zero mean across subjects), as well as dummy variables coding for the average connectivity at each site. For each site, a "contrast" vector was coded to measure the difference in average connectivity between this site and the grand average of functional connectivity combining all other sites. A p-value was generated for each connection to quantify the probability that the observed effect using this contrast was significantly different from zero (Worsley and Friston, 1995). The number of false discovery was also controlled (q = 0.05) using a BenjaminiHochberg false discovery rate (FDR) procedure (Benjamini and Hochberg, 1995). To quantify the severity of inter-site bias, we derived Cohen's d effect size measure for each connection: |β c |/σ, with β c being the weight associated with the contrast. The standard deviation from the noiseσ was calculated asσ = e 2 /(N − K), e being the residuals from the GLM, N the sample size and K the number of covariates in the model. As secondary analyses, t-tests were also implemented in the GLM to validate that age, sex as well as residual motion made significant contributions to the model.

Simulations
Data generation process. We implemented Monte-Carlo simulations to assess the detection sensitivity of group differences in rs-fMRI connectivity. The simulations were based on the 1000 FCP sample, with 8 sites totaling 345 subjects. The multisite simulations were sampled from 148 subjects, available across S = 7 sites. The monosite simulations were sampled from 195 subjects available at S = 1 site (Cambridge). For each simulation, a subset of subjects of a given size N was selected randomly and stratified by site. For each site, a ratio W of the selected subjects were randomly assigned to a so-called "patient" group. We focus our analysis on connections showing a fair-to-good test-retest reliability based on a previous study reporting 11 connections likely impacted by Alzheimer's disease, see Orban et al. (2015) for details. For each connection, a "pathology" effect was added to the connectivity measures of the subjects belonging to the "patient" group. This additive shift in connectivity for "patients" was selected as to achieve a specified effect size, defined below.
Effect size (Cohen's d). The Cohen's d was used to quantify the effect size. For a group comparison, Cohen's d is defined as the difference µ between the means of the two groups, divided by the standard deviation of the measures within each group, here assumed to be equal. For a given connection between brain regions i and j, let y i,j be the functional connectivity measure for a particular subject of the 1000 FCP sample. If the subject was assigned to the "patient" group in a particular simulation, an effect was added to generate a simulated connectivity measure y * i,j equals to y i,j + µ. For a specified effect size d, the parameter µ was set to d × s i,j , where s i,j is the standard deviation of connectivity between region i and j. The parameter s i,j was estimated as the standard deviation of connectivity measures across subjects in the mono-site sample (Cambridge), without any "pathological" effect simulated.
GLM tests. In order to detect changes between the simulated groups at each pair of connection, a GLM was estimated from the simulated data, using age, sex and residual motion as confounds (corrected to have a zero mean across subjects). To account for site-specific biases, S − 1 dummy variables (binary vectors coding for each site) were added to the model, with S being the total number of sites used in the study, in addition to an intercept accounting for the global average. Finally, one dummy variable coded for the "patient" group. The regression coefficients of the linear model were estimated with ordinary least squares, and a t-test, with associated p-value, was calculated for the coefficient of the "patient" variable. A significant pathology effect was detected if the p value was smaller than a prescribed α level. The α level needs to be adjusted for multiple comparisons (in our case 11 connections, but this would depend on the number of connections selected in a particular study), which can be done in an adaptive manner using FDR. We tested different typical values for α in {0.001, 0.01, 0.05}. For each simulation sample b and each connection, we derived a p-value p ( * b) , and the effect was deemed detected if p * b was lesser than α. The sensitivity of the test for a particular connection was evaluated by the frequency of positive detections over all simulation samples.
Prediction accuracy. In addition to mass univariate GLM tests, we also investigated a linear SVM (Cortes and Vapnik, 1995) using a Monte Carlo simulation framework similar to the one described above. For SVM simulations, all possible connections between the 100 brain regions were used simultaneously to predict the presence of the simulated pathology in a given subject. For a participant assigned to the "patient" group, a "pathology" effect was only simulated in a set percentage of connections, which were randomly selected. The proportion of connections with a non-null effect was denoted as π 1 . For a given simulation at sample size N , the SVM model was trained on N subjects selected randomly and stratified by site. The accuracy of the model was evaluated on a separate simulation implemented with the remaining subjects, unused during training. For example, for a multisite simulation with N = 80 subjects for SVM training, the model accuracy was estimated on 68 subjects: 148 (available subjects) minus 80 (subjects in the training set). During training, a 10-fold cross-validation was used to optimize the hyper-parameters of the SVM independently for each simulation. The mean and standard deviation of accuracy scores across all samples were derived for each simulation scenario.
Simulation experiments. All the simulation parameters have been summarized below: • Sample size N .
• Patient allocation ratio W .
• Number of sites S.
• The type of detection method, either GLM or SVM.
• For GLM tests, the false-positive rate α.
• The effect size d.
For a given set of simulation parameters, we generated B = 10 3 Monte-Carlo samples to estimate either the sensitivity (for GLM test) or the accuracy (for SVM prediction) of the method. For all experiments, we investigated effect sizes d ∈ {0, 2} with a step of 0.01 and α ∈ {0.001, 0.01, 0.05}. The number of site(s) was S = 1 for the monosite analysis and S = 7 for the multisite analysis. We implemented the following experiments: • (E 1 ) Test the impact of the sample size on GLM N ∈ {40, 80, 120}, with a fixed allocation ratio W = 0.5.
• (E 2 ) Test the impact of the allocation ratio on GLM W ∈ {0.5, 0.3, 0.15} for a fixed sample size N = 120.
• (E 3 ) Test the impact of multisite correction and affected connection volume (π 1 ) on the prediction accuracy. For the prediction scenario, we used a range of π 1 ∈ {0.1, 1, 5%}, and two sample sizes N ∈ {80, 120} subjects for training, with model accuracy estimated on N = 68 and N = 28, respectively. The second column shows the significant differences between the average functional connectivity maps of one site versus all the others. Panel C shows the number of sites with a significant inter-site difference for each brain region.

Inter-site bias in fMRI connectivity
Site bias in the default-mode network. We first focused on the connections associated with a seed region located in the posterior cingulate cortex, a key node of the default-mode network (DMN), which is one of the most widely studied resting-state network (Greicius et al., 2004). The connections were based on the Cambridge 100 parcellation, and were represented as a connectivity map, (Figure 1). Figure 1A shows the posterior cingulate cortex connectivity map, averaged across all subjects and all sites. The key regions of the DMN are easily identifiable, and include the posterior cingulate cortex, precuneus, inferior parietal lobule, anterior cingulate cortex, medial pre-frontal cortex (dorsal, anterior and ventral), superior frontal gyri and the medial temporal lobe (Damoiseaux et al., 2006;Dansereau et al., 2014;Yan et al., 2013b). Using a GLM, the average connectivity map of the DMN was then extracted for each site, Figure 1B. Qualitatively, the DMN maps were consistent across sites, as expected based on the literature. We then tested for the significance of the site bias, i.e. the difference in average connectivity at a given site and the average connectivity at all remaining sites. The statistical maps were corrected for multiple comparisons across the brain with FDR at q ≤ 0.05 (Benjamini and Hochberg, 1995). A significant bias for at least one connection could be identified for every site, without exception, Figure 1B. Figure 1C shows how reproducible were the significant biases in connectivity across the brain and sites. The identified significant connections locations were quite variable across sites, most of them being identified at less than three sites.
Site bias across the connectome. In order to extend these observations outside of the DMN, we derived the entire connectome using the Cambridge 100 parcellation. Figure 2A shows the average connectome, pooling all subjects and sites together. The regions have been re-ordered based on a hierarchical clustering with Ward criterion. A network structure is clearly visible as squares of high connectivity on the diagonal of the connectome (as outlined by black lines). Each diagonal square corresponds to the intra-network connectivity for a partition into 7 networks (Figure 2A). These 7 networks were consistent with the major resting-state networks reported using a cluster analysis in previous works (e.g. van den Heuvel et al., 2008;Bellec et al., 2010a;Yeo et al., 2011;Power et al., 2011): the DMN, visual, sensorimotor, dorsal and ventral attentional networks, mesolimbic and cerebellar networks were identified ( Figure 2B). Figure 2C shows how this large-scale connectome organization varied from site to site. The average connectivity per site as well as significant differences with the average of the remaining sites (q ≤ 0.05) is shown in Figure 2C. Visually, consistent with our previous observations in the DMN, the organization of the average connectome into large-scale resting-state networks was preserved across all sites. Some significant site effects were still detected in the connectivity both within each network, as well as between networks. By counting the number of sites showing a significant effect for each pair of region, it was apparent that significant site effects were quite variable in their localization and spread across the full connectome ( Figure 2D). Concerning the association with the other confounding variables in the model (sex, age and motion) many connections were found to be significantly associated with motion, see Supplementary Material Figure S3, although very few connections were found to be significantly associated with the sex and age, see Supplementary Material Figure S4 and S5. This  Site bias vs. within-site variations across subjects. We measured the amplitude of inter-site bias, represented as violin plots across connections using either the absolute difference in average connectivity ( Figure 3A,C) or Cohen's d effect size measures ( Figure 3B,D). The violin plots include either every pairs of connections from the BASC Cambridge parcellation ( Figure 3A,B), or only the 11 connections selected for Monte-Carlo simulations ( Figure 3C,D). For absolute differences, the distributions were mostly consistent across sites, with a median around 0.06, 5% percentile near 0 and 95% percentiles in the 0.08-0.1 range. For Cohen's d, the distributions were also consistent across sites, with a median around 0.33, 5% percentile near 0 and 95% percentiles in the 0.4-0.6 range. , either for a monosite (S = 1, in red) or a multisite (S = 7, in blue) sample, when testing differences between two groups with a GLM and a false-positive rate α = 0.001. The plain curves are the average statistical power across 11 connections, and the shaded area represent ±1 standard deviation across connections. In panel A, the patient allocation ratio is fixed (W = 50%) and three different sample sizes have been tested, N ∈ {40, 80, 120} (Experiment (E 1 )). In panel B, the sample size is fixed (N = 120) and three different patient allocation ratios have been tested W ∈ {15%, 30%, 50%} (Experiment (E 2 )).
This effect size would be deemed small-to-moderate, which suggests that the impact of additive inter-site bias on statistical tests will be limited. Similar findings were observed across all possible connections, or across the 11 pairs of connections selected in the simulation study. The region-to-region maps of within-site, across-subjects standard deviations are presented for the DMN in Supplementary Figure S1, and for the entire connectome in the Supplementary Material S2.

Multisite Monte-Carlo simulations
Statistical power and effect size. Figure 4A shows the relationship between effect size and a GLM detection power in experiment (E 1 ), i.e. for a fixed allocation ratio (W = 50%) and three different sample sizes, N ∈ {40, 80, 120}. The average and std of detection power was plotted across the 11 selected connections. The variations of statistical power across connections were very small for monosite simulations, as the effect size was adjusted based on the standard deviation of each connection within that sample. As expected, the sensitivity increased with sample size, quite markedly. In multisite simulations (S = 7), for a large effect size (d = 1), the detection power was 20% with 40 subjects , 80% with 80 subjects and 95% with 120 subjects. The sensitivity was larger with a single site than multisite sample, yet the difference between the two decreased as sample size increased. With N = 40 and d = 1, the detection power was close to 30% for a single site sample, compared to 20% for the multisite sample. With N = 120 and d = 1, the difference in sensitivity was only of a few percent. The same trend was apparent for all tested effect sizes as well as for α ∈ {0.01, 0.05} (not shown).
Figure 5: Effect size detectable at 80% sensitivity as a function of sample size, for different false-positive rate α ∈ {0.05, 0.01, 0.001} (experiment (E 1 )). All simulations used a balanced patient allocation ratio W = 50%. The monosite performance is shown in red and the multisite in blue. The dotted black line shows the detectable effect size for a classical parametric t-test.
Statistical power and group allocation ratio. Figure 4B shows the relationship between effect size and a GLM detection power in experiment (E 2 ), i.e. for a fixed sample size (N = 120) and three different patient allocation ratio, W ∈ {15%, 30%, 50%}. Overall, we found that the detection power increased with W . For example, with d = 1, the detection power was 65% for W = 15%, and increased to 90% with W = 30%, and finally 95% for W = 50%. The impact of W was observed in both monosite and multisite samples, with an optimal allocation ratio of W = 50% for both. This observation was also made for α ∈ {0.01, 0.05} (not shown).
Detectable effect size, as a function of sample size. An alternative summary of experiment (E 1 ) is to represent the effect size that can be detected with 80% sensitivity, as a function of sample size for monosite and multisite configurations, see Figure 5. As a reference, we computed the same curve for parametric t-test comparisons, under assumptions of normality. As expected, the detectable effect size for parametric t-tests closely followed the monosite estimation. For a small sample size (N = 40), the detectable effect size was notably larger in multisite configurations than in a monosite configuration (difference of about 0.25 in Cohen's d for α = 0.001). However, the difference decreased for large sample sizes to become smaller than 0.1 with N = 120 and α = 0.001. The lowest detectable effect size for a sensitivity of 80% at α = 0.05 was about d = 0.8, achieved in a monosite configuration with N = 120. At this sample size, the difference between single and multisite configurations was marginal, with only a few percent's of difference in detectable effect sizes. Figure 6: Prediction accuracy of patient vs. controls as a function of effect size. Three simulation settings are presented on each plot: monosite (red curve), multisite with regression of site effects (S = 7, blue curve), and multisite without regression of site effects (S = 7, black curve). Accuracy was estimated over B = 10 3 simulation samples with a patient allocation ratio W = 50% and 3 volumes of affected connections π 1 = 0.1% (left column), π 1 = 1% (middle column) and π 1 = 5% (right column). Two sample sizes were tested: N = 120 randomly selected subjects for training, with the remaining N = 28 to estimate accuracy (first row), and N = 80 randomly selected subjects for training, with the remaining N = 68 to estimate accuracy (second row).
Prediction accuracy. In experiment (E 3 ), we examined the impact of effect size and the volume of affected connections on prediction accuracy in a SVM, see Figure 6. The volume of changes π 1 had a major impact on prediction accuracy. At π 1 = 0.1% (around 5 connections) the accuracy level was at chance level across all tested effect sizes, ( Figure 6A). With π 1 = 1%, accuracy slightly increased, but effect sizes larger than d = 2 were still required to reach over 80% accuracy ( Figure 6B). With π 1 = 5%, 95% accuracy was achieved at the same effect size (about d = 1.5) for monosite and multisite simulations, although the accuracy in multisite simulations was notably lower than for monosite simulations across most effect sizes ( Figure 6C). The relationship between effect size and accuracy followed a sigmoidal curve in both settings, yet a sharper, and latter transition between very low and very high accuracy was observed in multisite simulations. Interestingly, correcting for site effects by regressing out the dummy variable before running the SVM classifier had no impact on accuracy levels. The sample size (N = 80 vs N = 120 for training) did have a moderate effect on prediction accuracy: for π 1 = 5% and d = 1 and monosite simulations, accuracy was about 85% with N = 120 ( Figure 6C) and 75% with N = 80 ( Figure 6F).

Discussion and conclusions
Inter-site bias in rs-fMRI connectivity. Typical resting-state networks, such as the DMN, the attentional, visual and sensorimotor networks, were reliably found across sites. This was strongly expected given the relative consistency of their distribution across individuals, studies, preprocessing approaches or even methods used to extract networks (e.g. Damoiseaux et al., 2006;van den Heuvel et al., 2008;Bellec et al., 2010b;Yeo et al., 2011;Power et al., 2011). We however found that significant differences in average connectivity existed between sites (i.e. site biases), as previously reported by Yan et al. (2013b). This connectivity bias may undermine the generalization of the results derived at a single site. The inter-subject (intra-site) standard deviation of the connections was found to be more than twice as large as the inter-site absolute bias, on average across brain connections. This effect size measured in Cohen's d would be deemed small-to-moderate, which suggests that the impact of additive inter-site bias on statistical tests will be limited. This is a reassuring finding supporting the feasibility of statistical tests pooling fMRI data across multiple sites. Previous studies (Sutton et al., 2008;Brown et al., 2011) had reported inter-site variance up to 10 times smaller than inter-subject variability, but these studies had much more homogeneous scanning environments than ours and also used different fMRI outcome measures. In our case, we still investigated only 3T scanners, mostly Siemens, and inter-site bias may be larger when considering other manufacturers or field strengths.
Statistical power and multisite rs-fMRI. After accounting for site-related additive bias in a GLM, the multisite simulation pooling 7 sites together showed detection power close to that of a monosite simulation with equivalent sample size. The difference was noticeable for small sample size (total N = 40), and became very small for a sample size N = 120. Another observation was that, for a given detection power, the lowest effect size that we were able to detect was more variable across connections for a low sample size. Taken together, these observations suggest to use sample size larger than 100 subjects for GLM multisite studies. A limitation of this conclusion is that it is likely dependent on the number of sites pooled in the study and the actual number of subjects in each of those sites.
Statistical power and sample size. For a medium effect size, e.g. d = 0.5, the sensitivity was low (below 20%), even for monosite simulations with N = 120 subjects. This sobering result supports the current trend in the literature to pool multiple data samples to increase sample size, at the cost of decreased homogeneity. We also found that resting-state studies based on 40 subjects or less, even at a single site, are seriously underpowered, except for extremely large effect sizes (Cohen's d greater than 1.5). Finally, unbalanced patient allocation ratio in site samples greatly reduces sensitivity, even in monosite studies. Balanced datasets, i.e. with equal numbers of patients and controls at each site, should therefore be favored.
Prediction. Comparing the monosite and the multisite accuracy curves reveals a substantial drop in accuracy from monosite to multisite across a broad range of effect sizes. However, it should be noted that classifiers trained across multiple data sources will likely generalize better to new observations, which is likely a critical feature in most applications and reflects the true potential clinical utility of this type of technique. Our conclusions are consistent with the work of Nielsen et al. (2013), which compares the prediction of a clinical diagnosis of autism in monosite vs. multisite settings. The authors concluded that the prediction accuracy for the multisite sample was significantly smaller than for the monosite sample. A somewhat surprising observation in our analysis was that linear correction for site-specific biases did not improve accuracy of prediction using SVM. The SVM model seems to learn features that are invariant across sites, maybe focusing on connections with the smallest site bias, or looking at differences between connections similarly impacted by a site bias. Finally, an important conclusion of our simulations was that the volume of brain connections affected by a disease impacts as much accuracy as the effect size per connection. This suggests that feature reduction and/or selection is a very important step to improve sensitivity to small effect sizes.
Beyond additive bias. An important limitation to our study is that we only investigated the impact of additive bias in brain connectivity across sites. Areas of future works include interactions between site effects and pathology, possibly in the form of polynomial and non-linear interactions. We hope that, in the future, fMRI data acquired on clinical cohorts at tens of sites will become available, which will enable researchers to test empirically the presence of such interaction effects.
Other types of multisite data. Another limitation of our study is that we only investigated multisite data featuring roughly equal sample sizes with fairly balanced patient allocation ratios at each site. Multisite studies including a very large number of sites with sometimes only a few subjects per site are however quite common, e.g. the Alzheimer's disease neuroimaging initiative (ADNI) (Mueller et al., 2005) and many pharmaceutical clinical trials at phase II and III 12 . In this type of design, the multisite effect may play a much more pronounced role than in our simulations as it cannot be modeled in the GLM, and will become an intrinsic added source of inter-subject variance. Unfortunately, this type of design could no be tested with the current dataset due to the limited number of sites available. This represents an important avenue of future work.
Underlying causes of the site bias. Not all sites seemed to be equally biased, with sites like Berlin or Saint-Louis showing a small number of connections significantly different then the grand average connectivity matrix, while sites like Baltimore, Queensland and Oxford showed much more biased connectivity measures. These differences may not be statistically significant, or they may reflect real differences due to protocol or scanner characteristics at these sites. Multiple causes may be interacting together to produce the site bias, as reported by Yan et al. (2013b), although some of these sources of variance could be better controlled like the scanner parameters, paired with the use of a phantom to promote more homogeneous configurations across sites Glover et al., 2012). Even in standardized experiments, it should be noted that differences in scanner protocols remain (Brown et al., 2011). A much larger multisite sample with systematically varying parameters could enable a data-driven identification of the critical parameters impacting site bias. The various releases made by the INDI initiative may fill that gap in the literature in the future, as the scanner protocols are much better described in recent releases, such as CoRR (Zuo et al., 2014), than they were in the initial FCP release.

Acknowledgments
Parts of this work were presented at the 2013 annual meeting of the Organization for Human Brain Mapping, as well as the 2013 Alzheimer's Association International Conference (AAIC) (Dansereau et al., 2013). The authors are grateful to the members of the 1000 functional connectome consortium for publicly releasing their datasets. The computational resources used to perform the data analysis were provided by ComputeCanada 13 and CLUMEQ 14 , which is funded in part by NSERC (MRS), FQRNT, and McGill University. This project was funded by NSERC grant number RN000028 and the Canadian Consortium on Neurodegeneration in Aging (CCNA), through a grant from the Canadian Institute of Health Research and funding from several partners including SANOFI-ADVENTIS R&D. PB is supported by a salary award from "Fonds de recherche du Québec -Santé" and the Courtois Foundation.