Why Model Batch and Condition Together?
For a given gene, measured expression can vary with both biological condition and technical batch. If the goal is to estimate a condition difference while accounting for batch, a model can include indicators for both factors. The key is that the coefficient for condition then describes a comparison holding batch fixed, under the model’s assumptions. It is not automatically a causal effect, and the data must contain enough overlap between batches and conditions to support that comparison.
This tutorial focuses on additive models: the condition contrast is assumed to be the same in every batch on the model’s chosen scale. The earlier material on batch-condition confounding establishes why the observed design matters for identifying that contrast. Here, the emphasis is on writing the model and interpreting its coefficients.
Reference-Coded Linear Models
Suppose there are \(B\) batches and two conditions, control and treated. For cell or sample \(i\) and gene \(g\), let \(Y_{ig}\) be an expression outcome on a scale for which a linear model is appropriate, such as log-transformed normalized expression. Choose one batch as the reference, and let control be the reference condition. An additive model is
Here, \(I(\cdot)\) equals 1 when its statement is true and 0 otherwise. The coefficient \(\beta_{0g}\) is the mean expression for gene \(g\) in the reference batch under control, according to the model. The coefficient \(\gamma_{bg}\) is the difference between batch \(b\) and the reference batch under the additive model. The condition coefficient \(\delta_g\) is the treated-minus-control difference in every batch. All differences are on the scale of \(Y\).
The linear model is commonly used with continuous transformed outcomes. Its standard assumptions for classical inference are a correctly specified conditional mean, independent errors at the observational unit used for inference, errors with mean zero and constant variance, and normally distributed errors for exact small-sample \(t\)- and \(F\)-based inference. A full-column-rank design matrix is also required for unique coefficient estimates. For single-cell data, treating cells from the same biological sample as independent replicates can understate uncertainty; the analysis must respect the actual experimental units.
Worked Example: Interpreting a Linear Condition Coefficient
Consider hypothetical mean log-normalized expression values for one gene. There are two batches, A and B, with control and treated samples represented in both. The group means are:
| Batch | Control | Treated |
|---|---|---|
| A | 2.0 | 5.0 |
| B | 4.0 | 7.0 |
Use batch A and control as references. The control mean in batch A gives \(\hat{\beta}_{0g}=2.0\). The batch B control mean is 4.0, so \(\hat{\gamma}_{Bg}=4.0-2.0=2.0\). In batch A, the treated-minus-control difference is \(5.0-2.0=3.0\); in batch B it is \(7.0-4.0=3.0\). Thus the data follow the additive pattern, with \(\hat{\delta}_g=3.0\).
For a treated sample in batch B, the fitted value is \(\hat{Y}_{ig}=2.0+2.0+3.0=7.0\), matching the table. The condition coefficient means that treated expression is 3.0 log-expression units higher than control expression at either batch, under this model. It does not mean that treatment caused the difference.
These four means illustrate coefficient interpretation; by themselves they do not establish the error distribution, independence, or uncertainty needed for a valid inferential conclusion.
Generalized Linear Models for Counts
Raw gene counts are nonnegative integers, and their variability typically depends on their mean. A linear model with constant-variance errors is therefore not generally a suitable model for raw counts. A generalized linear model (GLM) specifies a distribution for the outcome and relates its mean to predictors through a link function.
For gene \(g\), let \(Y_{ig}\) be the count in observation \(i\), and let \(s_i\) be a known positive size factor representing relative sequencing exposure. A negative-binomial GLM with a log link can be written as
The term \(\log(s_i)\) is an offset: it adjusts the expected count for sequencing exposure, with its coefficient fixed at 1. The negative-binomial dispersion \(\phi_g\) allows variance to exceed the mean. As in the linear model, batch A and control are reference levels. The intercept and batch coefficients describe log expected counts per unit of size factor. The condition coefficient \(\delta_g\) is an adjusted log rate ratio; exponentiating it gives the treated-to-control expected rate ratio at a fixed batch.
This model assumes the conditional mean and log link are appropriate, the offset adequately represents exposure, and counts are conditionally independent at the modeled observational unit with the stated negative-binomial variance structure. In practice, biological replication and cell-level dependence matter: many cells from one sample do not automatically provide the same information as independently sampled biological replicates. The count model also does not turn an observational condition contrast into a causal effect.
Worked Example: Interpreting a Negative-Binomial GLM Coefficient
Suppose a hypothetical fitted model for one gene has \(\hat{\delta}_g=\log(1.5)\), with control and batch A as references. The coefficient is on the log expected-count scale after accounting for the size-factor offset and batch indicators.
Exponentiating gives the adjusted rate ratio:
Thus, at a fixed batch and size factor, the model’s expected count per unit of exposure for treated observations is 1.5 times that for control observations—a 50% higher expected rate. This interpretation depends on the log-link model and the common condition effect assumption. It is not a claim that every treated cell has 50% more counts, nor does it establish a causal effect.
For illustration, if the fitted control expected count at a specified batch and size factor is 20, the corresponding treated expected count is \(20\times1.5=30\). These values illustrate the coefficient’s meaning; a fitted coefficient alone does not provide its standard error or a test of whether the contrast differs from zero.
Why the Additive Form Matters
The model’s additivity is a substantive assumption, not merely a convenient notation. It says that the condition contrast is shared across batches on the selected scale. For the linear model, the treated-minus-control difference is constant. For the log-link GLM, the treated-to-control rate ratio is constant. A condition effect that varies by batch is not represented by this additive form.
One short derivation shows how the condition coefficient becomes an adjusted contrast. Write the linear model’s conditional mean in the reference batch for control as \(\mu_{A,C}=\beta_{0g}\), and for treated as \(\mu_{A,T}=\beta_{0g}+\delta_g\). Subtracting gives
For another batch \(b\), its term \(\gamma_{bg}\) appears in both means and cancels in the same subtraction. Therefore the model implies the same difference \(\delta_g\) in every batch. Under a log link, subtracting the two linear predictors instead gives \(\log(\mu_{b,T}/\mu_{b,C})=\delta_g\), so the implied rate ratio is \(\exp(\delta_g)\) in each batch.
Worked Example: Checking Whether an Additive Linear Fit Captures the Pattern
Suppose hypothetical replicated samples produce these mean log-expression values, with the same number of observations in each batch-condition group:
| Batch | Control | Treated |
|---|---|---|
| A | 2.0 | 4.0 |
| B | 5.0 | 10.0 |
The condition differences are \(4.0-2.0=2.0\) in A and \(10.0-5.0=5.0\) in B. They are not equal, so the common-difference assumption does not describe these group means exactly.
For equal group sizes, the additive least-squares estimates of the condition and batch differences are the averages of the corresponding within-factor differences: \(\hat{\delta}_g=(2.0+5.0)/2=3.5\) and \(\hat{\gamma}_{Bg}=[(5.0-2.0)+(10.0-4.0)]/2=4.5\). The intercept is the grand mean minus half of each effect: \(\hat{\beta}_{0g}=5.25-4.5/2-3.5/2=1.25\). The resulting fitted means are 1.25, 4.75, 5.75, and 9.25 for A-control, A-treated, B-control, and B-treated, respectively. The residuals are \(+0.75,-0.75,-0.75,+0.75\).
The systematic residual pattern is evidence that the additive model is inadequate for these group means. It does not, by itself, identify the reason; measurement problems, sampling variation, or a condition contrast that differs by batch are possibilities to investigate. The key diagnostic is that the model-implied common contrast fails to match the observed pattern.
Diagnostics, Alternatives, and Limitations
Before interpreting a fitted coefficient, inspect whether the model is being asked to estimate a comparison supported by the design, and whether its assumptions are plausible. Useful checks include:
- Cross-tabulate batch and condition. Confirm that the data include observations for the comparisons the model is intended to make. If some combinations are absent, adjusted estimates may rely on stronger assumptions or may not be uniquely determined.
- Inspect residuals for a linear model. Look for unequal spread, strong skewness, influential observations, and systematic patterns by batch or condition. These can undermine standard errors or signal a poor mean model.
- Check count-model fit. Compare observed and fitted counts across groups, assess whether variability is compatible with the assumed dispersion, and look for batches or samples with persistent residual patterns.
- Compare condition contrasts across batches. Large differences among within-batch contrasts warn that the common-effect assumption may be questionable. A single adjusted coefficient can conceal that heterogeneity.
- Verify the analysis unit and offset. Confirm that observations used for inference are appropriately independent and that the size factor represents the exposure adjustment intended by the analysis.
The main alternatives depend on the outcome and scientific target. A linear model is a reasonable choice for a suitably transformed, approximately continuous expression outcome when its residual assumptions are plausible. A negative-binomial GLM is a common choice for overdispersed counts with an exposure offset. A Poisson GLM may be considered when its mean-variance assumptions are adequate, but it can understate uncertainty when counts are overdispersed. Neither changing the outcome scale nor choosing a different distribution removes the need for batch-condition overlap or makes the additive assumption true.
A model with batch and condition indicators adjusts for measured batch differences under its specified form; it does not guarantee that all technical effects have been removed. It also cannot resolve unmeasured differences between groups or justify causal language on its own. The estimates describe model-based, conditional associations. Their uncertainty and scientific interpretation depend on the design, the experimental unit, and the assumptions stated above.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- In a linear model with batch A and control as references, what does \(\gamma_{Bg}\) represent?
- A negative-binomial GLM has a condition coefficient of \(\log(0.8)\). What is the adjusted rate ratio, and how should it be interpreted?
- In a linear additive model, the treated-minus-control differences are 1.5 units in batch A and 1.5 units in batch B. Is this pattern consistent with the model’s common-effect assumption?
- Why is \(\log(s_i)\) included as an offset in the count model, rather than estimated as another coefficient?
- Why can a condition coefficient from a model that includes batch still not be called a causal effect automatically?
Solutions. (1) It is the modeled difference in expression between batch B and reference batch A under the additive model, holding condition fixed, in units of the linear outcome. (2) \(\exp(\log(0.8))=0.8\); at a fixed batch and exposure, the expected treated count rate is 0.8 times the control rate, or 20% lower. (3) Yes. The equal differences match the linear model’s assumption of one common condition contrast. (4) The offset represents known relative exposure and fixes its coefficient at 1, so expected counts scale proportionally with the size factor. (5) Adjustment for batch does not establish random assignment, eliminate unmeasured confounding, or otherwise provide all assumptions needed for causal inference.