Diagnosing Excess Zeros Against A Count Model

Why Diagnose Zeros Against a Count Model?

A large fraction of zero counts is not, by itself, evidence that a count model fails. When a gene has a low expected count, a Poisson or negative binomial model can predict many zeros even without any separate zero-generating process. The useful diagnostic question is therefore: given the fitted count model and the expected expression levels, are there more or fewer zeros than the model predicts?

This tutorial develops summaries for answering that question. It focuses on observed versus predicted zero frequencies, residuals within groups of similar fitted means, and simulation to assess how much variation is plausible. The earlier material on structural zeros and sampling zeros distinguishes possible interpretations of a zero; the summaries here assess fit to a specified count distribution, not the biological cause of any discrepancy.

Model-Based Expected Zeros and Their Uncertainty

Suppose a gene is measured in \(n\) cells, with count \(Y_i\) for cell \(i\). Under a Poisson model, assume that, conditional on the fitted mean \(\mu_i\), \(Y_i\) has a Poisson distribution with mean \(\mu_i\). Under a negative binomial model, use the mean-dispersion parameterization in which the conditional mean is \(\mu_i\) and the variance is \(\mu_i+\phi\mu_i^2\), where \(\phi\geq 0\). For the formulas below, the fitted means and any covariates or exposures are treated as fixed. The counts are assumed conditionally independent across cells. These assumptions matter: shared cell-level effects or other dependence can make the actual variability of the total number of zeros larger than the formulas suggest.

Definition: Let \(Z_i=1\) if \(Y_i=0\), and \(Z_i=0\) otherwise. Write \(q_i=P(Y_i=0)\) under the model. The observed number of zeros is \(O_0=\sum_{i=1}^n Z_i\), and its model-based expected value is \(E_0=\sum_{i=1}^n q_i\). Under a Poisson model, \(q_i=e^{-\mu_i}\). Under a negative binomial model with mean \(\mu_i\) and dispersion \(\phi>0\), \(q_i=(1+\phi\mu_i)^{-1/\phi}\).

If the observations are conditionally independent, each zero indicator is Bernoulli with probability \(q_i\). Therefore its variance is \(q_i(1-q_i)\), and variances add across independent cells:

E(O_0)=E_0=\sum_{i=1}^n q_i, \qquad \operatorname{Var}(O_0)=\sum_{i=1}^n q_i(1-q_i).

A basic standardized residual is the difference between observed and expected zeros divided by the model-based standard deviation:

R_0=\frac{O_0-E_0}{\sqrt{\sum_{i=1}^n q_i(1-q_i)}}.

A positive residual means more zeros were observed than expected; a negative residual means fewer. The denominator describes variation under the stated model with its parameters held fixed. It does not automatically account for uncertainty from estimating the means or dispersion, and a value such as \(R_0=2\) should not be treated as an exact normal-theory test without checking its calibration.

Worked Example: A Poisson Zero-Count Residual

Suppose, hypothetically, that a gene has counts in 100 independent cells and a fitted Poisson mean of \(0.2\) counts per cell. The model predicts a zero probability of \(q=e^{-0.2}=0.818731\). Assume 89 of the 100 observed counts are zero.

The expected number of zeros is \(E_0=100(0.818731)=81.8731\). Since all cells have the same predicted zero probability, the variance is \(100q(1-q)\):

\operatorname{Var}(O_0)=100(0.818731)(0.181269)=14.8411, \qquad R_0=\frac{89-81.8731}{\sqrt{14.8411}} \approx \frac{7.1269}{3.8524}=1.850.

The observed count is about 1.85 model-based standard deviations above its expectation. This is a diagnostic signal to examine, not proof that the Poisson model is false. The calculation assumes independent counts and treats the fitted mean as fixed; parameter estimation or dependence could change the appropriate reference distribution.

Stratifying by Fitted Mean

A single total can hide where a discrepancy occurs. Cells with different fitted means have different zero probabilities, so a useful diagnostic groups cells into ranges of similar fitted means. For a stratum \(b\), let \(I_b\) be its set of cells. Define \(O_b=\sum_{i\in I_b}Z_i\), \(E_b=\sum_{i\in I_b}q_i\), and \(V_b=\sum_{i\in I_b}q_i(1-q_i)\). Its standardized residual is \(R_b=(O_b-E_b)/\sqrt{V_b}\), provided \(V_b>0\). If a stratum is sufficiently large and no single cell dominates its variance, this residual can be a useful approximate scale for the departure; it is not guaranteed to follow a standard normal distribution.

The strata should be chosen to make the pattern interpretable, for example by intervals of fitted mean or predicted zero probability. They should not be selected after inspecting which grouping produces the largest residual. When expected counts or variances are small, a standardized residual can be unstable; show the observed and expected counts alongside it, and rely more on simulation than on a normal approximation.

Worked Example: A Pooled Summary Hides a Mean-Stratum Pattern

Consider a hypothetical set of 100 independent cells divided into two groups of 50. Under the fitted Poisson model, the first group has mean \(\mu=-\log(0.8)=0.223144\), so its predicted zero probability is \(0.8\). The second has mean \(\mu=-\log(0.2)=1.60944\), giving a zero probability of \(0.2\). Suppose 38 cells in the first group and 20 in the second are observed to be zero.

For the first group, \(E_1=50(0.8)=40\) and \(V_1=50(0.8)(0.2)=8\). For the second, \(E_2=50(0.2)=10\) and \(V_2=50(0.2)(0.8)=8\). The stratum residuals are

R_1=\frac{38-40}{\sqrt{8}}=-0.7071, \qquad R_2=\frac{20-10}{\sqrt{8}}=3.536.

Across both groups, there are 58 observed zeros and 50 expected zeros. The total variance is \(8+8=16\), so \(R_0=(58-50)/\sqrt{16}=2.000\). The pooled summary signals an excess, but the stratified results show that it is concentrated in the higher-mean group; the lower-mean group has slightly fewer zeros than predicted. This pattern can guide checks of mean specification or data quality. It does not identify the reason for the discrepancy.

Simulation When Parameters Are Estimated

In practice, the means and dispersion are usually estimated from the same data being checked. Plugging estimates into \(E_0\) and \(V_0\) gives a convenient summary, but it treats those estimates as known. That can understate or otherwise misrepresent uncertainty, especially when dispersion is estimated imprecisely or the model is complex. A parametric simulation provides a way to calibrate a diagnostic while repeating the fitting steps that produced the estimates.

Fit the specified count model.
Record the fitted means, dispersion if applicable, and the observed zero-count statistic or mean-stratum residuals.

Simulate data under the fitted model.
Generate counts using the fitted model while preserving the design features that the model conditions on, such as the number of cells and any fixed exposures.

Refit and recalculate.
For every simulated data set, repeat the parameter-estimation procedure and calculate the same diagnostic statistic. Refit when parameter estimation is part of the original analysis.

Compare the observed statistic with its simulated reference distribution.
For a two-sided discrepancy, count simulations whose absolute statistic is at least as large as the observed absolute statistic. Check whether the observed value is unusual relative to this distribution.

For \(B\) simulations and \(K\) simulated absolute statistics at least as large as the observed one, a common finite-simulation estimate of the two-sided tail probability is \((K+1)/(B+1)\). This avoids reporting a simulated tail probability of exactly zero. It is a Monte Carlo estimate: its precision depends on \(B\). It is not the probability that the model is true, nor does it resolve whether a departure is caused by omitted variation, dependence, a mean-model error, or a distinct zero-generating mechanism.

Worked Example: Simulated Calibration Under a Negative Binomial Model

Suppose, hypothetically, that 40 independent cell counts are fitted with a negative binomial model having common fitted mean \(\mu=0.5\) counts per cell and dispersion \(\phi=0.5\). The predicted probability of zero is \(q=(1+0.5(0.5))^{-1/0.5}=1.25^{-2}=0.64\). Suppose 32 counts are zero.

The plug-in expectation is \(E_0=40(0.64)=25.6\). Under the fixed fitted parameters, \(V_0=40(0.64)(0.36)=9.216\), so

R_0=\frac{32-25.6}{\sqrt{9.216}} =\frac{6.4}{3.0358} \approx 2.109.

Now suppose a hypothetical simulation procedure generates 2,000 data sets from the fitted model, refits the mean and dispersion in each, and recalculates the same standardized zero residual. If 71 simulated absolute residuals are at least 2.109, the finite-simulation tail estimate is \((71+1)/(2000+1)=0.03598\). Its approximate Monte Carlo standard error is \(\sqrt{0.03598(1-0.03598)/2001}=0.0042\). This estimate describes how often a statistic at least this extreme occurred under the specified simulation-and-refitting procedure. It is not a probability that the model is correct. The example assumes the simulated model and refitting procedure adequately represent the analysis; a misspecified simulation model cannot provide a trustworthy calibration.

Diagnostics, Alternatives, and Limitations

Interpret zero residuals alongside the model and data structure, rather than as a stand-alone verdict. In particular, check whether fitted means are sensible across cells, whether the dispersion estimate is stable, and whether residual patterns change across mean strata. Compare the observed and predicted zero counts in each stratum; a single extreme residual in a stratum with very small expected variance can be misleading. Also check whether cells share sources of variation that violate conditional independence. If they do, the sum of Bernoulli variances is not the total variance: covariance terms between zero indicators also contribute.

The main alternatives answer related but different diagnostic questions. A pooled zero fraction is easy to report but ignores how predicted zero probabilities vary across cells. Mean-stratified summaries retain some of that structure and can reveal localized departures, but their results depend on stratum definitions and can be noisy in small groups. Simulation can incorporate parameter estimation and the exact fitting procedure, but it depends on the assumed model and may be computationally demanding. A residual plot or a broader distributional check can reveal nonzero-count discrepancies that a zero-only summary misses. No one summary replaces checking the model’s assumptions and the full pattern of counts.

Do not compare raw zero fractions without accounting for means. A high zero fraction may be expected when fitted means are low.
Do not interpret a plug-in residual as an exact test automatically. Estimated parameters, sparse strata, and dependence can change its calibration.
Do not treat an unusual zero count as a diagnosis of cause. It signals a mismatch with the specified model; other aspects of the model or data may be responsible.
Do not ignore design or data processing. The fitted means, offsets, cell selection, and any filtering must correspond to the data being checked.

Key takeaway: Compare observed zeros with the sum of the model-predicted zero probabilities, examine residuals across fitted-mean strata, and use simulation with refitting when parameter uncertainty matters. Interpret any excess as a model-checking signal, not as proof of a particular biological or technical mechanism.

Check Your Understanding

Try each question first, then use the solutions to check your reasoning.

Under a Poisson model, 20 independent cells each have fitted mean \(0.1\). What is the expected number of zeros?
Why can a pooled zero fraction be misleading when cells have different fitted means?
A stratum has observed zero count 14, expected zero count 10, and model-based variance 4. Calculate its standardized residual and interpret its sign.
In a simulation with 999 replicates, 19 absolute simulated residuals are at least as large as the observed absolute residual. Calculate the finite-simulation two-sided tail estimate using the plus-one formula.

Solutions. (1) The Poisson zero probability is \(e^{-0.1}=0.9048\), so the expected number is \(20(0.9048)=18.10\) zeros. (2) Different fitted means imply different zero probabilities; a pooled fraction can obscure whether departures occur in low- or high-mean cells, or whether opposing stratum patterns cancel. (3) The residual is \((14-10)/\sqrt{4}=4/2=2.000\). It is positive, so there are more observed zeros than expected in that stratum; its magnitude alone does not establish a model failure. (4) The estimate is \((19+1)/(999+1)=20/1000=0.0200\). It is a simulation-based tail estimate under the specified procedure, with Monte Carlo uncertainty, not the probability that the null model is true.