Why Allow Treatment Effects to Vary Across Donors?
In a paired pseudobulk experiment, each donor may contribute samples in control and treated conditions. A model with a common treatment effect assumes that, after accounting for modeled variation, donors share the same treatment effect. That may be too restrictive: biological response can differ across donors. A random condition slope represents such donor-to-donor variation in response.
This is a model for pseudobulk counts, not for treating individual cells as independent experimental units. Each count is an aggregate for a donor, condition, and cell type (and, if applicable, a specified sample or time point). The donor remains the unit that links repeated samples and supports inference about between-donor variation.
Negative-Binomial Model With a Random Slope
Let \(Y_{ij}\) be the pseudobulk count for a gene in sample \(j\) from donor \(i\), and let \(x_{ij}\) indicate condition. A common parameterization uses a negative-binomial distribution with conditional mean \(\mu_{ij}\) and dispersion \(\phi>0\):
The log-link model is
Here \(L_{ij}\) is a library-size or normalization factor, so \(\log(L_{ij})\) is an offset; \(\mathbf{z}_{ij}\) contains any other specified fixed-effect covariates; and \(\boldsymbol{\gamma}\) contains their coefficients. The donor effects are modeled jointly as
The covariance matrix must be positive semidefinite: \(\tau_{00}\geq0\), \(\tau_{11}\geq0\), and \(\tau_{00}\tau_{11}-\tau_{01}^{2}\geq0\). The variance \(\tau_{11}\) describes between-donor variation in condition slopes; \(\tau_{01}\) allows baseline expression and treatment response to be associated. Setting \(\tau_{11}=0\) removes slope variation, while a random-intercept-only model also has no donor-specific condition effect.
The model assumes donors are independent draws from the specified random-effects distribution, and that observations are independent conditional on their donor effects and modeled covariates. The negative-binomial variance describes remaining count variation conditional on those effects. The log link makes conditional means positive. The specified offset must adequately represent the relevant library-size scaling, and the condition and covariate coding must match the design. These are modeling assumptions, not consequences of fitting the model.
Worked Example: Writing the Paired Count Model
Suppose a hypothetical study has a control and treated pseudobulk sample from each donor. For one gene, a donor's control count is 80 and treated count is 150; the corresponding normalization factors are 0.8 and 1.2. There are no additional covariates. Write the model for these two observations.
Use \(x=0\) for control and \(x=1\) for treatment. The control mean for donor \(i\) is \(\log(\mu_{i0})=\log(0.8)+\beta_0+b_{0i}\), and the treatment mean is \(\log(\mu_{i1})=\log(1.2)+\beta_0+\beta_1+b_{0i}+b_{1i}\). Each count follows the negative-binomial model with variance \(\mu+\phi\mu^2\), conditional on the donor effects. The observed counts, 80 and 150, are possible realizations under the model; they do not themselves determine the random effects or establish a treatment response.
The donor-specific log mean contrast after removing the offset difference is \(\beta_1+b_{1i}\). The raw count ratio \(150/80\) is not that contrast: it also reflects the different normalization factors and count variability. The paired design supplies within-donor condition information, which is essential for distinguishing donor-specific response from baseline differences.
How Random Slopes Create Covariance
Repeated samples from one donor share \(b_{0i}\) and \(b_{1i}\). Even though the counts are conditionally independent given those effects, their marginal distributions are dependent. The following derivation makes the source and size of that dependence explicit. Consider two distinct samples \(a\) and \(c\) from donor \(i\), and temporarily absorb the offset and fixed predictors into \(q_a\) and \(q_c\). Their conditional means are \(\mu_a=\exp(q_a+b_{0i}+x_ab_{1i})\) and \(\mu_c=\exp(q_c+b_{0i}+x_cb_{1i})\).
Write \(v_x=\tau_{00}+2x\tau_{01}+x^2\tau_{11}\), the variance of \(b_{0i}+xb_{1i}\). For a zero-mean normal random variable \(U\) with variance \(v\), \(E[\exp(U)]=\exp(v/2)\). Thus \(E[Y_a]=E[\mu_a]=\exp(q_a+v_{x_a}/2)\), and similarly for \(Y_c\). The covariance between the two random linear predictors is
Conditional independence gives \(E[Y_aY_c\mid b_{0i},b_{1i}]=\mu_a\mu_c\). Applying the normal exponential moment to the sum of the two random predictors gives \(E[\mu_a\mu_c]=E[\mu_a]E[\mu_c]\exp(c_{ac})\). Therefore:
This is the count-scale covariance under the stated lognormal random-effects model, for two distinct observations from the same donor. The conditional negative-binomial variances contribute to each count's marginal variance, but not to this cross-observation covariance when the counts are conditionally independent. If the observations are from different donors, their covariance is zero under the donor-independence assumption.
Worked Example: Calculating Cross-Condition Covariance
Consider hypothetical fitted random-effect variances \(\tau_{00}=0.16\), \(\tau_{01}=-0.04\), and \(\tau_{11}=0.09\). For one donor, suppose the fixed predictors and offsets imply \(q_0=\log(100)\) for control and \(q_1=\log(120)\) for treatment. Calculate the model-implied marginal means and covariance between one control and one treated count.
For control, \(x=0\), so \(v_0=0.16\). Its marginal mean is \(E[Y_0]=100\exp(0.16/2)=100\exp(0.08)=108.3\) counts. For treatment, \(x=1\), so \(v_1=0.16+2(-0.04)+0.09=0.17\). Its marginal mean is \(E[Y_1]=120\exp(0.17/2)=120\exp(0.085)=130.6\) counts.
For the cross-condition pair, the latent-predictor covariance is \(c_{01}=0.16+(0+1)(-0.04)+(0)(1)(0.09)=0.12\). Hence
The positive covariance means that, under this fitted model, donors with random effects that elevate one condition's mean tend also to have elevated means in the other condition. The shared intercept contributes to that dependence; the random slope and intercept-slope covariance affect its size. This is a model-implied covariance, not a causal effect or a claim that every donor has the same covariance in observed counts.
When Can Slope Variation Be Estimated?
The design must contain information that separates donor-specific treatment response from donor baseline and count noise. With control and treatment observations within each of many donors, the repeated-condition pattern can inform the slope variance. A single condition observation per donor generally cannot separate a random intercept from a random slope: the observed donor effect is only \(b_{0i}+x_ib_{1i}\), and the data do not independently reveal both components. Merely having many donors does not repair a design that lacks within-donor condition variation.
Even a paired design does not guarantee precise estimation. Information about \(\tau_{11}\) depends on the number of independent donors, the balance and completeness of their condition observations, count depth, dispersion, and the magnitude of slope variation relative to sampling noise. There is no universal minimum donor count that makes the variance estimable well. More cells aggregated into a sample may improve its count precision, but they do not create additional independent donors.
Worked Example: Comparing Designs for Slope Estimation
Compare two hypothetical studies of a binary condition. Study A has 30 donors, each measured in control and treatment. Study B has 60 donors, with 30 observed only in control and 30 only in treatment. Which design directly supports estimation of donor-specific condition-slope variation?
In Study A, each donor contributes both \(x=0\) and \(x=1\). The donor's shared intercept appears in both observations, while the slope effect enters only the treatment linear predictor. Across independent donors, the paired pattern can therefore inform variation in the within-donor condition contrast, subject to count noise and the model assumptions.
In Study B, each donor contributes only one condition. A control donor contributes \(b_{0i}\); a treated donor contributes \(b_{0i}+b_{1i}\). The observed between-donor variation in the treated group combines intercept and slope variation, and the design has no within-donor contrast to separate them. Although Study B has twice as many donors, it does not directly identify \(\tau_{11}\) separately from baseline heterogeneity without additional design information or restrictive assumptions. Study A is the more informative design for this target.
Fitting, Diagnostics, and Alternatives
Random-effect variances are commonly estimated by likelihood-based methods, with approximations or numerical integration used for the non-Gaussian count likelihood. The resulting estimates should be interpreted alongside the design and diagnostics. A variance estimate at or near zero may indicate little detectable slope variation, but it can also reflect weak information or difficulty distinguishing the slope variance from other components. A near-boundary estimate, a singular fit, convergence warnings, or large changes under modest model or data perturbations call for caution. The earlier material on convergence and singularity in mixed models and on donor-sample influence diagnostics gives relevant checks.
Inspect whether donors contribute both conditions, whether condition labels and offsets are correct, and whether a small number of donors dominate the estimated variation. Also check the count-model fit: a random slope does not compensate automatically for a poor mean-variance specification, problematic fitted means, or influential counts. Uncertainty in a variance component can be substantial, especially with few donors; a point estimate alone is not evidence of a precisely measured biological spread of responses.
| Approach | What it represents | Key limitation |
|---|---|---|
| Random intercept only | Donor baseline heterogeneity with a shared condition coefficient | Imposes no donor-to-donor slope variation |
| Random intercept and condition slope | Donor baseline heterogeneity and a distribution of condition effects | Needs repeated condition information across enough independent donors; can be weakly estimated |
| Donor-specific fixed condition effects | A separate condition contrast for each donor | Uses many parameters and does not model a population distribution of slopes in the same way |
A random-intercept-only model may be preferable when the design cannot support slope variation, but it changes the assumption being made; it is not proof that donor effects are identical. Conversely, adding a random slope just because samples are repeated can produce unstable estimates. The choice should reflect the scientific target and the information in the design. A robust standard-error adjustment addresses a different issue and does not, by itself, estimate donor-specific treatment-effect variation.
Common Pitfall
Wrong reasoning: “The model has many pseudobulk samples, so it can estimate a random slope even though each donor appears in only one condition.” This confuses the number of rows with the independent within-donor information needed to separate baseline and response variation.
Correction: Count independent donors and inspect the condition pattern within each donor. A donor observed in only one condition provides no within-donor contrast. More cells or more rows cannot substitute for condition replication within donors when the target is donor-to-donor variation in treatment response.
Check Your Understanding
Answer each question, then compare with the solutions.
- In the model above, what does \(\tau_{11}\) represent, and what restriction removes random-slope variation?
- For two distinct samples from the same donor, why can their counts be correlated even when they are conditionally independent?
- For a control-treatment pair, write \(c_{01}\) in terms of the three random-effect covariance parameters.
- Why is a study with one condition per donor generally insufficient for separating slope variance from intercept variance?
- Name two design or fitting features that should prompt caution about a random-slope estimate.
Solutions. (1) \(\tau_{11}\) is the between-donor variance of condition-slope deviations \(b_{1i}\); setting \(\tau_{11}=0\) removes slope variation. (2) The observations share the same donor random effects, so their conditional means co-vary across the random-effects distribution. (3) With \(x_0=0\) and \(x_1=1\), \(c_{01}=\tau_{00}+\tau_{01}\). (4) Each donor contributes only \(b_{0i}+x_ib_{1i}\), not separate within-donor information on the intercept and slope. (5) Examples include few independent donors, incomplete within-donor condition pairs, a near-zero variance estimate, a singular fit, convergence warnings, or high sensitivity to individual donors.