Estimating Negative Binomial Dispersion

Why Estimate Dispersion?

For a gene, cell counts vary partly because of sampling and partly because the count distribution may be more variable than a Poisson model allows. The negative binomial model represents this extra variation through a dispersion parameter. Estimating that parameter determines whether a gene’s observed variability is consistent with the Poisson limit or supports a more variable count distribution. It is an estimate under a model, not a direct measurement of a single biological mechanism.

The earlier material on mean-dispersion relationships across genes describes how expected variability can depend on mean expression. Here the focus is the likelihood for one gene: what data enter it, how its dispersion estimate is obtained, and what it means when the maximum lies at the Poisson boundary.

Model and Likelihood

Suppose $Y_i$ is the count for one gene in cell $i$, for $i=1,\ldots,n$. Use the negative binomial parameterization with mean $\mu_i$ and dispersion $\alpha$, where $\alpha\geq 0$. Conditional on the means, assume the cell counts are independent and have the same dispersion. The mean may be common across cells, or specified by a mean model such as $\mu_i=s_i\lambda$, where $s_i$ is a known cell-specific exposure or size factor and $\lambda$ is a gene-level mean parameter. The mean model and any offsets must be specified adequately; the likelihood does not correct a misspecified mean model.

Definition: For $\alpha>0$, a negative binomial count with mean $\mu_i$ and dispersion $\alpha$ has probability mass function $$ P(Y_i=y_i)= \frac{\Gamma(y_i+1/\alpha)}{\Gamma(1/\alpha)y_i!} \left(\frac{1}{1+\alpha\mu_i}\right)^{1/\alpha} \left(\frac{\alpha\mu_i}{1+\alpha\mu_i}\right)^{y_i}, \qquad y_i=0,1,2,\ldots. $$ Its variance is $\mu_i+\alpha\mu_i^2$. The case $\alpha=0$ is defined by the limiting Poisson model, $Y_i\sim\operatorname{Poisson}(\mu_i)$.

The likelihood is the joint probability of the observed counts, viewed as a function of unknown parameters. For fixed means, independence makes it a product of the cell-specific probabilities. In applications, the means may also be unknown; maximum likelihood estimation then chooses the mean parameters and $\alpha$ together to maximize this product. Often the mean parameters are first maximized for each fixed $\alpha$, producing a profile likelihood in $\alpha$.

L(\alpha,\boldsymbol{\mu};\mathbf{y}) =\prod_{i=1}^{n}P(Y_i=y_i\mid \mu_i,\alpha), \qquad \ell(\alpha,\boldsymbol{\mu};\mathbf{y})=\log L(\alpha,\boldsymbol{\mu};\mathbf{y}).

The log likelihood is used because it turns the product into a sum and has the same maximizer as the likelihood. When the means vary by cell, maximizing usually requires numerical methods: derivatives of the gamma-function terms involve the digamma function. This is a computational detail, not a reason to omit the model assumptions. The counts should be conditionally independent under the fitted model, the mean structure should account for systematic differences between cells, and a shared $\alpha$ should be a reasonable summary of the gene’s residual variation.

Using the Size Parameter

It is sometimes simpler to express the same model using the size parameter $r=1/\alpha$. For a common mean $\mu$, the probability mass function becomes

P(Y=y)= \frac{\Gamma(y+r)}{\Gamma(r)y!} \left(\frac{r}{r+\mu}\right)^r \left(\frac{\mu}{r+\mu}\right)^y.

As $r$ increases without bound, $\alpha$ approaches zero and the distribution approaches a Poisson distribution. Thus, maximizing over $r>0$ can have its maximum only in the limit $r\to\infty$; in dispersion notation this is the boundary estimate $\hat{\alpha}=0$. A finite estimate of $r$ corresponds to positive dispersion. The two conventions point in opposite directions: large $r$ means little extra variation, while large $\alpha$ means more.

Worked Example: Evaluating a Candidate Likelihood

Consider the hypothetical counts $0,1,2$ for a gene in three cells. For illustration, take a common mean $\mu=1$ and candidate size $r=2$, so $\alpha=1/r=0.5$. The assumed variance is $1+0.5(1)^2=1.5$ counts squared. Under this model, $P(Y=0)=4/9$, $P(Y=1)=8/27$, and $P(Y=2)=4/27$.

Because the observations are assumed independent, the likelihood at this candidate parameter is

L(0.5;0,1,2) =\frac{4}{9}\times\frac{8}{27}\times\frac{4}{27} =\frac{128}{6561} \approx 0.01951.

This calculation evaluates the likelihood at one parameter value; it does not establish that $0.5$ is the maximum-likelihood estimate. To estimate dispersion, compare this likelihood with values across the allowed parameter space, accounting for estimation of the mean if it is not fixed. That distinction prevents a common error: treating a plausible candidate value as an estimate without maximizing.

A Boundary Estimate at the Poisson Limit

A boundary estimate occurs when the likelihood increases as dispersion decreases toward zero, so no positive value of $\alpha$ gives a larger likelihood. This is not a numerical failure by definition. It can be the correct constrained maximum when the observed counts do not show evidence of extra-Poisson variation, although a small sample may also provide little information for distinguishing small positive dispersion from zero.

Worked Example: Four Counts at the Poisson Boundary

Suppose a hypothetical gene has counts $1,1,1,1$, and use a common mean estimated from these data. The sample mean is $\bar{x}=(1+1+1+1)/4=1$. For size $r$, the negative binomial probability of one count at mean one is

P(Y=1)=\frac{r^{r+1}}{(r+1)^{r+1}}.

The likelihood and log likelihood for the four observations are therefore

L(r)=\left(\frac{r^{r+1}}{(r+1)^{r+1}}\right)^4, \qquad \ell(r)=4(r+1)\left[\log r-\log(r+1)\right].

Differentiate with respect to $r$:

\ell'(r) =4\left[\log\left(\frac{r}{r+1}\right)+\frac{1}{r}\right] =4\left[\frac{1}{r}-\log\left(1+\frac{1}{r}\right)\right].

For every $x>0$, $\log(1+x)<x$. Taking $x=1/r$ shows that $\ell'(r)>0$ for every finite $r>0$. Thus the likelihood keeps increasing as $r$ increases; its supremum is the Poisson limit $r\to\infty$, or $\hat{\alpha}=0$. This is a boundary maximum, not a finite size estimate. The arithmetic also shows why the factor $r^{r+1}$ matters: omitting $r^r$ changes the likelihood and its derivative, and can reverse the conclusion.

A Finite Dispersion Estimate

In the following example, a large count alongside several zeros creates more variation than the Poisson boundary can accommodate. The calculation uses a common mean estimated from the counts. For a common-mean negative binomial model, maximizing the likelihood over the mean at fixed $r$ gives $\hat{\mu}=\bar{x}$; this follows because the score for the mean is proportional to $\sum_i y_i-n\mu$, which is zero at the sample mean.

Worked Example: Maximizing for an Overdispersed Count Vector

Consider the hypothetical counts $0,0,0,4$. Their sample mean is $\bar{x}=1$. Let $r=1/\alpha$. At this mean, the likelihood is $P(Y=0)^3P(Y=4)$, where

P(Y=0)=\left(\frac{r}{r+1}\right)^r, \qquad P(Y=4)=\frac{r(r+1)(r+2)(r+3)}{24} \frac{r^r}{(r+1)^{r+4}}.

Taking logarithms and differentiating gives the score

\ell(r) =4r\log\left(\frac{r}{r+1}\right)-4\log(r+1) +\sum_{j=0}^{3}\log(r+j)-\log(24), \qquad \ell'(r)=4\log\left(\frac{r}{r+1}\right) +\sum_{j=0}^{3}\frac{1}{r+j}.

This score has no convenient closed-form root. Numerically, it is positive at $r=0.15$ (approximately $0.171$) and negative at $r=0.20$ (approximately $-0.567$). Numerical maximization gives a root near $\hat{r}=0.1582$, and checking the likelihood over the positive parameter range identifies it as the maximum rather than a minimum. Therefore

\hat{\alpha}=\frac{1}{\hat{r}} \approx\frac{1}{0.1582} \approx 6.321.

The fitted variance at mean one is $1+6.321(1)^2=7.321$ counts squared. The large dispersion estimate reflects how unusual the observed spread is relative to the fitted mean under this small sample; it does not establish why that spread occurred. The estimate is also likely to be imprecise with only four cells.

Optimization, Diagnostics, and Alternatives

A practical fitting procedure optimizes the profile log likelihood over $\alpha\geq0$, including the Poisson boundary. A numerical routine that only accepts strictly positive dispersion must approximate that boundary with a very small value; it should not silently label the approximation a positive finite estimate. Likewise, optimization over $r$ should allow the Poisson limit to be considered, even though it is not a finite value of $r$.

Specify the mean model.
Choose the common mean or cell-specific means, including justified offsets or covariates. Confirm that the fitted mean scale is appropriate for the counts being modeled.

Form the gene-level likelihood.
Use the negative binomial probabilities for the cells and sum their log probabilities. Treat dispersion as shared across these observations unless the model explicitly says otherwise.

Maximize, including the boundary.
Profile out mean parameters as appropriate, then compare the maximum over positive dispersion with the Poisson limit at $\alpha=0$.

Check the fitted model.
Compare observed variability and count frequencies with what the fitted model predicts. Look for systematic residual patterns, groups of cells with different means, or dependence not represented by the model.

The maximum-likelihood estimate is most meaningful when the assumed count distribution and mean structure are reasonable. A wrong mean model can make residual variation appear larger and inflate estimated dispersion. Dependence among cells, unmodeled cell groups, or a small number of informative nonzero counts can also make the estimate unstable. A fitted dispersion is not itself a diagnostic of model adequacy: inspect the fitted distribution against the data and consider whether one shared dispersion is plausible.

A moment-based estimate instead compares the observed variance with the model relation $\operatorname{Var}(Y)=\mu+\alpha\mu^2$, giving the rough estimate $(s^2-\bar{x})/\bar{x}^2$ for a common mean when $\bar{x}>0$. It is easy to calculate, but can be negative when the observed variance is below the mean and can be noisy in small samples. Constraining it to zero creates a boundary estimate, but does not make it a likelihood estimate. Maximum likelihood uses the full count probabilities, including the locations of zeros and large counts, rather than only the first two sample moments. The Poisson model is the limiting alternative with no extra dispersion; it is not a substitute for checking whether the mean structure is correct.

Common Pitfalls

Confusing $r$ with $\alpha$. Since $r=1/\alpha$, a large $r$ means dispersion close to zero, not large overdispersion.
Ignoring the boundary. A constrained maximum at $\alpha=0$ is a legitimate result. A software limit at a tiny positive number may be a numerical approximation to that boundary.
Using the likelihood for one candidate as the estimate. An evaluated likelihood is not a maximizer. Compare values across the parameter space or use a verified optimization routine.
Interpreting dispersion as a cause. Positive dispersion describes variation beyond the Poisson variance under the fitted model; it does not identify whether that variation is technical, biological, or due to misspecification.
Trusting an unstable estimate. Few nonzero counts or an inadequate mean model can yield highly uncertain estimates. Check fit and report the modeling assumptions.

Key takeaway: Estimate a gene’s negative binomial dispersion by maximizing its count likelihood under an explicit mean model. Include the Poisson boundary, interpret $r$ and $\alpha$ consistently, and treat diagnostics and model assumptions as part of the estimation.

Check Your Understanding

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

If $\alpha=0.25$, what is $r$, and what is the variance when $\mu=2$?
For four hypothetical counts all equal to one with common mean one, why is the dispersion estimate at the Poisson boundary?
In the count vector $0,0,0,4$, what is the estimated common mean, and what does a finite dispersion estimate describe?
Name two assumptions that matter when estimating a shared gene-level dispersion across cells.
How does a moment estimate of dispersion differ from a maximum-likelihood estimate?

Solutions. (1) $r=1/0.25=4$, and the variance is $2+0.25(2^2)=3$ counts squared. (2) Its log likelihood as a function of $r$ has derivative $4[1/r-\log(1+1/r)]>0$, so it increases toward the Poisson limit as $r\to\infty$, equivalently $\alpha\to0$. (3) The common mean estimate is $(0+0+0+4)/4=1$. A finite estimate describes extra variation relative to the Poisson variance under the assumed model; it does not identify the source of that variation. (4) Examples are conditional independence of cell counts, an adequate mean model including relevant offsets, and a shared dispersion that reasonably summarizes the cells. (5) A moment estimate uses sample mean and variance, while maximum likelihood maximizes the full negative binomial count likelihood; the estimates need not agree, particularly in small samples.