Cell-Level Negative Binomial Mixed Models

Why Model Donors When the Data Are Cells?

A single-cell dataset may contain thousands of cells, but cells collected from the same donor share biological and technical features. Their gene counts are therefore not generally independent observations. If a condition is assigned to donors, the independent replication for comparing conditions comes from donors, not from the number of cells sequenced. A model that treats all cells as independent can mistake a large cell count for a large number of independent experimental units.

A cell-level negative binomial mixed model retains each cell’s count while representing shared donor-to-donor differences with a random effect. The introductory material on Pseudobulk Versus Cell-Level Mixed Models contrasts the broad data summaries; here the focus is the cell-level model’s specification and how its shared donor effect induces dependence.

Model Definition and Assumptions

For one gene at a time, let $Y_{ijg}$ be the observed count for gene $g$ in cell $j$ from donor $i$. Let $s_{ij}$ be a known, positive library-size or normalization factor for that cell, and let $\mathbf{x}_{ij}$ be a vector of observed covariates, including the condition indicator when appropriate. The donor random intercept for gene $g$ is $b_{ig}$.

Definition: A cell-level negative binomial mixed model with a donor random intercept specifies $$ Y_{ijg}\mid b_{ig},\mathbf{x}_{ij} \sim \operatorname{NB}(\mu_{ijg},\phi_g), \qquad \log(\mu_{ijg})=\log(s_{ij})+\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}_g+b_{ig}, \qquad b_{ig}\sim N(0,\tau_g^2). $$ Here the negative binomial parameterization has conditional mean $\mu_{ijg}$ and conditional variance $\mu_{ijg}+\phi_g\mu_{ijg}^2$, with dispersion $\phi_g\geq 0$. Conditional on the donor effects and modeled covariates, cell counts are assumed independent. The donor effects are assumed independent across donors and independent of the covariates.

The logarithm of the size factor is an offset: its coefficient is fixed at one, so a doubling of $s_{ij}$ doubles the modeled mean when the other terms are held fixed. The coefficient vector $\boldsymbol{\beta}_g$ describes gene-specific fixed effects. The random intercept shifts the log mean for every modeled cell from donor $i$ by the same donor-specific amount.

The model is conditional on $b_{ig}$. It assumes that, once the modeled covariates and donor effect are given, remaining cell-level variation follows the negative binomial variance relationship and counts from different cells are conditionally independent. It also assumes the random intercepts follow the stated normal distribution and do not depend on condition or other modeled covariates. These assumptions require a meaningful offset, a suitable fixed-effect design, and enough independent donors to estimate population-level effects. A large number of cells from only a few donors does not substitute for donor replication.

How a Shared Random Intercept Creates Correlation

Consider two different cells $j$ and $k$ from the same donor for the same gene. Write $a_j=s_{ij}\exp(\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}_g)$ and $a_k=s_{ik}\exp(\mathbf{x}_{ik}^{\mathsf T}\boldsymbol{\beta}_g)$, so their conditional means are $a_j\exp(b_{ig})$ and $a_k\exp(b_{ig})$. Conditional independence does not imply marginal independence: after the shared random effect is integrated out, both counts vary together with $b_{ig}$.

\begin{aligned} \operatorname{Cov}(Y_{ijg},Y_{ikg}) &=\operatorname{Cov}\bigl(E[Y_{ijg}\mid b_{ig}],E[Y_{ikg}\mid b_{ig}]\bigr)\\ &=a_j a_k\left\{E[\exp(2b_{ig})]-E[\exp(b_{ig})]^2\right\}\\ &=a_j a_k\left\{\exp(2\tau_g^2)-\exp(\tau_g^2)\right\}\\ &=m_jm_k\left\{\exp(\tau_g^2)-1\right\}, \end{aligned}

In the last line, $m_j=E[Y_{ijg}]=a_j\exp(\tau_g^2/2)$ and $m_k=a_k\exp(\tau_g^2/2)$, using the moment-generating function of a normal random variable. Thus, when $\tau_g^2>0$, the shared intercept gives positive covariance between counts from the same donor. Counts from different donors have covariance zero under the model’s independence assumptions. The negative binomial dispersion describes additional conditional count variability; it does not, by itself, represent the shared donor effect.

Worked Example: Covariance Between Two Cells

Suppose that for a hypothetical gene, two cells from one donor have conditional means $a_j=4$ and $a_k=8$ when the donor effect is zero, and the random-intercept variance is $\tau_g^2=0.25$. Find the marginal covariance induced by the shared donor effect.

The two marginal means are $m_j=4\exp(0.25/2)\approx 4.533$ counts and $m_k=8\exp(0.25/2)\approx 9.065$ counts. Therefore,

\begin{aligned} \operatorname{Cov}(Y_{ijg},Y_{ikg}) &=m_jm_k\{\exp(0.25)-1\}\\ &\approx (4.5327)(9.0653)(1.284025-1)\\ &\approx 11.67\ \text{counts}^2. \end{aligned}

Equivalently, the covariance is $4(8)\{\exp(0.50)-\exp(0.25)\}\approx 11.67\ \text{counts}^2$. This is covariance, not a correlation coefficient; a correlation would also require the marginal variances of both counts.

Interpreting and Calculating the Conditional Mean

The log link ensures positive expected counts. For two conditions coded $x_i=0$ and $x_i=1$, a simple fixed-effect component is $\beta_{0g}+\beta_{1g}x_i$. With other covariates held fixed and equal size factors, $\exp(\beta_{1g})$ is the condition-associated ratio of conditional mean counts. The same ratio applies to the marginal means under this random-intercept model because the shared random-effect distribution contributes the same multiplicative factor to each condition’s expectation. This is an adjusted association, not evidence that condition causally changes expression.

Worked Example: Computing a Conditional Fitted Mean

For a hypothetical gene and cell, suppose the intercept is $1.000$, the condition coefficient is $0.500$, the log size-factor offset is $0.100$, and the fitted donor effect is $0.0986$. The cell is in the condition coded $x_i=1$. Its linear predictor is

\eta_{ijg}=0.100+1.000+(0.500)(1)+0.0986=1.6986.

Exponentiating gives the conditional fitted mean:

\mu_{ijg}=\exp(1.6986)\approx 5.4663\ \text{counts}.

Thus the fitted expected count for this cell, conditional on its donor effect and modeled covariates, is approximately $5.466$ counts. If the gene’s dispersion is $\phi_g=0.200$, the conditional variance is

\operatorname{Var}(Y_{ijg}\mid b_{ig}) =\mu_{ijg}+\phi_g\mu_{ijg}^2 \approx 5.4663+0.200(5.4663)^2 \approx 11.44\ \text{counts}^2.

This variance is conditional on the specified donor effect. It is not the marginal variance across donors, which also includes variation in $b_{ig}$.

Worked Example: Interpreting a Condition Coefficient

Suppose a hypothetical fitted model has condition coefficient $\widehat{\beta}_{1g}=0.4055$ for a gene, with condition coded 0 or 1. Compare the modeled conditional means for two otherwise identical cells from the same donor, one in each condition.

Their log-mean difference is $(1-0)\widehat{\beta}_{1g}=0.4055$. The size factor, other covariates, and shared donor intercept cancel in the ratio, so

\frac{\mu_{ijg}(x_i=1)}{\mu_{ijg}(x_i=0)} =\exp(0.4055)\approx 1.500.

The model therefore estimates a conditional mean count about $1.500$ times as large in condition 1 as in condition 0, for equal offsets and covariates. This is a 50.0% higher modeled mean count, not a statement that every cell has 50.0% more counts. It is also not a test of whether the difference is statistically distinguishable from zero.

Specification, Comparisons, and Failure Modes

A random intercept is appropriate when cells from a donor may share an overall gene-specific expression level beyond what the fixed covariates explain. Its inclusion preserves cell-level observations while explicitly modeling that shared donor contribution. The design must still identify the fixed effects: for example, a condition perfectly confounded with batch cannot be separated from batch merely by adding a donor effect.

Compared with the main alternatives, the modeling choices have different consequences:

Approach	What it represents	Important distinction
Cell-level negative binomial mixed model	Cell counts, count overdispersion, and a shared donor intercept	Uses cells for conditional information while retaining donor-level dependence.
Pseudobulk negative binomial model	Counts aggregated within donor and sample, with donors as replicates	Summarizes cells before fitting; does not estimate a separate cell-level mean.
Cell-level negative binomial model without donor effects	Cell counts and count overdispersion, but no shared donor term	Can treat correlated cells as independent and understate uncertainty for donor-level comparisons.

A random intercept is not a universal dependence model. It represents a shared shift in log mean, not every possible pattern of correlation between cells, samples, or conditions. The negative binomial variance also does not repair a wrong experimental unit or omitted design structure. A model fit should be questioned if it fails to converge, estimates a donor variance at a boundary unexpectedly, or shows systematic disagreement between observed and fitted count patterns. Such signs motivate checking the model specification and data structure rather than simply reporting its estimates.

A common mistake is to argue, “There are 2,000 cells in each condition, so the comparison has 2,000 independent observations per condition.” That reasoning ignores the shared donor effect when those cells come from a smaller set of donors. The correction is to model donor dependence and base condition-level replication on the donors or independent samples that received the conditions. More cells can improve characterization of a donor’s cells, but cannot create additional independent donors.

Define the count and nesting.
Identify the gene-specific cell count and which donor each cell belongs to; specify the condition and covariates at the level where they were measured.

Specify the mean and dispersion.
Use a log link with a fixed-coefficient offset for the cell’s size factor, fixed effects for the intended comparisons, and a negative binomial conditional variance.

Add the donor random intercept.
State its distribution and shared role across cells from the same donor, then verify that the design and number of donors can support the intended inference.

Check Your Understanding

Answer each question, then compare with the solutions below.

Under the stated model, what is the conditional variance of a cell count with mean $8$ and dispersion $0.25$?
Why can two distinct cells from the same donor be dependent even though their counts are conditionally independent given the donor effect?
A condition coefficient is $-0.2231$. What is the modeled condition 1 to condition 0 mean ratio, with equal offsets and covariates?
Does adding 500 more cells from each of four donors create 500 additional independent donor replicates per condition? Explain.

Solutions. (1) The conditional variance is $8+0.25(8^2)=8+16=24\ \text{counts}^2$. (2) Both counts depend on the same random intercept; integrating it out produces positive covariance when its variance is positive. (3) The ratio is $\exp(-0.2231)\approx 0.8000$, so the modeled conditional mean is about 20.0% lower in condition 1. (4) No. The added cells are still nested within the same four donors; they may improve cell-level information but do not add independent donor replicates.