Deriving Linear Adjustment for Known Batch Effects

Why Derive the Adjustment?

A batch adjustment changes expression values to remove estimated systematic differences between processing groups. The key is to define what those differences mean. If an observed expression value reflects both biology and batch, subtracting a batch effect estimated without accounting for biology can remove part of the biological signal too. A linear model makes the intended separation explicit: it estimates batch shifts conditional on the biological variables in the model.

The derivation below applies to one feature at a time, such as one gene measured across $n$ samples on a continuous expression scale. The model can be fit separately for each feature; batch coefficients will generally differ across features.

Definition: For one feature, an additive batch model is $$ y_i = x_i^\mathsf{T}\beta + z_i^\mathsf{T}\gamma + \varepsilon_i,\qquad i=1,\ldots,n. $$ Here $y_i$ is the measured expression value for sample $i$, $x_i$ contains an intercept and the biological covariates to retain, and $z_i$ contains indicators for its batch. The coefficient vector $\beta$ describes the modeled biological effects; $\gamma$ describes batch shifts relative to a chosen reference coding. The error $\varepsilon_i$ is the part not explained by these terms.

Least-Squares Estimates and Their Conditions

Let $X$ be the $n$-row matrix of biological covariates and $Z$ the batch-indicator matrix. With $y$ the vector of expression values, ordinary least squares chooses $\hat{\beta}$ and $\hat{\gamma}$ to minimize the residual sum of squares:

\operatorname{RSS}(\beta,\gamma) = (y-X\beta-Z\gamma)^\mathsf{T}(y-X\beta-Z\gamma).

Write $D=[X\;Z]$ and $\theta=(\beta^\mathsf{T},\gamma^\mathsf{T})^\mathsf{T}$. If the columns of $D$ are linearly independent, the least-squares solution is unique and satisfies the normal equations $D^\mathsf{T}D\hat{\theta}=D^\mathsf{T}y$. Thus:

\hat{\theta}=(D^\mathsf{T}D)^{-1}D^\mathsf{T}y.

The full-rank condition matters: the observed design must contain enough separation of biology and batch to estimate their effects separately. This is the identifiability issue in batch-condition confounding. A reference batch is commonly coded with all-zero indicators, while the remaining batch indicators are included as columns of $Z$; the intercept in $X$ then represents the reference batch level. Alternative coding is possible, but a constraint is needed so that an intercept and all batch indicators do not duplicate one another.

The batch coefficients are estimated conditional on the biological covariates. To see the role of that adjustment, assume $X$ has full column rank and define $M_X=I-X(X^\mathsf{T}X)^{-1}X^\mathsf{T}$, the matrix that removes the part explained by $X$. The least-squares batch estimate can be written as:

\hat{\gamma} = (Z^\mathsf{T}M_X Z)^{-1}Z^\mathsf{T}M_X y,

provided $Z^\mathsf{T}M_X Z$ is invertible. This follows by first minimizing the residual sum of squares over $\beta$ for a fixed $\gamma$: the remaining residual is $M_X(y-Z\gamma)$. Minimizing its squared length over $\gamma$ gives the normal equation $Z^\mathsf{T}M_XZ\hat{\gamma}=Z^\mathsf{T}M_Xy$, and hence the displayed expression. In words, the fit compares batch differences after accounting for the modeled biology.

The additive model assumes that a batch contributes a shift on the chosen expression scale, not a batch-specific change in the biological response. For least-squares coefficient estimates to be unbiased conditional on the design, the errors must have mean zero given the modeled covariates. Independent, equal-variance errors are useful for standard classical inference, but are not required merely to calculate the least-squares fit. The model also presumes that the chosen biological terms and batch indicators adequately describe the systematic mean pattern.

From Coefficients to Adjusted Values

After fitting the model, the fitted batch shift for sample $i$ is $z_i^\mathsf{T}\hat{\gamma}$. To express all samples on the reference-batch scale, subtract that fitted shift:

y_{i,\mathrm{ref}}^{\mathrm{adj}} = y_i-z_i^\mathsf{T}\hat{\gamma}.

Since the fitted value is $x_i^\mathsf{T}\hat{\beta}+z_i^\mathsf{T}\hat{\gamma}$, the adjusted observation can also be written $y_{i,\mathrm{ref}}^{\mathrm{adj}}=x_i^\mathsf{T}\hat{\beta}+\hat{\varepsilon}_i$. Thus the modeled biological contribution remains, along with the residual; the fitted batch contribution is removed. Reference adjustment changes the overall location to the reference-batch scale.

A centered adjustment instead subtracts each fitted shift after subtracting the sample-average fitted shift. If $\bar d=n^{-1}\sum_i z_i^\mathsf{T}\hat{\gamma}$, then:

y_{i,\mathrm{cent}}^{\mathrm{adj}} = y_i-z_i^\mathsf{T}\hat{\gamma}+\bar d.

This removes differences between batches while restoring the original overall sample mean. Reference and centered adjustments therefore differ by a common additive constant, not by their estimated batch contrasts. When samples are unequally distributed across batches, $\bar d$ is weighted by the observed sample counts.

Worked Example: Estimating a Batch Shift While Retaining Condition

Suppose a hypothetical feature is measured in two conditions, A and B, in each of two batches. The values are: A in batch 1, 10; A in batch 2, 15; B in batch 1, 20; and B in batch 2, 25. Use batch 1 and condition A as reference levels. The model is $y=\mu+\tau I(B)+\gamma I(\mathrm{batch\ 2})+\varepsilon$.

The observed pattern is exactly additive: the condition B values exceed the corresponding condition A values by 10, and batch 2 exceeds batch 1 by 5 in each condition. The least-squares estimates are $\hat{\mu}=10$, $\hat{\tau}=10$, and $\hat{\gamma}=5$. For example, the fitted value for condition B in batch 2 is $10+10+5=25$, matching the observed value.

For reference-batch adjustment, subtract 5 from each batch 2 measurement and subtract 0 from each batch 1 measurement. In sample order A-batch 1, A-batch 2, B-batch 1, B-batch 2, the original values are $10,15,20,25$, so the adjusted values are $10-0=10$, $15-5=10$, $20-0=20$, and $25-5=20$. The condition difference remains $20-10=10$.

The batch shifts are 0 and 5, whose sample-average shift is $(0+5+0+5)/4=2.5$. Centering therefore subtracts each shift and adds 2.5: $10-0+2.5=12.5$, $15-5+2.5=12.5$, $20-0+2.5=22.5$, and $25-5+2.5=22.5$. The condition contrast is still 10. The reference-adjusted values are on the batch 1 scale; the centered values preserve the overall mean of 17.5.

Worked Example: A Continuous Biological Covariate

Consider four hypothetical samples with a dose $x$ and expression $y$: batch 1 has $(x,y)=(0,2)$ and $(1,5)$; batch 2 has $(0,6)$ and $(1,9)$. Fit $y=\beta_0+\beta_1x+\gamma I(\mathrm{batch\ 2})+\varepsilon$.

Within each batch, increasing dose from 0 to 1 raises expression by $5-2=3$ in batch 1 and $9-6=3$ in batch 2, so the fitted dose slope is $\hat{\beta}_1=3$. At dose 0, batch 1 has value 2 and batch 2 has value 6, giving $\hat{\beta}_0=2$ and $\hat{\gamma}=6-2=4$. These estimates reproduce all four observations: for example, batch 2 at dose 1 has fitted value $2+3(1)+4=9$.

Reference adjustment subtracts 4 from the batch 2 values. The values in the original order become $2,5,6-4=2,9-4=5$: both batches now share the same fitted dose pattern, 2 at dose 0 and 5 at dose 1. Centering uses the average fitted shift $(0+0+4+4)/4=2$, giving $4,7,4,7$. In either adjustment, the fitted dose contrast remains 3 expression units per dose unit. This preservation relies on including dose in the model and on the additive fit being appropriate.

Worked Example: Unequal Batch Sizes and Centering

For a hypothetical feature with no additional biological covariate, suppose batch 1 has values 10 and 12, while batch 2 has values 14, 16, and 18. Under an intercept-plus-batch model, the batch means are $(10+12)/2=11$ and $(14+16+18)/3=16$. The fitted batch 2 shift relative to batch 1 is $16-11=5$.

Reference adjustment leaves batch 1 unchanged and subtracts 5 from each batch 2 value. The adjusted values are $10,12,9,11,13$; the mean in each batch is now 11.

For centered adjustment, the sample-average fitted batch shift is $\bar d=(2(0)+3(5))/5=3$. Subtract the shift and add 3: batch 1 becomes $10+3=13$ and $12+3=15$; batch 2 becomes $14-5+3=12$, $16-5+3=14$, and $18-5+3=16$. Each batch mean is 14, the original overall mean $(10+12+14+16+18)/5=14$. Using the unweighted average of the two batch shifts would give the wrong centering constant here because the batch sizes differ.

Diagnostics, Alternatives, and Limits

Adjustment is not evidence that the model has separated biology and batch correctly. Before using the adjusted values, check the design and the fit.

Check overlap and rank. Inspect which biological groups occur in each batch and verify that the combined design matrix has full column rank. If each condition occurs in only one batch, the data cannot determine which difference is biological and which is a batch shift. Least squares cannot resolve that ambiguity by itself.
Inspect residuals by batch and biology. After fitting, residuals should not show a systematic batch pattern or a clear batch-specific biological trend. Persistent shifts suggest incomplete adjustment; different slopes across batches suggest that a single additive shift may be inadequate.
Check influential observations and stability. A few unusual samples can move a fitted shift substantially, especially when a batch is small. Compare estimates and residual plots with and without clearly influential observations, while investigating rather than automatically deleting them.
Use an appropriate expression scale. A constant additive shift is a statement about the scale on which the model is fit. A shift on a log-expression scale does not represent the same correction as a constant shift on a raw-count scale. Do not apply this continuous-response derivation blindly to raw sequencing counts.

The central alternative for many analyses is not to create adjusted values, but to keep the original expression values and include batch indicators alongside biological covariates in the downstream model. This directly estimates biological contrasts conditional on batch and avoids treating estimated adjusted values as if they were measured without uncertainty. If adjusted values are needed for visualization or another clearly specified purpose, state whether reference or centered adjustment was used.

A batch-only regression followed by subtraction is simpler, but it does not protect biology when batch and biology are associated: the batch-only fit may attribute part of a biological contrast to batch. Joint modeling protects only the effects represented in the biological design, and only when the design supports their separate estimation. An additive model also cannot represent batch-by-biology interactions. More elaborate batch methods may address different structures or estimation goals, but they do not remove the need to examine identifiability and biological preservation.

Common pitfall: “Batch 2 is higher, so subtract its average difference from every batch 2 sample.” This is unsafe if the batches have different biological composition: that average difference can include biology. Instead, fit batch and the biological variables together when the design permits, then subtract the fitted batch component. If biology and batch are confounded, no adjustment formula can recover their separate effects from those data alone.

Key takeaway: Least-squares adjustment removes the fitted batch term from an additive model. It preserves the biological terms included in that model, but only if the design identifies them and the additive assumptions are reasonable. Centering changes the overall level, not the fitted batch contrasts.

Check Your Understanding

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

In an additive model with a reference batch, what does the batch coefficient represent?
Three samples have fitted batch shifts 0, 2, and 2. What is the centered adjustment constant, and how does it enter the formula?
Why can a batch-only adjustment distort a biological contrast when batch and condition are associated?
Suppose a fit leaves a consistent residual elevation in one batch. What does this suggest, and name one diagnostic step.

Solutions. (1) It is the estimated expression shift for that batch relative to the reference, conditional on the modeled biological covariates. (2) The average fitted shift is $(0+2+2)/3=4/3$; subtract each sample’s fitted shift and add $4/3$. (3) The batch-only estimate can absorb biological differences that occur unevenly across batches, so subtracting it may remove part of the biological contrast. (4) The additive model may not have removed the batch pattern, or may be misspecified; plot residuals by batch and check for omitted structure or influential samples.

Deriving Linear Adjustment for Known Batch Effects

What You'll Learn