ComBat (Johnson et al., 2007) has become a popular way of removing scanner or site (“batch”) effects from group-level imaging data since Fortin and colleagues showed it works well for diffusion MRI and cortical thickness (Fortin et al., 2017, 2018). The usual call estimates each batch location and scale parameters from a sample and applies these parameters to that same sample to remove such unwanted effects. In the context of cross-validation, however, one hopes to avoid leakage of information between train and test data. Applying ComBat to the whole sample before splitting would create dependencies between train and test. A better alternative is to estimate the parameters from the training data and apply to the test data.

Here two MATLAB functions are provided. One that fits ComBat on a training fold and returns the parameters, and another that takes these parameters (without new fitting) and applies them to new data. These are combatfit.m and combatapply.m. If you are familiar with ComBat and with cross-validation, you don’t need to continue reading. Otherwise, some additional background is provided below.

The ComBat model

ComBat is a location-and-scale model. For feature $g$, sample $j$ in batch $i$ (Johnson et al., 2007):

$Y_{ijg} = \alpha_g + X_{ij}\beta_g + \gamma_{ig} + \delta_{ig}\varepsilon_{ijg}$

The covariates $X$ (age, sex, etc) are there so that biological variation are preserved and not absorbed into the scanner/site terms. Data are standardised using a grand mean and a pooled residual variance, batch-specific $\gamma$ and $\delta$ are estimated then shrunk with empirical Bayes across features. Batch effects are removed before the data are returned to the original scale. Biological effects are removed before the scale step and added back afterwards such that they are not distorted by $\delta$ or removed.

Here, combatfit.m does that fit and, unlike the original combat.m, it also returns the model parameters, i.e., grand mean, regression coefficients, pooled variance, the empirical-Bayes location and scale (gamma_star, delta_star), and the covariate means. The harmonised training data are also returned, and should match the original combat.m (although a few changes make combatfit.m faster for large datasets). Then, combatapply.m uses those parameters on new data. It does not re-estimate parameters.

Usage

Differently than original ComBat, combatfit.m and combatapply.m take the input data array organized as $N$ by $P$ (e.g., subjects by features), which is the usual layout for the general linear model (GLM), rather than features by samples. Do not include an intercept in mod; the batch dummy variables already span it.

A typical fold looks like:

% Fit ComBat on the training fold:
[yTrainH, gm, Bh, vp, gs, ds, mm] = ...
    combatfit(yTrain, batchTrain, modTrain, true);
% Apply the same parameters to the test fold:
yTestH = combatapply(yTest, batchTest, modTest, ...
    batchTrain, gm, Bh, vp, gs, ds, mm);

The fourth argument to combatapply.m must be the original training batch vector. The rows of gamma_star and delta_star follow unique(batch) (sorted), and that order has to be reconstructable.

For the detailed syntax type at the MATLAB/Octave prompt: help combatfit and help combatapply, or simply check the source code.

Covariates

If mod is used in training, the same columns must be supplied for the test set. New covariates are centred at the training means, then the training coefficients are applied, which is the out-of-sample analogue of the original model.

Moreover, mod should contain only covariates that must be preserved, and not the data that are being predicted. If the aim is to predicting diagnosis, putting that same variable in mod on the test fold leaks the label into the harmonised features. In other words, protect nuisance biology (age, sex, and similar), not the outcome.

Sites/scanners never seen

For a scanner or site to be harmonized in the test set, that scanner or site must be present in the training set. Thus, this pair of functions will not work in cross-validation schemes that leave an entire site or scanner out.

Availability

The functions run in MATLAB or Octave, and can be downloaded from GitHub: combatfit.m and combatapply.m. The original ComBat can be found in its GitHub repository.

References

• Fortin JP, Cullen N, Sheline YI, Taylor WD, Aselcioglu I, Cook PA, Adams P, Cooper C, Fava M, McGrath PJ, McInnis M, Phillips ML, Trivedi MH, Weissman MM, Shinohara RT. Harmonization of cortical thickness measurements across scanners and sites. Neuroimage. 2018;167:104–120.
• Fortin JP, Parker D, Tunç B, Watanabe T, Elliott MA, Ruparel K, Roalf DR, Satterthwaite TD, Gur RC, Gur RE, Schultz RT, Verma R, Shinohara RT. Harmonization of multi-site diffusion tensor imaging data. Neuroimage. 2017;161:149–170.
• Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 2007;8(1):118–127.