ComBat Harmonization With and Without Empirical Bayes Estimation for Resting-State Functional Connectivity in Pediatric Mild Traumatic Brain Injury: a CARE4Kids Study

ComBat Harmonization With and Without Empirical Bayes Estimation for Resting-State Functional Connectivity in Pediatric Mild Traumatic Brain Injury: a CARE4Kids Study | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article ComBat Harmonization With and Without Empirical Bayes Estimation for Resting-State Functional Connectivity in Pediatric Mild Traumatic Brain Injury: a CARE4Kids Study Adrian I Onicas¹, Finian Keleher, Kevin C Bickart, Christine L Mac Donald, and 6 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-9418750/v1 This work is licensed under a CC BY 4.0 License Status: Under Review Version 1 posted 4 You are reading this latest preprint version Abstract ComBat is commonly used to mitigate batch effects in multi-site neuroimaging studies. Although empirical bayes estimation is a key component of ComBat, its impact on resting-state functional connectivity measures, including site effect reductions and preservation of original within-site variability remains underexplored. Empirical bayes estimation was evaluated on functional connectivity in adolescents with mild TBI. Resting-state fMRI data were collected from 144 adolescents across six sites participating in the CARE4Kids study. Functional connectivity was computed over resting state networks defined using Seitzman parcellation. ComBat with and without empirical bayes estimation was assessed by quantifying: (1) site variability based on one-way ANOVA, (2) within-site consistency before and after harmonization (intraclass correlation coefficient), (3) age correlations, and (4) site effects of the principal components. All networks demonstrated significant moderate-to-large site effects (all p < 0.001) before harmonization. Harmonization without empirical bayes yielded low residual site effect sizes, whereas empirical bayes harmonization yielded low-to-medium site effects. Within-site consistency between the data before and after harmonization was consistently excellent when using empirical bayes estimation. Harmonization without empirical bayes was excellent across sites with higher sample sizes (n ≥ 13), but moderate-to-excellent for smaller sites (n ≤ 12). Both approaches improved age associations across connections involving the default mode and dorsal attention networks and reduced site effects in the first principal component ( p > 0.05). ComBat with or without empirical bayes effectively removes site variability while preserving age in functional connectivity measures, however, empirical bayes estimation better preserves the original within-site variability. multi-site neuroimaging resting-state fMRI functional connectivity ComBat harmonization empirical bayes mild traumatic brain injury Figures Figure 1 Figure 2 Figure 3 Figure 4 1 Introduction Neuroimaging research in pediatric concussion has been constrained by reliance on small, heterogeneous samples, which have contributed to variability and inconsistencies in reported outcomes [ 1 , 2 ]. Multi-site neuroimaging studies have the potential to address these challenges by involving larger sample sizes with enhanced statistical power, which may improve reproducibility, and enable the identification of neurobiological alterations associated with concussion in developing brains [ 3 , 4 ]. Site-specific variations in imaging protocols and scanner characteristics can introduce batch effects, potentially confounding biological findings and hindering the generalizability of connectivity analyses derived from resting-state functional MRI (rs-fMRI) [ 5 ]. Standardization efforts, such as processing pipelines that enforce uniform preprocessing across sites, help mitigate some sources of variability, but additional statistical harmonization techniques are necessary to account for residual scanner-related effects [ 6 ]. Correcting for site-related variability is necessary to ensure reliable functional connectivity measures, particularly in heterogeneous cohorts such as individuals with traumatic brain injury (TBI). Various harmonization methods have been developed, with ComBat emerging as a widely adopted approach in neuroimaging research [ 7 ]. Originally developed for genomic data [ 8 ], ComBat has been successfully adapted for neuroimaging applications, demonstrating its efficacy in accounting for site effects while preserving biological variability associated to age [ 7 , 9 ]. A key feature of the ComBat algorithm is the Empirical Bayes (EB) estimation, which pools information across features to improve estimation of batch effect parameters [ 8 ]. Specifically, ComBat uses EB to estimate scanner-specific location and scale parameters for each feature in the dataset. Location adjustment involves correcting for differences in the mean of features across different batches, addressing additive batch effects. Scale adjustment focuses on correcting for differences in the variance or spread of features across batches, targeting multiplicative batch effects. Subsequently, ComBat leverages the EB framework to pool information across all features, thereby refining the estimation of these parameters across different sites [ 7 ]. ​​When the EB=FALSE option is used, ComBat functions as a univariate location and scale (L/S) adjustment model, performing L/S adjustments independently for each feature [ 7 ]. The rationale of EB estimation is to stabilize site effect estimates by leveraging information across features, which can be particularly beneficial for batches with low variability or sparse data. ComBat documentation recommends using EB when the number of features is higher than the sample size, and using feature-specific adjustments (i.e., EB=FALSE) otherwise [ 10 ]. While a higher number of features does not imply that sufficient commonality in the batch effects is met to facilitate effective information pooling, average within- and between-network connectivity measures, commonly employed to study brain function after mild TBI [ 11 , 12 ], tend to have similar mean and variance across features. This suggests EB estimation could be appropriate even when the number of features is smaller than the sample size. The influence of EB estimation on ComBat harmonization likely depends on the characteristics of the data [ 8 , 9 ], and EB estimation has been widely applied and tested for features derived using structural neuroimaging in previous ComBat studies [ 7 , 9 , 13 ]. However, its use in functional connectivity remains limited, especially in special populations with heterogeneous data such as pediatric mild TBI. Previous work has raised concerns that multisite harmonization could alter original variability present in structural neuroimaging features, potentially removing meaningful within-site variability [ 14 ]. While the adaptation of ComBat for neuroimaging allows explicit inclusion of additional covariates to preserve biological variability, incorporating complex variables that contribute to site differences unrelated to batch effects can be challenging. Variability between sites may reflect genuine differences in participant characteristics beyond site effects, which may mirror more complex biological, psychological, or demographic characteristics, or heterogeneity of concussion. Assessing within-site consistency can help identify harmonization approaches that better preserve meaningful variability without compromising the removal of site effects. The current study addresses this gap by evaluating ComBat harmonization for within- and between-network functional connectivity measures in a pediatric mild TBI cohort. Specifically, the objective of this study is to evaluate the impact of using the EB estimation in ComBat harmonization on resting-state functional connectivity measures in adolescents with mild traumatic brain injury. We aim to assess how the use of EB affects: (1) the effectiveness of ComBat in reducing site-related variability for individual connectivity measures, (2) the preservation of original within-site variability, and (3) the ability to address multivariate site effects, as evaluated by examining the presence of site-related variance in the principal components from within- and between-network functional connectivity measures. We expect that harmonization without EB will more effectively remove site effects and preserve the original variability of individual features, whereas EB harmonization is anticipated to be more efficient at addressing multivariate site effects in the principal components, by leveraging mean and variance estimates across features derived from functional connectivity. 2 Methods 2.1 Study Design and Procedure Data were drawn from the CARE4Kids study [ 3 ], which employs a longitudinal design to investigate persistent post-concussive symptoms in adolescents. Participants were recruited from emergency departments, primary care clinics, concussion specialty clinics, athletic trainers, and through community outreach. Inclusion criteria include age range between 11 to 17.99 years, a concussion diagnosis by a healthcare provider, and persistence of at least one post-concussion symptom worse than pre-injury baseline within 35 days post-injury. Exclusion criteria include severe developmental delays, significant neurological disorders (with the exception of migraines), severe psychiatric illness, substance abuse, or a history of moderate/severe traumatic brain injury. The study protocol was approved by a central institutional review board at the University of Utah. All participants and their families were provided with detailed study information, ensuring voluntary participation and confidentiality. This multi-center study includes extensive biomarker, neuropsychological, and clinical data collection and was conducted across six consortium sites: University of California, Los Angeles (UCLA), University of Washington/Seattle Children's Research Institute (SEAT), Children's National Medical Center (CNMC), University of Texas Southwestern Medical Center (UTSW), University of Rochester (ROCH), and Atrium Health Wake Forest Baptist Hospital (WAKE). Three data collection points (ROCH, SEAT, and UCLA) were encoded with different sites due to a software update that may introduce variations in image acquisition and reconstruction [ 15 , 16 ]. This resulted in a total of 9 sites for the harmonization analysis (see Table 1 ). Table 1 Descriptive statistics across sites. The final dataset included 144 adolescents [54.9% female; mean (SD) age = 14.96 (1.85) years, mean (SD) time since injury = 21.94 (8.36) days]. DSI – days since injury. Kruskal-Wallis tests were used for continuous variables, and Chi-square tests were used for categorical variables. * p < 0.05 ** p < 0.01, *** p < 0.001 Site n*** (%) Sex [ n (% female)] Age [Mean (SD)] DSI** [Mean (SD)] CNMC 7 (4.9%) 6 (85.7%) 15.26 (2.15) 22.57 (8.92) ROCH1 17 (11.8%) 7 (41.2%) 14.37 (2.15) 24.24 (9.07) ROCH2 5 (3.5%) 4 (80%) 13.14 (2.03) 26 (8.12) SEAT1 12 (8.3%) 8 (66.7%) 14.64 (1.83) 22.25 (7.29) SEAT2 13 (9%) 10 (76.9%) 15.53 (1.96) 18.92 (9.6) UCLA1 8 (5.6%) 3 (37.5%) 15.39 (1.66) 29.25 (5.34) UCLA2 21 (14.6%) 10 (47.6%) 15.1 (1.71) 25.48 (8.27) UTSW 40 (27.8%) 19 (47.5%) 14.98 (1.65) 17.88 (6.46) WAKE 21 (14.6%) 12 (57.1%) 15.24 (1.93) 22 (8.38) 2.2 Image acquisition and preprocessing Structural T1-weighted (T1w) and resting-state functional brain MRI data were acquired at 3 tesla and were processed using the high-performance computing resources provided by the Center for High-Performance Computing (CHPC) at the University of Utah. Data from the baseline assessment (7 to 35 days post-injury) has been used for the purposes of this analysis. Two rs-fMRI runs were obtained for 5.11 min each with a repetition time of 800 ms (383 volumes per run), voxel size of 2.4 mm³, 60 slices, and echo time (TE) of 30 ms. The T1w imaging was performed with a voxel size of 1 mm³, 176–225 slices. Batches SEAT1 and SEAT2 had differing voxel sizes of 2.25 x 2.25 x 2.4 mm. The images preprocessed analyzed using HALFpipe (Harmonized AnaLysis of Functional MRI pipeline) version 1.2.2 [ 17 ], which is based on fMRIPrep version 20.2.7 [ 18 ], following standardized protocols to harmonize analysis and quality control across multiple sites (see http://enigma.ini.usc.edu/protocols/functional-protocols/ ). All denoising steps were performed using the eXtensible Connectivity Pipeline-DCAN (XCP-D) software after resampling the images to standard space [ 19 – 21 ]. Denoising included despiking (AFNI’s 3dDespike), the removal of six motion parameters, their temporal derivatives, and the quadratic expansion of the motion parameters and their temporal derivatives. The six translation and rotation head motion traces were band-stop filtered to remove signals between 12.0 and 20.0 breaths-per-minute using a fourth-order notch filter [ 22 ]. Physiological signals, including white matter (WM), cerebrospinal fluid (CSF), and global signal (GS), were also regressed out. A bandpass filter was applied to retain signal fluctuations within the frequency range of 0.01–0.1 Hz, which corresponds to low-frequency resting-state activity. The same filter was applied to the confound regressors. During quality control, one participant was excluded due to inconsistent acquisition (shorter scan length), five participants were excluded due to poor image processing (brain extraction and/or image registration), and an additional five participants were excluded due to motion artifacts (average framewise displacement > 0.25 mm). Additional 2 participants were excluded due to incomplete demographic information (age, sex or time since injury). The final dataset included 144 adolescents [54.9% female; mean (SD) age = 14.96 (1.85) years, mean (SD) time since injury = 21.94 (8.36) days]. 2.3 Connectome construction and harmonization Time series were extracted using the 300-region Seitzman parcellation [ 23 ], which includes 239 cortical, 34 subcortical, and 27 cerebellar regions, with subcortical and cerebellar regions functionally assigned to cortical networks using the winner-take-all partitioning method based on resting-state connectivity. 9 ROIs (3%) were excluded due to signal dropout. Average within- and between-network FC (i.e., Pearson correlation between ROI pairs) was calculated for the following resting state networks[ 24 ]: Auditory, Cingulo-Opercular, Default Mode, Dorsal Attention, Fronto-Parietal, Salience, Somatomotor (Dorsal and Lateral), Ventral Attention, and Visual (see Supplemental Table 1 for the final set of included ROIs and corresponding functional brain network). The extracted within- and between-network FC values were harmonized for scanner differences using neuroCombat function v1.0.13 in RStudio [ 7 , 10 ], with age at injury and sex included as covariates. Two ComBat versions were applied for data harmonization, with and without the EB estimation. 2.4 Data analysis Statistical analyses were conducted using R studio v4.4.1. To assess the performance of each harmonization approach (i.e., ComBat harmonization with and without empirical Bayes), the effect of site was evaluated through separate one-way ANOVA models for each within- and between-network connectivity measure. Site effects not reaching the statistical threshold ( p > 0.05) were interpreted as evidence that site-related variability had been successfully removed. No multiple comparison correction was applied across measures, which provided a more conservative evaluation for the presence of site effects. The within-scanner consistency was evaluated across connectivity measures before (unharmonized) and after each harmonization approach (before harmonization versus harmonization with EB, and before harmonization versus harmonization without EB) using the intraclass correlation coefficient (ICC). ICC values were interpreted as follows: <0.50 indicating poor consistency, 0.50–0.75 moderate, 0.75–0.90 good, and ≥ 0.90 excellent consistency [ 25 ]. Effective harmonization was expected to reduce site-related effects while preserving the within-scanner variability observed prior to harmonization [ 14 ]. To evaluate whether ComBat harmonization preserves biological variability associated with age, Pearson correlation was used to assess the relationship between age and the harmonized functional connectivity values before and after each harmonization approach. Principal component analysis (PCA) was performed using the prcomp function in R, which computes the singular value decomposition of the centered and scaled data matrix to extract orthogonal components representing uncorrelated sources of variance. The first two principal components (PC1 and PC2), accounting for the largest proportion of variance within each harmonization approach, were retained for further analysis. Consistent with the ANOVA procedures described above, these components were subsequently tested to evaluate potential residual site effects before and after the harmonization procedure. 3 Results 3.1 Presence of site effects All within- and between-network connectivity measures demonstrated significant site effects before harmonization (Fig. 1 and Supplementary Table 2 ). Effect sizes ranged from moderate to large, with connectivity between the salience and the somatomotor lateral network showing the highest effect ( F (8, 135) = 15.67, p < 0.001, Cohen’s f = 0.964), and within the dorsal attention network the lowest [ F (8, 135) = 2.10, p < 0.05, Cohen’s f = 0.353]. Following ComBat harmonization, no site effects remained significant, and effect sizes were low for both harmonization approaches, although harmonization with EB demonstrated larger effect sizes [mean Cohen’s f (SD) = 0.12 (0.04)] compared to harmonization without EB [mean Cohen’s f (SD) = 0.03 (0.02)]. 3.2 Within-site consistency Harmonization with EB demonstrated excellent within-site consistency across all resting state networks and sites [mean (SD) = 0.99 (0.01)]. Moderate to excellent consistency was observed following harmonization without EB (Fig. 2 ; Supplementary Fig. 1 ), ranging from 0.61 to 1 [mean (SD) = 0.97 (0.04)]. 3.3 Correlation between functional connectivity and age Age-correlation improved following harmonization with and without EB compared to the data before harmonization (Fig. 3 , A). The proportion of resting state networks significantly correlated with age increased from 0.02 before harmonization, to 0.29 after either of the harmonization approaches, with patterns of increased correlations present especially among the within and between-network connectivity involving default mode and dorsal attention networks (Fig. 3 , B). 3.4 Presence of site effects in the principal components The first principal component identified before harmonization explained 54.2% of the variance and showed significant site effects ( F (8, 135) = 17.70, p 0.05, Cohen’s f = 0.34). Following harmonization with EB, no site effects were present in the first ( F (8, 135) = 0.02, p > 0.05, Cohen’s f = 0.03) or the second principal component ( F (8, 135) = 0.33, p > 0.05, Cohen’s f = 0.14). Site effects were also below significance level for the first ( F (8, 135) = 0.02, p > 0.05, Cohen’s f = 0.04) and second principal components ( F (8, 135) = 0.01, p > 0.05, Cohen’s f = 0.02) following harmonization without EB. Both harmonization approaches explained less variability based on the first component (38.3% after harmonization with EB and 38.9% after harmonization without EB), but more variability based on the second principal component (15.4.3% after harmonization with EB and 14.2% after harmonization without EB) compared to the data before harmonization (Fig. 4 ). 4 Discussion The increasing availability of large-scale fMRI datasets in pediatric concussion offers unprecedented opportunities to examine brain function during critical developmental periods. In this multi-site study of concussed adolescents, we investigated the impact of using Empirical Bayes (EB) estimation within the ComBat harmonization framework on within and between-network resting-state functional connectivity measures. To our knowledge, this is the first study to apply ComBat harmonization specifically to within and between-network functional connectivity measures, and to assess the effect of EB estimation for neuroimaging data. ComBat harmonization was effective in reducing site effects regardless of whether EB estimation was used. This is consistent with past work showing that ComBat provides efficient means of data harmonization for neuroimaging features [ 7 , 9 , 13 , 14 , 26 , 27 ]. However, subtle differences were observed depending on the version of ComBat applied. When using EB estimation, the effect size of the residual site variability was slightly elevated compared to ComBat without EB, although the persistence of site effects did not reach statistical significance following either harmonization approach. Possibly, EB option’s attempt to balance the variability across features with variability across sites retains subtle residual site-related variability. However, since the residual variability associated with site did not reach statistical significance, it is unlikely that the use of EB estimation could meaningfully bias group-level statistical analyses when pooling information across features. Data harmonization approaches can encounter validation challenges. In the absence of an established ground truth for functional connectivity, we evaluated the extent to which harmonization preserved original data variability by looking at whether site-specific variability is maintained following harmonization. The current results are consistent with prior work examining how ComBat harmonization alters within-site variability in structural connectivity [ 14 ], showing moderate-to-excellent ICC when harmonizing graph theory derived from structural imaging without EB. In the current study, within-site consistency was consistently excellent when harmonizing with EB, indicative of better preservation of the original variability in functional brain networks even across sites with smaller sample sizes (n ≤ 12). This is relevant in developmental cohorts and individuals with TBI where interindividual variability can be pronounced, and suggests that EB estimation may be more appropriate when using within- and between-network connectivity measures. It has been indicated that EB estimation would be suitable when the number of features is high relative to the sample size [ 9 ]. However, our results suggest that EB estimation can successfully remove site effects, while also maintaining the original variability to a greater extent than harmonization without EB. This effect was prominent across sites with low sample size (see Fig. 2 , B and Supplementary Fig. 2 ), enforcing the idea that EB estimation can obtain better estimates for small sample sizes [ 7 , 8 ]. Averages of functional connectivity measures across brain networks tend to distribute similarly and may show relatively higher correlation compared to other features derived from structural imaging, thus EB harmonization might still be appropriate when the number of features is relatively smaller. Future studies should address how the similarity between features can affect the efficiency of harmonizing with and without EB estimation for different neuroimaging measures and modalities, and evaluate the persistence of site effects versus within-site consistency in the context of similarity across features. An overall improvement in the correlation between functional brain networks and age was present in the harmonized networks regardless of whether EB estimation was used. This finding is consistent with prior studies using ComBat harmonization, which have similarly reported enhanced age correlations after adjusting for site-related variability [ 9 , 27 ]. Such increases are typically interpreted as reflecting improved biological signal fidelity in the harmonized data, whereby true neural variability is more accurately preserved as site-related effects are reduced. Notably, we found that age-related correlations were predominantly strengthened across connections within and between the default mode network and the dorsal attention network. This aligns with developmental neuroimaging research demonstrating that functional connectivity within and between these networks systematically changes with age across childhood and adolescence [ 28 , 29 ], indicative of brain maturation that correspond to cognitive development and age-related shifts in attention and self-referential processing [ 30 , 31 ]. Our results thus reinforce the notion that harmonization enhances the ability to detect meaningful age-related neurodevelopmental trajectories by reducing confounding site effects and preserving biologically relevant connectivity patterns. The first principal component derived based on within- and between-network connectivity measures demonstrated site effects to a higher extent than the second component ( Fig. 4 ) , suggesting that site effects constitute a primary source of variability in multisite functional connectivity data. This may obscure meaningful biological variability for example by reducing age-related associations across sites. Site effects were absent from the principal components when PCA was applied on harmonized features following either harmonization approach, consistent with a previous study showing reduced site-dependent clustering following ComBat harmonization in structural MRI [ 13 ]. The presence of site effects in the first principal component before harmonization raises the possibility that sources of multisite variability could be present in the (multivariate) relationships across features. The EB estimation should be beneficial for the removal of site effects within principal components by pulling information across features. However, harmonization with EB estimation did not provide a clear advantage in controlling site-related artifacts in the principal components (see Fig. 4 ). This evaluation may be particularly relevant for studies addressing research questions that involve multivariate group-analysis, where the relative weighting of different features can be important when interpreting the results. Future studies might address whether the use of EB estimation provides advantages regarding how the relative weighting of different features is preserved. In addition, the proportion of explained variance decreased for the first component, but increased for the second component following either ComBat harmonization with or without EB. This potentially reflects more subtle biological variability that was previously overshadowed by batch effects and might have impacted the second component even though the first component already contained most of the variance associated with site effects. It also suggests that applying ComBat harmonization before, rather than after PCA can provide more meaningful components. Several limitations should be considered. Although ComBat harmonization effectively reduced site effects in the current sample, the generalizability of these results to other cohorts, scanner models, or acquisition protocols remains to be confirmed. We focused on within- and between-network connectivity measures, however, other functional connectivity estimates commonly used in TBI research should be evaluated in future studies. In addition, average within- and between-network connectivity features can be correlated, and features derived using other imaging modalities may show different correlation structure and may distribute differently post-harmonization depending on the use of EB estimation. Future studies using other rs-fMRI features or different imaging modalities should consider undergoing independent analyses to confirm the generalizability of the current findings. The efficiency of ComBat harmonization and the impact of EB estimation in datasets spanning broader age ranges, injury severities, and post-injury time points should be considered. Our results suggest a modest trade-off between fully removing site effects and preserving within-site, biologically meaningful variability. Both harmonization approaches provided effective means for removing site effects, with each offering distinct advantages. Standard ComBat with EB estimation better preserved original variability within sites, while ComBat without EB provided slightly stronger harmonization by more aggressively removing site effects, but at the cost of reducing within-site consistency. Future work should continue to explore strategies that balance these goals, particularly in large-scale pediatric neuroimaging studies where developmental variability is a key concern. Declarations Availability of data and materials Data analysis code used in this study is available is the following link: https://anonymous.4open.science/r/C4K_combat_analysis-3D38. Due to data-sharing restrictions, raw data are not publicly available and can be accessed only through the Federal Interagency Traumatic Brain Injury Research (FITBIR) Informatics System (https://fitbir.nih.gov). Qualified researchers may request access to the raw data via FITBIR, subject to appropriate data use agreements and approvals. Authors' contributions Conceptualization: Adrian I. Onicas, Finian Keleher, Kevin C. Bickart, Christine L Mac Donald, Anne Brown, Emily L. Dennis Data Curation: Adrian I. Onicas, Emily L. Dennis Formal Analysis: Adrian I. Onicas, Emily L. Dennis Funding Acquisition: Frederick P. Rivara, Gerard A. Gioia, Christopher C. Giza Investigation: Anne Brown Methodology: Adrian I. Onicas, Finian Keleher, Kevin C. Bickart, Christine L Mac Donald, Anne Brown, Emily L. Dennis Project Administration: Frederick P. Rivara, Gerard A. Gioia, Christopher C. Giza Resources: Lawrence Cook, Emily L. Dennis Software: Adrian I. Onicas, Emily L. Dennis Supervision: Frederick P. Rivara, Gerard A. Gioia, Christopher C. Giza, Emily L. Dennis Validation: Adrian I. Onicas, Finian Keleher, Emily L. Dennis Visualization: Adrian I. Onicas, Emily L. Dennis Writing – Original Draft: Adrian I. Onicas, Emily L. Dennis Writing – Review & Editing: All authors Funding U54NS121688: Christopher C. Giza; RO1NS122184: Emily L. Dennis; R61NS120249: Emily L. Dennis Declaration of Competing Interests Adrian I Onicas reports no disclosures relevant to the manuscript; Finian Keleher reports no disclosures relevant to the manuscript; Kevin C Bickart reports no disclosures relevant to the manuscript; Christine L Mac Donald Christine L Mac Donald has received research funding from the National Institutes of Health (NIH) and the Department of Defense (DOD) unrelated to the current study. Anne Brown reports no disclosures relevant to the manuscript; Lawrence Cook reports no disclosures relevant to the manuscript; Frederick P Rivara reports no disclosures relevant to the manuscript; Gerard A Gioia discloses the following relationships: (1) Grants or contracts: NINDS, CDC Injury Prevention Ctr., NIMH, Johnston Donation; (2) Royalties or licenses: Psychological Assessment Resources, Inc.; (3) Consulting fees: Zogenix Inc., Encoded Therapeutics; (4) Payment for expert testimony: BLANKROME*, Akerman, LLP*; (5) Support for attending meetings or travel: AAN*, SNS; (6) Memberships on panels/organizations/advisory committees related to concussion: USA Football, USA Lacrosse, Positive Coaching Alliance. Relationships marked with an asterisk (*) occurred more than five years ago. Christopher C Giza reports no disclosures relevant to the manuscript; Emily L Dennis reports no disclosures relevant to the manuscript; Acknowledgements This work was conducted on behalf of the Concussion Assessment, Research, and Education for Kids (CARE4Kids) consortium (Giza et al., 2024). We thank all members for their contributions and support. The authors would like to thank the adolescents and families who participated in the CARE4Kids study. References Eierud C, Craddock RC, Fletcher S, et al (2014) Neuroimaging after mild traumatic brain injury: Review and meta-analysis. Neuroimage Clin 4:283–294. https://doi.org/10.1016/j.nicl.2013.12.009 Kashou AW, Frees DM, Kang K, et al (2024) Drivers of resting-state fMRI heterogeneity in traumatic brain injury across injury characteristics and imaging methods: a systematic review and semiquantitative analysis. Front Neurol 15:1487796. https://doi.org/10.3389/fneur.2024.1487796 Giza CC, Gioia G, Cook LJ, et al (2024) CARE4Kids Study: Endophenotypes of Persistent Post-Concussive Symptoms in Adolescents: Study Rationale and Protocol. J Neurotrauma 41:171–185. https://doi.org/10.1089/neu.2023.0073 Yeates KO, Beauchamp M, Craig W, et al (2017) Advancing Concussion Assessment in Pediatrics (A-CAP): a prospective, concurrent cohort, longitudinal study of mild traumatic brain injury in children: study protocol. BMJ open 7:e017012 Noble S, Scheinost D, Finn ES, et al (2017) Multisite Reliability of MR-Based Functional Connectivity. Neuroimage 146:959–970. https://doi.org/10.1016/j.neuroimage.2016.10.020 Jovicich J, Minati L, Marizzoni M, et al (2016) Longitudinal reproducibility of default-mode network connectivity in healthy elderly participants: A multicentric resting-state fMRI study. Neuroimage 124:442–454. https://doi.org/10.1016/j.neuroimage.2015.07.010 Fortin J-P, Parker D, Tunç B, et al (2017) Harmonization of multi-site diffusion tensor imaging data. Neuroimage 161:149–170 Johnson WE, Li C, Rabinovic A (2007) Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 8:118–127 Fortin J-P, Cullen N, Sheline YI, et al (2018) Harmonization of cortical thickness measurements across scanners and sites. Neuroimage 167:104–120 Jfortin1 (2021) Harmonization of multi-site imaging data with ComBat. In: GitHub. https://github.com/Jfortin1/ComBatHarmonization. Accessed 6 Oct 2025 Iyer KK, Barlow KM, Brooks B, et al (2019) Relating brain connectivity with persistent symptoms in pediatric concussion. Annals of Clinical and Translational Neurology 6:954–961. https://doi.org/10.1002/acn3.764 Onicas A, Deighton S, Yeates KO, et al (2025) Brain Network Functional Connectivity in Children With a Concussion. Neurology 104:e213502. https://doi.org/10.1212/WNL.0000000000213502 Tassi E, Bianchi AM, Calesella F, et al (2024) Assessment of ComBat Harmonization Performance on Structural Magnetic Resonance Imaging Measurements. Hum Brain Mapp 45:e70085. https://doi.org/10.1002/hbm.70085 Onicas AI, Ware AL, Harris AD, et al (2022) Multisite Harmonization of Structural DTI Networks in Children: An A-CAP Study. Frontiers in Neurology 13: Medawar E, Thieleking R, Manuilova I, et al (2021) Estimating the effect of a scanner upgrade on measures of grey matter structure for longitudinal designs. PLOS ONE 16:e0239021. https://doi.org/10.1371/journal.pone.0239021 Potvin O, Khademi A, Chouinard I, et al (2019) Measurement Variability Following MRI System Upgrade. Front Neurol 10:726. https://doi.org/10.3389/fneur.2019.00726 Waller L, Erk S, Pozzi E, et al (2022) ENIGMA HALFpipe: Interactive, reproducible, and efficient analysis for resting-state and task-based fMRI data. Hum Brain Mapp 43:2727–2742. https://doi.org/10.1002/hbm.25829 Esteban O, Markiewicz CJ, Blair RW, et al (2019) fMRIPrep: a robust preprocessing pipeline for functional MRI. Nat Methods 16:111–116. https://doi.org/10.1038/s41592-018-0235-4 Ciric R, Wolf DH, Power JD, et al (2017) Benchmarking of participant-level confound regression strategies for the control of motion artifact in studies of functional connectivity. NeuroImage 154:174–187. https://doi.org/10.1016/j.neuroimage.2017.03.020 Mehta K, Salo T, Madison TJ, et al (2024) XCP-D: A robust pipeline for the post-processing of fMRI data. Imaging Neuroscience 2:1–26. https://doi.org/10.1162/imag_a_00257 Satterthwaite TD, Elliott MA, Gerraty RT, et al (2013) An improved framework for confound regression and filtering for control of motion artifact in the preprocessing of resting-state functional connectivity data. Neuroimage 64:240–256. https://doi.org/10.1016/j.neuroimage.2012.08.052 Fair DA, Miranda-Dominguez O, Snyder AZ, et al (2020) Correction of respiratory artifacts in MRI head motion estimates. Neuroimage 208:116400. https://doi.org/10.1016/j.neuroimage.2019.116400 Seitzman BA, Gratton C, Marek S, et al (2020) A set of functionally-defined brain regions with improved representation of the subcortex and cerebellum. Neuroimage 206:116290. https://doi.org/10.1016/j.neuroimage.2019.116290 Thomas Yeo BT, Krienen FM, Sepulcre J, et al (2011) The organization of the human cerebral cortex estimated by intrinsic functional connectivity. J Neurophysiol 106:1125–1165. https://doi.org/10.1152/jn.00338.2011 Koo TK, Li MY (2016) A guideline of selecting and reporting intraclass correlation coefficients for reliability research. Journal of chiropractic medicine 15:155–163 Pomponio R, Erus G, Habes M, et al (2020) Harmonization of large MRI datasets for the analysis of brain imaging patterns throughout the lifespan. NeuroImage 208:116450 Yu M, Linn KA, Cook PA, et al (2018) Statistical harmonization corrects site effects in functional connectivity measurements from multi-site fMRI data. Human brain mapping 39:4213–4227 Sherman LE, Rudie JD, Pfeifer JH, et al (2014) Development of the default mode and central executive networks across early adolescence: a longitudinal study. Dev Cogn Neurosci 10:148–159. https://doi.org/10.1016/j.dcn.2014.08.002 Fair DA, Cohen AL, Dosenbach NUF, et al (2008) The maturing architecture of the brain’s default network. Proceedings of the National Academy of Sciences 105:4028–4032. https://doi.org/10.1073/pnas.0800376105 Sanders AFP, Harms MP, Kandala S, et al (2023) Age-related differences in resting-state functional connectivity from childhood to adolescence. Cereb Cortex 33:6928–6942. https://doi.org/10.1093/cercor/bhad011 Satterthwaite TD, Wolf DH, Ruparel K, et al (2013) Heterogeneous impact of motion on fundamental patterns of developmental changes in functional connectivity during youth. Neuroimage 83:45–57. https://doi.org/10.1016/j.neuroimage.2013.06.045 Additional Declarations Competing interest reported. ● Adrian I Onicas reports no disclosures relevant to the manuscript; ● Finian Keleher reports no disclosures relevant to the manuscript; ● Kevin C Bickart reports no disclosures relevant to the manuscript; ● Christine L Mac Donald Christine L Mac Donald has received research funding from the National Institutes of Health (NIH) and the Department of Defense (DOD) unrelated to the current study. ● Anne Brown reports no disclosures relevant to the manuscript; ● Lawrence Cook reports no disclosures relevant to the manuscript; ● Frederick P Rivara reports no disclosures relevant to the manuscript; ● Gerard A Gioia discloses the following relationships: (1) Grants or contracts: NINDS, CDC Injury Prevention Ctr., NIMH, Johnston Donation; (2) Royalties or licenses: Psychological Assessment Resources, Inc.; (3) Consulting fees: Zogenix Inc., Encoded Therapeutics; (4) Payment for expert testimony: BLANKROME*, Akerman, LLP*; (5) Support for attending meetings or travel: AAN*, SNS; (6) Memberships on panels/organizations/advisory committees related to concussion: USA Football, USA Lacrosse, Positive Coaching Alliance. Relationships marked with an asterisk (*) occurred more than five years ago. ● Christopher C Giza reports no disclosures relevant to the manuscript; ● Emily L Dennis reports no disclosures relevant to the manuscript; Supplementary Files SupplementaryMaterial.docx Cite Share Download PDF Status: Under Review Version 1 posted Reviewers invited by journal 27 Apr, 2026 Editor assigned by journal 22 Apr, 2026 Submission checks completed at journal 22 Apr, 2026 First submitted to journal 14 Apr, 2026 You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-9418750","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":630377355,"identity":"4875f9b5-9827-4ba3-ac63-21d3598531f7","order_by":0,"name":"Adrian I Onicas¹","email":"","orcid":"","institution":"University of Utah School of Medicine","correspondingAuthor":false,"prefix":"","firstName":"Adrian","middleName":"I","lastName":"Onicas¹","suffix":""},{"id":630377356,"identity":"3b523ef3-c708-4ea1-addc-7ad0e65436d4","order_by":1,"name":"Finian Keleher","email":"","orcid":"","institution":"University of Utah School of Medicine","correspondingAuthor":false,"prefix":"","firstName":"Finian","middleName":"","lastName":"Keleher","suffix":""},{"id":630377357,"identity":"024260f2-1e65-416a-adc9-9bb782e091d3","order_by":2,"name":"Kevin C Bickart","email":"","orcid":"","institution":"University of California, Los Angeles","correspondingAuthor":false,"prefix":"","firstName":"Kevin","middleName":"C","lastName":"Bickart","suffix":""},{"id":630377358,"identity":"c412660f-cbfd-46e5-9d40-a4e6ce5c3a16","order_by":3,"name":"Christine L Mac Donald","email":"","orcid":"","institution":"University of Washington School of Medicine","correspondingAuthor":false,"prefix":"","firstName":"Christine","middleName":"L Mac","lastName":"Donald","suffix":""},{"id":630377359,"identity":"01c420d7-afef-40af-8377-b6993e015170","order_by":4,"name":"Anne Brown","email":"","orcid":"","institution":"University of California, Los Angeles","correspondingAuthor":false,"prefix":"","firstName":"Anne","middleName":"","lastName":"Brown","suffix":""},{"id":630377360,"identity":"bcee7ecc-ac49-4465-94c1-6eeee63e3cf2","order_by":5,"name":"Lawrence Cook","email":"","orcid":"","institution":"University of Utah School of Medicine","correspondingAuthor":false,"prefix":"","firstName":"Lawrence","middleName":"","lastName":"Cook","suffix":""},{"id":630377361,"identity":"665ad983-a374-41f5-891e-8a6d6190026c","order_by":6,"name":"Frederick P Rivara","email":"","orcid":"","institution":"University of Washington and Seattle Children's Research Institute","correspondingAuthor":false,"prefix":"","firstName":"Frederick","middleName":"P","lastName":"Rivara","suffix":""},{"id":630377362,"identity":"3935c8cb-258a-4294-b137-67ca4acddc3f","order_by":7,"name":"Gerard A Gioia","email":"","orcid":"","institution":"Children's National Hospital and George Washington University School of Medicine","correspondingAuthor":false,"prefix":"","firstName":"Gerard","middleName":"A","lastName":"Gioia","suffix":""},{"id":630377363,"identity":"63a16c77-1a16-415d-9874-6ae8f8c18fdd","order_by":8,"name":"Christopher C Giza","email":"","orcid":"","institution":"University of California, Los Angeles","correspondingAuthor":false,"prefix":"","firstName":"Christopher","middleName":"C","lastName":"Giza","suffix":""},{"id":630377364,"identity":"be3659a0-fada-4bcb-ae1a-5c72252ea21b","order_by":9,"name":"Emily L Dennis","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABCklEQVRIiWNgGAWjYFACxsYDDAzMIEYDw8cGiBiUi1NLA1wL48wGiAYCWhgYoFqAJC8xWvjFDjcc5mGwljOXbm7+bLujTp5BIv/g5wIGG9kNB7BrkZydCNKSbmw552CbdO6Zw4YNEsnM0jMY0oxxaTG4DdZyOHHDjcQ25ty2AwkMEskM0mARHFrsoVrqgVqaP1u21YG0MP/mYfiPU4uBNERLgsGNxAZpxjZmkBY2oC0HcGqRANpycI5BuiHIYZK9bYcN23gem1nzGCQbz8ShhX92+sMHbyqs5Q1upD/+8LOtTp6fPfHxbZ4KO9k+HFpAgInHAInHBnEwbuUgwPgDv/woGAWjYBSMdAAA9AZdi0DC9ZAAAAAASUVORK5CYII=","orcid":"","institution":"University of Utah School of Medicine","correspondingAuthor":true,"prefix":"","firstName":"Emily","middleName":"L","lastName":"Dennis","suffix":""}],"badges":[],"createdAt":"2026-04-14 18:23:12","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-9418750/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-9418750/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":108630379,"identity":"2fce1b68-6f49-4c3f-9d89-86a5df33270a","added_by":"auto","created_at":"2026-05-06 16:36:30","extension":"jpg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":252068,"visible":true,"origin":"","legend":"\u003cp\u003ePresence of site effects before and after harmonization (Cohen’s \u003cem\u003ef\u003c/em\u003e values).\u003c/p\u003e","description":"","filename":"image1.jpg","url":"https://assets-eu.researchsquare.com/files/rs-9418750/v1/ca6f6792fe3ead41b5fdbf8b.jpg"},{"id":108630381,"identity":"f9dd2072-6cc0-473f-9706-798ed9bc782a","added_by":"auto","created_at":"2026-05-06 16:36:30","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":231552,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003e(A) \u003c/strong\u003eWithin-site consistency (ICC) between functional connectivity measures before and after harmonization. Higher values mean that within-site variability is preserved following harmonization; \u003cstrong\u003e(B) \u003c/strong\u003eWithin-site consistency across functional connectivity measures by site.\u003c/p\u003e","description":"","filename":"image2.png","url":"https://assets-eu.researchsquare.com/files/rs-9418750/v1/5b7540478774bfa3ebae47f0.png"},{"id":108805135,"identity":"0700648a-0729-40e9-9fb9-bd0b60dc5de5","added_by":"auto","created_at":"2026-05-08 15:24:56","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":672220,"visible":true,"origin":"","legend":"\u003cp\u003eAssociation between resting state networks and age.\u003c/p\u003e","description":"","filename":"image3.png","url":"https://assets-eu.researchsquare.com/files/rs-9418750/v1/554f46aaea24f2af60d0abe9.png"},{"id":108805156,"identity":"ceae4d4d-5f58-4918-8fc0-eff3f5a011dc","added_by":"auto","created_at":"2026-05-08 15:25:01","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":486041,"visible":true,"origin":"","legend":"\u003cp\u003eDistribution of the first two principal components before and after harmonization.\u003c/p\u003e","description":"","filename":"image4.png","url":"https://assets-eu.researchsquare.com/files/rs-9418750/v1/787e97ddc89542f88d269ff1.png"},{"id":108809363,"identity":"05eff8a9-f9f2-4660-9dc3-7da9f9e07563","added_by":"auto","created_at":"2026-05-08 15:52:43","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":1545376,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-9418750/v1/670c639b-a78a-4bae-9f58-366ce12379f5.pdf"},{"id":108630378,"identity":"ca3eb8c7-867b-4a4a-b15a-c770d1a78c33","added_by":"auto","created_at":"2026-05-06 16:36:30","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":623018,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryMaterial.docx","url":"https://assets-eu.researchsquare.com/files/rs-9418750/v1/c2b386c7b4f531435e3cad48.docx"}],"financialInterests":"Competing interest reported. ●\tAdrian I Onicas reports no disclosures relevant to the manuscript;\n●\tFinian Keleher reports no disclosures relevant to the manuscript;\n●\tKevin C Bickart reports no disclosures relevant to the manuscript;\n●\tChristine L Mac Donald Christine L Mac Donald has received research funding from the National Institutes of Health (NIH) and the Department of Defense (DOD) unrelated to the current study.\n●\tAnne Brown reports no disclosures relevant to the manuscript;\n●\tLawrence Cook reports no disclosures relevant to the manuscript;\n●\tFrederick P Rivara reports no disclosures relevant to the manuscript;\n●\tGerard A Gioia discloses the following relationships: (1) Grants or contracts: NINDS, CDC Injury Prevention Ctr., NIMH, Johnston Donation; (2) Royalties or licenses: Psychological Assessment Resources, Inc.; (3) Consulting fees: Zogenix Inc., Encoded Therapeutics; (4) Payment for expert testimony: BLANKROME*, Akerman, LLP*; (5) Support for attending meetings or travel: AAN*, SNS; (6) Memberships on panels/organizations/advisory committees related to concussion: USA Football, USA Lacrosse, Positive Coaching Alliance. Relationships marked with an asterisk (*) occurred more than five years ago.\n●\tChristopher C Giza reports no disclosures relevant to the manuscript;\n●\tEmily L Dennis reports no disclosures relevant to the manuscript;","formattedTitle":"ComBat Harmonization With and Without Empirical Bayes Estimation for Resting-State Functional Connectivity in Pediatric Mild Traumatic Brain Injury: a CARE4Kids Study","fulltext":[{"header":"1 Introduction","content":"\u003cp\u003eNeuroimaging research in pediatric concussion has been constrained by reliance on small, heterogeneous samples, which have contributed to variability and inconsistencies in reported outcomes [\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e, \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e]. Multi-site neuroimaging studies have the potential to address these challenges by involving larger sample sizes with enhanced statistical power, which may improve reproducibility, and enable the identification of neurobiological alterations associated with concussion in developing brains [\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e, \u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e]. Site-specific variations in imaging protocols and scanner characteristics can introduce batch effects, potentially confounding biological findings and hindering the generalizability of connectivity analyses derived from resting-state functional MRI (rs-fMRI) [\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e]. Standardization efforts, such as processing pipelines that enforce uniform preprocessing across sites, help mitigate some sources of variability, but additional statistical harmonization techniques are necessary to account for residual scanner-related effects [\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e]. Correcting for site-related variability is necessary to ensure reliable functional connectivity measures, particularly in heterogeneous cohorts such as individuals with traumatic brain injury (TBI).\u003c/p\u003e \u003cp\u003eVarious harmonization methods have been developed, with ComBat emerging as a widely adopted approach in neuroimaging research [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e]. Originally developed for genomic data [\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e], ComBat has been successfully adapted for neuroimaging applications, demonstrating its efficacy in accounting for site effects while preserving biological variability associated to age [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e, \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e]. A key feature of the ComBat algorithm is the Empirical Bayes (EB) estimation, which pools information across features to improve estimation of batch effect parameters [\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e]. Specifically, ComBat uses EB to estimate scanner-specific location and scale parameters for each feature in the dataset. Location adjustment involves correcting for differences in the mean of features across different batches, addressing additive batch effects. Scale adjustment focuses on correcting for differences in the variance or spread of features across batches, targeting multiplicative batch effects. Subsequently, ComBat leverages the EB framework to pool information across all features, thereby refining the estimation of these parameters across different sites [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e]. ​​When the EB=FALSE option is used, ComBat functions as a univariate location and scale (L/S) adjustment model, performing L/S adjustments independently for each feature [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eThe rationale of EB estimation is to stabilize site effect estimates by leveraging information across features, which can be particularly beneficial for batches with low variability or sparse data. ComBat documentation recommends using EB when the number of features is higher than the sample size, and using feature-specific adjustments (i.e., EB=FALSE) otherwise [\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e]. While a higher number of features does not imply that sufficient commonality in the batch effects is met to facilitate effective information pooling, average within- and between-network connectivity measures, commonly employed to study brain function after mild TBI [\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e, \u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e], tend to have similar mean and variance across features. This suggests EB estimation could be appropriate even when the number of features is smaller than the sample size. The influence of EB estimation on ComBat harmonization likely depends on the characteristics of the data [\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e, \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e], and EB estimation has been widely applied and tested for features derived using structural neuroimaging in previous ComBat studies [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e, \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e, \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e]. However, its use in functional connectivity remains limited, especially in special populations with heterogeneous data such as pediatric mild TBI. Previous work has raised concerns that multisite harmonization could alter original variability present in structural neuroimaging features, potentially removing meaningful within-site variability [\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e]. While the adaptation of ComBat for neuroimaging allows explicit inclusion of additional covariates to preserve biological variability, incorporating complex variables that contribute to site differences unrelated to batch effects can be challenging. Variability between sites may reflect genuine differences in participant characteristics beyond site effects, which may mirror more complex biological, psychological, or demographic characteristics, or heterogeneity of concussion. Assessing within-site consistency can help identify harmonization approaches that better preserve meaningful variability without compromising the removal of site effects.\u003c/p\u003e \u003cp\u003eThe current study addresses this gap by evaluating ComBat harmonization for within- and between-network functional connectivity measures in a pediatric mild TBI cohort. Specifically, the objective of this study is to evaluate the impact of using the EB estimation in ComBat harmonization on resting-state functional connectivity measures in adolescents with mild traumatic brain injury. We aim to assess how the use of EB affects: (1) the effectiveness of ComBat in reducing site-related variability for individual connectivity measures, (2) the preservation of original within-site variability, and (3) the ability to address multivariate site effects, as evaluated by examining the presence of site-related variance in the principal components from within- and between-network functional connectivity measures. We expect that harmonization without EB will more effectively remove site effects and preserve the original variability of individual features, whereas EB harmonization is anticipated to be more efficient at addressing multivariate site effects in the principal components, by leveraging mean and variance estimates across features derived from functional connectivity.\u003c/p\u003e"},{"header":"2 Methods","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003e2.1 Study Design and Procedure\u003c/h2\u003e \u003cp\u003eData were drawn from the CARE4Kids study [\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e], which employs a longitudinal design to investigate persistent post-concussive symptoms in adolescents. Participants were recruited from emergency departments, primary care clinics, concussion specialty clinics, athletic trainers, and through community outreach. Inclusion criteria include age range between 11 to 17.99 years, a concussion diagnosis by a healthcare provider, and persistence of at least one post-concussion symptom worse than pre-injury baseline within 35 days post-injury. Exclusion criteria include severe developmental delays, significant neurological disorders (with the exception of migraines), severe psychiatric illness, substance abuse, or a history of moderate/severe traumatic brain injury. The study protocol was approved by a central institutional review board at the University of Utah. All participants and their families were provided with detailed study information, ensuring voluntary participation and confidentiality. This multi-center study includes extensive biomarker, neuropsychological, and clinical data collection and was conducted across six consortium sites: University of California, Los Angeles (UCLA), University of Washington/Seattle Children's Research Institute (SEAT), Children's National Medical Center (CNMC), University of Texas Southwestern Medical Center (UTSW), University of Rochester (ROCH), and Atrium Health Wake Forest Baptist Hospital (WAKE). Three data collection points (ROCH, SEAT, and UCLA) were encoded with different sites due to a software update that may introduce variations in image acquisition and reconstruction [\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e, \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e]. This resulted in a total of 9 sites for the harmonization analysis (see Table\u0026nbsp;\u003cspan refid=\"Tab1\" class=\"InternalRef\"\u003e1\u003c/span\u003e).\u003c/p\u003e \u003cp\u003e \u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab1\" border=\"1\"\u003e \u003ccaption language=\"En\"\u003e \u003cdiv class=\"CaptionNumber\"\u003eTable 1\u003c/div\u003e \u003cdiv class=\"CaptionContent\"\u003e \u003cp\u003eDescriptive statistics across sites. The final dataset included 144 adolescents [54.9% female; mean (SD) age\u0026thinsp;=\u0026thinsp;14.96 (1.85) years, mean (SD) time since injury\u0026thinsp;=\u0026thinsp;21.94 (8.36) days]. DSI \u0026ndash; days since injury. Kruskal-Wallis tests were used for continuous variables, and Chi-square tests were used for categorical variables. * p\u0026thinsp;\u0026lt;\u0026thinsp;0.05 ** p\u0026thinsp;\u0026lt;\u0026thinsp;0.01, *** p\u0026thinsp;\u0026lt;\u0026thinsp;0.001\u003c/p\u003e \u003c/div\u003e \u003c/caption\u003e \u003ccolgroup cols=\"5\"\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c5\" colnum=\"5\"\u003e\u003c/div\u003e \u003cthead\u003e \u003ctr\u003e \u003cth align=\"left\" colname=\"c1\"\u003e \u003cp\u003eSite\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c2\"\u003e \u003cp\u003e\u003cem\u003en***\u003c/em\u003e (%)\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c3\"\u003e \u003cp\u003eSex [\u003cem\u003en\u003c/em\u003e (% female)]\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c4\"\u003e \u003cp\u003eAge [Mean (SD)]\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c5\"\u003e \u003cp\u003eDSI** [Mean (SD)]\u003c/p\u003e \u003c/th\u003e \u003c/tr\u003e \u003c/thead\u003e \u003ctbody\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eCNMC\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e7 (4.9%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e6 (85.7%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e15.26 (2.15)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e22.57 (8.92)\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eROCH1\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e17 (11.8%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e7 (41.2%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e14.37 (2.15)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e24.24 (9.07)\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eROCH2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e5 (3.5%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e4 (80%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e13.14 (2.03)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e26 (8.12)\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eSEAT1\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e12 (8.3%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e8 (66.7%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e14.64 (1.83)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e22.25 (7.29)\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eSEAT2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e13 (9%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e10 (76.9%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e15.53 (1.96)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e18.92 (9.6)\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eUCLA1\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e8 (5.6%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e3 (37.5%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e15.39 (1.66)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e29.25 (5.34)\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eUCLA2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e21 (14.6%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e10 (47.6%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e15.1 (1.71)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e25.48 (8.27)\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eUTSW\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e40 (27.8%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e19 (47.5%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e14.98 (1.65)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e17.88 (6.46)\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eWAKE\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003e21 (14.6%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e12 (57.1%)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e15.24 (1.93)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c5\"\u003e \u003cp\u003e22 (8.38)\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003c/tbody\u003e \u003c/colgroup\u003e \u003c/table\u003e\u003c/div\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec4\" class=\"Section2\"\u003e \u003ch2\u003e2.2 Image acquisition and preprocessing\u003c/h2\u003e \u003cp\u003eStructural T1-weighted (T1w) and resting-state functional brain MRI data were acquired at 3 tesla and were processed using the high-performance computing resources provided by the Center for High-Performance Computing (CHPC) at the University of Utah. Data from the baseline assessment (7 to 35 days post-injury) has been used for the purposes of this analysis. Two rs-fMRI runs were obtained for 5.11 min each with a repetition time of 800 ms (383 volumes per run), voxel size of 2.4 mm\u0026sup3;, 60 slices, and echo time (TE) of 30 ms. The T1w imaging was performed with a voxel size of 1 mm\u0026sup3;, 176\u0026ndash;225 slices. Batches SEAT1 and SEAT2 had differing voxel sizes of 2.25 x 2.25 x 2.4 mm. The images preprocessed analyzed using HALFpipe (Harmonized AnaLysis of Functional MRI pipeline) version 1.2.2 [\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e], which is based on fMRIPrep version 20.2.7 [\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e], following standardized protocols to harmonize analysis and quality control across multiple sites (see \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://enigma.ini.usc.edu/protocols/functional-protocols/\u003c/span\u003e\u003cspan address=\"http://enigma.ini.usc.edu/protocols/functional-protocols/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003e).\u003c/span\u003e\u003c/p\u003e \u003cp\u003eAll denoising steps were performed using the eXtensible Connectivity Pipeline-DCAN (XCP-D) software after resampling the images to standard space [\u003cspan additionalcitationids=\"CR20\" citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e]. Denoising included despiking (AFNI\u0026rsquo;s 3dDespike), the removal of six motion parameters, their temporal derivatives, and the quadratic expansion of the motion parameters and their temporal derivatives. The six translation and rotation head motion traces were band-stop filtered to remove signals between 12.0 and 20.0 breaths-per-minute using a fourth-order notch filter [\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e]. Physiological signals, including white matter (WM), cerebrospinal fluid (CSF), and global signal (GS), were also regressed out. A bandpass filter was applied to retain signal fluctuations within the frequency range of 0.01\u0026ndash;0.1 Hz, which corresponds to low-frequency resting-state activity. The same filter was applied to the confound regressors.\u003c/p\u003e \u003cp\u003eDuring quality control, one participant was excluded due to inconsistent acquisition (shorter scan length), five participants were excluded due to poor image processing (brain extraction and/or image registration), and an additional five participants were excluded due to motion artifacts (average framewise displacement\u0026thinsp;\u0026gt;\u0026thinsp;0.25 mm). Additional 2 participants were excluded due to incomplete demographic information (age, sex or time since injury). The final dataset included 144 adolescents [54.9% female; mean (SD) age\u0026thinsp;=\u0026thinsp;14.96 (1.85) years, mean (SD) time since injury\u0026thinsp;=\u0026thinsp;21.94 (8.36) days].\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec5\" class=\"Section2\"\u003e \u003ch2\u003e2.3 Connectome construction and harmonization\u003c/h2\u003e \u003cp\u003eTime series were extracted using the 300-region Seitzman parcellation [\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e], which includes 239 cortical, 34 subcortical, and 27 cerebellar regions, with subcortical and cerebellar regions functionally assigned to cortical networks using the winner-take-all partitioning method based on resting-state connectivity.\u003c/p\u003e \u003cp\u003e9 ROIs (3%) were excluded due to signal dropout. Average within- and between-network FC (i.e., Pearson correlation between ROI pairs) was calculated for the following resting state networks[\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e]: Auditory, Cingulo-Opercular, Default Mode, Dorsal Attention, Fronto-Parietal, Salience, Somatomotor (Dorsal and Lateral), Ventral Attention, and Visual (see \u003cb\u003eSupplemental Table\u0026nbsp;1\u003c/b\u003e for the final set of included ROIs and corresponding functional brain network). The extracted within- and between-network FC values were harmonized for scanner differences using neuroCombat function v1.0.13 in RStudio [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e, \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e], with age at injury and sex included as covariates. Two ComBat versions were applied for data harmonization, with and without the EB estimation.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec6\" class=\"Section2\"\u003e \u003ch2\u003e2.4 Data analysis\u003c/h2\u003e \u003cp\u003eStatistical analyses were conducted using R studio v4.4.1. To assess the performance of each harmonization approach (i.e., ComBat harmonization with and without empirical Bayes), the effect of site was evaluated through separate one-way ANOVA models for each within- and between-network connectivity measure. Site effects not reaching the statistical threshold (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026gt;\u0026thinsp;0.05) were interpreted as evidence that site-related variability had been successfully removed. No multiple comparison correction was applied across measures, which provided a more conservative evaluation for the presence of site effects.\u003c/p\u003e \u003cp\u003eThe within-scanner consistency was evaluated across connectivity measures before (unharmonized) and after each harmonization approach (before harmonization versus harmonization with EB, and before harmonization versus harmonization without EB) using the intraclass correlation coefficient (ICC). ICC values were interpreted as follows: \u0026lt;0.50 indicating poor consistency, 0.50\u0026ndash;0.75 moderate, 0.75\u0026ndash;0.90 good, and \u0026ge;\u0026thinsp;0.90 excellent consistency [\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e]. Effective harmonization was expected to reduce site-related effects while preserving the within-scanner variability observed prior to harmonization [\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eTo evaluate whether ComBat harmonization preserves biological variability associated with age, Pearson correlation was used to assess the relationship between age and the harmonized functional connectivity values before and after each harmonization approach.\u003c/p\u003e \u003cp\u003ePrincipal component analysis (PCA) was performed using the prcomp function in R, which computes the singular value decomposition of the centered and scaled data matrix to extract orthogonal components representing uncorrelated sources of variance. The first two principal components (PC1 and PC2), accounting for the largest proportion of variance within each harmonization approach, were retained for further analysis. Consistent with the ANOVA procedures described above, these components were subsequently tested to evaluate potential residual site effects before and after the harmonization procedure.\u003c/p\u003e \u003c/div\u003e"},{"header":"3 Results","content":"\u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003e3.1 Presence of site effects\u003c/h2\u003e \u003cp\u003eAll within- and between-network connectivity measures demonstrated significant site effects before harmonization (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e and \u003cb\u003eSupplementary Table\u0026nbsp;2\u003c/b\u003e). Effect sizes ranged from moderate to large, with connectivity between the salience and the somatomotor lateral network showing the highest effect (\u003cem\u003eF\u003c/em\u003e(8, 135)\u0026thinsp;=\u0026thinsp;15.67, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.001, Cohen\u0026rsquo;s \u003cem\u003ef\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.964), and within the dorsal attention network the lowest [\u003cem\u003eF\u003c/em\u003e (8, 135)\u0026thinsp;=\u0026thinsp;2.10, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05, Cohen\u0026rsquo;s \u003cem\u003ef\u0026thinsp;=\u003c/em\u003e\u0026thinsp;0.353]. Following ComBat harmonization, no site effects remained significant, and effect sizes were low for both harmonization approaches, although harmonization with EB demonstrated larger effect sizes [mean Cohen\u0026rsquo;s \u003cem\u003ef\u003c/em\u003e (SD)\u0026thinsp;=\u0026thinsp;0.12 (0.04)] compared to harmonization without EB [mean Cohen\u0026rsquo;s \u003cem\u003ef\u003c/em\u003e (SD)\u0026thinsp;=\u0026thinsp;0.03 (0.02)].\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec9\" class=\"Section2\"\u003e \u003ch2\u003e3.2 Within-site consistency\u003c/h2\u003e \u003cp\u003eHarmonization with EB demonstrated excellent within-site consistency across all resting state networks and sites [mean (SD)\u0026thinsp;=\u0026thinsp;0.99 (0.01)]. Moderate to excellent consistency was observed following harmonization without EB (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e; \u003cb\u003eSupplementary Fig.\u0026nbsp;1\u003c/b\u003e), ranging from 0.61 to 1 [mean (SD)\u0026thinsp;=\u0026thinsp;0.97 (0.04)].\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec10\" class=\"Section2\"\u003e \u003ch2\u003e3.3 Correlation between functional connectivity and age\u003c/h2\u003e \u003cp\u003eAge-correlation improved following harmonization with and without EB compared to the data before harmonization (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e, A). The proportion of resting state networks significantly correlated with age increased from 0.02 before harmonization, to 0.29 after either of the harmonization approaches, with patterns of increased correlations present especially among the within and between-network connectivity involving default mode and dorsal attention networks (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e, B).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003e3.4 Presence of site effects in the principal components\u003c/h2\u003e \u003cp\u003eThe first principal component identified before harmonization explained 54.2% of the variance and showed significant site effects (\u003cem\u003eF\u003c/em\u003e(8, 135)\u0026thinsp;=\u0026thinsp;17.70, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.001, Cohen\u0026rsquo;s \u003cem\u003ef\u003c/em\u003e\u0026thinsp;=\u0026thinsp;1.02). The second component explained 11.7% of the variance and showed no significant site effects (\u003cem\u003eF\u003c/em\u003e(8, 135)\u0026thinsp;=\u0026thinsp;2.00, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026gt;\u0026thinsp;0.05, Cohen\u0026rsquo;s \u003cem\u003ef\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.34). Following harmonization with EB, no site effects were present in the first (\u003cem\u003eF\u003c/em\u003e(8, 135)\u0026thinsp;=\u0026thinsp;0.02, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026gt;\u0026thinsp;0.05, Cohen\u0026rsquo;s \u003cem\u003ef\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.03) or the second principal component (\u003cem\u003eF\u003c/em\u003e(8, 135)\u0026thinsp;=\u0026thinsp;0.33, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026gt;\u0026thinsp;0.05, Cohen\u0026rsquo;s \u003cem\u003ef\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.14). Site effects were also below significance level for the first (\u003cem\u003eF\u003c/em\u003e(8, 135)\u0026thinsp;=\u0026thinsp;0.02, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026gt;\u0026thinsp;0.05, Cohen\u0026rsquo;s \u003cem\u003ef\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.04) and second principal components (\u003cem\u003eF\u003c/em\u003e(8, 135)\u0026thinsp;=\u0026thinsp;0.01, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026gt;\u0026thinsp;0.05, Cohen\u0026rsquo;s \u003cem\u003ef\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.02) following harmonization without EB. Both harmonization approaches explained less variability based on the first component (38.3% after harmonization with EB and 38.9% after harmonization without EB), but more variability based on the second principal component (15.4.3% after harmonization with EB and 14.2% after harmonization without EB) compared to the data before harmonization (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e"},{"header":"4 Discussion","content":"\u003cp\u003eThe increasing availability of large-scale fMRI datasets in pediatric concussion offers unprecedented opportunities to examine brain function during critical developmental periods. In this multi-site study of concussed adolescents, we investigated the impact of using Empirical Bayes (EB) estimation within the ComBat harmonization framework on within and between-network resting-state functional connectivity measures. To our knowledge, this is the first study to apply ComBat harmonization specifically to within and between-network functional connectivity measures, and to assess the effect of EB estimation for neuroimaging data.\u003c/p\u003e \u003cp\u003eComBat harmonization was effective in reducing site effects regardless of whether EB estimation was used. This is consistent with past work showing that ComBat provides efficient means of data harmonization for neuroimaging features [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e, \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e, \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e, \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e, \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e, \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e]. However, subtle differences were observed depending on the version of ComBat applied. When using EB estimation, the effect size of the residual site variability was slightly elevated compared to ComBat without EB, although the persistence of site effects did not reach statistical significance following either harmonization approach. Possibly, EB option\u0026rsquo;s attempt to balance the variability across features with variability across sites retains subtle residual site-related variability. However, since the residual variability associated with site did not reach statistical significance, it is unlikely that the use of EB estimation could meaningfully bias group-level statistical analyses when pooling information across features.\u003c/p\u003e \u003cp\u003eData harmonization approaches can encounter validation challenges. In the absence of an established ground truth for functional connectivity, we evaluated the extent to which harmonization preserved original data variability by looking at whether site-specific variability is maintained following harmonization. The current results are consistent with prior work examining how ComBat harmonization alters within-site variability in structural connectivity [\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e], showing moderate-to-excellent ICC when harmonizing graph theory derived from structural imaging without EB. In the current study, within-site consistency was consistently excellent when harmonizing with EB, indicative of better preservation of the original variability in functional brain networks even across sites with smaller sample sizes (n\u0026thinsp;\u0026le;\u0026thinsp;12). This is relevant in developmental cohorts and individuals with TBI where interindividual variability can be pronounced, and suggests that EB estimation may be more appropriate when using within- and between-network connectivity measures. It has been indicated that EB estimation would be suitable when the number of features is high relative to the sample size [\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e]. However, our results suggest that EB estimation can successfully remove site effects, while also maintaining the original variability to a greater extent than harmonization without EB. This effect was prominent across sites with low sample size (see Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e, B \u003cb\u003eand Supplementary Fig.\u0026nbsp;2\u003c/b\u003e), enforcing the idea that EB estimation can obtain better estimates for small sample sizes [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e, \u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e]. Averages of functional connectivity measures across brain networks tend to distribute similarly and may show relatively higher correlation compared to other features derived from structural imaging, thus EB harmonization might still be appropriate when the number of features is relatively smaller. Future studies should address how the similarity between features can affect the efficiency of harmonizing with and without EB estimation for different neuroimaging measures and modalities, and evaluate the persistence of site effects versus within-site consistency in the context of similarity across features.\u003c/p\u003e \u003cp\u003eAn overall improvement in the correlation between functional brain networks and age was present in the harmonized networks regardless of whether EB estimation was used. This finding is consistent with prior studies using ComBat harmonization, which have similarly reported enhanced age correlations after adjusting for site-related variability [\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e, \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e]. Such increases are typically interpreted as reflecting improved biological signal fidelity in the harmonized data, whereby true neural variability is more accurately preserved as site-related effects are reduced. Notably, we found that age-related correlations were predominantly strengthened across connections within and between the default mode network and the dorsal attention network. This aligns with developmental neuroimaging research demonstrating that functional connectivity within and between these networks systematically changes with age across childhood and adolescence [\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e, \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e], indicative of brain maturation that correspond to cognitive development and age-related shifts in attention and self-referential processing [\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e, \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e]. Our results thus reinforce the notion that harmonization enhances the ability to detect meaningful age-related neurodevelopmental trajectories by reducing confounding site effects and preserving biologically relevant connectivity patterns.\u003c/p\u003e \u003cp\u003eThe first principal component derived based on within- and between-network connectivity measures demonstrated site effects to a higher extent than the second component \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e\u003cb\u003e)\u003c/b\u003e, suggesting that site effects constitute a primary source of variability in multisite functional connectivity data. This may obscure meaningful biological variability for example by reducing age-related associations across sites. Site effects were absent from the principal components when PCA was applied on harmonized features following either harmonization approach, consistent with a previous study showing reduced site-dependent clustering following ComBat harmonization in structural MRI [\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e]. The presence of site effects in the first principal component before harmonization raises the possibility that sources of multisite variability could be present in the (multivariate) relationships across features. The EB estimation should be beneficial for the removal of site effects within principal components by pulling information across features. However, harmonization with EB estimation did not provide a clear advantage in controlling site-related artifacts in the principal components (see Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e). This evaluation may be particularly relevant for studies addressing research questions that involve multivariate group-analysis, where the relative weighting of different features can be important when interpreting the results. Future studies might address whether the use of EB estimation provides advantages regarding how the relative weighting of different features is preserved. In addition, the proportion of explained variance decreased for the first component, but increased for the second component following either ComBat harmonization with or without EB. This potentially reflects more subtle biological variability that was previously overshadowed by batch effects and might have impacted the second component even though the first component already contained most of the variance associated with site effects. It also suggests that applying ComBat harmonization before, rather than after PCA can provide more meaningful components.\u003c/p\u003e \u003cp\u003eSeveral limitations should be considered. Although ComBat harmonization effectively reduced site effects in the current sample, the generalizability of these results to other cohorts, scanner models, or acquisition protocols remains to be confirmed. We focused on within- and between-network connectivity measures, however, other functional connectivity estimates commonly used in TBI research should be evaluated in future studies. In addition, average within- and between-network connectivity features can be correlated, and features derived using other imaging modalities may show different correlation structure and may distribute differently post-harmonization depending on the use of EB estimation. Future studies using other rs-fMRI features or different imaging modalities should consider undergoing independent analyses to confirm the generalizability of the current findings. The efficiency of ComBat harmonization and the impact of EB estimation in datasets spanning broader age ranges, injury severities, and post-injury time points should be considered.\u003c/p\u003e \u003cp\u003eOur results suggest a modest trade-off between fully removing site effects and preserving within-site, biologically meaningful variability. Both harmonization approaches provided effective means for removing site effects, with each offering distinct advantages. Standard ComBat with EB estimation better preserved original variability within sites, while ComBat without EB provided slightly stronger harmonization by more aggressively removing site effects, but at the cost of reducing within-site consistency. Future work should continue to explore strategies that balance these goals, particularly in large-scale pediatric neuroimaging studies where developmental variability is a key concern.\u003c/p\u003e"},{"header":"Declarations","content":"\u003ch2\u003eAvailability of data and materials\u003c/h2\u003e\n\u003cp\u003eData analysis code used in this study is available is the following link: https://anonymous.4open.science/r/C4K_combat_analysis-3D38. Due to data-sharing restrictions, raw data are not publicly available and can be accessed only through the Federal Interagency Traumatic Brain Injury Research (FITBIR) Informatics System (https://fitbir.nih.gov). Qualified researchers may request access to the raw data via FITBIR, subject to appropriate data use agreements and approvals.\u003c/p\u003e\n\u003ch2\u003eAuthors\u0026apos; contributions\u003c/h2\u003e\n\u003cul class=\"decimal_type\"\u003e\n \u003cli\u003eConceptualization: Adrian I. Onicas, Finian Keleher, Kevin C. Bickart, Christine L Mac Donald, Anne Brown, Emily L. Dennis\u003c/li\u003e\n \u003cli\u003eData Curation: Adrian I. Onicas, Emily L. Dennis\u003c/li\u003e\n \u003cli\u003eFormal Analysis: Adrian I. Onicas, Emily L. Dennis\u003c/li\u003e\n \u003cli\u003eFunding Acquisition: Frederick P. Rivara, Gerard A. Gioia, Christopher C. Giza\u003c/li\u003e\n \u003cli\u003eInvestigation: Anne Brown\u003c/li\u003e\n \u003cli\u003eMethodology: Adrian I. Onicas, Finian Keleher, Kevin C. Bickart, Christine L Mac Donald, Anne Brown, Emily L. Dennis\u003c/li\u003e\n \u003cli\u003eProject Administration: Frederick P. Rivara, Gerard A. Gioia, Christopher C. Giza\u003c/li\u003e\n \u003cli\u003eResources: Lawrence Cook, Emily L. Dennis\u003c/li\u003e\n \u003cli\u003eSoftware: Adrian I. Onicas, Emily L. Dennis\u003c/li\u003e\n \u003cli\u003eSupervision: Frederick P. Rivara, Gerard A. Gioia, Christopher C. Giza, Emily L. Dennis\u003c/li\u003e\n \u003cli\u003eValidation: Adrian I. Onicas, Finian Keleher, Emily L. Dennis\u003c/li\u003e\n \u003cli\u003eVisualization: Adrian I. Onicas, Emily L. Dennis\u003c/li\u003e\n \u003cli\u003eWriting \u0026ndash; Original Draft: Adrian I. Onicas, Emily L. Dennis\u003c/li\u003e\n \u003cli\u003eWriting \u0026ndash; Review \u0026amp; Editing: All authors\u003c/li\u003e\n\u003c/ul\u003e\n\u003ch2\u003eFunding\u003c/h2\u003e\n\u003cp\u003eU54NS121688: Christopher C. Giza; RO1NS122184: Emily L. Dennis; R61NS120249: Emily L. Dennis\u003c/p\u003e\n\u003ch2\u003eDeclaration of Competing Interests\u003c/h2\u003e\n\u003cul class=\"decimal_type\"\u003e\n \u003cli\u003eAdrian I Onicas reports no disclosures relevant to the manuscript;\u003c/li\u003e\n \u003cli\u003eFinian \u0026nbsp; Keleher reports no disclosures relevant to the manuscript;\u003c/li\u003e\n \u003cli\u003eKevin C Bickart reports no disclosures relevant to the manuscript;\u003c/li\u003e\n \u003cli\u003eChristine L Mac Donald Christine L Mac Donald has received research funding from the National Institutes of Health (NIH) and the Department of Defense (DOD) unrelated to the current study.\u003c/li\u003e\n \u003cli\u003eAnne Brown reports no disclosures relevant to the manuscript;\u003c/li\u003e\n \u003cli\u003eLawrence \u0026nbsp; Cook reports no disclosures relevant to the manuscript;\u003c/li\u003e\n \u003cli\u003eFrederick P Rivara reports no disclosures relevant to the manuscript;\u003c/li\u003e\n \u003cli\u003eGerard A Gioia discloses the following relationships: (1) Grants or contracts: NINDS, CDC Injury Prevention Ctr., NIMH, Johnston Donation; (2) Royalties or licenses: Psychological Assessment Resources, Inc.; (3) Consulting fees: Zogenix Inc., Encoded Therapeutics; (4) Payment for expert testimony: BLANKROME*, Akerman, LLP*; (5) Support for attending meetings or travel: AAN*, SNS; (6) Memberships on panels/organizations/advisory committees related to concussion: USA Football, USA Lacrosse, Positive Coaching Alliance. Relationships marked with an asterisk (*) occurred more than five years ago.\u003c/li\u003e\n \u003cli\u003eChristopher C Giza reports no disclosures relevant to the manuscript;\u003c/li\u003e\n \u003cli\u003eEmily L Dennis reports no disclosures relevant to the manuscript;\u003c/li\u003e\n\u003c/ul\u003e\n\u003ch2\u003eAcknowledgements\u003c/h2\u003e\n\u003cp\u003eThis work was conducted on behalf of the Concussion Assessment, Research, and Education for Kids (CARE4Kids) consortium (Giza et al., 2024). We thank all members for their contributions and support. The authors would like to thank the adolescents and families who participated in the CARE4Kids study.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eEierud C, Craddock RC, Fletcher S, et al (2014) Neuroimaging after mild traumatic brain injury: Review and meta-analysis. Neuroimage Clin 4:283\u0026ndash;294. https://doi.org/10.1016/j.nicl.2013.12.009 \u003c/li\u003e\n\u003cli\u003eKashou AW, Frees DM, Kang K, et al (2024) Drivers of resting-state fMRI heterogeneity in traumatic brain injury across injury characteristics and imaging methods: a systematic review and semiquantitative analysis. Front Neurol 15:1487796. https://doi.org/10.3389/fneur.2024.1487796 \u003c/li\u003e\n\u003cli\u003eGiza CC, Gioia G, Cook LJ, et al (2024) CARE4Kids Study: Endophenotypes of Persistent Post-Concussive Symptoms in Adolescents: Study Rationale and Protocol. J Neurotrauma 41:171\u0026ndash;185. https://doi.org/10.1089/neu.2023.0073 \u003c/li\u003e\n\u003cli\u003eYeates KO, Beauchamp M, Craig W, et al (2017) Advancing Concussion Assessment in Pediatrics (A-CAP): a prospective, concurrent cohort, longitudinal study of mild traumatic brain injury in children: study protocol. BMJ open 7:e017012 \u003c/li\u003e\n\u003cli\u003eNoble S, Scheinost D, Finn ES, et al (2017) Multisite Reliability of MR-Based Functional Connectivity. Neuroimage 146:959\u0026ndash;970. https://doi.org/10.1016/j.neuroimage.2016.10.020 \u003c/li\u003e\n\u003cli\u003eJovicich J, Minati L, Marizzoni M, et al (2016) Longitudinal reproducibility of default-mode network connectivity in healthy elderly participants: A multicentric resting-state fMRI study. Neuroimage 124:442\u0026ndash;454. https://doi.org/10.1016/j.neuroimage.2015.07.010 \u003c/li\u003e\n\u003cli\u003eFortin J-P, Parker D, Tun\u0026ccedil; B, et al (2017) Harmonization of multi-site diffusion tensor imaging data. Neuroimage 161:149\u0026ndash;170 \u003c/li\u003e\n\u003cli\u003eJohnson WE, Li C, Rabinovic A (2007) Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 8:118\u0026ndash;127 \u003c/li\u003e\n\u003cli\u003eFortin J-P, Cullen N, Sheline YI, et al (2018) Harmonization of cortical thickness measurements across scanners and sites. Neuroimage 167:104\u0026ndash;120 \u003c/li\u003e\n\u003cli\u003eJfortin1 (2021) Harmonization of multi-site imaging data with ComBat. In: GitHub. https://github.com/Jfortin1/ComBatHarmonization. Accessed 6 Oct 2025 \u003c/li\u003e\n\u003cli\u003eIyer KK, Barlow KM, Brooks B, et al (2019) Relating brain connectivity with persistent symptoms in pediatric concussion. Annals of Clinical and Translational Neurology 6:954\u0026ndash;961. https://doi.org/10.1002/acn3.764 \u003c/li\u003e\n\u003cli\u003eOnicas A, Deighton S, Yeates KO, et al (2025) Brain Network Functional Connectivity in Children With a Concussion. Neurology 104:e213502. https://doi.org/10.1212/WNL.0000000000213502 \u003c/li\u003e\n\u003cli\u003eTassi E, Bianchi AM, Calesella F, et al (2024) Assessment of ComBat Harmonization Performance on Structural Magnetic Resonance Imaging Measurements. Hum Brain Mapp 45:e70085. https://doi.org/10.1002/hbm.70085 \u003c/li\u003e\n\u003cli\u003eOnicas AI, Ware AL, Harris AD, et al (2022) Multisite Harmonization of Structural DTI Networks in Children: An A-CAP Study. Frontiers in Neurology 13: \u003c/li\u003e\n\u003cli\u003eMedawar E, Thieleking R, Manuilova I, et al (2021) Estimating the effect of a scanner upgrade on measures of grey matter structure for longitudinal designs. PLOS ONE 16:e0239021. https://doi.org/10.1371/journal.pone.0239021 \u003c/li\u003e\n\u003cli\u003ePotvin O, Khademi A, Chouinard I, et al (2019) Measurement Variability Following MRI System Upgrade. Front Neurol 10:726. https://doi.org/10.3389/fneur.2019.00726 \u003c/li\u003e\n\u003cli\u003eWaller L, Erk S, Pozzi E, et al (2022) ENIGMA HALFpipe: Interactive, reproducible, and efficient analysis for resting-state and task-based fMRI data. Hum Brain Mapp 43:2727\u0026ndash;2742. https://doi.org/10.1002/hbm.25829 \u003c/li\u003e\n\u003cli\u003eEsteban O, Markiewicz CJ, Blair RW, et al (2019) fMRIPrep: a robust preprocessing pipeline for functional MRI. Nat Methods 16:111\u0026ndash;116. https://doi.org/10.1038/s41592-018-0235-4 \u003c/li\u003e\n\u003cli\u003eCiric R, Wolf DH, Power JD, et al (2017) Benchmarking of participant-level confound regression strategies for the control of motion artifact in studies of functional connectivity. NeuroImage 154:174\u0026ndash;187. https://doi.org/10.1016/j.neuroimage.2017.03.020 \u003c/li\u003e\n\u003cli\u003eMehta K, Salo T, Madison TJ, et al (2024) XCP-D: A robust pipeline for the post-processing of fMRI data. Imaging Neuroscience 2:1\u0026ndash;26. https://doi.org/10.1162/imag_a_00257 \u003c/li\u003e\n\u003cli\u003eSatterthwaite TD, Elliott MA, Gerraty RT, et al (2013) An improved framework for confound regression and filtering for control of motion artifact in the preprocessing of resting-state functional connectivity data. Neuroimage 64:240\u0026ndash;256. https://doi.org/10.1016/j.neuroimage.2012.08.052 \u003c/li\u003e\n\u003cli\u003eFair DA, Miranda-Dominguez O, Snyder AZ, et al (2020) Correction of respiratory artifacts in MRI head motion estimates. Neuroimage 208:116400. https://doi.org/10.1016/j.neuroimage.2019.116400 \u003c/li\u003e\n\u003cli\u003eSeitzman BA, Gratton C, Marek S, et al (2020) A set of functionally-defined brain regions with improved representation of the subcortex and cerebellum. Neuroimage 206:116290. https://doi.org/10.1016/j.neuroimage.2019.116290 \u003c/li\u003e\n\u003cli\u003eThomas Yeo BT, Krienen FM, Sepulcre J, et al (2011) The organization of the human cerebral cortex estimated by intrinsic functional connectivity. J Neurophysiol 106:1125\u0026ndash;1165. https://doi.org/10.1152/jn.00338.2011 \u003c/li\u003e\n\u003cli\u003eKoo TK, Li MY (2016) A guideline of selecting and reporting intraclass correlation coefficients for reliability research. Journal of chiropractic medicine 15:155\u0026ndash;163 \u003c/li\u003e\n\u003cli\u003ePomponio R, Erus G, Habes M, et al (2020) Harmonization of large MRI datasets for the analysis of brain imaging patterns throughout the lifespan. NeuroImage 208:116450 \u003c/li\u003e\n\u003cli\u003eYu M, Linn KA, Cook PA, et al (2018) Statistical harmonization corrects site effects in functional connectivity measurements from multi-site fMRI data. Human brain mapping 39:4213\u0026ndash;4227 \u003c/li\u003e\n\u003cli\u003eSherman LE, Rudie JD, Pfeifer JH, et al (2014) Development of the default mode and central executive networks across early adolescence: a longitudinal study. Dev Cogn Neurosci 10:148\u0026ndash;159. https://doi.org/10.1016/j.dcn.2014.08.002 \u003c/li\u003e\n\u003cli\u003eFair DA, Cohen AL, Dosenbach NUF, et al (2008) The maturing architecture of the brain\u0026rsquo;s default network. Proceedings of the National Academy of Sciences 105:4028\u0026ndash;4032. https://doi.org/10.1073/pnas.0800376105 \u003c/li\u003e\n\u003cli\u003eSanders AFP, Harms MP, Kandala S, et al (2023) Age-related differences in resting-state functional connectivity from childhood to adolescence. Cereb Cortex 33:6928\u0026ndash;6942. https://doi.org/10.1093/cercor/bhad011 \u003c/li\u003e\n\u003cli\u003eSatterthwaite TD, Wolf DH, Ruparel K, et al (2013) Heterogeneous impact of motion on fundamental patterns of developmental changes in functional connectivity during youth. Neuroimage 83:45\u0026ndash;57. https://doi.org/10.1016/j.neuroimage.2013.06.045 \u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"brain-informatics","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"brai","sideBox":"Learn more about [Brain Informatics](https://braininformatics.springeropen.com/about)","snPcode":"40708","submissionUrl":"https://submission.nature.com/new-submission/40708/3","title":"Brain Informatics","twitterHandle":"@SpringerOpen","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"BMC/SO AJ","inReviewEnabled":true,"inReviewRevisionsEnabled":true},"keywords":"multi-site neuroimaging, resting-state fMRI, functional connectivity, ComBat harmonization, empirical bayes, mild traumatic brain injury","lastPublishedDoi":"10.21203/rs.3.rs-9418750/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-9418750/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eComBat is commonly used to mitigate batch effects in multi-site neuroimaging studies. Although empirical bayes estimation is a key component of ComBat, its impact on resting-state functional connectivity measures, including site effect reductions and preservation of original within-site variability remains underexplored. Empirical bayes estimation was evaluated on functional connectivity in adolescents with mild TBI. Resting-state fMRI data were collected from 144 adolescents across six sites participating in the CARE4Kids study. Functional connectivity was computed over resting state networks defined using Seitzman parcellation. ComBat with and without empirical bayes estimation was assessed by quantifying: (1) site variability based on one-way ANOVA, (2) within-site consistency before and after harmonization (intraclass correlation coefficient), (3) age correlations, and (4) site effects of the principal components. All networks demonstrated significant moderate-to-large site effects (all \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.001) before harmonization. Harmonization without empirical bayes yielded low residual site effect sizes, whereas empirical bayes harmonization yielded low-to-medium site effects. Within-site consistency between the data before and after harmonization was consistently excellent when using empirical bayes estimation. Harmonization without empirical bayes was excellent across sites with higher sample sizes (n\u0026thinsp;\u0026ge;\u0026thinsp;13), but moderate-to-excellent for smaller sites (n\u0026thinsp;\u0026le;\u0026thinsp;12). Both approaches improved age associations across connections involving the default mode and dorsal attention networks and reduced site effects in the first principal component (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026gt;\u0026thinsp;0.05). ComBat with or without empirical bayes effectively removes site variability while preserving age in functional connectivity measures, however, empirical bayes estimation better preserves the original within-site variability.\u003c/p\u003e","manuscriptTitle":"ComBat Harmonization With and Without Empirical Bayes Estimation for Resting-State Functional Connectivity in Pediatric Mild Traumatic Brain Injury: a CARE4Kids Study","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2026-05-06 16:36:26","doi":"10.21203/rs.3.rs-9418750/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"reviewersInvited","content":"","date":"2026-04-27T09:18:57+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2026-04-22T08:35:29+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2026-04-22T08:34:59+00:00","index":"","fulltext":""},{"type":"submitted","content":"Brain Informatics","date":"2026-04-14T18:05:44+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"brain-informatics","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"brai","sideBox":"Learn more about [Brain Informatics](https://braininformatics.springeropen.com/about)","snPcode":"40708","submissionUrl":"https://submission.nature.com/new-submission/40708/3","title":"Brain Informatics","twitterHandle":"@SpringerOpen","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"BMC/SO AJ","inReviewEnabled":true,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"d0bf5850-2521-49ea-aeeb-d54924681dc2","owner":[],"postedDate":"May 6th, 2026","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"under-review","subjectAreas":[],"tags":[],"updatedAt":"2026-05-06T16:36:26+00:00","versionOfRecord":[],"versionCreatedAt":"2026-05-06 16:36:26","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-9418750","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-9418750","identity":"rs-9418750","version":["v1"]},"buildId":"XKTyCvWXoU3ODBz1xrDgd","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}