MM-ComBat and MM-CovBat: Multivariate Frameworks for Joint Harmonization of Multi-Metric Neuroimaging Data

Aggregating neuroimaging data across sites and studies is increasingly common, yet site- and scanner-related batch effects can obscure meaningful biological variation and introduce spurious associations. Although ComBat and its extensions are widely used, they are primarily designed for single-metric harmonization. In practice, neuroimaging studies often involve multiple biologically coupled metrics (e.g., cortical thickness, surface area, and gray-matter volume) measured across multiple features (e.g., regional values), with shared covariance structure both within and across metrics. Applying single-metric ComBat independently to each metric ignores these cross-metric dependencies. Using data from the NIH Acute to Chronic Pain Signatures (A2CPS) program, we show that batch effects occur not only in means and variances but also in covariance across cortical regions and metrics — relationships that single-metric ComBat does not fully remove. We propose MM-ComBat, a multivariate extension of ComBat that jointly harmonizes multiple metrics by borrowing strength across them, better capturing cross-metric dependence while also improving covariance estimation across features. Because joint harmonization whitens residual covariance toward a standardized baseline, it risks distorting biologically meaningful cross-metric structure when batch effects are moderate. We therefore introduce two complementary formulations: a baseline formulation suited to batch-dominated settings, and a target-covariance formulation that remaps adjusted covariances toward an estimated shared biological structure rather than fully whitening them. Both empirical Bayes (EB) and Bayesian Markov Chain Monte Carlo (MCMC) implementations of MM-ComBat effectively reduce batch effects. In our experiments, EB is more robust to measurement error, whereas MCMC more accurately recovers cross-metric correlations when priors are well specified. Recognizing that batch effects can also affect feature-level covariance, CovBat was recently introduced as an extension of ComBat that harmonizes both first- and second-order moments across sites. We extend CovBat to the multivariate framework as MM-CovBat, which performs a second-stage latent-space harmonization to directly address covariance-related batch effects across features and metrics. Simulations confirm that MM-ComBat improves correlation recovery and better preserves biological effects in the mean structure relative to single-metric ComBat, particularly for moderate-to-strong effects, and that MM-CovBat further improves separation of true biological variation from batch effects when independence assumptions are violated. Together, these methods provide a flexible and unified framework for harmonizing complex, multi-metric neuroimaging data in large-scale, multi-site studies.