Checking Fitted Means and Excess Zeros

Why Check the Zeros a Model Predicts?

A pseudobulk count can be zero because a gene has low expected expression and a count of zero is plausible under ordinary sampling variation. A zero is not, by itself, evidence of a special “dropout” process. The diagnostic question is more specific: given the fitted mean, dispersion, design, and dependence structure, does the model predict about as many zeros as are observed?

For a gene, a negative binomial model specifies a distribution for each pseudobulk count, conditional on modeled covariates and any random effects. Its fitted mean determines the probability of a zero, while its dispersion determines how much probability lies near zero and in the upper tail. Thus an excess of observed zeros can reflect an inadequate mean model, an inadequate dispersion model, dependence that was not represented, or a process beyond the fitted negative binomial model. The zero comparison helps identify lack of fit; it does not by itself identify the cause.

Definition: For pseudobulk sample \(i\), let \(Y_i\) be the gene count and let \(\hat{\mu}_i\) be its fitted mean. Under the NB2 parameterization used here, \(Y_i\) has conditional variance \(\mu_i+\mu_i^2/\theta\), where \(\theta>0\) is the dispersion parameter. The fitted probability of a zero is \[ \hat{p}_{0i}=\left(\frac{\hat{\theta}}{\hat{\theta}+\hat{\mu}_i}\right)^{\hat{\theta}}. \] The observed zero count is \(Z_{\mathrm{obs}}=\sum_{i=1}^n I(Y_i=0)\), and its fitted expectation, conditional on the fitted probabilities, is \(\sum_{i=1}^n\hat{p}_{0i}\).

The mean model commonly uses a log link, so fitted means may vary with condition, batch, sample-level covariates, and a library-size offset. A model with a random effect also has conditional means that depend on the random effect. The zero probability must be calculated from the mean appropriate to the question: conditional on a specified random effect for a conditional check, or averaged over the random-effect distribution for a population-level check.

Why the Negative Binomial Predicts Zeros

With mean \(\mu\) and dispersion \(\theta\), the NB2 probability mass function for a nonnegative integer \(y\) is

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

Setting \(y=0\) makes the factorial and gamma ratio equal to one, and the final factor is raised to zero. Therefore

P(Y=0)=\left(\frac{\theta}{\theta+\mu}\right)^\theta.

This expression makes the role of fitted means explicit: at a fixed dispersion, increasing \(\mu\) decreases the zero probability. Dispersion matters too. For example, as \(\theta\) grows without bound, the NB2 model approaches a Poisson distribution and its zero probability approaches \(e^{-\mu}\). Consequently, it is generally wrong to predict zeros by applying the zero-probability formula to the average fitted mean when fitted means differ across samples. Calculate each sample’s probability first, then add them.

Worked Example: A Single Fitted Mean

Suppose a hypothetical gene has \(n=12\) independent pseudobulk counts. For this illustration, the fitted model gives the same mean, \(\hat{\mu}_i=2\) counts, for every sample, with fitted dispersion \(\hat{\theta}=1\). Five of the twelve observed counts are zero. What zero count does the model predict?

For each sample, \[ \hat{p}_{0i}=\left(\frac{1}{1+2}\right)^1=\frac{1}{3}=0.3333. \] The expected number of zeros is therefore \[ E(Z_{\mathrm{obs}})=\sum_{i=1}^{12}\hat{p}_{0i} =12\left(\frac{1}{3}\right)=4. \] The observed fraction is \(5/12=0.4167\), compared with a fitted expected fraction of \(4/12=0.3333\). This difference is a diagnostic signal to examine, not proof that the model is wrong: counts fluctuate, and the model’s parameters were estimated rather than known.

Allowing Fitted Means to Differ

In a real design, a gene may have distinct fitted means across conditions, covariate values, or library sizes. Under independent counts and fixed fitted probabilities, the zero indicators \(I(Y_i=0)\) are Bernoulli variables with probabilities \(\hat{p}_{0i}\). Their sum has a Poisson-binomial distribution, with

E(Z_{\mathrm{obs}})=\sum_{i=1}^n\hat{p}_{0i}, \qquad \operatorname{Var}(Z_{\mathrm{obs}})=\sum_{i=1}^n\hat{p}_{0i}(1-\hat{p}_{0i}).

The variance formula follows because independent Bernoulli variables have variances \(p(1-p)\), and the variance of an independent sum is the sum of the variances. These formulas give a useful reference when probabilities are treated as fixed. They do not fully account for parameter estimation, and they are not valid as written when zero indicators are dependent.

Worked Example: Unequal Means and an Approximate Check

Consider twelve hypothetical independent pseudobulk samples, four each with fitted means \(1\), \(2\), and \(4\) counts, and a common fitted dispersion \(\hat{\theta}=2\). Suppose six observed counts are zero. Calculate the expected and standard deviation of the zero count under the fixed fitted probabilities.

For the four samples with fitted mean 1, \[ \hat{p}_0=\left(\frac{2}{2+1}\right)^2=\frac{4}{9}=0.4444. \] For mean 2, \(\hat{p}_0=(2/4)^2=0.2500\); for mean 4, \(\hat{p}_0=(2/6)^2=1/9=0.1111\). Thus, \[ E(Z_{\mathrm{obs}})=4\left(\frac{4}{9}+\frac{1}{4}+\frac{1}{9}\right) =3.222. \] The fixed-probability variance is \[ 4\left[\frac{4}{9}\left(1-\frac{4}{9}\right) +\frac{1}{4}\left(1-\frac{1}{4}\right) +\frac{1}{9}\left(1-\frac{1}{9}\right)\right] =2.133, \] so the standard deviation is \(\sqrt{2.133}=1.460\). The observed count is \(6\), or \( (6-3.222)/1.460=1.903\) standard deviations above the fixed-probability expectation. This is a descriptive approximation, not a calibrated formal test: fitted parameters, dependence, and model selection can change the reference distribution.

Using the overall mean of \( (4(1)+4(2)+4(4))/12=2.333\) instead would erase the sample-to-sample differences in fitted means and yield a different zero prediction. The correct calculation sums the twelve probabilities associated with their own fitted means.

Mixed Models and Design-Aware Calibration

For repeated pseudobulk measures, zero indicators from the same donor may be dependent. A mixed model represents such dependence through random effects, but a conditional zero probability is not automatically the population-average probability. If the model has random effect \(b_i\) and conditional mean \(\mu_i(b_i)\), the conditional probability is \(P(Y_i=0\mid b_i)\). A marginal probability averages this quantity over the random-effect distribution:

P(Y_i=0)=\int \left(\frac{\theta}{\theta+\mu_i(b)}\right)^\theta \,dF(b),

where \(F\) is the specified distribution of the random effect. In general, this integral is not equal to the zero probability evaluated at the marginal mean. For a check of the observed study, simulation from the fitted mixed model is often more direct than relying on a simple variance formula.

Define the check.
Specify the gene, set of pseudobulk samples, and whether the question concerns conditional or marginal zero frequencies. Use the same retained samples and offsets as the fitted analysis.

Compute the observed and fitted quantities.
Count observed zeros and calculate model-predicted zero probabilities using each sample’s fitted mean and the fitted dispersion.

Examine where discrepancies occur.
Compare by fitted-mean range and relevant design groups. A discrepancy concentrated in one condition or at a particular range of library sizes may indicate a mean-model problem rather than a uniform excess-zero process.

Calibrate the comparison if needed.
Simulate counts under the fitted model while preserving the design, offsets, and random-effects structure. For a formal reference that reflects parameter estimation, refit the model to each simulated data set and recompute the diagnostic.

Worked Example: An Omitted Covariate Can Mimic Excess Zeros

Suppose ten hypothetical independent samples are split evenly between two known groups. A reduced model omits group and fits a common mean of \(1\) count with \(\hat{\theta}=1\). There are seven observed zeros. The reduced model predicts \[ 10\left(\frac{1}{1+1}\right)=5 \] zeros. Now consider a model that includes group and gives fitted means \(0.25\) for five samples and \(1.75\) for the other five, still with \(\hat{\theta}=1\). Its predicted zero count is \[ 5\left(\frac{1}{1+0.25}\right)+5\left(\frac{1}{1+1.75}\right) =5(0.8000)+5(0.3636)=5.818. \] The two models have the same average fitted mean of \(1\), but their predicted zero counts differ. The group-adjusted model predicts more zeros because the low-mean samples have a high probability of being zero. Seven observed zeros may therefore point toward a missing mean structure; the zero count alone cannot distinguish this from other kinds of lack of fit.

The calculation does not establish that the group-adjusted model is correct. It shows why zero diagnostics should be read alongside fitted means and design information, rather than interpreted as a direct test for a separate zero-generating mechanism.

Interpreting Discrepancies and Choosing Alternatives

A practical diagnostic compares observed zeros with the distribution of zero counts generated under the fitted model, not just with a single expected value. In a parametric simulation, generate counts under the fitted negative binomial model using the original design and offsets; for a mixed model, generate random effects according to its fitted structure. The simulated zero counts form a reference distribution. If the model was fitted to the same data, refitting each simulated data set helps reflect the fact that the parameters were estimated. A discrepancy is evidence of model miscalibration for this diagnostic, not a probability that the model is true or false.

Patterns matter. Too many zeros mainly among low fitted means may be ordinary sampling variation; too many among high fitted means are more concerning. A group-specific discrepancy may suggest an omitted interaction, covariate, or incorrect mean trend. Broad excess variation can implicate the dispersion model, while clustered discrepancies within donors can signal inadequately represented dependence. Check that zero counts were not created by filtering or preprocessing choices inconsistent with the fitted sample set.

The main alternatives are not interchangeable fixes. A Poisson model imposes variance equal to the mean and can predict too few zeros when counts are overdispersed. A negative binomial model allows variance above the mean and already accommodates many zeros. A zero-inflated model adds a separate probability of a structural zero to a count distribution; a hurdle model separates the zero-versus-positive process from the positive-count distribution. Either extension needs a defensible interpretation, adequate data, and evidence of improved fit or predictive calibration. It should not be selected merely because the observed zero fraction exceeds a rough plug-in expectation.

Common pitfall: “The observed zero percentage exceeds the percentage calculated from the average fitted mean, so the gene needs a zero-inflated model.” This reasoning ignores sample-specific fitted means, parameter uncertainty, and dependence. Compute probabilities from each fitted mean, examine the pattern of discrepancies, and use design-aware simulation before considering a more elaborate zero model.

Check Your Understanding

Answer each question, then compare with the solutions below.

For a negative binomial count with \(\mu=3\) and \(\theta=2\), calculate the probability of zero.
Four independent samples have fitted zero probabilities \(0.1, 0.2, 0.5,\) and \(0.8\). What are the expected zero count and its fixed-probability variance?
Why can a mixed model’s marginal probability of zero differ from the probability evaluated at its marginal mean?
What is one reason that a large observed-minus-expected zero count does not, by itself, establish a need for a zero-inflated model?

Solutions. (1) \(P(Y=0)=(2/(2+3))^2=(2/5)^2=0.1600\). (2) The expectation is \(0.1+0.2+0.5+0.8=1.6\). The variance is \(0.1(0.9)+0.2(0.8)+0.5(0.5)+0.8(0.2)=0.09+0.16+0.25+0.16=0.6600\). (3) The zero probability is a nonlinear function of the conditional mean, so averaging that probability over random effects generally differs from evaluating it at an averaged mean. (4) The discrepancy could arise from estimated parameters, unmodeled covariates, dependence, or an inadequate dispersion or mean model; the zero count does not identify a unique cause.