Gene-Specific Batch Effects in Mixed Models

Why Batch Effects Need to Be Gene-Specific

A technical batch shift does not have to affect every gene equally. A processing difference might produce a large shift in one gene, a small shift in another, and almost none in a third. A model with one common batch effect across all genes cannot represent that pattern. Gene-specific batch terms allow the magnitude of batch variation to differ by gene, while cell-state terms help distinguish technical shifts from genuine differences in the mix of cells being measured.

Here, expression $y_{ig}$ denotes the normalized, typically log-scale expression value for cell $i$ and gene $g$. Let $b(i)$ identify the batch containing cell $i$, and let $s(i)$ identify its cell state. The model below includes state-specific average expression for each gene and a gene-specific random intercept for batch. Its batch effect is common across states for a given gene. It does not model a batch shift that changes with cell state; such a pattern requires a different model.

Definition: A gene-specific batch random-intercept model for gene $g$ is $$ y_{ig}=\beta_{0g}+\boldsymbol{x}_{s(i)}^\mathsf{T}\boldsymbol{\beta}_{g}+u_{g,b(i)}+\varepsilon_{ig}. $$ Here, $\beta_{0g}$ is the gene-specific intercept, $\boldsymbol{x}_{s(i)}$ encodes the cell state, and $\boldsymbol{\beta}_{g}$ contains the corresponding state effects. The batch term $u_{g,b}$ is a random effect for batch $b$ and gene $g$. The residual $\varepsilon_{ig}$ is cell-level variation not explained by the included terms.

A standard version of the model assumes $u_{g,b}\sim N(0,\tau_g^2)$ and $\varepsilon_{ig}\sim N(0,\sigma_g^2)$, with batch effects independent across batches, residuals independent across cells conditional on the model, and batch effects independent of residuals. In particular, $\tau_g^2$ is allowed to differ across genes. The zero-mean assumption makes the fixed intercept the expected expression for the reference state, averaged over the batch-effect distribution. It does not say that every observed batch has exactly zero effect.

These assumptions require care in single-cell data. If cells from the same biological sample have residual dependence beyond the batch effect, the independent-residual assumption is not adequate. Also, state effects and batch effects can only be separated when the observed design contains information to distinguish them—for example, when relevant states are represented in more than one batch. A model cannot recover a state comparison that the data do not support.

What the Gene-Specific Batch Term Means

For a particular gene, $u_{g,b}$ is the departure of batch $b$ from the gene’s mean after accounting for the modeled cell-state differences. The variance $\tau_g^2$ describes how much these batch departures vary across batches for that gene. It is not a measure of biological variability between states, and it does not imply that the same batch has the same effect on every gene.

The model is usually fitted gene by gene, or through a fitting procedure that retains a separate batch variance for each gene. With few batches, the variance $\tau_g^2$ may be estimated imprecisely. A small estimated variance is not proof that the gene has no technical sensitivity; it may reflect limited information. Conversely, a large estimate based on very few batches can be unstable.

For interpretation, keep the expression scale in view. If $y_{ig}$ is log base 2 normalized expression, a batch offset of $0.5$ is a half-unit shift on that scale. It corresponds to a multiplicative factor of $2^{0.5}$, about 1.414, on the underlying linear scale, subject to the interpretation of the normalization and transformation used.

Worked Example: Shrinking a Batch Offset for One Gene

Suppose a hypothetical gene is measured in three batches, with four cells per batch. After accounting for the fitted cell-state means, assume the average residuals are $1.2$, $-0.6$, and $-0.6$ log2 expression units. Suppose the fitted variance components are $\tau_g^2=0.64$ and $\sigma_g^2=1.00$. For this calculation, treat the fixed state means as known.

The residual mean in a batch estimates its random effect, but it is noisy. Under the normal model, the conditional mean of the batch effect given its residual mean is a shrunken version of that mean:

\widehat{u}_{g,b} = \frac{\tau_g^2}{\tau_g^2+\sigma_g^2/n_b}\,\bar{r}_{g,b},

where $n_b$ is the number of cells in the batch and $\bar{r}_{g,b}$ is its average residual after state adjustment. The shrinkage factor is

\frac{0.64}{0.64+1.00/4} = \frac{0.64}{0.89} \approx 0.7191.

Thus the three predicted batch effects are $0.7191(1.2)=0.863$, $0.7191(-0.6)=-0.431$, and $0.7191(-0.6)=-0.431$ log2 expression units. The unshrunk residual means sum to zero, and the shrunken effects also sum to approximately zero (up to rounding). These are model-based, shrunken predictions, not direct observations of fixed batch constants.

Why Shrinkage Occurs

The shrinkage formula follows from combining the random-effect distribution with the sampling uncertainty in a batch mean. Conditional on the fixed state means, write $\bar{r}_{g,b}=u_{g,b}+\bar{\varepsilon}_{g,b}$. The assumptions give $u_{g,b}\sim N(0,\tau_g^2)$ and $\bar{\varepsilon}_{g,b}\sim N(0,\sigma_g^2/n_b)$, independently. The conditional distribution of $u_{g,b}$ given $\bar{r}_{g,b}$ is normal. Its variance is the reciprocal of the sum of prior and data precisions, and its mean is that variance times the data precision multiplied by $\bar{r}_{g,b}$:

\operatorname{Var}(u_{g,b}\mid\bar{r}_{g,b}) = \left(\frac{1}{\tau_g^2}+\frac{n_b}{\sigma_g^2}\right)^{-1}, \qquad E(u_{g,b}\mid\bar{r}_{g,b}) = \frac{\tau_g^2}{\tau_g^2+\sigma_g^2/n_b}\bar{r}_{g,b}.

The conditional mean is the best predictor under squared-error loss, and is commonly called the best linear unbiased predictor (BLUP) when the variance parameters are treated as known. In practice, the variance parameters are estimated, so fitted predictions also reflect estimation uncertainty that this simple derivation omits. The formula nevertheless shows the key behavior: a noisy batch mean is pulled more strongly toward zero, while a batch mean supported by many cells is pulled less strongly.

Worked Example: Comparing Batch Variation Across Genes

Consider two hypothetical genes fitted with the same design. Their estimated batch variances are $\widehat{\tau}_A^2=0.25$ and $\widehat{\tau}_B^2=1.44$, in squared log2 expression units. The corresponding estimated standard deviations of batch effects are

\widehat{\tau}_A=\sqrt{0.25}=0.500, \qquad \widehat{\tau}_B=\sqrt{1.44}=1.20.

Under the fitted model, gene B has greater between-batch variability in its adjusted expression offset: its estimated batch-effect standard deviation is $1.20$ log2 units, compared with $0.500$ for gene A. The estimated variance for B is $1.44/0.25=5.76$ times that for A. This comparison concerns the spread of batch effects, not the average expression level of the genes, and it does not by itself establish that a particular batch shifts gene B upward or downward.

The comparison is meaningful only if the models use comparable expression scales and account for relevant cell-state composition. A larger variance can otherwise reflect differences in state representation or an inadequate mean model rather than a purely technical difference.

Cell-State Composition and Identifiability

State adjustment matters because batches can contain different proportions of cell states. If expression differs strongly by state, a batch with more cells of a high-expression state can have a higher raw average even when there is no within-state batch shift. The following example uses state means that are identical across batches to show how composition alone can create an apparent difference.

Worked Example: A Raw Batch Difference Caused by State Composition

Suppose a hypothetical gene has mean expression $2$ log2 units in state A and $6$ log2 units in state B, with no batch effect within either state. Batch 1 contains eight state-A cells and two state-B cells; batch 2 contains two state-A cells and eight state-B cells.

The raw mean for batch 1 is

\frac{8(2)+2(6)}{10} = \frac{28}{10} = 2.8.

The raw mean for batch 2 is

\frac{2(2)+8(6)}{10} = \frac{52}{10} = 5.2.

The unadjusted difference is $5.2-2.8=2.4$ log2 units, despite the stipulated absence of a batch shift within either state. After including state in the model, the state-A means are both $2$, and the state-B means are both $6$. The adjusted batch effect is therefore $0$ in this hypothetical example. This conclusion relies on having both states represented in both batches; if a state occurs in only one batch, its state effect and that batch’s effect may be inseparable.

When state and batch are perfectly confounded, the design matrix does not contain enough independent information to estimate their separate effects. For example, if all state-A cells occur in batch 1 and all state-B cells occur in batch 2, an observed difference could be due to state, batch, or both. Adding a random-effects distribution does not create information absent from the design. Estimates in such cases depend heavily on modeling assumptions and should not be described as empirically separated technical and biological effects.

Model Checks, Alternatives, and Limits

A fitted variance component should be assessed alongside the design and the residual pattern. Useful checks include examining whether each state is represented across batches, whether the number of batches is adequate to estimate a batch variance, and whether state-adjusted residuals still show systematic batch differences. A pattern in which a gene’s batch residuals differ in direction or size across states signals that the common-within-gene batch offset may be inadequate. This model does not resolve that pattern by allowing state-specific batch slopes.

Other warning signs include a variance estimate at or near zero, a very large estimate driven by one batch, unstable estimates across reasonable fitting choices, and residual dependence among cells from the same biological sample. A large cell count cannot compensate for having very few independent batches when the target is between-batch variability. Cell-level precision and the number of independent batch units are different sources of information.

The principal alternatives answer different questions:

Approach	What it allows	Main trade-off
Fixed batch terms	A separate unconstrained mean offset for each observed batch and gene	Flexible for the observed batches, but does not model a distribution of batch effects for generalization to new batches
Gene-specific random intercepts	Batch offsets that vary across batches, with a separate variance for each gene	Provides partial pooling, but relies on the random-effect distribution and can be weakly informed with few batches
One pooled batch variance	Batch variation governed by a common variance across genes	More stable with limited information, but may conceal strong gene-to-gene differences

Fixed batch effects and random batch effects should not be treated as interchangeable labels for the same analysis. Fixed terms estimate offsets for the batches observed; random terms assume those batches are draws from a distribution and estimate its variability. Choosing between them depends on the inferential target, the number and nature of batches, and the intended use of the fitted values.

A gene-specific variance is also not a guarantee of accurate correction. The model can still be misspecified, poorly identified, or influenced by unmodeled sample structure. Nor does a fitted batch association establish that the batch caused the expression difference. The measured offsets describe conditional patterns under the model; attributing them to a technical mechanism requires appropriate study design and supporting evidence.

Key takeaway: A gene-specific batch random intercept lets the amount of between-batch variation differ by gene, while state terms account for differences in cell-state composition. Its interpretation depends on adequate state-batch overlap, enough independent batches, and a credible residual and random-effect model.

Check Your Understanding

Answer each question using the model assumptions and interpretations stated above.

If $\widehat{\tau}_g^2=0.81$, what is the estimated standard deviation of batch effects, and what are its units?
For a batch with $\bar{r}_{g,b}=0.8$, $n_b=4$, $\tau_g^2=0.36$, and $\sigma_g^2=1$, calculate the shrunken prediction.
Why can different cell-state proportions create a raw batch difference even when there is no within-state batch effect?
What is the main interpretation difference between fixed batch terms and random batch terms?
If each cell state appears in only one batch, can state and batch effects be cleanly separated from these data? Explain.

Solutions. (1) The estimated standard deviation is $\sqrt{0.81}=0.900$, in the same units as expression, such as log2 expression units. (2) The shrinkage factor is $0.36/(0.36+1/4)=0.36/0.61\approx0.5902$; the predicted effect is $0.5902(0.8)\approx0.472$ expression units. (3) A batch with more cells from a state whose mean expression is higher will have a higher weighted raw average, even if the means within each state are identical across batches. (4) Fixed terms estimate offsets for the observed batches without treating them as draws from a batch-effect distribution; random terms model offsets as draws from a distribution and estimate its variance. (5) No. State and batch are confounded, so the observed differences cannot uniquely identify separate state and batch effects without additional information or assumptions.