Mitigating site effects in covariance for machine learning in neuroimaging data

1The need for larger samples in human subjects research have led to a growing number of multi‐site studies that aggregate data across multiple locations. This trend is especially prevalent in neuroimaging research where the reliability and generalizabilty of findings from the conventional single‐site studies are often limited by the ability to recruit and study sufficiently large and representative samples from the population. Many consortia have been formed to address such issues (Mueller et al.,; Sudlow et al.,; Trivedi et al.,; Van Essen et al., ). The larger samples obtained through these efforts promote greater power to detect significant associations as well as better generalizability of results. However, these study designs also introduce heterogeneity in acquisition and processing that, if not appropriately addressed, may impact study findings.

2Several researchers have determined that variability driven by site differences, often called site effects, reduce the reliability of derived measurements, and can introduce bias. Neuroimaging measurements have been repeatedly shown to be affected by scanner manufacturer, model, magnetic field strength, head coil, voxel size, and acquisition parameters (Han et al.,; Kruggel, Jessica, Tugan Muftuler, & Initiative,; Reig et al.,; Wonderlick et al., ). Yet even in scanners of the exact same model and manufacturer, differences still exist for certain neuroimaging biomarkers (Takao, Hayashi, & Ohtomo, ).

3Until recently, neuroimaging analyses primarily involved mass univariate testing which treats features as independent and does not leverage covariance between features. Under this paradigm, the impact of site effects is through changes in the mean and variance of measurements. Increasingly, researchers have used sets of neuroimaging features as inputs into prediction algorithms or state‐of‐the‐art machine learning (ML) methods. This approach has become a powerful tool for leveraging both functional and structural neuroimaging for research into pain (Smith, López‐Solà, McMahon, Pedler, & Sterling,; Wager et al., ), neural representations (Haxby, Connolly, & Swaroop Guntupalli, ), and psychiatric illnesses (Koutsouleris et al., ). One of the major benefits of ML is that it leverages the joint distribution and correlation structure among multivariate brain features in order to better characterize a phenotype of interest (Gregorutti, Michel, & Saint‐Pierre,; O'Toole et al., ). As a result, site effects on the covariance of measurements are likely to impact findings substantially. In fact, a recent investigation showed that an ML algorithm was able to detect scanner with high accuracy and that the detection of sex depended heavily on the scanners included in the training and test data (Glocker, Robinson, Castro, Dou, & Konukoglu, ).

4The major statistical harmonization techniques employed in neuroimaging have generally corrected for differences across sites in mean and variance, but not covariance (Fortin et al.,; Fortin, Sweeney, Muschelli, Crainiceanu, & Shinohara,; Rao, Monteiro, & Mourao‐Miranda,; Yamashita et al., ). Increasingly, the ComBat model (Johnson, Cheng, & Rabinovic, ) has become a popular harmonization technique in neuroimaging and has been successfully applied to structural and functional measures (Bartlett et al.,; Fortin et al.,,; Marek et al.,; Yu et al., ). However, this model does not address potential site effects in covariance.

5Recently, another stream of data‐driven harmonization methods have aimed to apply generative adversarial networks (GANs) or distance‐based methods to unify distributions of measurements across sites. However, the GAN‐based harmonization methods have only been tested for harmonization of images and lack options to retain clinical associations of interest (Gao, Liu, Wang, Shi, & Jinhua,; Nguyen, Morris, Harris, Korgoankar, & Ramos,; Zhong et al., ). A recent distance‐based method is applicable to derived measurements and has been tested in classification of Alzheimer's disease (AD) using support vector machines (Zhou et al., ). However, the method in (Zhou et al., ) has not been tested for detection of site via ML and also requires several conditions which may not hold in studies with sufficiently heterogeneous sites or major differences in subject demographics across sites.

6In this article, we examine whether site effects influence ML results. In particular, we study the cortical thickness measurements derived from images acquired by the Alzheimer's Disease Neuroimaging Initiative (ADNI) and demonstrate the existence of site effects in covariance of structural imaging measures. We then propose a novel harmonization method called Correcting Covariance Batch Effects (CovBat) that removes site effects in mean, variance, and covariance. We apply CovBat and show that within‐site correlation matrices are successfully harmonized. Furthermore, we find that ML methods are unable to detect Siemens scanners after our proposed harmonization is applied, and that the CovBat‐harmonized data retain accurate prediction of disease group. We also assess the performance of the proposed method in simulated data, and again find that the method mitigates site effects and maintains detection of meaningful associations. Our results demonstrate the need to consider covariance in harmonization methods, and suggest a novel procedure that can be applied to better harmonize data from multi‐site imaging studies.

7The data for this article consist of baseline scans from ADNI (http://adni.loni.usc.edu/ which are processed using the ANTs longitudinal single‐subject template pipeline (Tustison et al., ) with code available on GitHub (https://github.com/ntustison/CrossLong). Informed consent was obtained by all participants in the ADNI study. Institutional review boards approved the study at all of the contributing institutions.

8We briefly summarize the steps involved. First, we download raw T1‐weighted images from the ADNI‐1 database, which were acquired using MPRAGE for Siemens and Philips scanners and a works‐in‐progress version of MPRAGE on GE scanners (Jack et al., ). We choose ADNI‐1 to highlight a severe site effect situation driven by greater variability in scanner properties, including magnetic field strength, which are standardized in later ADNI releases. For each subject, we estimate a single‐subject template from all the image timepoints. After rigid spatial normalization to this single‐subject template, each normalized timepoint image is then processed using the single image cortical thickness pipeline consisting of brain extraction (B. Avants, Klein, Tustison, Woo, & Gee, ), denoising (Manjón, Coupé, Luis, Louis Collins, & Robles, ), N4 bias correction (N. J. Tustison et al., ), Atropos n‐tissue segmentation (B. B. Avants, Tustison, Wu, Cook, & Gee, ), and registration‐based cortical thickness estimation (Das, Avants, Grossman, & Gee, ). For our analyses, we use the 62 baseline cortical thickness values as defined by the Desikan–Killiany–Tourville atlas (Klein & Tourville, ). The sample covariance matrix for these cortical thicknesses in the largest site is shown with labels in Figure S1.

9We define site based on information contained within the Digital Imaging and Communications in Medicine (DICOM) headers for each scan. Specifically, subjects are considered to be acquired at the same site if they share the scanner location, scanner manufacturer, scanner model, head coil, and magnetic field strength. In total, this definition yields 142 distinct sites of which 78 had less than three subjects and were removed from analyses. The final sample consists of 505 subjects across 64 sites, with 213 subjects imaged on scanners manufactured by Siemens, 70 by Philips, and 222 by GE. The sample has a mean age of 75.3 (SD 6.70) and is comprised of 278 (55%) males, 115 (22.8%) AD patients, 239 (47.3%) late mild cognitive impairment (LMCI), and 151 (29.9%) cognitively normal (CN) individuals.

10The ADNI sample demographics are considerably different across sites, which precludes application of certain harmonization methods. For example, (Zhou et al., ) relies on “nontrivial overlap” of the potential confounders across sites and proposed a subsampling approach that performs distributional shifts on subsamples of data matched by the discrete stratum of the confounders. Given that our data are sufficiently heterogeneous, it is challenging to form bins matched by age, sex, and diagnosis status to ensure that each site has at least one individual in each bin. This prevents protection of age effects in applying the harmonization method proposed by (Zhou et al., ).

11We first review ComBat (Fortin et al.,,; Johnson et al., ) for harmonization of neuroimaging measures. ComBat seeks to remove the mean and variance site effects of the data in an empirical Bayes framework. Let yij=yij1yij2…yijpT, i=1,2,…,M, j=1,2,…,ni denote the p×1 vectors of observed data where i indexes site, j indexes subjects within sites, ni is the number of subjects acquired on site i, and p is the number of features. Our goal is to harmonize these vectors across the M sites. ComBat assumes that the features indexed by v followwhere αv is the intercept, xij is the vector of covariates, βv is the vector of regression coefficients, γiv is the mean site effect, and δiv is the variance site effect. The errors eijv are assumed to independently follow eijv∼N0σv2. ComBat first finds least‐squares estimates α^v and β^v for each feature. To estimate the site effects using empirical Bayes, ComBat assumes that the γiv follow independent normal distributions and the δiv follow independent inverse gamma distributions. The hyperparameters are then estimated via method of moments using data across all features. The empirical Bayes point estimates γiv* and δiv* are then obtained as the means of the posterior distributions. Finally, ComBat residualizes with respect to these estimates to obtain harmonized data

12To address potential covariance site effects, we build on the existing ComBat framework. We again assume that the features follow

13However, the error vectors eij=eij1eij2…eijpT∼N0∑i may be spatially correlated and differ in covariance across site. Analogous to how ComBat modifies observations to bring each within‐site variance to the pooled variance across sites, our proposed method modifies principal component (PC) scores to shift each within‐site covariance to the pooled covariance structure. We achieve this by approximating within‐site covariance structures using the PCs and PC scores obtained across all observations. We propose the CovBat algorithm, which accounts for the joint distribution of ComBat‐adjusted observations as follows:

14Step 1. We first perform ComBat to remove the mean and variance shifts in the marginal distributions of the cortical thickness measures. Then, we additionally residualize with respect to the intercept and covariates to obtain ComBat‐adjusted residuals denoted eijComBat=eij1ComBateij2ComBat…eijpComBatT where p is the number of features. We then define these residuals using notation from Section 2.2 aswhere i=1,2,…,M, j=1,2,…,ni, M is the number of sites, and ni is the number of subjects acquired at site i. α^v, xijT, β^v, γiv*, and δiv* are defined in Section 2.2.

15Step 2. The eijComBat are assumed to have mean 0; their covariance matrices which we denote by ∑i, however, may differ across sites. We first perform principal components analysis (PCA) on the full data residuals and represent the full data covariance matrix as ∑=∑k=1qλkϕkϕkT where the rank q=min∑i=1Mnip, λk are the eigenvalues of ∑, and ϕk are the PCs obtained as the eigenvectors of ∑. In practice, PCA is performed on the sample covariance matrix ∑^ and we obtain estimated eigenvalues λ^k and eigenvectors ϕ^k. The ComBat‐adjusted residuals can then be expressed as eijComBat=∑k=1qξijkϕ^k where ξijk are the principal component scores.

16We then aim to bring each within‐site covariance matrix ∑i to the pooled covariance across sites. Since our goal is to recover covariance structures resembling ∑, we approximate the within‐site covariance matrices as ∑^i=∑k=1qλ^ikϕ^kϕ^kT where λ^ik are within‐site eigenvalues estimated as the sample variance of the principal component scores λ^ik=∑j=1niξijk−∑j=1niξijk/ni2/ni−1 and ϕ^k are estimated from the full data covariance. This model assumes that the covariance site effect is contained within the variances of the PC scores with the principal components estimated from the full data. This assumption may not hold in some cases, but harmonization of these PC score variances will bring the within‐site covariance matrices closer to the pooled covariance. This is analogous to how ComBat brings site‐specific variance closer to the variance estimated using observations across all sites.

17Step 3. Thus, we posit:where ϵijk∼N0τk2, τk is the error standard deviation, and μik, ρik are the center and scale parameters corresponding to principal components k=1,2,…K where K≤q is a tuning parameter chosen to capture the desired proportion of the variation in the observations. If K is chosen such that K=q, all principal components are harmonized. Note that this is similar to the ComBat model, applied to each of the k principal component scores instead of the original measures. We can then estimate each of the K pairs of center and scale parameters by finding the values that bring each site's mean and variance in scores to the pooled mean and variance, which we denote μ^ik) and ρ^ik. We then remove the site effect in the scores via ξijkCovBat=ξijk−μ^ik/ρ^ik.

18Step 4. We obtain CovBat‐adjusted residuals eijCovBat=eij1CovBateij2CovBat…eijpCovBatT by projecting the adjusted scores back into the residual space via,

19We then add the intercepts and covariates effects estimated in Step 1 to obtain CovBat‐adjusted observations

20Let yij, i=1,2,3, j=1,2,…,ni be vectors of length p representing simulated cortical thickness values for three sites, each with ni observations. The yij are generated using the following model:where xij is a simulated diagnosis variable drawn from a Bernoulli distribution with probability 0.25, α is the first p/2 elements in each hemisphere from the sample mean vector of Scanner B observations in the ADNI data, β is the vector of simulated diagnosis effects on the mean, and eij is the vector of error terms. We simulate mean and variance site effects based on the assumptions of ComBat and CovBat. The mean site effects γi=γi1γi2…γipT are vectors drawn from independent and identically distributed (i.i.d.) normal distributions with mean zero and standard deviation 0.1. The variance site effects δi=δi1δi2…δipT are vectors drawn from site‐specific inverse gamma distributions with chosen parameters. For our simulations, we chose to distinguish the site‐specific scaling factors by assuming δ1v∼i.i.dInverse Gamma46,50, δ2v∼i.i.dInverse Gamma51,50, and δ3v∼i.i.dInverse Gamma56,50 for v=1,2,…,p.

21We first assess whether CovBat can recover the underlying covariance structure when the covariance site effects are captured by its PC directions. We refer to this simulation setting as the Simple Covariance Effects simulation. For p=62, we set the underlying covariance S as the sample correlation matrix of cortical thickness observations in the ADNI data, with eigendecomposition S=∑k=162λ^kϕ^kϕ^kT. We then generate error terms eij that contain site‐specific shifts in the first eigenvalue viawhere ci controls the severity of the covariance shift. For our simulations, we set c1=−1/2, c2=0, and c3=1/2 so that the pooled covariance structure is equivalent to S. We choose βi=−0.5 for p/4 regions of interest in both the left and right hemispheres to associate the simulated diagnosis with decreases in mean simulated cortical thicknesses.

22We also investigate how the rank of the covariance effect influences CovBat harmonization results. We modify the rank of the simple covariance effect by varying K in the generation of error termswhere ci takes the same values as previously, c1=−1/2, c2=0, and c3=1/2. We simulate datasets while choosing K as 2, 6, and 12 PCs and evaluate how harmonization influences detection of site via ML.

23To evaluate how CovBat performs when the covariance site effects are not easily captured by the principal components and when the simulated diagnosis may affect covariance, we modify the simulation to incorporate high‐rank covariance shifts due to site via Ωi and due to simulated diagnosis as Ψ. For p=62, the error terms eij are now generated viawhere ∑ij=S+xijΨ+Ωi, S is the sample correlation matrix of cortical thickness observations in the ADNI data, Ψ is a chosen diagnosis‐driven covariance shift matrix, and Ωi are site‐specific covariance shift matrices. The matrices Dij=diagdij where dijk=1/∑ijkk for k=1,2,…,62 ensure that these covariance effects do not modify the marginal variances of eij. To constrain the covariance matrices to be positive definite, we set the negative eigenvalues of ∑ij equal to a small constant, 10−12. For p<62, we instead generate eij from the p×p submatrices of Dij∑ijDij constructed from the rows and columns corresponding to the first p/2 features in each hemisphere. In this simulation scenario, S is no longer the pooled covariance structure since the covariance site effects Ωi can take any form and these site‐specific covariance structures do not necessarily combine across sites to resemble S. Instead of focusing on recovery of an underlying structure, we evaluate site effects throughout these simulations via ML.

24We design four simulation experiments to test how CovBat performs in multiple settings with varied covariance effects. In our ComBat simulation, we generate data without any covariance site effect from a model that resembles the ComBat model. In the Diagnosis Affects Mean simulation, we then introduce a covariance site effect to assess if CovBat can outperform ComBat in harmonization of covariance. We next introduce a simulated diagnosis effect on covariance in the Diagnosis Affects Covariance simulation, which better illustrates how the detection of the simulated diagnosis is affected by related covariance effects. Finally, we design a Covariance Only simulation, which has no site or diagnosis effects in mean or variance, but still contains site and diagnosis effects in covariance. This final simulation illustrates how effects on covariance can influence ML results and be addressed through our proposed harmonization method.

25In the ComBat simulation, we impose site effects in mean and variance while having the simulated diagnosis affect only the mean of the observations. We choose βi=−0.5 for p/4 regions of interest in both the left and right hemispheres to impose that about half of the ROIs are negatively associated with the simulated diagnosis. We also choose the Ωi and Ψ to be 62×62 zero matrices to ensure that the covariance does not depend on site or the simulated diagnosis.

26In the Diagnosis Affects Mean simulation, we again impose that the simulated diagnosis only affects the mean of the measurements but also introduce a site effect in covariance. We keep the same β and Ψ as in the ComBat simulation. However, we choose Ωi to be distinct 62×62 high‐rank matrices to distinguish covariance structures across sites. These Ωi are constructed by downsampling three distinct images to 62 × 62 pixels and then scaling the values so that the diagonal is a vector of ones. These matrices thus have different eigenvectors from each other and also from the sample covariance matrix S, which simulates complex site covariance structures which do not have the same eigenvectors.

27In the Diagnosis Affects Covariance simulation, we assume that the simulated diagnosis affects not only mean, but also covariance. We choose the diagnosis effect on covariance to be proportional to a site's covariance shift. This scenario represents a situation where detection of the diagnosis using ML could be highly influenced by the presence of site effects. We use the same β value as in the ComBat and diagnosis affects mean simulations but choose Ψ to be related to Ω3 to force confounding of Site 3 and diagnosis effects on covariance. To achieve this, we set Ψ=−3/4Ω3.

28In the Covariance Only simulation, we assume that both site and the diagnosis influence the covariance, not the mean or variance. We fix γ=0 and δ=1 to remove site effects in mean and variance while also using the same Ωi as in the Diagnosis Affects Mean and Diagnosis Affects Covariance simulations. Furthermore, we modify the diagnosis effect by setting β=0 while keeping Ψ=−3/4Ω3 for the diagnosis effect in covariance.

29In our simple covariance effects simulation with rank one covariance effects, we generate 1,000 datasets across varying sample sizes and number of features to evaluate recovery of the underlying covariance structure. For each site, we calculate the average Frobenius distance across datasets between each sample within‐site covariance matrix and the true covariance S. We then report the average across sites before and after harmonization, where CovBat harmonization is performed on PCs that explain 95% of the variation.

30For the other simulation settings, we generate 1,000 datasets for several choices of within‐site sample size and number of features and perform experiments to evaluate detection of site and disease using ML. For each dataset, we (a) randomly split the sample into 50% training data and 50% validation data, (b) train a random forests algorithm to recognize if the observations either are from Site 1 or have the simulated diagnosis, and (c) assess predictive performance on the testing data. Our random forests algorithm is implemented via the randomForest package on CRAN with default parameters using R version 3.6.1. Lower AUC for detection of Site 1 and higher AUC for detection of the simulated diagnosis indicate improved harmonization. CovBat harmonization is performed on PCs that explain 95% of the variation. For prediction of Site 1, we avoid the possibility that site could be detected through the simulated diagnosis by using linear models to residualize out disease from each simulated cortical thickness value. For simulated datasets where the training set does not contain observations with the simulated diagnosis, the random forest algorithm cannot be trained so we generate another dataset to replace it. We repeat these ML experiments while varying the number of PCs included in the CovBat harmonization. For CovBat including PCs that explain 95% of the variation, we also perform MANOVA for testing associations with site and simulated diagnosis and report the rejection rate at a type I error rate of 0.05 across the 1,000 datasets.

31We first perform simulations to assess whether either harmonization method can recover the underlying covariance structure in our simple simulation setting and harmonize covariance matrices generally. For ComBat and CovBat, we include our simulated diagnosis status as a covariate. We apply our CovBat method using the number of PCs that explain 95% of the variation. In the Simple Covariance Effects simulation, Figure 1 shows that CovBat outperforms ComBat in recovery of the true covariance structure (denoted S in Section 2.5) across all parameters considered. Remaining deviation from the true covariance can be explained by error in covariance estimation; even with 10,000 samples per site and 62 features the distance of the pooled covariance estimate from the true covariance is still 3.14. In Table 1, we observe in the same simulation setting with a within‐site sample size of 250 and 62 simulated features that CovBat performs best in harmonizing within‐site covariance matrices. That result is replicated with a more complex covariance site effect in the Diagnosis Affects Mean simulation as shown in Table 2.

32We then evaluate CovBat through two ML experiments across all simulation settings considered. For our main analyses, we simulate 62 cortical thicknesses for 250 subjects per site. We begin by examining a ComBat simulation, where we generate data from the original ComBat model. That is, we impose mean and variance site effects while also simulating a diagnosis that has an effect on the mean. Figure 2a,b show that CovBat performs almost identically to ComBat in this scenario, showing that our method performs competitively in the absence of covariance effects. In the Diagnosis Affects Mean simulation, Figure 2c,d shows that CovBat also substantially reduces the chance of detecting site and performs similarly to ComBat for detection of the simulated diagnosis. These simulations demonstrate that CovBat performs at least as well or better than ComBat when the simulated diagnosis affects only the mean of the observations.

33We then incorporate a diagnosis effect on covariance in our Diagnosis Affects Covariance simulation and show that CovBat reduces detection of site (Figure 3a) and maintains the association with the simulated diagnosis (Figure 3b). In order to further emphasize the importance of covariance effects, we investigate a Covariance Only simulation where both site and diagnosis effects exist only in the covariance of observations, but not the mean or variance. In unharmonized data, we observe high mean AUC values for detection of Site 1 and detection of the simulated diagnosis, both of which are essentially unaffected after implementing ComBat. After CovBat though, we see substantial improvements on both metrics (Figure 3c,d). While CovBat performs well in both simulations, we find that CovBat does not entirely remove the severe covariance site effect. Regardless, we observe that CovBat offers notable improvements over ComBat owing to its harmonization of covariance.

34We then conduct additional analyses to assess the robustness of CovBat to reductions in sample size per site and number of features. In Table S1, we show that the random forests largely hold for simulations without a diagnosis effect on covariance; however, in the absence of a covariance site effect CovBat performs slightly worse overall and in situations with small sample size per site (ni=25) and larger number of features (p≥48) we observe that CovBat can inflate detection of site. Table S2 shows that similar results hold in simulations with a diagnosis effect on covariance with CovBat showing good performance overall but poor performance with small sample sizes and large number of features. We also show in Table S3 through the Covariance Only simulation that site and diagnosis can both be detected even without affecting the mean of observations, as demonstrated by high AUCs for detection across all settings. To assess if our findings may be tied to the ML paradigm, we additionally perform MANOVA for site and diagnosis status across Diagnosis Affects Mean and Diagnosis Affects Covariance simulations and show in Tables S4 and S5 that associations with site are reduced after CovBat across all scenarios while associations with diagnosis are preserved. We repeat these MANOVA analyses in the Simple Covariance Effects simulation and show in Table S6 that CovBat performs better than ComBat at controlling the rejection rate across nearly all parameters considered. However, CovBat only successfully controls the rejection rate at 5% in simulations with less features, small sample size, and large sample size. The controlled rejection rate in small samples however could be explained by the low power of MANOVA in small sample size relative to the number of features (Stevens, ).

35To better inform the choice of PCs in CovBat, we evaluate the simulation results obtained across varying number of PCs. For simple covariance effects with varying low rank structures, Figure S2 shows that harmonizing smaller numbers of PCs yields suboptimal results for reducing the chance of detecting sites. Across varying ranks of simulated site effects contained in the covariance structures (2–12 PCs), simulation results show that the median AUC for detecting site achieves the lowest value when CovBat includes PCs that explain around 90% of total variation for a sample size of 100. With a larger sample size of 250, the best performance was achieved when CovBat includes PCs that explain 95% of the total variation. Across both sample sizes, the U‐shape curves in median AUC for detecting sites indicate that including excess numbers of PCs in CovBat hurts the performance due to overfitting and the optimal number of PCs increases with sample size. For the high rank covariance effect in our Diagnosis Affects Covariance, we observe in Figure S3 that for a sample size of 250 increasing the number of PCs lowers the AUC for detection of site across the whole range of considered PCs. However, looking at a sample size of 50, we observe that AUC for site detection increases as we select PCs that explain between roughly 95% and 100% of the variation. For detection of simulated diagnosis, Figure S4 shows that CovBat largely maintains detection of diagnosis across all numbers of PCs included and sample sizes, with minor increases in median AUC as more PCs are included.

36We evaluate the performance of harmonization in the ADNI data by comparing three different approaches. First, we test the original unharmonized data in our subsequent experiments. Second, we produce ComBat‐harmonized data by performing ComBat while including age, sex, and diagnosis status as covariates. Third, we obtain CovBat‐harmonized data by running CovBat with varying number of PCs, while also including age, sex, and diagnosis status. For our primary analyses, we run CovBat while including 37 PCs, which cumulatively explain 95% of variation. We include additional analyses evaluating how our results vary with the choice of PCs.

37Using the full ADNI sample, we evaluate whether the harmonization procedures affect the results of ML using the neuroimaging measures as patterns for a prediction algorithm. We achieve this through a Monte Carlo split‐sample experiment where we (i) randomly split the subjects into 50% training set and 50% validation set, (ii) train a random forests algorithm to detect either scanner manufacturer or a binary clinical covariate, and (iii) assess predictive performance on the validation set via AUC. Our random forests algorithm is implemented via the randomForest package on CRAN with default parameters using R version 3.6.1. We train separate models for unharmonized, ComBat‐harmonized, and CovBat‐harmonized data where both harmonization methods are performed including age, sex, and diagnosis status as covariates. We perform steps (i)–(iii) 100 times for each dataset and report the mean AUC along with standard deviation. For these experiments, lower AUC for detection of scanner manufacturer and higher AUC for detection of clinical covariates would indicate improved harmonization. For prediction of scanner manufacturer, we avoid the possibility that scanner could be detected through the covariates age, sex, and disease status by residualizing out these variables from each cortical thickness value via linear models. To assess performance across choices of PCs in CovBat, we repeat our ML experiments across varying numbers of PCs explaining between 44 and 100% of the variation.

38We also evaluate the harmonization methods using multivariate analysis of variance (MANOVA). MANOVA tests for differences in mean across groups in multivariate data, but is known to be sensitive to differences in covariance across groups. We perform MANOVA across scanner manufacturer, sex, and diagnosis status using Pillai's trace, which is known to be more robust to inhomogeneity in covariance than other alternative test statistics (Olson, ). We report p‐values for these associations before and after harmonization.

39We first examine the empirical covariance of the ADNI data before and after harmonization. To evaluate site differences, we investigate the three largest ADNI sites. Site A consists of 23 subjects acquired on a Siemens Symphony 1.5T scanner while Sites B and C each consist of 20 subjects acquired on GE Signa Excite 1.5T scanners. See Table 3 for demographic details. To avoid influence of site demographics on the covariance matrices, we residualize the cortical thickness measures across the three sites jointly on age, sex, and diagnosis status in each dataset using a linear model. Figure 4 shows the covariance matrices for each site using the residualized cortical measures both before and after harmonization (ROI labels are shown in Figure S1). The differences between the unharmonized covariance matrices are striking. Especially notable are the increased positive covariance across most pairs of cortical regions in Site A and the weakened correlation between the right and left hemispheres in Site C, visible as the diagonal line in the top‐left and bottom‐right quadrants. Visually, the covariance differences remain similar after applying ComBat. These inter‐site differences are considerably mitigated after CovBat.

40We also provide quantitative comparisons for pairwise distances across sites before and after harmonization in Table 4. A tuning parameter of the CovBat model is the desired proportion of variance explained in the dimension reduction space, which we select at 95% (37 PCs). To ensure that our results do not depend strongly on the choice of tuning parameter, we also report the minimum and maximum of the pairwise Frobenius norms after applying CovBat with percent variation explained ranging from 44% (2 PCs) to 100% (62 PCs). We report the results of this sensitivity analysis in parentheses. We find that ComBat adjustment can modestly harmonize the covariance matrices but CovBat adjustment shows large reductions in the between‐site distances across a range of tuning parameter choices.

41To evaluate the potential impact of site effects in covariance using ML, we conduct a Monte‐Carlo split‐sample experiment for prediction of scanner manufacturer labels using all 213 ADNI subjects before and after harmonization with existing methods. We train using data harmonized with the state‐of‐the‐art ComBat method and our proposed method, CovBat. Figure 5a shows that Siemens sites are easily identifiable based on unharmonized cortical thickness measurements (median area‐under‐the‐curve [AUC] 0.89, IQR 0.87–0.90), which is consistent with recent findings (Glocker et al., ). We also note that scanner manufacturer is still detected after ComBat is applied (0.66, 0.64–0.68). After CovBat, the ML method's performance for differentiating between sites is close to chance (0.46, 0.44–0.48). CovBat's performance depends on the number of PCs included in the model, but Figure S5 shows that the performance gain for each PC becomes negligible around the number of PCs that explain 95% of the variation. DeLong's test results shown in Figure S2 suggest that these AUC values for site detection are significantly different between ComBat and CovBat. Using MANOVA, Table 5 shows that the association with scanner manufacturer is statistically significant in unharmonized and ComBat‐adjusted data but is eliminated in CovBat‐adjusted data.

42It is well‐known that cortical thickness differs substantially by sex and AD status (Lerch et al.,; Sowell et al., ). To assess whether CovBat maintains biological associations of interest, we perform two ML experiments using the full ADNI data to classify healthy versus AD and to differentiate patients by sex. Figure 5b,c show that both of these biological associations are retained after either harmonization method. For AD classification, the median AUC increases from 0.74 (IQR 0.72–0.75) in unharmonized data to 0.78 (0.75–0.79) in ComBat‐adjusted data to 0.79 (0.77–0.81) in CovBat‐adjusted. Similarly, the median AUC for detection of sex increased from 0.66 (0.64–0.68) to 0.69 (0.67–0.71) to 0.70 (0.67–0.72). For detection of both AD and male sex, DeLong's test results plotted in Figure S6 support that the AUCs are not significantly different between ComBat and CovBat. These findings suggest that CovBat not only provides thorough removal of site effects, but also maintains clinical associations. Figure S4 shows that CovBat retains these associations across varying number of PCs included in the model. Appendix A1 shows that similar results hold for prediction of age, where both ComBat and CovBat reduce root‐mean‐square error for prediction of age compared to the unharmonized data. Appendix A2 shows that these results for both detection of site and biological associations largely hold even when CovBat is trained on a subset of the data and all sites are included in both the training and validation sets. MANOVA results in Table 5 show that the significant associations with diagnosis and sex are retained after either harmonization method.

43The growing number of multi‐site studies across diverse fields has spurred the development of harmonization methods that are general, but also account for field‐specific challenges. In neuroimaging research, the rise of ML in neuroimaging has established an unmet need for harmonization of covariance. We demonstrate that strong site effects in covariance exist, influence downstream ML experiments, and remain after performing the state‐of‐the‐art harmonization. We then propose a novel method and demonstrate that it is effective in removing site differences in covariance and retaining the detection of biological associations via ML. Simulation studies show similar ML results and demonstrate that CovBat performs well across a variety of settings and sample sizes.

44In ADNI data, we show that substantial differences exist in the covariance structures of cortical thickness observations and can be mitigated through our proposed method. We furthermore show that ML can detect these site effects, whether through ML or conventional multivariate analyses. These results mirror recent studies that predict scanner from neuroimaging features with high accuracy (Glocker et al., ) and a recent study demonstrating that ComBat is insufficient to prevent detection of Siemens‐manufactured scanners in a large multi‐site dataset (Nielson et al., ). We then demonstrate that CovBat can almost entirely prevent site detection in the ADNI dataset. To ensure generalizability of these results, implementation of CovBat in other multi‐site studies of varying experimental designs should be pursued in the future.

45In simulation, CovBat shows generally strong performance in removing site effects in medium and large sample sizes across varying number of features and complexity of the site effect. We demonstrate that CovBat almost fully removes covariance site effects when they exist in the principal component directions, but deviations from the true covariance still remain due to error in covariance estimation, overcorrection of PCs without site difference, and the remaining site effects in the marginal variances. Additionally, ML results show that considerable site effects may remain in more complex scenarios. Caution should be taken in attempting to address covariance site effects in smaller samples with many features. We show potential increases in site detection in these situations, which are potentially the result of poor covariance estimation in high‐dimensional settings. Through investigating the performance of CovBat across varying number of PCs being harmonized, we conclude that the chance of detecting sites could increase when including excessive number of PCs in cases with simple site covariance effects or small sample size relative to the number of features. We generally recommend selecting number of PCs that cumulatively explain 90–95% of the total variation depending on the sample size. While we do not observe these limitations through MANOVA, we also acknowledge that MANOVA may be underpowered in our scenarios with low sample size and high dimensionality as shown in previous studies (Stevens, ).

46Our proposed method harmonizes covariance across sites by removing mean and variance shifts in the principal components space, which we show to be effective in addressing the covariance effects we observe. This idea resembles spectral models, which also relate covariates to the eigendecomposition of covariance matrices (Boik, ). Our method assumes that the ideal covariance structure exists in the eigenspace of the full data covariance matrix. As we show through our simulations, in some cases this model may be insufficient to remove site effects, which do not resemble the covariance structure of the full data. Potential extensions could incorporate methods that model site effects as separate low‐rank structures (Hoff & Niu, ) or identify projections most related to site (Zhao, Wang, Mostofsky, Caffo, & Luo, ). However, implementation of these methods in a harmonization framework may not be as straightforward as our proposed method.

47A limitation of our methodology is that CovBat is a covariate‐assisted harmonization method similar to ComBat and requires specification of covariates to protect in the data. Associations with covariates not included in the harmonization step can certainly be removed alongside site effects. Furthermore, nonlinear covariate associations may not be adequately captured by the linear model and future extensions of CovBat could consider incorporating previous work on general additive models in ComBat to capture complex covariate effects (Pomponio et al., ). Recent articles have identified situations where spurious associations can be introduced via ComBat (Nygaard, Rødland, & Hovig,; Zindler, Frieling, Neyazi, Bleich, & Friedel, ). While we do not observe CovBat introducing false positives in our investigation, care must be taken in implementing CovBat protecting for the outcome of interest especially in unbalanced study designs. We reiterate previous advice that analyses should be performed with and without harmonization and the analysis design be made very clear to ensure that results can be interpreted properly (Zindler et al., ).

48Our study demonstrates that site effects can exist in the covariance of structural neuroimaging data and can be mitigated via our proposed methodology. Future studies should determine how scanner properties can influence the covariance structure of the data and if other multi‐site multivariate neuroimaging studies contain similar effects. Further methodological work could utilize other covariance modeling strategies in order to address more complex site effects. Since our method operates on general multivariate data, our findings extend directly to functional, metabolic, and other imaging modalities. However, our method does not currently handle time series from functional imaging, which could be the subject of future investigation. Further studies should also determine the extent to which multivariate statistical and ML studies of genomic data are susceptible to the biases documented.

49All of the postprocessing analysis was performed in the R statistical software (V3.6.1). CovBat is available for both R and Python (https://github.com/andy1764/CovBat_Harmonization). Reference implementations for ComBat are available in R and Matlab (https://github.com/Jfortin1/ComBatHarmonization) and in Python (https://github.com/ncullen93/neuroCombat). Our R implementation of CovBat runs in reasonable time, even in large samples. We generate observations with 62 features with varying sample sizes. On a MacBook Pro (16‐in., 2019) with a 2.3 GHz 8‐Core Intel Core i9 and 32 GB 2667 MHz DDR4 memory, CovBat runs in 0.014 s for 100 samples, 0.090 s for 1,000 samples, 1.088 s for 10,000 samples, and 12.235 seconds for 100,000 samples.

50The authors declare no potential conflict of interest.

51Russell T. Shinohara and Haochang Shou contributed equally to this work. Andrew A. Chen: formal analysis, conceptualization, methodology, software, investigation, writing—original draft. Joanne C. Beer: data curation, writing—review & editing. Nicholas J. Tustison: data curation, Writing—review & editing. Philip A. Cook: data curation, writing—review & editing. Russell T. Shinohara: supervision, conceptualization, methodology, writing—review & editing. Haochang Shou: supervision, conceptualization, methodology, writing—review & editing.