Dispersion Estimation Across Biological Samples

Why Dispersion Is Estimated Across Samples

For a gene in a cell type, pseudobulk analysis gives one count for each biological sample by aggregating the cells from that sample. The negative-binomial model allows those counts to vary more than a Poisson model permits. Dispersion quantifies that extra variation relative to the modeled mean. It matters because uncertainty in estimated expression differences depends not only on the expected counts, but also on how much biological samples vary around those expectations.

The relevant replication is the set of independent biological samples, not the set of cells pooled into each pseudobulk. More cells can make a sample’s aggregate count more informative about that sample, but they do not create additional independent observations of between-sample variation. Dispersion is therefore estimated from the pattern of counts across biological samples, after accounting for the design and expected mean.

Definition: For gene $g$ and biological sample $i$, let $Y_{gi}$ be the pseudobulk count and $\mu_{gi}=E(Y_{gi})$ its modeled mean. The gene-wise dispersion $\alpha_g$ is the nonnegative parameter in the negative-binomial variance model $$ \operatorname{Var}(Y_{gi}\mid X_i)=\mu_{gi}+\alpha_g\mu_{gi}^{2}. $$ The term $\mu_{gi}$ is the Poisson sampling component; $\alpha_g\mu_{gi}^{2}$ is the additional variance allowed by the model. Dispersion is dimensionless. A value of zero gives the Poisson variance, while a larger value gives greater variance at the same mean.

Model and Assumptions

In a pseudobulk negative-binomial generalized linear model, the mean for a gene is related to sample-level predictors through a log link. With library-size offset $o_i$ and design row $x_i$, a common specification is:

\log(\mu_{gi})=o_i+x_i^\mathsf{T}\beta_g, \qquad \operatorname{Var}(Y_{gi}\mid X_i)=\mu_{gi}+\alpha_g\mu_{gi}^{2}.

Here $\beta_g$ contains the gene’s model coefficients, and $o_i$ is the offset established from the pseudobulk library size. The offset adjusts the expected count for differences in library scale; it is not a biological replicate and does not directly estimate dispersion. For a fixed gene, the model uses one dispersion parameter across its samples, while the means may differ with condition, covariates, or offsets.

The standard likelihood approach assumes that, conditional on the design and offsets, the pseudobulk counts follow the specified negative-binomial distribution and are independent across independent biological samples. It also requires a mean model that adequately represents the design: if important sample-level structure is omitted, apparent extra variation can be attributed to dispersion. Repeated measurements or other dependence between samples require a model that represents that dependence; treating correlated observations as independent can distort dispersion estimates and uncertainty.

The negative binomial is often parameterized by a size parameter $r=1/\alpha_g$. Under this parameterization its probability mass function is:

P(Y_{gi}=y)= \frac{\Gamma(y+r)}{\Gamma(r)y!} \left(\frac{r}{r+\mu_{gi}}\right)^r \left(\frac{\mu_{gi}}{r+\mu_{gi}}\right)^y, \qquad y=0,1,2,\ldots .

For $\alpha_g>0$, a likelihood estimate chooses the dispersion (and mean-model coefficients) that best fits the observed counts under this distribution. In practice, the coefficients are re-estimated or profiled for candidate dispersion values. The estimate is driven by how far the observed sample counts vary from the fitted means, relative to the variation already expected from their count magnitudes.

Why the Variance Has This Form

A short derivation shows what the dispersion parameter represents. Suppose a count is Poisson conditional on a latent rate $\Lambda$, and that rate varies across samples according to a gamma distribution with mean $\mu$ and variance $\alpha\mu^2$. Then the law of total variance gives:

\begin{aligned} \operatorname{Var}(Y) &=E\{\operatorname{Var}(Y\mid\Lambda)\} +\operatorname{Var}\{E(Y\mid\Lambda)\}\\ &=E(\Lambda)+\operatorname{Var}(\Lambda)\\ &=\mu+\alpha\mu^2. \end{aligned}

Thus the negative-binomial variance can be understood as ordinary count noise plus variation in the latent rate. This interpretation motivates the model, but it does not mean that a fitted $\alpha_g$ isolates a single biological mechanism. It summarizes extra variation among the modeled pseudobulk observations, which can reflect genuine sample heterogeneity as well as unmodeled structure or model inadequacy.

A Moment Estimate and Its Limits

A simple method-of-moments calculation illustrates the information used to estimate dispersion. If $n$ independent samples have the same mean $\mu$, the variance relationship suggests estimating $\alpha_g$ by comparing the observed sample variance $s^2$ with the Poisson variance $\bar{x}$:

\widehat{\alpha}_{\mathrm{MOM}} = \max\left(0,\frac{s^2-\bar{x}}{\bar{x}^{\,2}}\right).

The truncation at zero reflects the model restriction that dispersion cannot be negative. This is an illustrative estimator, not a universal replacement for likelihood estimation: real designs have different fitted means, offsets, and estimated coefficients. A likelihood fit uses those means directly. In either approach, with few biological samples the observed variability is noisy, so a gene-wise estimate can be unstable.

Worked Example: Estimating Dispersion From Replicate Counts

Suppose four independent hypothetical biological samples have pseudobulk counts 80, 100, 120, and 140 for a gene. Assume they are from one condition, have equal offsets, and share a common mean. These assumptions make the simple moment calculation appropriate for illustration.

The mean is $\bar{x}=(80+100+120+140)/4=110$. The sample variance is:

s^2=\frac{(80-110)^2+(100-110)^2+(120-110)^2+(140-110)^2}{4-1} =\frac{2000}{3}=666.7.

Therefore:

\widehat{\alpha}_{\mathrm{MOM}} = \frac{666.7-110}{110^2} = \frac{556.7}{12100} = 0.0460.

The estimate is positive, so the observed replicate variance exceeds the Poisson variance of 110 counts. Under this fitted moment model, the variance at mean 110 is $110+0.0460(110^2)\approx666.6$ count-squared units. The calculation does not establish why the samples vary, and with only four samples the estimate is uncertain.

Offsets and Fitted Means Matter

Raw count variation is not itself a dispersion estimate. Samples with different library sizes can have different expected counts even when their expression rate is the same. Dispersion estimation must compare each observation with its model-based mean, which incorporates the offset and design. A calculation that treats unequal expected counts as though they shared one mean can mistake the mean pattern for extra variation.

Worked Example: Avoiding a False Dispersion Signal From Library Size

Consider four hypothetical samples with relative library exposures $0.5,1,1.5,$ and $2$. Their observed counts are $50,100,150,$ and $200$, respectively. Suppose the fitted expression rate is 100 counts per unit of exposure. The model means are therefore $50,100,150,$ and $200$, matching the observations.

The raw-count mean is 125 and the raw sample variance is:

s^2= \frac{(-75)^2+(-25)^2+25^2+75^2}{3} = \frac{12500}{3} = 4166.7.

Using the equal-mean moment formula on these raw counts would give $(4166.7-125)/125^2=0.2587$. That number is not a valid dispersion estimate here: the samples do not share a common expected count. The large raw variance is explained by the different exposures.

With the offset-based fitted means, all residuals $Y_{gi}-\mu_{gi}$ are zero in this constructed example. There is no observed excess residual variation to support positive dispersion; a constrained moment calculation would be at the zero boundary. Real datasets will not generally have exact fits, but the principle remains: estimate dispersion from variation around fitted, offset-adjusted means, not from raw counts alone.

Biological Samples, Cells, and Information

Within a pseudobulk sample, aggregating more cells can increase the count and may reduce the influence of cell-level sampling noise on the aggregate. However, dispersion across samples is learned from differences among independent sample-level observations. If a study has four biological samples, it has four independent count observations for a gene in a given condition, whether each aggregate contains 500 cells or 50,000 cells.

The mean–variance formula also explains why a larger count does not automatically imply a more precise biological comparison. As $\mu$ grows, the Poisson component grows as $\mu$, but the extra component grows as $\alpha_g\mu^2$. When dispersion is appreciable, variation across samples can remain substantial even for large pseudobulk counts. Adding cells cannot be treated as adding independent biological replicates or as a way to divide the sample-level variance by the number of cells.

Worked Example: More Cells Do Not Create More Replicates

Suppose a gene has mean pseudobulk count $\mu=100$ and dispersion $\alpha_g=0.04$ across independent samples. The modeled variance per sample is:

$$ 100+0.04(100^2)=500. $$

For four independent biological samples with this mean and variance, the standard error of their average count under this simple model is $\sqrt{500/4}=11.18$ counts. If each sample’s aggregate instead contains ten times as many cells, the number of independent sample observations is still four. The calculation does not become $\sqrt{500/40}$: there are not 40 biological replicates.

More cells may change the counts, their expected means, and the contribution of cell-level sampling to the aggregate. Those changes should be represented through the data and mean model. They do not change the number of independent biological samples available to estimate between-sample variation.

Diagnostics, Alternatives, and Failure Modes

A fitted dispersion should be assessed alongside the mean model, not interpreted as an isolated biological quantity. Useful checks include comparing observed counts with fitted means, examining residual patterns across samples and conditions, and checking whether a small number of samples dominate the variability. Strong residual structure by condition, batch, or another known sample-level factor can indicate that the mean model is incomplete. Extra zeros, outliers, or dependence may also make a simple negative-binomial model inadequate.

The main approaches have different roles:

Approach	How it uses variation	Main caution
Method of moments	Matches observed variability to the negative-binomial variance formula.	Simple calculations can be noisy and need adjustment for unequal means and fitted parameters.
Negative-binomial likelihood	Fits the count distribution while accounting for the design and offsets.	With few samples or weak information, a gene-wise estimate can be unstable.
Poisson model	Sets dispersion to zero, so variance equals the mean.	Can understate uncertainty when samples vary more than Poisson counts allow.
Quasi-likelihood approach	Uses a mean–variance relationship without requiring the full negative-binomial probability model.	Its assumptions and inferential procedure differ from a negative-binomial likelihood fit.

A common failure is to estimate dispersion from a single sample or from thousands of cells as if they were independent biological replicates. Across-sample dispersion cannot be learned from one independent sample, and cell-level counts do not substitute for sample-level replication. Another failure is interpreting a large estimate as proof of biological heterogeneity: omitted covariates, technical differences, outliers, or a misspecified mean can also produce excess residual variation.

Gene-wise estimates based on few samples can be erratic. Stabilizing information across genes is a separate step from defining or estimating the raw gene-wise dispersion; it does not remove the need for biological replication or correct model specification. If residual diagnostics show systematic departures, simply accepting a larger dispersion may conceal a poorly specified model.

Key takeaway: Gene-wise dispersion describes excess count variation around fitted pseudobulk means across biological samples. Estimate it using independent sample-level observations and a mean model that accounts for design and library-size offsets; the number of cells is not the number of biological replicates.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

For a gene with fitted mean 50 and dispersion 0.02, what variance does the negative-binomial model imply?
Four equal-exposure sample counts are 40, 50, 60, and 50. Calculate the simple method-of-moments dispersion estimate.
Why is raw-count variance not a valid dispersion estimate when samples have different library-size offsets?
A pseudobulk study has six independent biological samples, each containing many cells. How many independent sample-level observations contribute per gene, and why can’t the cells be counted as additional replicates?
What might a large estimated dispersion indicate besides genuine biological sample heterogeneity?

Solutions. (1) The variance is $50+0.02(50^2)=50+50=100$ count-squared units. (2) The mean is 50 and the sample variance is $\{(-10)^2+0^2+10^2+0^2\}/3=200/3=66.67$. Thus $\widehat{\alpha}_{\mathrm{MOM}}=(66.67-50)/50^2=0.00667$. (3) Different offsets imply different expected counts even at the same expression rate, so raw variation combines mean differences with residual variation. (4) There are six independent sample-level observations per gene. Cells contribute to each aggregate, but cells from the same biological sample are not independent biological replicates. (5) Omitted sample-level factors, technical differences, outliers, dependence, or a misspecified mean or count distribution can also produce excess variation.