Why Adjust for Sample-Level Covariates?
Pseudobulk analysis aggregates cell counts within each biological sample, often separately for each cell type, and models the resulting sample-level counts. Across samples, age, sex, or another measured characteristic may vary along with condition. If such a covariate is related to expression and is unevenly distributed between conditions, an unadjusted comparison can mix the condition difference with the covariate difference. Adjustment separates their modeled contributions, provided the design contains enough information to estimate them.
The unit contributing independent replication is the biological sample (or donor, with repeated samples modeled appropriately), not each cell used to form that sample's pseudobulk. Recording age or sex for thousands of cells from a few donors does not create thousands of independent covariate observations. The sample-level covariate has one value per sample, and inference about its relationship to expression is supported by the number and arrangement of biological samples.
Model and Interpretation
For gene \(g\) and biological sample \(i\), let \(Y_{gi}\) be the pseudobulk count and \(\mu_{gi}=E(Y_{gi})\). A common model is a negative binomial generalized linear model (GLM):
where \(\phi_g\) is the gene-specific dispersion. With a library-size offset, binary condition \(C_i\), centered age \(A_i\), and a binary sex indicator \(S_i\), one possible linear predictor is
Here \(s_i\) is the sample's size factor, the coefficient for the offset \(\log(s_i)\) is fixed at 1, and \(A_i\) is age minus a chosen centering value. The reference levels are \(C_i=0\) and \(S_i=0\). A design matrix records the intercept and these sample-level predictor columns; the coefficients are estimated from the biological samples. This uses the negative binomial and offset framework from earlier material on pseudobulk GLMs.
With no condition-by-covariate interaction, \(\beta_{Cg}\) is the condition log-fold change at any fixed age and sex. The corresponding fold change is \(\exp(\beta_{Cg})\). This common conditional contrast follows directly by comparing two samples that have the same covariates and offset:
That cancellation is the central interpretation: adjustment compares modeled means at common covariate values rather than comparing groups with whatever covariate distributions happened to be observed. It also exposes an assumption: this model assigns the same condition log-fold change across the included covariate values. If the condition effect varies with age, for example, a model without a condition-by-age interaction does not describe that variation.
Worked Model Interpretation
Worked Example: Testing an Adjusted Condition Effect
Suppose a hypothetical pseudobulk analysis fits the model above, adjusting for centered age and sex. For one gene, the fitted condition coefficient is \(\hat{\beta}_{Cg}=0.4055\), with standard error \(SE=0.1800\). The Wald test assesses \(H_0:\beta_{Cg}=0\) against a two-sided alternative, using the large-sample standard normal reference distribution.
The Wald statistic is
The two-sided p-value is \(2[1-\Phi(2.25278)]\approx 0.0243\). At \(\alpha=0.05\), this is evidence against a zero adjusted condition coefficient under the fitted model. The condition fold change is
Thus, for this gene, the fitted mean pseudobulk count is about 1.50 times as large in the \(C_i=1\) condition as in \(C_i=0\), for samples with the same modeled age, sex, and size factor. This is an estimated conditional association; it does not establish that condition caused the difference. The p-value is not the probability that the null hypothesis is true.
Worked Example: Predicting at a Common Covariate Profile
For a second hypothetical fitted model, suppose the estimated coefficients for one gene are \(\hat{\beta}_0=2.00\), \(\hat{\beta}_C=0.30\), \(\hat{\beta}_A=0.018\) per year, and \(\hat{\beta}_S=0.25\). Age is centered at 50 years, and \(S=1\) denotes the non-reference sex category. Compare the conditions at age 60, \(S=1\), and size factor \(s=1.2\).
At age 60, \(A=60-50=10\). For \(C=0\), the linear predictor excluding the offset is
The expected count is \(1.2\exp(2.43)\approx 13.64\). For \(C=1\), add the condition coefficient: \(2.43+0.30=2.73\), giving expected count \(1.2\exp(2.73)\approx18.40\). Their ratio is \(18.40/13.64\approx1.349\), subject to rounding; directly, \(\exp(0.30)\approx1.350\). The model therefore predicts about 35.0% higher mean count in condition 1 at this common covariate profile. Different common profiles yield different absolute means, but the same ratio under the no-interaction model.
Estimability and Design Checks
Adjustment is possible only when the observed design distinguishes the condition column from the covariate columns. A coefficient cannot be recovered by statistical software merely because it is included in a formula. Before interpreting estimates, inspect the sample-level design: how many biological samples are in each condition, whether covariate values overlap across conditions, and whether any predictor is an exact linear combination of the others.
Worked Example: Condition and Age Are Completely Confounded
Consider six hypothetical independent samples: three controls all aged 20 years and three treated samples all aged 40 years. Let \(C=0\) for control and \(C=1\) for treated, and let \(A\) be age centered at 20. For every sample, \(A=20C\). In a model with intercept, condition, and age, the age column is therefore exactly 20 times the condition column. The design matrix is rank deficient.
Consequently, the data cannot determine whether a difference in expression is associated with condition, age, or a combination of the two. No unique adjusted condition coefficient is estimable from this design; dropping one predictor would produce a coefficient, but it would answer a different, unadjusted question. More samples with age variation represented in both conditions could provide the overlap needed to separate the effects. The issue is lack of information in the design, not a small p-value or a software setting.
A related issue arises in repeated measures. If age is constant within each donor and the model includes a separate fixed effect for every donor, age is a linear combination of donor indicators and its main effect cannot also be estimated. The within-donor condition contrast may still be estimable if condition varies within donors, but the data do not identify a separate age coefficient in that fixed-effects model. A random-intercept model can represent donor dependence and estimate a between-donor covariate association under its distributional assumptions; it does not manufacture information absent from the sample design.
Choosing and Checking Covariates
Covariates should be selected because they are scientifically relevant to the comparison and measured at the correct level, not because they happen to yield a smaller p-value. A covariate associated with both condition assignment and expression may help address measured confounding. But adjusting for a mediator or a collider can change the scientific question or introduce bias. In observational data, measured adjustment alone cannot guarantee a causal interpretation.
For a continuous variable such as age, a linear term on the log-mean scale assumes a constant change in log expected count per year. Centering age changes the meaning of the intercept but not the fitted condition coefficient in this no-interaction model. If the age relationship is visibly nonlinear or scientifically expected to be nonlinear, a suitable nonlinear functional form may be needed. The number of biological samples must support the additional terms; adding many covariates to a small design can make estimates unstable.
For repeated observations from the same donor, use a model that represents the repeated-sample dependence, as in the earlier material on donor effects and covariance structures. Treating repeated pseudobulks as independent would understate uncertainty. Conversely, duplicating sample-level age or sex across cells and fitting a cell-level analysis does not solve a replication problem: cells from a donor share biological context and are not independent replicates of that donor's age or sex.
Compare an adjusted model with a clearly defined unadjusted model to understand how the estimate changes, but do not select the preferred result simply by which is more significant. Large changes can indicate covariate imbalance, a changed target contrast, or model sensitivity. Inspect the design matrix for rank deficiency and near-collinearity, verify covariate coding and units, examine whether condition groups have covariate overlap, and confirm that the intended offset and sample identifiers are correct. Model fit and influential-sample checks are also important; adjustment cannot repair an inappropriate count model or a result driven by one biological sample.
Alternatives and Limitations
An unadjusted condition model is simpler and can be appropriate when covariates are balanced or when the desired comparison is the observed, unadjusted contrast. It does not, however, separate condition from covariate differences. Stratified comparisons can make patterns across covariate groups visible, but they may be inefficient and require adequate samples in each stratum. A multivariable GLM estimates the conditional contrast in one model, but relies on its functional form and on sufficient covariate overlap.
In a paired or repeated design, donor fixed effects or a suitable random-effects/covariance model address within-donor structure; these are not interchangeable with simply adding age or sex. A covariate adjustment term does not itself model repeated-measure correlation. Likewise, an adjusted coefficient is conditional on the other variables in the model, and may differ from a population-average contrast when effects vary across covariate profiles.
Finally, pseudobulk adjustment operates on sample-level counts for a chosen cell population. It does not estimate a cell-specific covariate effect from each cell as an independent observation, and it does not explain why cell abundance or composition differs between samples. Its scope is the modeled gene-expression contrast for the defined pseudobulk observations.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- In a pseudobulk model, why does recording age for 8,000 cells from four donors not provide 8,000 independent age observations?
- A no-interaction model estimates a condition coefficient of \(-0.20\). What is its condition fold change, approximately, and how should it be interpreted?
- Suppose every treated sample is older than every control sample, but age still varies within each condition. Is the condition effect necessarily unidentifiable? What design feature should be checked?
- Age is constant within donor, and a model includes a fixed effect for each donor. Can the model also estimate a separate age main effect? Explain.
- Why can an adjusted condition coefficient not automatically be described as a causal effect?
Solutions. (1) Age is measured at donor or sample level and is shared by that donor's cells. The independent biological replication for age is four donors, not 8,000 cells. (2) The fold change is \(\exp(-0.20)\approx0.819\). Conditional on the modeled covariates, the mean is estimated to be about 0.819 times as large in the \(C=1\) condition, or about 18.1% lower. (3) Not necessarily. Exact confounding occurs when the design columns are linearly dependent. Check age overlap between conditions, rank, and whether the observed age ranges provide adequate information to separate the terms; weak overlap can still make estimates unstable. (4) No. The age column is determined by the donor indicators, so the separate age coefficient is not identifiable alongside donor fixed effects. (5) The model adjusts only for included covariates and depends on its assumptions. Unmeasured confounding, inappropriate adjustment, or model misspecification can remain, so conditional association alone does not establish causation.