Why Add a Donor Random Intercept?
When a donor contributes pseudobulk observations under more than one condition, those observations are repeated measurements, not independent biological replicates. A donor may have a consistently higher baseline expression level for a gene than another donor. A model that represents this donor-specific baseline can account for the resulting within-donor dependence while estimating an average condition effect across donors.
For a given gene and cell type, let \(Y_{di}\) be the summed pseudobulk count for donor \(d\) in observation \(i\). Let \(x_{di}\) indicate condition, and let \(s_{di}\) be the sample's library-size normalization factor. This tutorial uses a negative-binomial mixed model for counts. A related random-intercept model can be used for suitably transformed expression, but its residual distribution and interpretation must be stated separately.
For a two-condition comparison, code the reference condition as \(x_{di}=0\) and the comparison condition as \(x_{di}=1\). The coefficient \(\beta_1\) is the change in log conditional expected count, comparing conditions for the same donor and the same library-size factor. Thus \(\exp(\beta_1)\) is the conditional fold change. The library-size offset adjusts the expected count for sequencing depth; it is not an estimated condition effect.
The conditional-independence assumption is that, given the donor intercept and modeled covariates, the observations have the specified negative-binomial distributions and do not retain further dependence. The shared intercept accounts for one particular form of within-donor association: donors differ in their baseline level, but all observations from a donor share the same shift on the log-mean scale. The model does not, by itself, represent every possible repeated-measure covariance pattern.
Writing and Interpreting the Model
A donor intercept shifts the log mean for all of that donor's observations. If \(u_d\) is positive, that donor's conditional expected counts are higher than the fixed-effects part of the model predicts; if it is negative, they are lower. The condition coefficient is shared across donors. It describes a common condition contrast, rather than a separate condition effect for each donor.
A useful consequence is that the fixed condition coefficient is a within-donor contrast in a paired design, while the model also estimates a population distribution of donor baselines. This differs from fitting a separate unrestricted fixed coefficient for every donor. The random-intercept approach assumes donor effects are draws from a common distribution and independent of condition assignment, whereas a donor fixed-effects model conditions on the particular donors and allows each its own unrestricted intercept.
The coefficient is a conditional log fold change, but under the stated random-intercept model it also gives the marginal fold change when comparing conditions with the same offset and the same donor-effect distribution. To see this, integrate over \(u_d\):
The final equality uses \(u_d\sim N(0,\tau^2)\), for which \(E\{\exp(u_d)\}=\exp(\tau^2/2)\). At a common offset, the ratio of marginal expected counts for \(x=1\) versus \(x=0\) is therefore \(\exp(\beta_1)\): the shared multiplier cancels. This result relies on the random-intercept distribution being the same across conditions and independent of the modeled condition.
Worked Example: Specifying a Paired Count Model
Suppose three hypothetical donors each contribute a control and a treated pseudobulk sample for one gene in one cell type. The observed counts are 12 and 22 for donor A, 8 and 13 for donor B, and 15 and 20 for donor C, respectively. For each observation, the model uses its library-size factor \(s_{di}\). Write the model and explain the comparison it targets.
Set \(x_{di}=0\) for control and \(x_{di}=1\) for treated. The model for all six observations is
For donor A, the conditional expected control count is \(s_{A0}\exp(\beta_0+u_A)\); the treated expected count is \(s_{A1}\exp(\beta_0+\beta_1+u_A)\). Their ratio after accounting for the two offsets is
The same cancellation holds for donors B and C. Thus the condition contrast is not the raw ratio \(22/12\) for donor A, nor the average of the three raw count ratios. It is the model-based comparison adjusted for the library-size factors and shared donor baselines. The six observed counts are illustrative data; they do not, on their own, supply fitted values for \(\beta_1\), \(\phi\), or \(\tau^2\).
Worked Example: Interpreting a Hypothetical Condition Estimate
Suppose a fitted model gives \(\hat{\beta}_1=0.40\) with standard error \(SE=0.16\). For this illustration, use a normal-reference Wald test of \(H_0:\beta_1=0\) against \(H_a:\beta_1\ne 0\), and a 95% Wald interval. These are hypothetical model results, not results from the counts in the previous example.
The Wald statistic and approximate two-sided p-value are
At a 0.05 significance level, this test provides evidence against a zero condition coefficient under the model and the normal-reference approximation. It does not give the probability that the null hypothesis is true.
The estimated fold change is \(\exp(0.40)=1.492\), or about 49.2% higher expected expression in the treated condition than in control, at a common library-size factor. The 95% interval on the log scale is
Exponentiating both limits gives \((\exp(0.0864),\exp(0.7136))=(1.090,\;2.041)\). Under the model and the interval procedure, the data are compatible with a treated-to-control fold change from about 1.09 to 2.04. This interval describes uncertainty in the condition effect; it is not a 95% probability statement about a fixed parameter being inside this particular interval.
What Dependence Does the Intercept Represent?
The common donor term creates association because two observations from the same donor include the same \(u_d\). On a Gaussian transformed-expression scale, this can be seen directly with the simpler model \(Z_{di}=\beta_0+\beta_1x_{di}+u_d+\varepsilon_{di}\), where \(u_d\sim N(0,\tau^2)\), \(\varepsilon_{di}\sim N(0,\sigma^2)\), and residuals are independent of each other and of \(u_d\). For two distinct observations from the same donor, the shared term contributes to their covariance:
This calculation is for the stated Gaussian model, not an ICC formula to apply automatically to negative-binomial counts. In the count model, the shared random intercept likewise induces dependence after it is integrated out, but the marginal count distribution and correlation have a different form.
Worked Example: Checking the Implied Gaussian-Scale Correlation
For a hypothetical transformed-expression random-intercept model, suppose the fitted donor-intercept variance is \(\tau^2=0.36\) and the residual variance is \(\sigma^2=0.64\), in squared log-expression units. The implied correlation between two observations from the same donor is
The model therefore implies a within-donor correlation of 0.36 on this Gaussian scale. This is an implication of the assumed model, not proof that the covariance pattern is adequate for the data. The value should not be transferred to the negative-binomial count scale without a model-specific calculation.
Assumptions, Diagnostics, and Alternatives
A random-intercept analysis is most defensible when the design and data support its core assumptions. In particular, donor effects should be plausibly represented by a common distribution, be independent of condition and included covariates, and act as a constant shift on the log-mean scale. The condition effect is assumed common across donors unless interactions or random slopes are added. The count model also requires an appropriate negative-binomial mean-variance relationship and suitable library-size offsets.
Check the design before fitting: observations must be correctly linked to donors, the condition coding and reference level must be clear, and the fixed condition effect must be identifiable. In a paired design, donors should contribute information on both conditions for the within-donor comparison. If condition is constant within each donor, the condition contrast is between donors rather than a paired within-donor contrast; it then depends on comparisons across the donor groups and the assumption that their random effects are independent of condition.
Warning signs include convergence failures, estimated variance components at a boundary, implausibly large or unstable standard errors, strong residual patterns, and systematic donor-specific condition contrasts. These can indicate too little information to estimate the model, an unsuitable mean or variance specification, or a condition effect that varies among donors. A boundary estimate of \(\tau^2\) does not prove that all donors are identical; it can also reflect limited information. Diagnostics should assess both the count fit and whether the shared-intercept structure is a plausible representation of repeated observations.
The main alternatives answer related but not identical modeling questions:
- Donor fixed effects: give each observed donor an unrestricted intercept. This avoids modeling donor intercepts as draws from a normal distribution, but does not estimate a population distribution of donor baselines. Donor-level covariates that do not vary within donor cannot be separately estimated alongside unrestricted donor indicators.
- Independent pseudobulk model: omits donor dependence. Treating repeated observations as independent can understate uncertainty because it fails to account for shared donor baselines.
- Random intercept plus random slope: allows the condition effect to vary across donors as well as their baseline levels. This is a different, more parameter-intensive model and requires enough repeated information to estimate the additional variation.
A random intercept is not a universal solution for repeated measurements. It imposes one particular dependence structure, and it does not account for arbitrary time patterns, condition-specific residual correlations, or unmodeled technical effects. The variance component itself is a model parameter whose interpretation and inferential details require care; a small number of donors can make its estimation especially uncertain.
Check Your Understanding
Answer each question, then compare with the short solutions.
- In the model above, what is the treated-to-control conditional fold change if \(\beta_1=-0.25\)?
- Why does a shared donor intercept induce dependence between two observations from the same donor?
- Suppose each donor contributes only a treated sample or only a control sample. Is the condition coefficient a within-donor contrast?
- For the Gaussian random-intercept model, calculate the implied within-donor correlation when \(\tau^2=0.50\) and \(\sigma^2=1.50\).
Solutions. (1) \(\exp(-0.25)=0.7788\), so the treated expected count per common offset is about 0.779 times the control expected count, or about 22.1% lower. (2) Both observations contain the same random quantity \(u_d\), so their outcomes share a source of variation even if their conditional residuals are independent. (3) No. Condition varies between donors, so the contrast is between donor groups, not within donors; its interpretation relies on the model assumptions about donor effects and condition. (4) The correlation is \(\tau^2/(\tau^2+\sigma^2)=0.50/(0.50+1.50)=0.25\).