Mixed-Effects Models for Batch Integration

Why Model Batch and Biological State Together?

Cells from different batches can differ for technical reasons, such as processing conditions, as well as for biological reasons, such as differing proportions of cell states. If expression is compared across batches without accounting for state, a genuine difference in cell composition can look like a technical shift. Conversely, if a technical shift is left untreated, cells in the same biological state may appear farther apart than they should. A mixed-effects model provides a way to represent both sources of variation in one model.

Here the response is a transformed expression measurement for a gene in a cell, not a raw count. The model is useful only to the extent that this response and its assumptions are suitable for the intended analysis. It does not itself determine which cells are biologically equivalent, nor can it recover a biological contrast that the design does not distinguish from batch.

Definition: For gene $g$ in cell $i$, let $y_{ig}$ denote a transformed expression value, $s_i$ the cell's biological state, and $j(i)$ its batch. A random-intercept mixed-effects model is $$ y_{ig}=\mu_g+\beta_{g,s_i}+b_{j(i),g}+\varepsilon_{ig}, $$ where $\mu_g$ is a reference-state mean, $\beta_{g,s_i}$ is the state effect (zero for the reference state), $b_{j(i),g}$ is the batch effect, and $\varepsilon_{ig}$ is cell-level residual variation. In the usual Gaussian formulation, $b_{jg}\sim N(0,\tau_g^2)$ and $\varepsilon_{ig}\sim N(0,\sigma_g^2)$, with batch effects independent across batches and residuals independent conditional on the modeled terms.

The state terms are fixed effects: they describe the biological states whose differences are of direct interest or are to be preserved. The batch terms are random effects: batches are treated as draws from a population of possible processing batches, with a common mean of zero and variance $\tau_g^2$. A fitted model estimates the state effects and variance components, then obtains a batch-effect estimate for each observed batch. The model is written gene by gene here to make the response explicit; this introduction does not assume that every gene has the same batch shift.

Assumptions and Interpretation

The Gaussian model is a working model for transformed expression, often a log-scale normalized measurement. Its conditional mean should be adequately represented by the included state and batch terms, and residual variation should be reasonably described by the chosen distribution. A technical batch shift is represented as additive on this response scale. The normality assumption concerns the model terms and residuals, not the original count measurements.

Conditional independence of residuals is another assumption. Cells from the same batch share a random effect, so their expression measurements are not marginally independent under the model. Additional dependence, such as cells clustered by donor, may require additional structure; omitting important clustering can make uncertainty appear smaller than it is. In particular, a large number of cells from only a few batches does not provide the same information about batch-to-batch variability as a large number of independent batches.

The design must also distinguish state from batch. If every state occurs in multiple batches, the state contrasts can be estimated separately from batch shifts more credibly. If a state occurs in only one batch, its state effect and that batch's effect may be inseparable. State labels must be sufficiently accurate: misclassification can leave state differences in the apparent batch effects or remove real state structure.

Worked Example: Separating a State Difference From a Batch Shift

Suppose a hypothetical gene has the following mean transformed expression values in two batches. Both batches contain cells in a reference state and in a second state.

Batch	Reference-state mean	Second-state mean
A	4.8	6.8
B	5.4	7.4

The within-batch difference between the second state and the reference state is $6.8-4.8=2.0$ in batch A and $7.4-5.4=2.0$ in batch B. Thus the data support a state difference of 2.0 expression units, rather than a difference caused solely by unequal state composition. Both means in batch B are 0.6 units higher than the corresponding means in batch A, consistent with an additive batch shift.

With equally weighted observations in this simplified example, the pooled reference-state mean is $(4.8+5.4)/2=5.1$. The corresponding batch deviations are $4.8-5.1=-0.3$ for A and $5.4-5.1=0.3$ for B. The fitted additive description is therefore a reference mean of 5.1, a second-state effect of 2.0, and batch deviations of -0.3 and 0.3. In an actual mixed-effects fit, estimated batch effects may be shrunk toward zero, so these conditional mean deviations need not equal the final random-effect estimates.

This calculation relies on both states being represented in both batches and on comparable measurement scales. It illustrates the model's separation of terms; it does not prove that the entire observed batch difference is technical or that all genes have identical shifts.

Estimating and Applying Batch Effects

For a given gene, the estimated batch effect $\hat b_{jg}$ represents the model's estimate of the batch-specific deviation after accounting for biological state. A simple conditional estimate begins with residuals after subtracting the fitted state mean. A mixed-effects estimate also accounts for the fact that batch effects are modeled as draws from a distribution. This typically shrinks a noisy batch estimate toward the population mean of zero.

For a cell in batch $j$ and state $s$, a basic adjusted value is $y_{ig}-\hat b_{jg}$. Under the additive model, this removes the estimated batch contribution while retaining the fitted state contribution. It is an alignment on the chosen expression scale, not a claim that the adjusted value is the cell's error-free biological expression. Uncertainty in $\hat b_{jg}$, model misspecification, and genuine differences between batches may remain.

Why Shrinkage Occurs

The shrinkage can be derived in a simplified setting. Suppose the state terms have already been accounted for, and the mean residual in batch $j$, $\bar r_j$, estimates its batch effect with sampling variance $\sigma_g^2/n_j$. Under the random-effects model, $b_{jg}\sim N(0,\tau_g^2)$, and the mean residual is $b_{jg}$ plus sampling error. The conditional mean of the batch effect given this mean residual is

E(b_{jg}\mid \bar r_j) = \frac{\tau_g^2}{\tau_g^2+\sigma_g^2/n_j}\bar r_j.

To see this, the covariance of $b_{jg}$ with $\bar r_j$ is $\tau_g^2$, while the variance of $\bar r_j$ is $\tau_g^2+\sigma_g^2/n_j$. For jointly normal variables, the conditional mean is the covariance divided by the variance of the conditioning variable, multiplied by its observed value. The multiplier is between zero and one: a small batch, with a noisy mean residual, receives stronger shrinkage; a precise batch mean receives less. In practice, the variance components are estimated from the data, and the fitted model may be more complex than this derivation.

Worked Example: Calculating a Shrunk Batch Estimate

For a hypothetical gene, suppose the state-adjusted residual mean in one batch is $\bar r_j=0.30$ expression units, based on $n_j=16$ cells. Assume the fitted variance components are $\tau_g^2=0.16$ and $\sigma_g^2=0.64$, in squared expression units. Using the conditional-mean formula:

\hat b_{jg} = \frac{0.16}{0.16+0.64/16}(0.30) = \frac{0.16}{0.20}(0.30) = 0.80(0.30) = 0.24.

The estimated batch shift is $0.24$ expression units, rather than the raw residual mean of $0.30$. If a cell in that batch has observed expression $5.30$ and its fitted biological-state mean is $5.00$, subtracting the batch estimate gives $5.30-0.24=5.06$. This adjusted value remains 0.06 units above the state mean; the adjustment removes the estimated batch contribution, not all residual variation.

Identifiability, Diagnostics, and Alternatives

A model can return estimates even when its assumptions or design make those estimates unreliable. Begin by tabulating cells by state and batch. Look for absent state-batch combinations, very small batches, and state labels that are unevenly supported. Then inspect residuals after fitting: persistent shifts by batch suggest that the additive batch term is insufficient, while state-specific patterns across batches can indicate that a single shift does not describe the data. Strong outliers or markedly unequal residual spread also challenge the Gaussian working model. These checks identify warning signs; they do not by themselves establish a corrected model.

Worked Example: Recognizing a Confounded Design

Consider a hypothetical experiment in which batch A contains only reference-state cells with mean expression 5.0, while batch B contains only second-state cells with mean expression 8.0. The observed difference is $8.0-5.0=3.0$ units. There are no observations that compare the two states within either batch, and no observations that compare the batches within either state.

Under the model, the difference could be explained by a state effect of 3.0 and no batch shift, by a batch shift of 3.0 and no state effect, or by many combinations whose sum is 3.0. The data alone cannot decide among these explanations. Fitting software may impose constraints or return a numerical answer, but that does not supply the missing information. Cross-batch adjustment in this design would depend on untestable assumptions, so the estimated alignment should not be interpreted as a data-identified technical correction.

Mixed-effects adjustment is not the only integration strategy. Simple batch centering is transparent but does not, by itself, account for differences in biological-state composition; it may shift genuine biology if batches contain different cell mixtures. MNN and anchor-based approaches use relationships between cells or local neighborhoods to align shared states. They can be useful when suitable shared populations and reliable correspondences exist, but their results depend on those correspondences. A mixed-effects model instead makes the state and batch terms explicit and produces interpretable batch estimates under its assumptions. None of these approaches can safely align a state that is completely confounded with batch without additional information.

Common pitfall: “The average expression differs between batches, so subtracting each batch mean will remove only technical variation.” This reasoning ignores state composition. If one batch has more cells in a high-expression state, its overall mean can be higher for biological reasons. Estimate batch differences conditional on state, check that state and batch are distinguishable in the design, and treat the correction as model-dependent rather than as proof that all between-batch differences were technical.

The model is most useful when there are enough batches to estimate between-batch variation, relevant states are represented across batches, and a Gaussian additive description is reasonable for the transformed response. It is less trustworthy with severe confounding, very few batches, poor state labels, or systematic batch-specific patterns remaining after adjustment. For inferential conclusions, uncertainty and the experimental unit matter: many cells from one batch do not turn that batch into many independent batches.

Check Your Understanding

Answer each question, then compare with the solutions.

In the model, what does the batch effect represent after biological state is included?
Using $\tau_g^2=0.09$, $\sigma_g^2=0.36$, $n_j=9$, and $\bar r_j=0.20$, calculate the shrunk batch estimate.
Why is it not possible to separate state from batch when each state occurs in only one distinct batch?
Name one diagnostic that can reveal a poor additive batch model and one limitation of simple batch centering.

Solutions. (1) It is the estimated batch-specific deviation in expression after accounting for the modeled state effect. (2) The multiplier is $0.09/(0.09+0.36/9)=0.09/0.13\approx0.6923$, so $\hat b_{jg}\approx0.6923(0.20)=0.1385$ expression units. (3) The observed difference can be attributed to either state or batch, because the design contains no within-batch state comparison or within-state batch comparison. (4) Residual patterns that persist by batch or state can signal a poor additive model; simple centering can remove genuine biological differences when state composition varies across batches.