Negative-Binomial Deviance for Feature Ranking

Why Rank Genes by Negative-Binomial Deviance?

A gene’s raw count variation across cells is difficult to interpret on its own. Counts tend to increase with cell library size, and their variability also changes with their mean. The earlier material on mean-variance relationships and regularized dispersion addresses these features of count data. Negative-binomial deviance offers another way to rank genes: it measures how much better a cell-specific fit describes a gene’s counts than a null fit in which the gene has one common expression rate after accounting for library size.

The resulting score is a measure of departure from that null model, not a test result by itself. A large score can make a gene a useful candidate feature, but it does not establish that the gene is biologically important or identify the cause of its variation.

Model and Definition of Deviance

For gene $g$ in cell $i$, let $Y_{ig}$ be the observed count and let $s_i>0$ be a fixed library-size exposure factor. Assume counts are conditionally independent across cells and follow a negative-binomial distribution with mean $\mu_{ig}$ and dispersion $\alpha_g>0$:

Y_{ig}\mid s_i \sim \operatorname{NB}(\mu_{ig},\alpha_g), \qquad \operatorname{Var}(Y_{ig}\mid s_i)=\mu_{ig}+\alpha_g\mu_{ig}^{2}.

Equivalently, the negative-binomial size parameter is $r_g=1/\alpha_g$. The exposure factors are treated as known for this fit, and the dispersion is held fixed while the means are fitted. It may be supplied by a gene-specific regularized estimate; deviance ranking does not itself remove the need to choose a defensible dispersion. The model also assumes the offsets account adequately for library-size differences and that the cells being compared are appropriate to analyze together.

The null fit assigns gene $g$ one rate $q_g$ shared across cells, with mean $\mu_{ig}=s_iq_g$. The saturated fit gives each cell its own mean, fitted to that cell’s observed count. Under this fit, $\hat{\mu}_{ig}=Y_{ig}$, including the boundary value zero for a zero count.

Definition: The negative-binomial deviance for gene $g$ is twice the log-likelihood difference between its saturated fit and its null fit. With $r_g=1/\alpha_g$, its cell-level contribution, evaluated at null mean $\hat{\mu}_{ig}=s_i\hat q_g$, is $$ d_{ig}=2\left[ Y_{ig}\log\left(\frac{Y_{ig}}{\hat{\mu}_{ig}}\right) -(Y_{ig}+r_g)\log\left(\frac{Y_{ig}+r_g}{\hat{\mu}_{ig}+r_g}\right) \right], \qquad D_g=\sum_i d_{ig}. $$ A term $0\log(0/a)$ is defined as zero. Thus a zero observed count has contribution $2r_g\log(1+\hat{\mu}_{ig}/r_g)$.

Fitting the Null Rate

Ignoring terms in the negative-binomial log likelihood that do not depend on the mean, the contribution from one cell is $Y_{ig}\log(\mu_{ig})-(Y_{ig}+r_g)\log(\mu_{ig}+r_g)$. Substituting $\mu_{ig}=s_iq_g$ and differentiating with respect to $q_g$ gives the null-fit score. For a positive fitted rate, the score equation is:

\frac{\partial \ell_g}{\partial q_g} =\frac{r_g}{q_g}\sum_i \frac{Y_{ig}-s_iq_g}{s_iq_g+r_g}=0, \qquad\text{so}\qquad \sum_i\frac{Y_{ig}-s_iq_g}{s_iq_g+r_g}=0.

This equation generally requires numerical solution when exposures differ. Its left-hand side is strictly decreasing in $q_g$: the derivative of its $i$th term is $-s_i(r_g+Y_{ig})/(s_iq_g+r_g)^2<0$. If at least one count is positive, the left-hand side is positive at $q_g=0$ and approaches a negative value as $q_g$ increases without bound. It therefore has a unique positive root, which gives the maximum-likelihood estimate. If all observed counts are zero, the score is negative for every $q_g>0$, and the maximum occurs at the boundary $\hat q_g=0$.

This boundary case matters in sparse data: the positive-root rule does not apply to an all-zero gene. Its null and saturated means are both zero, so its deviance is zero. Such a gene normally contributes no evidence for cell-to-cell variation and should not be treated as though a positive rate had been estimated.

Worked Example: The Null Fit With Equal Exposures

Suppose a hypothetical gene has counts $0,2,4,2$ in four cells, all with exposure $s_i=1$, and dispersion $\alpha_g=0.5$, so $r_g=2$. For equal exposures, the score equation reduces to $\sum_i(Y_{ig}-q_g)/(q_g+r_g)=0$. Its denominator is common across cells, so the fitted rate is the sample mean:

\hat q_g=\frac{0+2+4+2}{4}=2, \qquad \hat{\mu}_{ig}=1(2)=2\quad\text{for each cell}.

There are positive counts, so the fit is an interior solution. In contrast, if the four counts were all zero, the maximum-likelihood estimate would be $\hat q_g=0$, not a positive solution to the score equation. This example also illustrates why the null fit is a shared rate rather than a shared observed count: the fitted mean is multiplied by each cell’s exposure.

How the Deviance Formula Arises

Deviance compares the likelihood at the saturated means with the likelihood at the null means. For one cell, set its saturated mean to $Y_{ig}$ and its null mean to $\hat{\mu}_{ig}$. The terms in the log likelihood that depend on the mean yield the difference

\ell(Y_{ig};Y_{ig})-\ell(Y_{ig};\hat{\mu}_{ig}) = Y_{ig}\log\left(\frac{Y_{ig}}{\hat{\mu}_{ig}}\right) -(Y_{ig}+r_g)\log\left(\frac{Y_{ig}+r_g}{\hat{\mu}_{ig}+r_g}\right).

The other likelihood terms are identical for the two fits and cancel. Summing the differences over cells and multiplying by two gives $D_g$. Each contribution is nonnegative because the saturated mean maximizes that cell’s likelihood over its mean parameter. A contribution is zero when the observed count equals its null fitted mean; it grows as the count departs from that mean, with the size of the departure interpreted under the chosen dispersion.

This is a likelihood-based measure of fit improvement, not a direct measurement of variance. In particular, the dispersion and the null mean affect the scale of the contribution. A high score says that cell-specific means fit substantially better than a common exposure-adjusted rate under the specified negative-binomial model.

Worked Example: Computing a Gene’s Deviance

Continue with the hypothetical counts $0,2,4,2$, null mean $\hat{\mu}_{ig}=2$, and size $r_g=2$. The two cells with count 2 contribute zero. For a zero count, the contribution is

d(0,2)=2\left[0-2\log\left(\frac{2}{4}\right)\right] =4\log(2)=2.772588722\ldots

For the count of 4, the contribution is

d(4,2)=2\left[ 4\log\left(\frac{4}{2}\right) -6\log\left(\frac{6}{4}\right) \right] =2\left[4\log(2)-6\log(1.5)\right] =0.679596147\ldots

Therefore the total deviance is

D_g=2\left[2\{2\log(2)\}+ \{4\log(2)-6\log(1.5)\}\right] =6.224773592\ldots\approx 6.225.

The two zero-count cells and the cell with count 4 drive the score; the two counts equal to the fitted mean contribute nothing. The score describes departure from a fitted mean of 2 counts per unit exposure under dispersion $0.5$. It is not a probability, nor does this calculation alone establish that the gene is significant.

Worked Example: Similar Mean, Different Deviance

Compare the preceding gene with another hypothetical gene having counts $1,1,3,3$, the same unit exposures, and the same dispersion. Its null fit is also $\hat q_g=2$, since its counts sum to 8. For a count of 1 and a count of 3, the unmultiplied contributions inside the deviance formula are respectively

1\log(1/2)-3\log(3/4)=0.1698990368\ldots,

3\log(3/2)-5\log(5/4)=0.1006775677\ldots.

There are two of each count, so

D_g=2\left[2(0.1698990368\ldots)+2(0.1006775677\ldots)\right] =1.082306418\ldots\approx1.082.

The two genes have the same null mean, but the first has the larger deviance because its counts are more discrepant from that mean under the specified model. This comparison illustrates how the ranking captures distributional departure beyond a difference in average expression.

Worked Example: Unequal Library Exposures

Consider counts $1$ and $4$ in two hypothetical cells with exposures $s_1=1$ and $s_2=2$, and take $r_g=2$. The null means are $q_g$ and $2q_g$, so the score equation is

\frac{1-q_g}{q_g+2}+ \frac{4-2q_g}{2q_g+2}=0.

Multiplying by the positive denominators gives $(1-q_g)(q_g+1)+(4-2q_g)(q_g+2)=0$, or $9-3q_g^2=0$. The positive root and fitted means are

\hat q_g=\sqrt{3}\approx1.732, \qquad (\hat\mu_{1g},\hat\mu_{2g})=(\sqrt{3},2\sqrt{3}) \approx(1.732,3.464).

The simple pooled count-to-exposure ratio is $5/3\approx1.667$, which is not the negative-binomial maximum-likelihood estimate here. The dispersion-weighted score equation, rather than that ratio, defines the null fit when exposures differ. Ranking should use the fitted null means from the chosen model.

Using Rankings and Checking Model Fit

In practice, fit a null model for each retained gene, calculate its deviance, and rank genes from larger to smaller values. The ranking is conditional on the cell set, exposure factors, and dispersion values. It is useful when a gene’s counts depart from the expected negative-binomial pattern after accounting for library size, but it does not tell whether the departure is caused by a continuous gradient, a rare cell group, a technical artifact, or an omitted biological grouping.

Several diagnostics help reveal when the ranking is being driven by model failure rather than informative variation:

Check the offsets. Plot or summarize counts against exposure. Remaining library-size structure can inflate deviance because the null means are then systematically wrong.
Inspect the count pattern. A few extreme cells can dominate a gene’s score. Check whether these cells are plausible observations and whether a small subpopulation accounts for the departure.
Review dispersion choices. A dispersion that is too small can make ordinary count fluctuations appear surprising and inflate scores; an overly large dispersion can suppress real departures.
Look for groups in the analysis set. Different cell types, conditions, or batches can create high deviance even if within-group expression is stable. The score does not distinguish these explanations.
Check sparse and boundary cases. An all-zero gene has deviance zero. For genes with very few positive counts, one or two observations can have disproportionate influence.

Deviance ranking differs from regularized dispersion ranking in what it summarizes. Dispersion methods compare observed variability with a mean-variance trend and may regularize estimates across genes. Deviance instead measures the likelihood gain from allowing cell-specific means over a shared rate, given the selected negative-binomial dispersion and offsets. Neither score is automatically a formal hypothesis test: a ranked list does not supply calibrated p-values or control a false discovery rate.

A further limitation is that a high score is not necessarily a reason to select a gene. It may reflect a technical artifact or a difference between cell populations that is not relevant to the intended analysis. Conversely, a rare but important population may contribute too few cells to create a high overall score. Deviance is a screening criterion that should be interpreted alongside biological design and data quality, not as a substitute for them.

Key takeaway: Negative-binomial deviance is twice the log-likelihood gain of a cell-specific saturated fit over a library-adjusted common-rate null fit. A large value indicates greater departure from that null model under the chosen dispersion; it is not, by itself, a significance result or a biological interpretation.

Check Your Understanding

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

What does the null model assume about a gene’s mean across cells, and how do the exposure factors enter?
For a gene with all-zero counts, what is the null maximum-likelihood rate and why?
With $r=2$, find the deviance contribution of an observed count of zero when its null fitted mean is 2.
Two genes have the same fitted average but different deviances. What does the higher score indicate, and what does it not establish?
Name two diagnostics that could reveal why a gene has a high deviance.

Solutions. (1) The null model uses one rate $q_g$ for the gene, with cell mean $s_iq_g$, so expected counts scale with exposure. (2) The estimate is $\hat q_g=0$: for all-zero counts the likelihood is maximized at the boundary, and the score is negative for every positive rate. (3) The contribution is $2[0-2\log(2/4)]=4\log(2)\approx2.773$. (4) A higher score indicates a larger likelihood improvement from cell-specific means under the specified model; it does not establish a p-value, causation, or biological importance. (5) Examples include checking for residual library-size structure, inspecting extreme cells, reviewing the dispersion, and checking whether distinct cell groups are mixed in the analysis.