ComBat Harmonization Overview

• ComBat harmonization is an empirical Bayes framework that removes scanner- and protocol-related batch effects from imaging data.
• It preserves biological signals by modeling covariate effects alongside additive and multiplicative adjustments.
• The method has been extended for nonlinear covariates and federated settings, enhancing reproducibility in multi-site studies.

ComBat harmonization is an empirical Bayes location–scale framework for removing unwanted batch effects—typically site, scanner, protocol, or tracer differences—from imaging-derived features while preserving biological variability tied to covariates of interest. In neuroimaging, these batch effects are commonly expressed as additive shifts and multiplicative scale distortions in feature distributions, and ComBat estimates them feature-wise while borrowing strength across features to stabilize estimation under finite per-batch sample sizes. The method is widely used in diffusion MRI, structural MRI, PET, and radiomics, and has been extended to nonlinear covariates, normative-reference alignment, pathology-robust estimation, spatially informed image models, and federated deployment (McMaster et al., 31 Mar 2026, Jodoin et al., 19 May 2025, Silva et al., 19 Jan 2026).

1. Concept and scope

ComBat was originally developed for microarray gene expression data and later widely adopted in neuroimaging to remove batch effects while retaining biological signal (Xu et al., 2024). In its neuroimaging usage, a “batch” is typically a site, scanner, scanner model, sequence, protocol, or tracer, depending on which factor is regarded as the most proximal driver of technical variability (McMaster et al., 31 Mar 2026, Zhang et al., 16 Jun 2026). The preserved signal is encoded through covariates such as age, sex, handedness, diagnosis, intracranial volume, or other study-specific variables that should remain interpretable after harmonization (David et al., 18 Mar 2026, Silva et al., 19 Jan 2026).

The method is especially relevant in multi-site imaging, where acquisition heterogeneity violates the i.i.d. assumption that underlies many statistical and machine-learning procedures. Diffusion MRI studies describe scanner/protocol differences as inducing site-specific mean and dispersion changes in derived metrics such as fractional anisotropy, mean diffusivity, or bundle-level features (Jodoin et al., 19 May 2025). Structural MRI and PET studies report analogous batch structure in cortical thickness and SUVR maps, respectively (Reinhardt et al., 2 Apr 2026, Zhang et al., 16 Jun 2026). This places ComBat at the intersection of batch-effect correction, covariate adjustment, and reproducible multi-center inference.

A central design principle is that harmonization is not mere normalization. Classical global z-scoring or histogram matching does not explicitly separate biological covariate effects from scanner effects. ComBat does so by estimating the biological mean structure and then adjusting only the residual location and scale associated with batch membership (Sharifian et al., 27 Aug 2025).

2. Statistical formulation and empirical Bayes estimation

The standard ComBat model expresses an observed feature as the sum of a feature-specific intercept, preserved biological covariate effects, and batch-specific additive and multiplicative terms. In a common notation,

$$y_{ijv} = \alpha_v + \mathbf{x}_{ij}^{\top}\boldsymbol{\beta}_v + \gamma_{b(i)v} + \delta_{b(i)v}\,\varepsilon_{ijv}, \qquad \varepsilon_{ijv}\sim\mathcal{N}(0,\sigma_v^2).$$

Here, $y_{ijv}$ is feature $v$ for subject $i$, $\mathbf{x}_{ij}$ contains preserved covariates, $\gamma_{b(i)v}$ is the batch-specific additive effect, and $\delta_{b(i)v}$ is the batch-specific multiplicative effect (David et al., 18 Mar 2026, Silva et al., 19 Jan 2026).

Estimation proceeds by first fitting the biological mean structure and standardizing residuals. A typical standardization step is

$$z_{ijv} = \frac{y_{ijv} - \hat{\alpha}_v - \mathbf{x}_{ij}^{\top}\hat{\boldsymbol{\beta}}_v}{\hat{\sigma}_v}.$$

ComBat then estimates batch-specific location and scale parameters under empirical Bayes priors such as

$$\gamma_{bv}\sim \mathcal{N}(\mu_{\gamma,v},\tau_{\gamma,v}^2), \qquad \delta_{bv}^{2}\sim \mathrm{InvGamma}(a_{\delta,v},b_{\delta,v}),$$

or equivalent site-wise normal and inverse-gamma parameterizations (McMaster et al., 31 Mar 2026, Hoang et al., 2024). The essential mechanism is shrinkage: batch-specific means and variances are pulled toward pooled priors estimated across features, which stabilizes harmonization when some batches are modestly populated (Sharifian et al., 27 Aug 2025).

After posterior or posterior-like empirical Bayes estimates are obtained, harmonized values are reconstructed by removing the estimated batch shift, rescaling the residual variance, and restoring the preserved biological mean function. A standard adjustment is

$$y^{\mathrm{ComBat}}_{ijv} = \hat{\alpha}_v + \mathbf{x}_{ij}^{\top}\hat{\boldsymbol{\beta}}_v + \frac{ y_{ijv}-\hat{\alpha}_v-\mathbf{x}_{ij}^{\top}\hat{\boldsymbol{\beta}}_v-\hat{\gamma}_{b(i)v} }{ \hat{\delta}_{b(i)v} }.$$

This location–scale back-transformation defines the canonical ComBat adjustment across much of the literature (Zhang et al., 16 Jun 2026, Marzi et al., 2022).

Classical ComBat is feature-wise: it does not model dependence across features, graph topology, or spatial structure unless a later extension explicitly introduces such information (Xu et al., 2024, Reinhardt et al., 2 Apr 2026). That property is computationally convenient and partly explains its broad adoption, but it also defines several of its known failure modes.

3. Major variants and methodological extensions

A large methodological family now sits under the ComBat label. The most established extension is ComBat-GAM, which replaces a strictly linear covariate term with a generalized additive mean function, typically to capture nonlinear age effects. In infant diffusion MRI, this is described as essential because neurodevelopmental trajectories are steep and nonlinear even over weeks (McMaster et al., 31 Mar 2026). Other variants modify the target of harmonization, the prior structure, the robustness mechanism, or the deployment regime (McMaster et al., 31 Mar 2026, David et al., 18 Mar 2026, Girard et al., 6 Nov 2025, Reinhardt et al., 2 Apr 2026, Hoang et al., 2024, Silva et al., 19 Jan 2026).

Variant	Key modification	Reported use
ComBat-GAM	Smooth nonlinear covariates such as age	Infancy and lifespan modeling
Pairwise-ComBat	Align moving site directly to a fixed reference	Normative/reference harmonization
CovBat	Harmonize covariance across features	Multivariate feature consistency
Clinical-ComBAT	Site-wise reference mapping with polynomial basis and adaptive variance priors	Small or incrementally added clinical cohorts
Robust-ComBat	Exclude flagged outliers from site-effect estimation	Clinical dMRI with high pathology prevalence
Tensor-ComBat	Spatial Bayesian tensor model plus post-hoc ComBat adjustment	Voxel-level structural MRI

ComBat-GAM augments the model with a smooth term such as $y_{ijv}$0,

$y_{ijv}$1

thereby preserving nonlinear age-related variation while still using empirical Bayes shrinkage for batch effects (McMaster et al., 31 Mar 2026). Pairwise-ComBat instead aligns each moving site to a fixed normative reference rather than to a pooled multi-site mean, which is particularly attractive for normative modeling because the target does not move as sites are added (Jodoin et al., 19 May 2025, David et al., 18 Mar 2026).

Clinical-ComBAT retools ComBat for clinical deployment by harmonizing each site independently to a normative reference, replacing the linear covariate model with a polynomial basis and using variance priors adapted to small cohorts (Girard et al., 6 Nov 2025). Robust-ComBat addresses pathological contamination by inserting a subject-level outlier-identification stage; flagged subjects are excluded from empirical Bayes estimation of site parameters, but all subjects are subsequently harmonized using the filtered estimates (David et al., 18 Mar 2026).

Two other extensions address structural limitations of the original feature-wise model. Tensor-ComBat uses Bayesian tensor response regression with low-rank PARAFAC decomposition to estimate spatially distributed scanner effects and then applies a post-hoc ComBat adjustment, explicitly incorporating voxel spatial configuration rather than treating voxels as interchangeable (Reinhardt et al., 2 Apr 2026). Cluster ComBat and Fed-ComBat adapt the framework to distributed settings: the former clusters sites or samples so that harmonization parameters can be shared across clusters and reused for newly joined sites, whereas the latter replaces centralized fixed-effect estimation with federated optimization and local empirical Bayes correction (Hoang et al., 2024, Silva et al., 19 Jan 2026).

4. Reported applications and empirical behavior

In infant diffusion MRI from the HBCD release 1.1, ComBat-GAM was applied to bundle-level scalar metrics from 437 infants scanned across six unique scanner models. The study restricted analysis to subjects aged between 0 and 0.2 years, used scanner model as the batch variable, and preserved age through a smooth spline and sex linearly. After ComBat-GAM, the authors reported “zero statistically significant differences between the distributions from any scanner model following FDR correction,” and scanner-related effect sizes collapsed from predominantly medium-to-large before harmonization to almost entirely small after harmonization: 266/268 small effects, 2/268 medium, and 0/268 large (McMaster et al., 31 Mar 2026).

In diffusion MRI connectivity matrices, ComBat was applied directly to vectorized upper-triangular matrix entries and then graph measures were recomputed. In a two-site cohort of 168 age-matched, sex-matched normal subjects, ComBat effectively eliminated site effects for global efficiency and modularity and outperformed mean shift and CycleGAN on those measures. However, all methods exhibited poor performance when harmonizing average betweenness centrality derived from mean streamline length matrices, and ComBat nearly tripled the coefficient of variation for that measure (Xu et al., 2024). This establishes that successful edge-wise harmonization does not guarantee uniformly successful graph-level harmonization.

Voxel-wise structural MRI has also been used as a ComBat target. In cross-domain sex classification from 3D T1-weighted MRI, domain shift between IXI and OASIS3 reduced cross-domain performance to chance level, whereas ComBat harmonization of voxel intensities improved the cross-domain balanced accuracy of the best model from approximately 0.50 to 0.61 while maintaining within-domain performance (Sharifian et al., 27 Aug 2025). In a distinct structural MRI setting, Tensor-ComBat was evaluated on over 2100 ADNI-1 T1w-MRI scans and achieved greater scanner effect removal, improved biological prediction, and superior reproducibility relative to state-of-the-art approaches such as ComBat; mean voxel/scanner $y_{ijv}$2 was 0.75 for Tensor-ComBat versus 0.63 for ComBat, with fewer residual outliers (Reinhardt et al., 2 Apr 2026).

Radiomics studies have shown that the appropriate ComBat configuration can depend on biological subgrouping. In pulmonary nodule CT radiomics, harmonizing without distinction rendered an average 2.1% of features acquisition-independent, harmonizing with a covariate for subgroup membership raised that to 27.3%, and harmonizing benign, malignant, and screening subgroups separately raised it to 90.9%. Downstream LASSO-SVM malignancy prediction on screening scans also improved relative to collective harmonization when subgroup differences were preserved or handled separately (Huchthausen et al., 2024).

Classical ComBat can also underperform when source and target cohorts differ biologically. In cross-tracer tau PET harmonization, ComBat corrected the global mean shift but left region-dependent biases and subgroup misalignments. On held-out test data, ComBat had higher residual site separability than the proposed transport-based alternative, and tau-positivity sign mismatches occurred exclusively in the neg→pos direction: 33/503 left hemispheres and 26/503 right hemispheres (Zhang et al., 16 Jun 2026).

5. Assumptions, limitations, and recurrent controversies

Standard ComBat assumes that preserved covariate effects are correctly specified, typically linearly unless a GAM extension is used, and that batch effects can be represented as additive and multiplicative distortions around that biological mean structure (Jodoin et al., 19 May 2025, McMaster et al., 31 Mar 2026). It also relies on adequate overlap of covariate distributions across sites. Severe covariate–batch confounding can cause over-correction or under-correction, and infant diffusion MRI work explicitly recommends checking whether age or sex distributions differ systematically across scanners or sites (McMaster et al., 31 Mar 2026).

A major technical limitation is slope heterogeneity. In the dMRI best-practices study, standard ComBat did not correct multiplicative distortion of biological slopes, and severe misalignment was observed when slope scaling departed strongly from unity, including at $y_{ijv}$3 and $y_{ijv}$4. The same study reported large train–test gaps for very small moving-site sample sizes $y_{ijv}$5, markedly improved behavior around $y_{ijv}$6–32, and diminishing returns beyond about $y_{ijv}$7. It also found poor generalization for narrow age spans of 10–20 years, large errors for male-only versus female-only harmonization, and compression of pathological cases toward the normative range when mixed healthy-control and pathology data were used to train harmonization parameters (Jodoin et al., 19 May 2025).

Pathology contamination is a recurrent concern because classical ComBat presumes site-wise subject distributions that are approximately Gaussian and representative of the same underlying biological population. Robust-ComBat formalized this problem in clinical diffusion MRI and showed that including pathological cases in site-effect estimation induces significant distortions when harmonizing to a normative reference. Across scenarios with pathology prevalence up to 80%, an MLP-based robust outlier compensation scheme consistently outperformed conventional statistical baselines and lowered harmonization error across all ComBat variants considered (David et al., 18 Mar 2026).

Another controversy concerns data leakage in machine-learning pipelines. Harmonizing the entire dataset before splitting allows test data to influence empirical Bayes estimates, and a 36-site structural MRI study showed that this can falsely overestimate downstream performance; the remedy proposed there was a harmonizer transformer fitted only on training folds and applied to held-out data (Marzi et al., 2022). A related but distinct issue arises when site and target are statistically dependent. In such settings, including the target as a preserved covariate can require access to the true test label at harmonization time, creating test-target leakage. The PrettYharmonize study documented this problem for ComBat-based methods and proposed pretending the target label at application time to avoid leakage while retaining target-related signal (Nieto et al., 2024).

A further misconception is that ComBat is uniformly effective on any feature representation. The connectivity-matrix study showed that feature-wise adjustment of edges ignores dependence across edges and graph topology, which can degrade topology-sensitive metrics such as average betweenness centrality (Xu et al., 2024). The tau PET study similarly showed that linear location–scale correction may be inadequate when batch effects are region-dependent and nonlinear or when cohort subgroup composition differs markedly (Zhang et al., 16 Jun 2026).

6. Implementation, diagnostics, and deployment practice

Across applications, implementation begins with careful batch definition. Scanner model, site, tracer, or acquisition parameter can all serve as the batch variable, but several studies recommend modeling the most proximal driver of technical variability and avoiding overly fine batch splitting when sample sizes are small (McMaster et al., 31 Mar 2026, Huchthausen et al., 2024). Covariates to preserve are study-dependent but commonly include age and sex, with handedness, diagnosis, APOE, intracranial volume, or subgroup indicators added where appropriate (David et al., 18 Mar 2026, Silva et al., 19 Jan 2026). In infancy and lifespan settings, age is frequently modeled nonlinearly with GAM splines rather than linearly (McMaster et al., 31 Mar 2026, Silva et al., 19 Jan 2026).

The operational workflow is typically: extract a subject-by-feature matrix, fit the ComBat model on training data only, estimate empirical Bayes batch parameters, transform held-out data with the learned parameters, and then perform downstream inference or prediction. Train-only fitting is essential. The harmonizer-transformer formulation introduced for machine-learning pipelines encapsulates ComBat in a fit/transform interface so that $y_{ijv}$8, $y_{ijv}$9, $v$0, $v$1, and $v$2 are estimated only on training folds, thereby avoiding leakage (Marzi et al., 2022).

Diagnostics are modality-specific but highly consistent in spirit. Reported checks include pre/post density and quantile–quantile comparisons, ANOVA across batches with Benjamini–Hochberg FDR correction, effect-size summaries such as Cohen’s $v$3, Bhattacharyya distance for reference alignment, coefficient of variation, site-classifier balanced accuracy or $v$4, subgroup Wasserstein or KS distances, and inspection of estimated $v$5 and $v$6 parameters for plausibility (McMaster et al., 31 Mar 2026, David et al., 18 Mar 2026, Zhang et al., 16 Jun 2026). In reference-based workflows, the choice of reference matters, and Clinical-ComBAT recommends a large, diverse healthy-control reference site together with goodness-of-fit assessment based on Bhattacharyya distance (Girard et al., 6 Nov 2025).

Deployment requirements differ by variant. Pairwise and clinical reference-based models support harmonization to a fixed normative target, which is useful when reproducible z-scores are needed across newly added cohorts (Jodoin et al., 19 May 2025, Girard et al., 6 Nov 2025). Cluster ComBat and Fed-ComBat instead address decentralized settings by exchanging site-level sufficient statistics or model parameters rather than raw data, permitting harmonization of newly joined or unseen sites without full retraining (Hoang et al., 2024, Silva et al., 19 Jan 2026). Spatially structured images may benefit from Tensor-ComBat when scanner effects are spatially distributed and uncertainty quantification is required (Reinhardt et al., 2 Apr 2026).

Taken together, the literature presents ComBat as a flexible family of empirical Bayes harmonization methods rather than a single fixed algorithm. Its enduring value lies in an interpretable location–scale model that preserves designated biological effects, but its validity depends on covariate specification, demographic overlap, batch composition, leakage control, and the geometry of the feature space to which it is applied.