ComBat Empirical-Bayes Adjustment

Why ComBat Uses Information Across Features

A batch can change measurements in more than one way. It may shift the typical value of a feature, alter its variability, or do both. A feature-by-feature correction can estimate such changes, but estimates from small batches may be noisy. ComBat addresses this by estimating batch effects for each feature and shrinking those estimates toward batch-level distributions learned across features. The resulting correction is still feature-specific; the sharing occurs through the estimation of the batch-effect distributions.

The method requires a batch label for each sample and a design describing biological covariates whose fitted effects should be retained. It is intended for continuous measurements on features measured in a sufficiently comparable way across batches. Its adjustment is conditional on the supplied model and design: empirical-Bayes shrinkage does not solve a design that cannot distinguish batch from biology.

Definition: Let $Y_{gj}$ be the measurement for feature $g$ in sample $j$, and let $k(j)$ denote that sample’s batch, with $k=1,\ldots,K$. Let $m_{gj}=\alpha_g+\mathbf{x}_j^T\boldsymbol{\beta}_g$ be the fitted biological mean for feature $g$, where $\mathbf{x}_j$ contains the specified covariates and excludes an intercept. ComBat’s location-scale model is $$ Y_{gj}=m_{gj}+\sigma_g\left[\gamma_{g,k(j)}+\delta_{g,k(j)}\epsilon_{gj}\right],\qquad \epsilon_{gj}\sim N(0,1). $$ Here, $\sigma_g>0$ is a feature-level residual scale, $\gamma_{gk}$ is the standardized location effect for feature $g$ in batch $k$, and $\delta_{gk}>0$ is its standardized scale effect. The factor $\sigma_g$ multiplies both batch terms.

This placement of $\sigma_g$ matters. Once the biological mean is removed and the residual is divided by $\sigma_g$, the standardized value is

Z_{gj}=\frac{Y_{gj}-m_{gj}}{\sigma_g} =\gamma_{g,k(j)}+\delta_{g,k(j)}\epsilon_{gj}.

Thus the batch location and scale parameters are both on the standardized residual scale. Subtracting a location parameter defined on the original measurement scale, or dividing by $\sigma_g$ again after estimating a standardized location effect, would be inconsistent with this model.

Inputs and Conditions for Adjustment

A typical input is a feature-by-sample matrix, a batch label for every sample, and a design matrix of covariates to preserve. ComBat estimates a feature’s covariate mean and residual scale, standardizes its residuals, estimates batch effects, then returns adjusted measurements on the original feature scale. Missing values, features with zero or near-zero variation, and features not measured comparably across batches require careful handling; they cannot simply be treated as ordinary complete measurements in this model.

The central assumptions are that the specified biological mean is adequate, that standardized residuals within a feature and batch are reasonably represented by the location-scale model, and that the residual errors are independent across samples conditional on the model. The empirical-Bayes layer further assumes that, within each batch, feature-specific location effects share a distribution and feature-specific scale effects share a distribution. The feature effects need not be identical: the distributions provide a common basis for shrinkage.

The design must contain enough information to estimate biological covariate effects separately from batch. In practice, biological groups should be represented across batches sufficiently for the intended covariate effects to be estimable. The earlier material on rank and connectivity checks gives tools for assessing that support. A numerical output from ComBat does not establish that the design can separate biology from batch.

How Empirical-Bayes Shrinkage Works

After standardization, ComBat estimates a location and scale effect for each feature in each batch. It then estimates batch-specific hyperparameters from the collection of feature estimates. A common parametric specification is

\gamma_{gk}\sim N(\bar{\gamma}_k,\tau_k^2), \qquad \delta_{gk}^2\sim \operatorname{InvGamma}(a_k,b_k).

The first distribution describes the spread of standardized location effects across features in batch $k$. The second describes the spread of squared scale effects. Here, $\operatorname{InvGamma}(a,b)$ has density proportional to $u^{-(a+1)}\exp(-b/u)$ for $u>0$. Its mean is $b/(a-1)$ when $a>1$. The hyperparameters are estimated from the feature-wise effects, rather than chosen as known constants.

For a fixed scale effect $\delta_{gk}$, suppose the standardized observations in batch $k$ for one feature have mean $\bar{Z}_{gk}$ and count $n_k$. Their likelihood for $\gamma_{gk}$ is normal with mean $\gamma_{gk}$ and variance $\delta_{gk}^2/n_k$. Combining that likelihood with the normal prior gives the posterior mean below.

Justification: The likelihood precision is $n_k/\delta_{gk}^2$, and the prior precision is $1/\tau_k^2$. Normal-normal updating adds precisions, so $$ \widehat{\gamma}^{\,*}_{gk} = \frac{(n_k/\delta_{gk}^2)\bar{Z}_{gk}+(1/\tau_k^2)\bar{\gamma}_k} {n_k/\delta_{gk}^2+1/\tau_k^2} = \frac{n_k\tau_k^2\bar{Z}_{gk}+\delta_{gk}^2\bar{\gamma}_k} {n_k\tau_k^2+\delta_{gk}^2}. $$ This is a weighted average of the feature’s batch mean and the across-feature prior mean. The feature mean receives more weight when its batch has more samples or less residual variation; the prior mean receives more weight when the feature estimate is noisy.

The scale estimate is also shrunk using the batch-specific inverse-gamma distribution. For example, if $\delta_{gk}^2\sim\operatorname{InvGamma}(a_k,b_k)$ and the standardized observations $Z_{gj}$ in that batch are conditionally normal with mean $\gamma_{gk}$, then, conditional on $\gamma_{gk}$, the posterior for $\delta_{gk}^2$ has shape $a_k+n_k/2$ and scale $b_k+\frac12\sum_j(Z_{gj}-\gamma_{gk})^2$. Its posterior mean, when the shape exceeds 1, is the scale divided by the shape minus 1. Practical empirical-Bayes procedures estimate the hyperparameters from the features and use iterative or plug-in estimates of the feature effects. Details can vary by implementation, but the aim is the same: stabilize feature-wise effects using the across-feature distributions.

Applying the Adjustment

Let $\widehat{m}_{gj}$ and $\widehat{\sigma}_g$ be the estimated biological mean and feature scale, and let $\widehat{\gamma}^{\,*}_{gk}$ and $\widehat{\delta}^{\,*}_{gk}$ be the empirical-Bayes batch estimates. A ComBat-style adjustment removes the estimated location effect and rescales by the estimated scale effect, then adds back the fitted biological mean:

Y^{\mathrm{adj}}_{gj} = \widehat{m}_{gj} + \widehat{\sigma}_g \left( \frac{Z_{gj}-\widehat{\gamma}^{\,*}_{g,k(j)}} {\widehat{\delta}^{\,*}_{g,k(j)}} \right), \qquad Z_{gj}=\frac{Y_{gj}-\widehat{m}_{gj}}{\widehat{\sigma}_g}.

The fitted mean is added back so that covariate effects included in the model are retained in the adjusted data. If no biological covariates are specified beyond an intercept, the adjustment retains the feature-level fitted mean. “Retained” means retained according to the specified model; it does not guarantee that omitted biology is protected from adjustment.

Worked Example: Adjusting One Measurement

Suppose a hypothetical feature has fitted biological mean $\widehat{m}_{gj}=10$ units and estimated scale $\widehat{\sigma}_g=2$ units. For a sample in batch 2, let the observed value be $Y_{gj}=17.6$ units, with empirical-Bayes estimates $\widehat{\gamma}^{\,*}_{g2}=0.8$ and $\widehat{\delta}^{\,*}_{g2}=1.5$. These values are assumed to have been estimated from a design that supports the required covariate and batch effects.

First standardize the residual:

Z_{gj}=\frac{17.6-10}{2}=\frac{7.6}{2}=3.8.

Then remove the estimated standardized location and scale effects:

\frac{Z_{gj}-\widehat{\gamma}^{\,*}_{g2}} {\widehat{\delta}^{\,*}_{g2}} = \frac{3.8-0.8}{1.5} = \frac{3.0}{1.5} =2.

Finally, return to the original scale and add back the fitted mean:

Y^{\mathrm{adj}}_{gj}=10+2(2)=14\text{ units}.

Under the specified model, the standardized residual before correction corresponds to $\gamma+\delta\epsilon=0.8+1.5(2)=3.8$, matching the computed $Z_{gj}$. This check confirms that the location, scale, and feature-scale terms have been applied consistently.

Worked Example: Shrinking a Location Estimate

Consider a hypothetical feature in a batch with four samples whose standardized residuals are $1.0,\ 1.4,\ 0.8,\ 1.2$. Assume for this calculation that the scale variance is $\delta_{gk}^2=1$, and that across-feature estimates give a prior mean $\bar{\gamma}_k=0.2$ and prior variance $\tau_k^2=0.25$. The sample mean is

\bar{Z}_{gk}=\frac{1.0+1.4+0.8+1.2}{4}=\frac{4.4}{4}=1.1.

Using the posterior-mean formula,

\widehat{\gamma}^{\,*}_{gk} = \frac{4(0.25)(1.1)+(1)(0.2)} {4(0.25)+1} = \frac{1.1+0.2}{2} =0.65.

The raw batch mean of $1.1$ is pulled toward the across-feature mean of $0.2$, giving an estimated location effect of $0.65$. With $\widehat{\delta}^{\,*}_{gk}=1$, the adjusted standardized residuals are $0.35,\ 0.75,\ 0.15,\ 0.55$, each obtained by subtracting $0.65$. Their mean is $0.45$, not zero: shrinkage does not force every batch-feature mean to zero. It trades some residual batch offset for less variable estimates when the batch contains limited information.

Estimating the Across-Feature Scale Distribution

The parametric empirical-Bayes approach estimates hyperparameters by summarizing feature-wise estimates within each batch. For squared scale estimates $u_g=\widehat{\delta}_{gk}^2$, a moment-based estimate can match their mean and variance to the inverse-gamma distribution. If its mean is $m$ and variance is $v$, then

m=\frac{b}{a-1},\qquad v=\frac{b^2}{(a-1)^2(a-2)},\qquad a=2+\frac{m^2}{v},\qquad b=m(a-1).

These relationships require $a>2$ for a finite variance. They also illustrate why a very small estimated across-feature variance implies a concentrated prior. Moment estimates are themselves uncertain, especially when few features contribute, so the chosen prior family and the features included in its estimation matter.

Worked Example: Moment Estimates for Scale Shrinkage

Suppose, hypothetically, that four features in one batch have squared scale estimates $0.8,\ 1.0,\ 1.2,\ 1.0$. Their mean is $m=(0.8+1.0+1.2+1.0)/4=1.0$. Using the sample variance, the squared deviations sum to $0.04+0+0.04+0=0.08$, so $v=0.08/(4-1)=0.02667$.

The moment estimates for the inverse-gamma parameters are

\widehat{a}=2+\frac{1.0^2}{0.02667}\approx39.50, \qquad \widehat{b}=1.0(39.50-1)\approx38.50.

The implied prior mean is $38.50/(39.50-1)=1.00$, as required by the moment match. This is a distribution for squared scale effects, not for $\delta_{gk}$ itself. In a real application, these estimates would be based on many features, and their adequacy would need to be assessed rather than assumed from a small illustrative set.

Diagnostics, Alternatives, and Limitations

Before adjustment, check the batch counts, covariate representation across batches, feature scales, and missingness. After fitting, inspect standardized residual distributions and compare batch patterns before and after adjustment, both overall and within biological groups. A large number of extreme estimates, nearly zero scales, strong non-normality, or substantial residual batch structure can signal that the model or its inputs are unsuitable. Such checks diagnose potential problems; they do not prove that the corrected values are free of batch effects.

ComBat is useful when a shared, feature-wise location-scale representation is reasonable and there are many features from which to estimate batch-level hyperparameters. It differs from estimating and subtracting a separate batch coefficient for every feature without shrinkage: the latter can be unstable in small batches, while shrinkage can introduce bias if the across-feature distributions are poorly represented. A mixed-effects model offers another approach when batch is treated as a random effect and the modeling goal includes that structure; direct inclusion of batch in a downstream regression can instead retain the batch term rather than create a corrected matrix. These approaches answer different questions, and none removes the need for an estimable design and a defensible model.

ComBat’s output should not be read as proof that biological signal has been preserved or that technical variation has been completely removed. Results depend on which covariates are specified, the batch definitions, the parametric assumptions, and the quality of overlap in the design. A correction is an estimated transformation, not a substitute for the design information needed to distinguish biological from technical differences.

Wrong reasoning: “The data are standardized, so subtracting $\gamma$ and dividing by $\delta$ directly gives adjusted measurements.” Standardized values are residuals after removing the fitted biological mean and dividing by the feature scale.
Correction: Standardize as $Z=(Y-\widehat{m})/\widehat{\sigma}$, adjust $Z$, multiply by $\widehat{\sigma}$, and add back $\widehat{m}$.
Wrong reasoning: “Empirical Bayes makes a confounded batch design identifiable.” Shrinkage stabilizes estimates under its distributional assumptions; it cannot supply missing information to separate batch from biology.
Correction: Check design support first, and interpret the adjustment only for effects estimable under that design.

Key takeaway: ComBat standardizes feature residuals, estimates batch-specific location and scale effects, shrinks them using across-feature distributions, and then reconstructs adjusted measurements while adding back the specified biological mean. The location-scale model, design, and prior assumptions must all be coherent for that adjustment to be interpretable.

Check Your Understanding

Answer each question, then compare your reasoning with the brief solutions.

If $Y=15$, $\widehat{m}=9$, and $\widehat{\sigma}=3$, what is the standardized residual?
Why can the posterior location estimate differ from the raw batch mean?
In the ComBat model, what does the parameter $\delta_{gk}$ scale?
What design condition is needed to distinguish a batch effect from a biological covariate effect?
Why does the adjusted standardized batch mean in the shrinkage example not have to equal zero?

Solutions. (1) $Z=(15-9)/3=2$. (2) The posterior estimate combines the feature’s batch mean with the across-feature prior mean, with weights determined by their relative precision. (3) It scales the standardized residual error within feature $g$ and batch $k$, after the standardized location effect is accounted for. (4) The design must contain enough cross-batch representation to estimate the biological contrast separately from batch; perfect confounding does not provide that information. (5) Shrinkage pulls the estimated location toward the prior mean instead of forcing it to equal the raw batch mean, so subtracting the shrunk estimate need not center the observed batch residuals at zero.