A Worked Mixed-Model Integration

From a Model to an Integration Adjustment

A mixed-effects integration is not just a rule for subtracting batch means. It encodes assumptions about which expression differences are technical batch shifts, which differences are associated with cell state, and how much information each batch provides about its own shift. Working through a small example makes these choices visible: we can fit a model, estimate its variance components, calculate batch adjustments, and verify that the state-associated contrast remains in the adjusted values.

For clarity, the example uses one gene’s normalized, log-scale expression in a balanced design with three batches and two states. The numbers are hypothetical and are intended to show the calculations, not to represent a real study. A separate model would generally be fitted for each gene, with careful attention to how expression was normalized and which cells or samples are appropriate observational units.

Model: For cell $i$ in batch $j$, let $y_{ij}$ be log-scale expression and let $S_{ij}$ equal $-1/2$ for State A and $+1/2$ for State B. The random-intercept model is $$ y_{ij}=\beta_0+\beta_S S_{ij}+b_j+\epsilon_{ij}. $$ Here $\beta_0$ is the overall expression level, $\beta_S$ is the State B minus State A contrast, $b_j$ is the batch-specific random intercept, and $\epsilon_{ij}$ is the residual. The usual model assumes $b_j\sim N(0,\tau^2)$, $\epsilon_{ij}\sim N(0,\sigma^2)$, and independence between random effects and residuals, with residuals independent across cells conditional on the batch effect. The state term is fixed; the batch intercepts are modeled as draws from a common distribution.

The centered coding of state makes $\beta_0$ the overall midpoint of the two state means, while $\beta_S$ is their difference. A random intercept shifts both states in a batch by the same amount. Thus the model can represent a shared batch offset, but it cannot represent a batch effect that changes with state. The earlier material on random slopes covers that distinct model.

For this model to distinguish state-associated expression from batch, both states need to be represented across batches, or the design needs other information that identifies their effects. If State A occurs only in one batch and State B only in another, state and batch are confounded: their contributions cannot be separately determined from the expression values alone. The example below has both states in every batch.

A Small Hypothetical Data Set

Each batch contains two observations in each state. The entries are individual expression values; units are log-expression units. Within each state and batch, the two observations are one unit below and one unit above their cell-group mean.

Batch	State A values	State B values
1	8.500, 10.500	9.833, 11.833
2	10.000, 12.000	12.333, 14.333
3	5.500, 7.500	7.833, 9.833

Worked Example: Estimate the State Contrast and Batch Residual Means

First compute the mean expression for each batch-state group. The means are $9.500$ and $10.833$ in Batch 1, $11.000$ and $13.333$ in Batch 2, and $6.500$ and $8.833$ in Batch 3, for States A and B respectively.

The state contrasts are therefore $10.833-9.500=1.333$, $13.333-11.000=2.333$, and $8.833-6.500=2.333$ log-expression units. Because the design is balanced, the fixed state estimate is their average:

\hat{\beta}_S=\frac{1.333+2.333+2.333}{3}=2.000.

The midpoint of the two state means is $10.167$ in Batch 1, $12.167$ in Batch 2, and $7.667$ in Batch 3. The overall midpoint is their average, $10.000$, so $\hat{\beta}_0=10.000$. Subtracting this fitted overall midpoint from each batch midpoint gives the residual batch means $r_1=0.167$, $r_2=2.167$, and $r_3=-2.333$ log-expression units. Their sum is zero, as expected with an intercept.

These residual batch means are not yet the mixed-model adjustments. They are the observed batch deviations after fitting the common intercept and state contrast. A mixed model shrinks them toward zero because each is estimated with uncertainty and the batch effects are modeled as draws from a distribution centered at zero.

Estimating the Variance Components

The amount of shrinkage depends on the estimated batch variance $\tau^2$ relative to the residual variance $\sigma^2$. In this balanced example, the sums of squares separate into replicate variation within batch-state groups, variation in batch-specific state contrasts beyond the common state effect, and variation among batch midpoints. This separation also provides a transparent REML calculation.

Within each of the six batch-state groups, the two observations are one unit below and one unit above their group mean. Their sum of squared deviations is $1^2+(-1)^2=2$. Across six groups, that component is $6\times2=12$, with $6$ degrees of freedom.

The observed batch-specific state contrasts differ from the common contrast $2.000$ by $-2/3,\;1/3,\;1/3$. With state coded as $-1/2,+1/2$, the residuals around the fitted state effect are $-q_j/2$ and $+q_j/2$ for the two states in batch $j$, where $q_j$ is that batch’s contrast deviation. Each state has two observations, so the sum of squared residuals for batch $j$ is $q_j^2$. Thus this component is

\left(-\frac{2}{3}\right)^2+ \left(\frac{1}{3}\right)^2+ \left(\frac{1}{3}\right)^2 =\frac{2}{3},

with $2$ degrees of freedom. Together, the within-batch residual sum of squares is $12+2/3=38/3$, with $8$ degrees of freedom. The REML estimate of the residual variance is therefore

\hat{\sigma}^2=\frac{12+2/3}{8}=\frac{19}{12}=1.583.

The three residual batch means are $1/6,\;13/6,\;-7/3$. Their squared deviations sum to

\left(\frac{1}{6}\right)^2+ \left(\frac{13}{6}\right)^2+ \left(-\frac{7}{3}\right)^2 =\frac{61}{6}.

There are $2$ degrees of freedom among these three means. A batch mean averages four observations, so its residual-error variance is $\sigma^2/4$, while its batch-effect variance is $\tau^2$. The variance of a batch mean is consequently $\tau^2+\sigma^2/4$. Equating the observed mean square to this quantity gives the REML component estimate:

\hat{\tau}^2+\frac{\hat{\sigma}^2}{4} =\frac{61/6}{2}=\frac{61}{12}, \qquad \hat{\tau}^2=\frac{61}{12}-\frac{19/12}{4} =\frac{75}{16}=4.688.

Why the variance calculation works here: In this balanced design, the batch-midpoint variation is orthogonal to the within-batch residual variation. The within-batch mean square estimates $\sigma^2$, and the batch-midpoint mean square estimates $\tau^2+\sigma^2/4$. Subtracting the residual contribution yields the batch variance estimate. This simple decomposition is specific to the balanced setup; general mixed models typically estimate variance components numerically using REML rather than this hand calculation.

Worked Example: Check the REML Arithmetic

A frequent source of error is to report the batch-midpoint mean square, $61/12=5.083$, as though it were the batch variance. It includes residual variation in each four-cell batch mean. Here that residual contribution is $\hat{\sigma}^2/4=(19/12)/4=19/48=0.396$. Subtracting it gives

5.083-0.396=4.688=\hat{\tau}^2.

The calculation uses the stated data and model throughout. It does not treat the variance estimates as externally stipulated values: $\hat{\sigma}^2=19/12$ and $\hat{\tau}^2=75/16$ follow from the sums of squares shown above.

From Batch Residuals to Adjustments

For known variance components and Gaussian random effects, the best linear unbiased predictor (BLUP) of a batch effect is the conditional expectation of that effect given its observed residual batch mean. In this balanced example, write $r_j=b_j+\bar{\epsilon}_j$, where $\bar{\epsilon}_j$ is the average residual for the four observations in batch $j$. Then $\operatorname{Var}(r_j)=\tau^2+\sigma^2/4$ and $\operatorname{Cov}(b_j,r_j)=\tau^2$. For jointly normal variables, the conditional expectation is the covariance divided by the variance, multiplied by the observed value. Therefore,

E(b_j\mid r_j) =\frac{\operatorname{Cov}(b_j,r_j)}{\operatorname{Var}(r_j)}r_j =\frac{\tau^2}{\tau^2+\sigma^2/4}r_j.

In practice the variance components are estimated, so the resulting predictor is an empirical BLUP. Substituting the REML estimates gives a shrinkage factor of

\hat{\gamma} =\frac{75/16}{75/16+(19/12)/4} =\frac{225}{244}=0.9221.

Worked Example: Calculate the Shrunken Batch Adjustments

Multiply each residual batch mean by $225/244$. For Batch 1, the adjustment is $(225/244)(1/6)=0.154$ log-expression units. For Batch 2, it is $(225/244)(13/6)=1.998$. For Batch 3, it is $(225/244)(-7/3)=-2.152$. Rounding to three decimals, the adjustments sum to zero.

A common batch-adjusted value is the observed expression minus the estimated batch effect:

y_{ij}^{\mathrm{adj}}=y_{ij}-\hat{b}_j.

For example, the Batch 2 State A observation $10.000$ becomes $10.000-1.998=8.002$ log-expression units. The correction subtracts the same estimated shift from both states in that batch; it does not alter their within-batch difference.

Worked Example: Check the State Contrast After Correction

The average State A expression across the three batches is $(9.500+11.000+6.500)/3=9.000$. The average State B expression is $(10.833+13.333+8.833)/3=11.000$. Their difference is $2.000$ log-expression units.

Each batch adjustment is subtracted from both state means, so the difference within that batch is unchanged. For instance, in Batch 2 the adjusted means are $11.000-1.998=9.002$ for State A and $13.333-1.998=11.335$ for State B, giving a difference of $2.333$ log-expression units, just as before correction. Averaging the three within-batch contrasts still gives $2.000$ log-expression units. Thus the shared random-intercept correction removes part of the batch-level midpoint shifts while retaining the fitted state contrast.

This preservation is a property of the chosen adjustment and balanced comparisons, not evidence that the model has captured every biological difference. In particular, the example contains variation among batch-specific state contrasts, which a random-intercept-only model does not represent as a batch-by-state random slope.

Diagnostics, Alternatives, and Limits

A numerical fit is useful only if its assumptions are reasonably compatible with the data and design. Inspect residuals against fitted values, batch, and state; look for strong skew, unequal residual spread, outliers, or patterns that suggest an unmodeled batch-by-state interaction. Also inspect the estimated batch effects and variance components. In this example, the batch-specific state contrasts $1.333,\;2.333,\;2.333$ are not identical, so a modeller should ask whether that spread is sampling noise or evidence that a common state contrast is inadequate. With only three batches, that distinction is especially uncertain.

The model also assumes that the batch effects are drawn from a common, mean-zero distribution and that cells are conditionally independent after accounting for the modeled effects. These assumptions may be poor if batches have very different technical processes, if cells from the same biological sample remain correlated, or if a measured covariate is omitted. Cells are not automatically independent biological replicates; the experimental design and sampling hierarchy must guide the model. A batch effect can only be estimated separately from state when the design contains sufficient overlap.

The main alternatives answer different questions. Treating batch as a fixed effect estimates a separate coefficient for each observed batch and does not shrink those coefficients toward a population mean; it can be preferable when interest is confined to those specific batches, but uses more parameters. Complete pooling applies no batch adjustment. A random-slope model allows state-associated expression to vary by batch, but requires enough batches and information to estimate that additional variation. Anchor-based corrections use relationships among cells or cell populations rather than this explicit random-intercept distribution. None is universally best: the design, assumptions, and intended adjusted-value interpretation matter.

Do not interpret adjusted values as newly observed biological measurements or as proof that all technical variation has been removed. This model estimates a shared additive shift on the modeled expression scale. It cannot repair confounding, recover a state absent from a batch, or ensure that all genes share the same batch effect. The later benchmarking material addresses how to define and evaluate integration targets; here, the immediate check is whether the fitted model and the state contrasts make sense for the data at hand.

Common pitfall: Subtracting each raw batch mean in full treats it as a known technical effect. In a mixed model, the residual batch mean is noisy, so its estimated effect is shrunken by a factor determined by the batch and residual variances. Conversely, assuming that any retained state contrast proves biological preservation ignores possible confounding, unmodeled state-dependent batch shifts, and model misspecification.

Check Your Understanding

Try each question, then compare your reasoning with the solutions.

In the model, what does $\beta_S$ measure when state is coded $-1/2,+1/2$?
Why is the batch-midpoint mean square not itself the estimate of $\tau^2$ in the worked example?
If the estimated shrinkage factor were $0.50$ and a batch residual mean were $1.2$ log-expression units, what would its predicted batch adjustment be?
Why can a random-intercept correction preserve a state contrast within a batch but fail to describe batch-specific state shifts?

Solutions. (1) $\beta_S$ is the State B minus State A contrast in log-expression units, holding the shared batch shift fixed. (2) A batch midpoint averages four residuals, so its variance includes $\sigma^2/4$ in addition to $\tau^2$; that contribution must be removed. (3) The predicted adjustment is $0.50\times1.2=0.600$ log-expression units. (4) A random intercept adds the same value to both states in a batch, so subtraction cancels in their difference. It cannot represent a batch-by-state interaction, which requires a state-dependent batch term such as a random slope.