Simulating Confounding and Correction Failure

Why Simulate Confounding?

A correction can produce an expression matrix with visually aligned batches without recovering the biological condition effect. Simulation makes that distinction measurable: the analyst specifies the true condition and batch effects, generates data under a stated design, applies an adjustment, and compares the result with the known condition effect. The central question is not whether a correction returns numbers, but whether the study design contains enough information for those numbers to estimate the target.

The introductory material on batch-condition confounding establishes that overlap between batches and conditions matters. Here, simulation is used to make that principle concrete. We will compare perfect confounding, where the two effects cannot be separated, with partial and near confounding, where separation is possible but may be imprecise. The examples use one gene and a linear expression scale; they illustrate design and estimation, not a complete model for all single-cell counts.

A Generative Model for a Simulation

For a single gene, let \(Y_i\) denote the expression value for observation \(i\), let \(C_i=0\) for control and \(C_i=1\) for treated, and let \(b(i)\) identify its batch. A simple additive model is

Here, \(\mu\) is the control mean in a reference batch, \(\tau\) is the condition effect of interest, and \(\gamma_b\) is the effect of batch \(b\), with one batch effect set to zero as a reference. The error \(\varepsilon_i\) represents variation not explained by condition or batch.

For single-cell studies, independence at the cell level can be unrealistic because cells from the same biological sample are correlated. A simulation intended to evaluate biological-condition inference should therefore represent the relevant replication structure, such as independent samples or donors, rather than treating a large number of cells from one sample as independent biological replicates. The examples below isolate the confounding problem using independent observations so the calculations are transparent.

A useful simulation records the true \(\tau\) and \(\gamma_b\) before generating data. It then varies the allocation of observations across the batch-condition cells. Holding the true effects fixed while changing that allocation reveals how the design, rather than the correction algorithm alone, changes what can be estimated.

Perfect Confounding: The Effects Cannot Be Separated

This is a structural failure, not a matter of increasing the number of observations in the same two groups. More observations can make the combined difference more precise, but cannot supply the missing comparison. A simulation should label this case as non-identifiable rather than report a condition-effect estimate as though it had been recovered.

Partial Overlap: Adjustment Uses Within-Batch Comparisons

When both conditions occur within at least some batches, the condition effect can be estimated from within-batch contrasts under the additive model. For batch \(b\), let \(\bar{Y}_{1b}\) and \(\bar{Y}_{0b}\) be the treated and control sample means, and let \(n_{1b}\) and \(n_{0b}\) be their sample sizes. The batch contrast is \(D_b=\bar{Y}_{1b}-\bar{Y}_{0b}\). The batch effect cancels in this subtraction.

If the batches have independent errors with the same variance, the inverse-variance weight for \(D_b\) is proportional to \(w_b=n_{1b}n_{0b}/(n_{1b}+n_{0b})\). Combining the contrasts gives an adjusted estimate that gives more weight to batches with more precise within-batch comparisons.

Near Confounding: Identifiable Does Not Mean Precise

In a near-confounded design, both conditions occur within batches, but one of the groups needed for a within-batch contrast may be very small. The condition effect is then identifiable under the model, yet the estimate can be noisy. This distinction is important: a finite estimate from a fitted model does not show that the design has strong information about the target.

Justification of the Weighted Estimate

The weighting formula follows from the model for the batch contrasts. If \(D_b\) has mean \(\tau\) and variance \(v_b=\sigma^2(1/n_{1b}+1/n_{0b})\), consider a linear estimator \(\sum_b a_bD_b\) with \(\sum_ba_b=1\), so it is unbiased for \(\tau\). Its variance is \(\sum_ba_b^2v_b\), because the batches are independent. Minimizing this variance subject to the weights summing to 1 gives \(a_b\) proportional to \(1/v_b\), hence proportional to \(w_b=n_{1b}n_{0b}/(n_{1b}+n_{0b})\). Normalizing those weights yields the stated adjusted estimator. Its variance is \(\sigma^2/\sum_bw_b\).

This result relies on the stated common effect and variance assumptions. If condition effects differ by batch, the \(D_b\) do not share a common mean \(\tau\); a weighted average may still be a descriptive summary, but it is not automatically the single common effect in the additive model. If variances differ across batches, equal-variance weights need not be efficient. Simulations can test these departures by generating interaction effects or unequal error variances and checking how the chosen estimator behaves.

Comparing Adjustment Strategies and Diagnosing Failure

The main alternatives answer different questions. An unadjusted treated-control difference uses both between- and within-batch comparisons; if batch allocation differs by condition, it can mix condition and batch effects. A model that includes batch as a covariate estimates the condition effect from the design's available overlap under its modeling assumptions. Blind residualization removes variation associated with batch without preserving condition by design, and can remove genuine condition signal when condition and batch are associated. In perfect confounding, no method can identify the two effects from the observed data alone.

A useful simulation reports more than whether an algorithm completed. Across many repetitions, compare each estimate with the known \(\tau\). For estimates \(\hat{\tau}^{(r)}\) from \(R\) repetitions, empirical bias is the mean estimation error, and root mean squared error reflects both bias and variability.

For a non-identifiable design, these summaries should not be computed as if a unique adjusted condition estimate existed. For identifiable designs, also examine the distribution of estimates and interval coverage if uncertainty intervals are part of the evaluation. An estimator can have little average bias but still be too variable for useful inference.

• Check the design before interpreting the estimate. Tabulate counts by batch and condition. Empty cells indicate missing comparisons; very small cells warn of unstable contrasts.
• Check model behavior. Large residual patterns, markedly different within-batch contrasts, or poor calibration in repeated simulations can indicate that the additive, common-effect assumptions are inadequate.
• Separate algorithm failure from information failure. A correction that returns an estimate under perfect confounding has not created information to identify the condition effect.

A common mistake is to treat “the model produced an adjusted number” as evidence that confounding was resolved. That reasoning confuses computation with identification. The correction is to compare the fitted target with the simulation's known truth, inspect the design's within-batch overlap, and report non-identifiability or high uncertainty when the allocation cannot support a reliable estimate.