Increased power by harmonizing structural MRI site differences with the ComBat batch adjustment method in ENIGMA

1After the early reporting of ventricular enlargement in patients with schizophrenia (SCZ) using pneumoencephalography, there has been an exponential increase in the number of studies that use imaging techniques to detect brain differences in people with psychiatric disorders. This increase is most evident for studies using magnetic resonance imaging (MRI), probably due to its high resolution and its wide availability around the globe. However, most MRI studies have examined relatively small sample sizes, a limitation that may prevent the detection of true differences (type II errors), and because of the use of liberal thresholds, may even lead to increased detection of false differences (type I errors). Consequently, reports of unreliable, inconsistent and even contradictory results are not uncommon.

2Collaborative multi-site initiatives provide an opportunity to assemble larger and more diverse groups of subjects, leading to increased power and findings that may be more representative of the general population. Among these initiatives, the ENIGMA (Enhancing Neuro Imaging Genetics through Meta-Analysis; http://enigma.ini.usc.edu) Consortium stands out for including hundreds of groups worldwide and facilitating the sharing of tens of thousands of neuroimages. One great advantage of this consortium is the harmonization of the protocols to pre-process the MRI data, which has decreased the heterogeneity between the sites related to methodological factors. All sites apply the same pre-processing pipelines to obtain thickness and surface area estimates for cortical regions of interest (ROI) and volume estimates for subcortical ROIs; similar harmonized protocols are in use for standardized analysis of diffusion MRI, resting state fMRI and EEG data, as well as various kinds of omics data (GWAS and epigenetic data).

3However, even though all sites participating in an ENIGMA project apply the same pre-processing protocol, data from different sites still show relevant methodological heterogeneity due to systematic differences in MRI scanning devices and acquisition sequences. Also, prior studies have reported that the results of the FreeSurfer segmentation process, for morphometric analysis of MRI, can be affected even by using different FreeSurfer versions, workstations or operating systems. Most ENIGMA projects address this residual heterogeneity by random-effects meta-analysis (RE-Meta), but estimation and control of heterogeneity in site-aggregated meta-analyses may be suboptimal. It is worth noting that a few ENIGMA studies have analyzed shared individual data (rather than site-aggregated statistical data). These “mega-analyses” of individual data considered the “site” as a random factor within a linear mixed-effects model (ME-Mega), and in several cases examined so far, showed higher statistical power than RE-Meta (Table 1).

4Here, we tested whether ME-Mega may be further improved using a recently developed method to control for batch effects. Standard ME-Mega assumes that the error terms follow the same normal distribution at all sites, which is rarely the case as sites usually have different error variances. In addition, both RE-Mega and ME-Meta estimate the heterogeneity of each ROI independently, while it is likely that all ROIs share some heterogeneity. One method that overcomes these issues is ComBat, a batch adjustment method developed for genomics data. Fortin and colleagues have shown that ComBat mega-analysis (ComBat-Mega) outperformed other methods for removing the effects of site from cortical thickness data obtained using the ANTs cortical thickness pipeline from a moderately small number of different sites (≤7 sites). Specifically, ComBat decreased scan-related heterogeneity and increased statistical power and reproducibility. The current study examines whether this harmonization result can be extended to ENIGMA data obtained using a standardized FreeSurfer pipeline. Moreover, we did not know whether the use of a larger number of sites could minimize the advantages of ComBat-Mega as compared to ME-Mega. To answer these questions, we analyzed the main structural MRI data from the ENIGMA Schizophrenia Working Group using RE-Meta, ME-Mega and ComBat-Mega, and then compared the findings. The RE-Meta of these data have been already published; in those analyses, individuals with SCZ showed widespread thinner cortex and smaller surface area, as well as smaller hippocampus, amygdala, thalamus and accumbens volumes, and larger pallidum and lateral ventricle volumes.

5We hypothesized that ComBat-Mega would show improvements over RE-Meta and ME-Mega in detecting differences between groups of individuals with SCZ and healthy controls (CON), with standard errors of these effects scaling by method: Combat-Mega < ME-Mega < RE-Meta. We further provide the R code (http://enigma.ini.usc.edu/wp-content/uploads/combat_for_ENIGMA_sMRI/combat_for_ENIGMA_sMRI.R) for the application of ComBat harmonization for other ENIGMA mega-analyses or other multi-site structural imaging work even if the imaging data are partially missing in some ROIs (the original ComBat function did not accept missing data).

6Before detailing the collection of data and analyses in the present study, we will briefly explain the three methodological approaches. To exemplify the explanation, we will refer to a simple comparison of cortical thickness between groups of individuals with SCZ and CON, after covarying for effects of age and sex, but the concepts are applicable to other measures and statistical contrasts. We also conducted an alternative analysis covarying for age, sex, and intracranial volume (ICV).

7In the random-effects meta-analysis (RE-Meta), a linear model estimates the difference in cortical thickness between SCZ and CON for each ROI at each site, covarying for age and sex: where yr,i,j is the measurement of cortical thickness of the rth ROI from the jth individual of the ith site, αr,i is the estimate overall cortical thickness of the rth ROI from individuals of the ith site, Xi,j are the values of the variables (disorder, age, and sex) of the jth individual of the ith site, βr,i are the estimates of the coefficients of these variables for the rth ROI from individuals of the ith site, and εr,i,j is the error term for the rth ROI in the jth individual of the ith site.

8Estimates of coefficients of interest (e.g., βr,i,1, the difference between SCZ and CON are then pooled to obtain a single estimate for each ROI (βr,meta,1). A typical method to pool the coefficients is the weighted mean of the coefficient of each site: where wr,i the weight of ith site for the rth ROI, and is calculated as the inverse of the variance of βr,i,1, plus the heterogeneity for the rth ROI:

9Frequently, the analyst does not use the coefficients but effect sizes, such as Hedges’ g, but the concept is similar.

10Some problems of RE-Meta are that βr,i may be poorly estimated in sites with small sample sizes, or that may be poorly estimated in some scenarios.

11In the standard mixed effects mega-analysis (ME-Mega), a linear mixed-effects model is performed on shared individual subject data to estimate the overall difference in cortical thickness between SCZ and CON, for each ROI, covarying for age and sex. This analysis is conducted in a single step, with “site” included in the model as a random factor: where αr is the estimate overall cortical thickness of the rth ROI from all individuals, βr are the estimates of the coefficients of the variables for the rth ROI from all individuals, and γr,i are the additive effects of the ith site in the rth ROI.

12This approach benefits from a more robust estimation of αr and βr as it is based on the data from all sites, as well as from a more precise estimation of the heterogeneity. However, it still may have some minor issues. It assumes that the error terms follow the same normal distribution at all sites, which may seldom be the case. We acknowledge that it is possible to create linear mixed-effects models that consider a different variance for each site, but they involve the specification of variance structures for each statistical test, which may substantially complicate the analyses. In addition, the effects of site are estimated independently for each ROI, which may be suboptimal because the effects of site, even if different for each ROI, may still share some commonalities (e.g., an MRI device may yield a better signal contrast than another across the brain).

13As compared to ME-Mega, the Combat mega-analysis (ComBat-Mega) assumes that the error terms may follow varying normal distributions at different sites: where δr,i are the multiplicative effects of the ith site in the rth ROI.

14In addition, it assumes that the additive and multiplicative effects of the sites are not completely independent across ROIs but, rather, they share a common distribution. Such considerations prevent the use of standard linear models, but ComBat uses an empirical Bayes framework to estimate the distribution of the effects of site. Once estimated, it derives the additive error terms:

15These terms allow the derivation of harmonized data:

16These simpler data can then be analyzed with standard linear models to estimate the overall difference in cortical thickness between SCZ and CON groups, for each ROI.

17Fortin and colleagues modified the original “combat” function, in the “sva” package for R, so that it could be applied to imaging data. However, Fortin’s “combat” function may not be easily applicable to ENIGMA projects as it requires that the dataset has no missing data, which is seldom the case. In addition, it finds the harmonization parameters and applies them to the data within the same function, while some analyses - such as machine learning - require that the parameters are found in a training set and later applied to an independent test set (this is not the case here, but it might be the case in future studies). We further modified the “combat” function to allow for missing data and to separate the fitting and the application of the harmonization.

18First, we divided the function into two subfunctions: “combat_fit”, which finds the harmonization parameters, and “combat_apply”, which applies them to the same or to another set. The “combat_fit” function automatically imputes missing data so that the function can find the harmonization parameters without errors. These imputations are predictions based on linear models of the ROI values by the covariates, separately for each ROI and each site:

19The covariates are the variables introduced into the “combat_fit” function, which in the present study were the diagnosis, age, and sex. The “combat_fit” function also discards ROIs with no variance, which returned errors in the previous “combat” function. Importantly, these imputations are temporary and only aimed to avoid errors during the fitting of the parameters, they are not saved. To apply the parameters, the user must use the “combat_apply” function with the original data, and missing values are not imputed.

20The reader may download the adapted ComBat functions for R from http://enigma.ini.usc.edu/wp-content/uploads/combat_for_ENIGMA_sMRI/combat_for_ENIGMA_sMRI.R.

21The data for this paper includes the cortical thickness, surface area and subcortical volumes from 33 sites of the ENIGMA Schizophrenia Working Group who shared individual subject level FreeSurfer data for this project. The overall sample included 2897 individuals with a diagnosis of SCZ (mean age 34 years, 34% females) and 3141 CON (mean age 33 years, 49% females). For SCZ, the mean age of onset was 23 years and their Positive and Negative Syndrome Scale (PANSS) scores for total/positive/negative symptoms were 61/16/17, respectively. The researchers at each of the sites had collected the data after obtaining participants’ written informed consent, with protocols that had been approved by local institutional review boards. We provide a description of the overall sample in Table 2 and a description of the sample from each site in Supplementary Table S1.

22All sites had processed the data with FreeSurfer versions 4.0 to 5.3, except for version 5.2 which was found to produce low intra-class correlations compared to the other versions, and within site all patients and controls were processed using the same FreeSurfer version according to the ENIGMA protocols, which are available at http://enigma.usc.edu/protocols/imaging-protocols. For cortical ROIs, they involved the estimation of cortical vertex-wise statistics, the extraction of cortical thickness and surface area for 70 Desikan-Killiany (DK) atlas regions, and quality checks. For subcortical ROIs, they involved the estimation of subcortical volumes and quality checks.

23We conducted comparisons of MRI data between individuals with SCZ and CON to assess the statistical significance, power and familywise error rate (FWER) using RE-Meta, ME-Mega and ComBat-Mega. We formally tested whether ComBat-Mega increases the statistical significance and power of the differences between individuals with SCZ and CON by attenuating site-effects, using a permutation test and a small-subset strategy respectively. We also used the data of the permutation test to check the FWER.

24We conducted the RE-Meta in two steps. In the first step, we compared the values of each ROI between SCZ and CON via a standard linear model, with age and sex as covariates, separately for each site. We then converted the difference to a Hedges’ g and its variance for each site and ROI. In the second step, we conducted a random-effects meta-analysis of the Hedges’ g of each ROI with the “metafor” package for R, and we corrected the p-values for multiple comparisons with the Holm method.

25For ME-Mega, we compared the values of each ROI between SCZ and CON via a linear mixed-effects model, with age and sex as covariates and site as a random factor, with the “lme4” and “lmerTest” packages for R. We then divided the difference by the standard deviation (derived from the model) and corrected it for small-sample bias to obtain a Hedges’ g and its variance, and we corrected the p-values for multiple comparisons using the Holm method.

26Finally, for ComBat-Mega, we first removed the effects of site using the ComBat functions (modelling the effects of diagnosis, age, and sex), and then compared the values of each ROI (e.g., cortical thickness of the frontal pole) between SCZ and CON via a standard linear model, with age and sex as covariates. Note that the ComBat functions use covariates (e.g., age and sex) to better estimate the effects of site, but they do not remove the effects of these covariates; for this reason, we included these covariates in the subsequent linear model. As for ME-Mega, we converted the difference to a Hedges’ g and its variance, and we corrected the p-values for multiple comparisons with the Holm method. Note that we applied a single ComBat harmonization for different types of data (cortical thickness, cortical surface area, and subcortical volume) because we considered that they were related. We also conducted an alternative analysis with a separate harmonization for each type of data.

27To test whether ComBat-Mega had improved the statistical significance we used a permutation approach. We followed the Draper-Stoneman procedure, which according to results from a study comparing different algorithms, is one of the procedures that best controls the FWER and that can be safely applied here. Note that other algorithms such as Freedman Lane would produce different permuted data for RE-Meta, ME-Mega and ComBat-Mega, which would be problematic in our study because these unwanted differences could confound other potential differences between the methods. Specifically, we randomly permuted the diagnosis among the individuals within each site and repeated all comparison analysis 1000 times.

28To show the differences in statistical significance between methods expected by chance, we plotted the histogram of the median difference in the logit-transformed p-values between the methods across the permutations (Fig. 1). For example, in one permutation we randomly assigned study participants to patient or control status. We then compared these randomly assigned patients and controls using RE-Meta, ME-Mega and ComBat-Mega. We then calculated differences between logit-transformed p-values of the ComBat-Mega comparison and logit-transformed p-values of the RE-Meta (or ME-Mega) comparisons for each ROI. From these, we only saved the median between logit-transformed p-value difference. Note that this median difference should be very close to zero, given that participant assignment was random, and there should therefore be no patient-control group differences other than by chance. By conducting multiple of these permutations, we were able plot the histogram of the median differences expected by chance alone. Finally, we compared the median difference of the original analysis (with correctly assigned patient and control status) with the histogram of the median differences expected by chance. Only median differences were used in this analysis to simplify the test as doing so avoids the need to correct for multiple comparisons.

29We must note that without the logit (or other) transforms, the detection of differences in statistical significance would be too sensitive for large p-values and too little sensitive for small p-values. For example, if the (non-transformed) p-value using one approach was 0.6 and the (non-transformed) p-value using another approach was 0.4, the difference in p-values would be very large (0.6–0.4 = 0.2) even if the two p-values might be considered conceptually very similar, whereas if the (non-transformed) p-value using one approach was 0.003 and the (non-transformed) p-value using another approach was 0.001, the difference in p-values would be very small (0.003–0.001 = 0.002) even if one p-value is three times the size of the other. With the logit transform, the p-values of the first example would be 0.4 and −0.4, with a difference of 0.8, and the p-values of the second example would be −5.8 and −6.9, with a difference of 1.1.

30The use of a permutation test implied that both the estimated probability of obtaining the observed median difference in (logit-transform) p-values was discrete, i.e., it could only be 0.001, or 0.002, or 0.003, etcetera. However, we were only interested in assessing if this estimation was <0.05, for what this level of precision should not pose any problems.

31We also tested whether ComBat-Mega increases the statistical power using a small-subset strategy. Specifically, we repeated 500 times the analyses but including each time only a random sample of 10 sites. We then counted the number of times that these analyses using only 10 sites were able to detect differences between SCZ and CON. We only used ROIs in which the differences between SCZ and CON were strongly statistically significant in the main analyses using the 33 sites (FWER<0.001 for RE-Meta, for ME-Mega, and for ComBat-Mega), as we assumed that they have true differences. Finally, we conducted a Wilcoxon signed-ranked test to compare the statistical power across ROIs between ComBat-Mega and RE-Meta, as well as between ComBat-Mega and ME-Mega.

32We also used the permutation data created above to check whether the FWER for the three methods were appropriate, i.e., we counted the proportion of permutations in which at least one ROI had a Holm-corrected p-value < 0.05. Again, the use of a permutation test implied that the estimated FWER was discrete, but we were only interested in assessing whether it was <0.05.

33With ComBat-Mega, on average, individuals with a diagnosis of SCZ showed thinner cortex and smaller surface area in nearly all cortical ROIs (Table 3). The only exceptions were the bilateral pericalcarine fissures and right entorhinal cortex (where between-group differences in thickness did not reach statistical significance after correction for multiple comparisons) and the left isthmus of the cingulate and right temporal pole (where between-group differences in surface area did not reach statistical significance after correction for multiple comparisons). The SCZ group also showed, on average, smaller bilateral thalamus, hippocampus, amygdala, and right accumbens volumes, and larger bilateral lateral ventricle, putamen, and pallidum volumes. Smaller left accumbens and larger bilateral caudate volumes were not statistically significant after correction for multiple comparisons.

34Results were in the same direction for the RE-Meta and ME-Mega, though RE-Meta did not detect thinner cortex in three ROIs (bilateral rostral anterior cingulate and left caudal anterior cingulate) and smaller surface area in six ROIs (bilateral pericalcarine fissure, left posterior cingulate and temporal pole, and right isthmus cingulate and insula).

35The Hedges’ g estimates for the differences were similar across the different analytic methods, but their statistical significance was greater in ComBat-Mega as compared to RE-Meta and ME-Mega (Figs. 2 and 3). The difference in statistical significance was relatively minor when comparing ComBat-Mega to ME-Mega, whereas particularly relevant when comparing ComBat-Mega to RE-Meta (Fig. 3).

36The median difference between logit-transformed ComBat-Mega p-values and logit-transformed RE-Meta p-values in the original data was 13.9. This was substantially larger than any of the median differences in the permuted data (all < 0.61), indicating that the higher statistical significance of ComBat-Mega findings was unlikely due to chance (probability 0.001) (Fig. 4). For the comparison between ComBat-Mega and ME-Mega, the median difference was smaller (3.2), but still unlikely due to chance (all median differences in the permuted data <0.52, probability 0.001).

37Interestingly, a plot of the ComBat-Mega-related increase in statistical significance as a function of the intra-site variance/total variance ratio, showed that the increase in statistical significance was larger in those ROIs in which intra-site variance was only ~50–70% of total variance compared to those ROIs in which intra-site variance was ~90–100% of total variance (p < 0.001, Fig. 5).

38In the evaluation of statistical power using the small-subset strategy, the statistical power was higher for ComBat-Mega (statistical power = 83.5%) than for RE-Meta (statistical power = 53.7%; Wilcoxon p-value < 0.001) or ME-Mega (statistical power = 80.4%; Wilcoxon p-value < 0.001).

39The empirical FWER was ≤0.05 for all analytic methods (RE-Meta: 0.024; ME-Mega: 0.027; ComBat-Mega: 0.025).

40When we applied the ComBat harmonization separately for cortical thickness data, cortical surface area data and subcortical volume data, we found the same differences with nearly identical Hedges’ g (Supplementary Figure S1). The statistical significance was minimally lower (median difference between single ComBat logit-transformed p-values and separate ComBat logit-transformed p-values was 0.1), the statistical power in the small-subset strategy was 83.5%, and the empirical FWER was 0.026.

41When we covaried ComBat-Mega by age, sex and ICV, results were similar: The only differences were that the right frontal pole, isthmus of the cingulate and pericalcarine and left parahippocampal and temporal pole decreases in surface area were no longer statistically significant, whereas the left pericalcarine decrease in surface area and the bilateral caudate increases in volume reached statistical significance. Results were again in the same direction for the RE-Meta and ME-Mega, though RE-Meta did not detect statistically significant differences in 36 of the ROIs showing differences with ComBat-Mega, and ME-Mega did not detect smaller right accumbens volume (and detected smaller surface area in left parahippocampal and right pericalcarine but not in left paracentral and right entorhinal). The Hedges’ g estimates for the differences were again similar across the different analytic methods, but their statistical significance was again greater in ComBat-Mega as compared to RE-Meta and ME-Mega (Supplementary Figure S2).

42In this study, we analyzed individual subject level data pooled by the ENIGMA Schizophrenia Working Group using three methods to account for the effects of site: random-effects meta-analysis (RE-Meta), linear mixed-effects models (ME-Mega), and ComBat harmonization followed by standard linear models (ComBat-Mega). The results of the comparison between SCZ and CON using ComBat-Mega were similar to the studies already published by the ENIGMA Schizophrenia Working Group: SCZ showed a widespread thinner cortex and smaller surface area, smaller hippocampus, amygdala, thalamus and accumbens, and larger lateral ventricles, putamen and pallidum than CON. The results of the same comparison using RE-Meta and ME-Mega were in the same direction and had similar effect sizes, although with a lower statistical significance (i.e. wider confidence intervals, larger p-values), especially for RE-Meta. In other words, the use of ComBat increased the statistical significance (i.e., narrower confidence intervals, smaller p-values) of the differences between SCZ and CON. This was specially apparent in those ROIs in which intra-site variance was only ~50–70% of total variance. ComBat Mega also showed increased statistical power when we repeated the analyses with fewer sites. All approaches controlled well the FWER, even too strictly probably due to the use of the Holm method, which is more powerful than the Bonferroni method but still conservative. Findings were similar when covarying by ICV.

43Based on these findings, we recommend that ENIGMA mega-analysis projects consider applying the ComBat function to reduce the effects of site, followed by standard statistical analysis without including site as a fixed or random effect in the statistical model. To apply ComBat harmonization, we provide easy-to-use functions for R that work even if there are missing data and they can be trained with data from one set and then applied to data from another.

44We must note that we conducted these analyses with the three main types of data used in ENIGMA projects: thickness of cortical ROIs, surface area of cortical ROIs, and volumes of subcortical nuclei. However, some ENIGMA projects use other types of data, such as mean fractional anisotropy of white matter tracts, and we have not explored whether the application of ComBat would be beneficial for these projects. Two notions suggest that ComBat should be broadly beneficial. On the one hand, the ComBat algorithm is not specific for a given type of imaging data. Indeed, while it was developed for genomics data, we here successfully applied it to three types of ENIGMA imaging data. Moreover, Fortin and colleagues found that ComBat outperforms other harmonization methods for voxel-based fractional anisotropy and mean diffusivity, and Yu et al. found similar results for resting-state functional connectivity and network measures.

45While our findings suggest that ComBat harmonization will be useful for most ENIGMA mega-analyses and other multi-site structural imaging work, we suggest caution when combining different types of data. We conducted a single ComBat harmonization for different types of MRI data because we considered that thickness, area, and volume are related, as they are obtained from the same FreeSurfer output of the T1-weighted image and all measure amounts of gray matter. Indeed, an alternative analysis with separate ComBat harmonization for each type of data yielded nearly identical results. However, we do not know whether the application of a single ComBat harmonization on other combinations of data would behave similarly.

46Other popular approaches for pooling neuroimaging data are the voxel-based meta-analytic methods, such as Seed-based d Mapping (SDM) or Activation Likelihood Estimation (ALE). These methods can include imaging studies even if they only report the coordinates of the peaks of the clusters of statistical significance. Therefore, a great advantage of these methods is the exhaustive inclusion of studies. In addition, the analyses are conducted at the voxel level (rather than using ROIs). These methods traditionally tested whether the reported findings tended to converge in a few brain voxels, but novel methods are able to directly test whether there are differences - even if they are widespread and do not converge. In view of the results of the present study, one could wonder whether these voxel-based methods should also conduct ComBat mega-analysis instead of meta-analysis. However, to use ComBat they would need access to individual subject level data, which at present are often not available. Another aspect to consider is whether we need SDM or ALE meta-analyses after an ENIGMA ComBat mega-analysis is published. Here, we must remember that SDM and ALE are voxel-based and include virtually all published studies, whereas most ENIGMA studies are ROI-based and include only the data that authors agree to share. Therefore, these different approaches present interesting complementary information.

47Our study has some limitations. First, we already stated that we have not explored whether the application of ComBat would be beneficial for projects using other types of data, although several facts suggest that ComBat should be broadly beneficial. Second, we also acknowledged that we do not know whether the application of a single ComBat harmonization on other combinations of data would behave similarly. Third, our analysis is focused on the differences between SCZ and CON, whose distribution is roughly similar across sites. The effects of site and thus the importance of their removal might be larger for conditions with few cases in each site, where pooling data is more beneficial. Fourth, ComBat-Mega addresses some issues but not others, which still need to be investigated, such as site by nuisance confounds. For example, a site with poor quality data may also be a site with a mean age higher than other sites. Future studies addressing these questions could point to methods other than ComBat. Finally, there is a conceptual difference in the effects of site that are modeled in ComBat/ME-Mega and the effects of site that are modeled in RE-Meta. The former effects are in (individual) raw data and refer to site-specific constants that are added to or that multiply the measurements. The latter effects, conversely, are in (group) effect sizes, and are probably a mix of several factors such as site-specific constants that multiply the measurements, heterogeneity in the differences between SCZ and CON, or differences in precision between studies.

48To conclude, this paper provides evidence of the superiority of ComBat harmonization over standard mega-analyses and meta-analyses in reducing site-related heterogeneity and thus increase statistical power. We therefore recommend that ENIGMA mega-analysis projects and other multi-site structural imaging work consider applying the ComBat function, which we provide employing easy functions for R. The provided code works with missing data and allows for harmonization of a test set based on the training set (a requirement for machine learning and possibly replication studies). We hope that future ENIGMA mega-analysis projects will improve between-site harmonization using ComBat.