Covariate-Preserving Batch Correction

Why the protected covariates matter

Batch correction is not simply the removal of every pattern associated with batch. A biological variable can be unevenly represented across batches, so some measured differences may reflect both biology and technical shifts. A covariate-preserving adjustment specifies which biological effects the model should retain while estimating batch contrasts conditional on those covariates.

The protected covariates determine what the model treats as biological structure worth preserving. If treatment group is included, its fitted effect is retained; if age is included, the fitted age relationship is retained. This is a modelling choice, not proof that the retained component is entirely biological or that the removed component is entirely technical. The earlier material on estimable biological contrasts under confounding is especially relevant: adjustment cannot distinguish effects that the observed design does not identify.

Definition: For one molecular feature measured in $n$ samples, a covariate-preserving additive batch model is $$ y_g = Z\beta_g + C\gamma_g + \varepsilon_g. $$ Here $y_g$ is the $n$-vector of measurements for feature $g$; $Z$ is the design matrix for protected biological covariates, usually including an intercept; $C$ is a matrix of batch contrasts; $\beta_g$ and $\gamma_g$ are feature-specific coefficients; and $\varepsilon_g$ is the unexplained component. With a chosen reference batch, the corrected measurements are $y_g^* = y_g-C\hat{\gamma}_g$. They are expressed on the reference-batch scale while retaining the fitted effects represented by $Z\hat{\beta}_g$.

Model, assumptions, and the correction

For example, $Z$ might contain an intercept, treatment indicators, and age, while $C$ contains indicators for all batches other than a designated reference batch. A coefficient in $\gamma_g$ then represents an additive difference from that reference, conditional on the covariates in $Z$. A continuous covariate such as age can be represented by a linear term; if the intended relationship is nonlinear, the model must include appropriate terms, such as a spline basis, for that relationship to be protected.

The adjustment relies on a defensible design and mean model. In particular, the columns of $[Z\ C]$ must be linearly independent to estimate all included coefficients uniquely. The conditional mean should be adequately represented by the specified covariate and batch terms, and the additive model assumes the batch contrast does not vary across covariate values. For ordinary least squares point estimates, a standard condition is $E(\varepsilon_g\mid Z,C)=0$. Independence and constant variance are not required for the algebraic correction, but are assumptions of the usual classical standard-error calculations. If samples are dependent or variance differs, uncertainty calculations need to account for that structure.

Fit the model separately for each feature, or use a suitable joint method that implements the same design. In ordinary least squares, $\hat{\beta}_g$ and $\hat{\gamma}_g$ minimize the sum of squared residuals. Subtracting only the fitted batch component leaves the fitted covariate component and residuals:

y_g^* = y_g-C\hat{\gamma}_g = Z\hat{\beta}_g+\hat{\varepsilon}_g.

This identity gives the preservation property. Since the least-squares residual $\hat{\varepsilon}_g$ is orthogonal to every column of the fitted design, $Z^\mathsf{T}\hat{\varepsilon}_g=0$. Therefore, when $Z$ has full column rank, fitting $y_g^*$ on $Z$ gives

(Z^\mathsf{T}Z)^{-1}Z^\mathsf{T}y_g^* =(Z^\mathsf{T}Z)^{-1}Z^\mathsf{T} (Z\hat{\beta}_g+\hat{\varepsilon}_g) =\hat{\beta}_g.

The correction thus preserves the fitted covariate coefficients in this model. It does not guarantee that every downstream analysis will reproduce the same result: downstream models may use different covariates, transformations, weights, or sample subsets. Nor does the identity establish that the chosen model is scientifically appropriate.

Worked Example: Preserving a Treatment Contrast

Consider four hypothetical measurements of one feature. Treatment and batch are crossed, with one control and one treated sample in each of two batches. Measurements are in arbitrary normalized units.

Batch	Treatment	Measurement
A (reference)	Control	10
A (reference)	Treated	14
B	Control	13
B	Treated	17

Use an intercept, a treated indicator, and a batch-B indicator. The model is $y_i=\beta_0+\beta_1T_i+\gamma_B B_i+\varepsilon_i$. The control measurement in batch A gives $\hat{\beta}_0=10$. Within batch A, the treatment difference is $14-10=4$, so $\hat{\beta}_1=4$. The batch-B control is 13, which is 3 units above the batch-A control; hence $\hat{\gamma}_B=13-10=3$. The batch-B treated value checks the additive model: $10+4+3=17$, matching the observation.

Correct to batch A by subtracting $3B_i$. The corrected measurements are $10,14,13-3=10,$ and $17-3=14$ units. In each batch the treatment contrast is $14-10=4$ units, and the fitted treatment coefficient remains 4 units. The adjustment removes the estimated three-unit batch contrast, not the specified treatment difference.

Choosing what the corrected values represent

Reference coding gives a direct interpretation: corrected values are on the scale of the reference batch, conditional on the fitted covariates. A different target can be defined by centering the fitted batch effects around a stated average rather than using one reference. The target must be explicit because the corrected numerical values depend on it, even when covariate contrasts are retained. With unequal batch sizes, a sample-weighted average and an equally weighted average across batches need not be the same.

The protected set should match the scientific estimand and intended use. A categorical treatment indicator preserves its fitted group difference; a continuous covariate preserves its fitted slope or specified curve. Omitting a biological covariate means its variation is not separately protected by the model. Including unnecessary or poorly supported terms can also make the design unstable or rank deficient. Covariate preservation is therefore a property of a specified model and estimable design, not a universal guarantee about all biology.

Worked Example: Preserving an Age Trend

Suppose four hypothetical samples have ages $x=0,1,1,2$ years after a chosen origin. The first two are in reference batch A, and the last two are in batch B. Measurements, in arbitrary units, are respectively 5, 7, 10, and 12. Consider $y_i=\beta_0+\beta_1x_i+\gamma_BB_i+\varepsilon_i$.

The A observations give an intercept of 5 units and a slope of $(7-5)/(1-0)=2$ units per year. The batch-B observation at age 1 is 10 units, while the fitted reference-batch value at age 1 is $5+2(1)=7$, so the estimated batch contrast is $10-7=3$ units. The age-2 observation checks the model: $5+2(2)+3=12$. These four points therefore fit exactly, with $\hat{\beta}_0=5$, $\hat{\beta}_1=2$, and $\hat{\gamma}_B=3$.

Subtract 3 units from each batch-B measurement and zero from each batch-A measurement. The corrected values are 5, 7, 7, and 9 units. Across the corrected data, the fitted age slope is still $(9-5)/(2-0)=2$ units per year, and the age-1 samples both have value 7. The coefficient describes the modelled trend; it should not be interpreted as causal solely because it is preserved.

When the design cannot separate biology and batch

If batch and a protected covariate are perfectly confounded, the model cannot uniquely estimate their separate effects. Adjustment software cannot solve this by retaining the covariate name in the formula: with dependent design columns, multiple coefficient combinations give the same fitted values. Earlier material on partial confounding explains why even imperfect overlap can increase uncertainty; complete confounding is more severe because the separate contrast is not estimable.

Worked Example: A Confounded Treatment and Batch

Suppose two hypothetical control samples are measured only in batch A, with values 8 and 10 units, while two treated samples are measured only in batch B, with values 15 and 17 units. The batch means are $(8+10)/2=9$ and $(15+17)/2=16$ units; the observed treated-minus-control difference is $16-9=7$ units.

In a model with an intercept, treatment indicator, and batch-B indicator, treatment and batch-B have identical columns: both are 0 for the controls and 1 for the treated samples. Only their sum of coefficients is identified. For any number $a$, assigning a treatment coefficient $a$ and a batch coefficient $7-a$ gives the same group-mean difference of 7 units. The data therefore do not determine how much of that difference is treatment and how much is batch.

No unique covariate-preserving correction can be justified from these observations alone. The correct conclusion is that the treatment effect and batch effect cannot be separated in this design, not that a particular correction has successfully preserved biology.

Diagnostics, alternatives, and limitations

Before adjustment, inspect the design matrix for rank and confirm that the protected contrasts are estimable. Check how covariate values and groups are represented within each batch; a coefficient may technically be estimable but rely heavily on sparse overlap. After fitting, examine residuals against fitted values, batch, and protected covariates, and inspect whether substantial batch-associated patterns remain. Residual structure can indicate an inadequate mean model, an omitted interaction, outliers, or dependence; it does not by itself identify the cause. The earlier material on residual checks after batch adjustment provides a broader diagnostic framework.

A covariate-preserving adjustment is not the only way to handle batch. If the planned analysis can include batch as a covariate, modelling batch directly in that analysis avoids creating a corrected matrix and can keep the inferential model explicit. A simple batch-mean subtraction is easier but does not condition on protected biology, so it can remove biological differences when group representation varies by batch. ComBat is an alternative that uses empirical-Bayes estimation for batch adjustment; its assumptions and design choices still matter, and the protected-covariate principle remains relevant.

The additive feature-level model has important limits. It may not capture batch-by-biology interactions, nonlinear covariate relationships, or feature-specific variance changes. A fitted batch contrast can also be noisy, especially with few samples or weak overlap. Correction removes an estimated component, so it adds estimation uncertainty even though the corrected values are often treated as fixed in later analyses. The method does not establish that a retained effect is biologically meaningful or that removed variation was purely technical.

Common pitfall: “Treatment and batch are imbalanced, so subtract each batch mean and the treatment effect will remain.” That reasoning ignores the possibility that the batch means include different treatment mixtures. Instead, fit batch and the specified treatment covariate together, verify that their effects are estimable, and subtract only the fitted batch contrast. If the design is confounded, report that the effects cannot be separated rather than presenting an arbitrary allocation as a correction.

Specify the target.
Name the biological covariates to preserve and the batch scale to which corrected values should refer.

Check identifiability.
Build the joint design and verify that the protected contrasts and batch contrasts are estimable.

Fit and correct.
Estimate feature-specific coefficients under the stated model and subtract only its fitted batch component.

Inspect and report.
Check residual patterns and describe the covariates, reference batch, model assumptions, and limitations.

Check Your Understanding

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

In $y_g=Z\beta_g+C\gamma_g+\varepsilon_g$, what is subtracted to correct to the reference batch, and which modelled component is retained?
A fitted model estimates an age slope of 1.8 units per year and a nonreference-batch contrast of 2.5 units. What is the correction for a sample in that nonreference batch?
Why does a design with treatment indicators identical to batch indicators fail to separate treatment from batch?
Name one model failure that could make the additive correction inadequate, and one diagnostic that could reveal a problem.

Solutions. (1) Subtract $C\hat{\gamma}_g$; the fitted protected component $Z\hat{\beta}_g$ remains in the corrected measurements. (2) Subtract 2.5 units from that sample; the specified age effect remains in the model. (3) Identical columns are linearly dependent, so only the sum of the corresponding coefficients is identified, not their separate values. (4) A batch-by-age interaction can violate the additive model; for example, inspect residuals against age within each batch for systematic patterns.