Tutorials › Statistical Methods for Single-Cell RNA Sequencing › Simulation Checks For Dropout Models

Noise Models and Dropout · Tutorial 73 of 700

Simulation Checks For Dropout Models

Learn to build model-based simulation checks for count noise, quantify false alarms and detection, and recognize when a sample is too small to distinguish dropout mechanisms.

Intermediate 12 min read

What You'll Learn

  • Define Poisson, negative-binomial, and zero-inflated simulation mechanisms for gene counts.
  • Use simulated zero counts to check calibration under a sampling-only model.
  • Compare diagnostic false-positive rates and detection probabilities across alternatives.
  • Account for sample size, unequal zero probabilities, and parameter fitting in simulation checks.
  • Recognize why a diagnostic result does not by itself identify a dropout mechanism.

Why Check a Diagnostic by Simulation?

A count diagnostic can look convincing on paper and still behave poorly at the sample size available. A rule that flags an unusually large number of zeros may rarely flag data generated by its reference model, yet also have little ability to distinguish that model from a plausible alternative. Simulation makes both questions concrete: how often does the diagnostic flag data when the model is adequate, and how often does it detect a specified departure?

Here the focus is on checking diagnostics by generating hypothetical count datasets from specified mechanisms. This complements posterior predictive checks: the key practical requirement is to reproduce the entire analysis step being assessed, including parameter estimation when parameters would be estimated from the data. A simulation check evaluates a diagnostic under stated assumptions; it does not establish that the simulated mechanism is the true biological or technical explanation.

Definition: A simulation check generates repeated datasets from a specified count model, calculates a diagnostic statistic for each dataset, and compares the observed statistic with that simulated reference distribution. Its calibration concerns behavior under the specified model; its ability to detect a departure concerns performance under a specified alternative.

Specify the Count Mechanisms

For clarity, begin with one gene measured in \(n\) cells with equal exposure and independent counts. These are simplifying assumptions, not universal properties of a single-cell dataset. A sampling-only Poisson model has \(Y_i\sim\operatorname{Poisson}(\lambda)\), so each cell has probability \(e^{-\lambda}\) of a zero. An overdispersed negative-binomial model can be written with mean \(\mu\) and size \(r\), for which \(\operatorname{Var}(Y_i)=\mu+\mu^2/r\) and \(P(Y_i=0)=(r/(r+\mu))^r\). A zero-inflated Poisson model has a structural-zero probability \(\pi\); with probability \(1-\pi\), the count is drawn from \(\operatorname{Poisson}(\lambda)\). Its zero probability is \(\pi+(1-\pi)e^{-\lambda}\).

These formulas are useful for designing and checking simulations, but the simulated data should follow the full model, not just its zero probability. For example, a zero-inflated Poisson draw first selects the structural-zero state and otherwise draws a Poisson count. Nonzero counts, the sample mean, and the spread can therefore help distinguish models too. In real data, exposures may differ by cell, gene counts may be dependent, and parameters may be estimated. A credible check must represent those features when they matter to the diagnostic.

Assumptions for a simple zero-count benchmark: If cell counts are independent and have the same zero probability \(p_0\), then the number of zeros \(K\) has a \(\operatorname{Binomial}(n,p_0)\) distribution. If exposures differ, zero probabilities differ by cell, or cells are dependent, this binomial reference need not hold.

The justification is direct. Define \(I_i=1\) when cell \(i\) has a zero and \(I_i=0\) otherwise. Under the stated assumptions, the \(I_i\) are independent Bernoulli variables with success probability \(p_0\), and \(K=\sum_{i=1}^{n} I_i\). Therefore,

$$ P(K=k)=\binom{n}{k}p_0^k(1-p_0)^{n-k},\qquad k=0,\ldots,n. $$

If each cell instead has its own zero probability \(p_{0i}\), then \(K\) is a sum of independent Bernoulli variables with different success probabilities, called a Poisson-binomial count. If counts are dependent, even that distribution may be inappropriate. Simulation can accommodate such structures when they are included in the data-generating process.

Worked Example: Calibrating a Zero-Count Threshold

Suppose \(n=100\) independent cells have counts generated by a Poisson model with known \(\lambda=2\). Consider the diagnostic rule “flag if at least 20 cells have zero counts.” Under this null model, \(p_0=e^{-2}=0.135335\), so \(K\sim\operatorname{Binomial}(100,0.135335)\). The exact binomial tail is

$$ P(K\ge20)=\sum_{k=20}^{100}\binom{100}{k}(0.135335)^k(0.864665)^{100-k} \approx 0.04592. $$

Thus this particular threshold has an approximately 4.592% false-positive rate under the stated Poisson model and sample size. It is not exactly a 5% rule, and the calibration depends on the specified mean, independence, and equal exposure assumptions.

Now compare specified alternatives. Under a negative-binomial model with \(\mu=2\) and \(r=2\), the zero probability is \((2/(2+2))^2=0.25\). Hence \(K\sim\operatorname{Binomial}(100,0.25)\), and the probability of crossing the same threshold is \(P(K\ge20)\approx0.90047\), or about 90.05%. This is the detection probability for this specific alternative, not a universal power value.

For a zero-inflated Poisson alternative, set \(\pi=0.20\) and choose \(\lambda\) so the overall zero probability is \(0.3021\): \(0.20+0.80e^{-\lambda}=0.3021\), giving \(\lambda\approx2.059\). Then \(K\sim\operatorname{Binomial}(100,0.3021)\), and \(P(K\ge20)\) is approximately 0.99. The three rates show that a zero-count threshold can distinguish these specified alternatives at \(n=100\); they do not show that a high zero count uniquely identifies zero inflation.

Use Simulations to Measure Finite-Sample Performance

A simulation check is most informative when it mirrors the planned analysis. First specify the mechanism, sample size, exposures, and any dependence. Then generate many datasets, apply the diagnostic to each, and summarize how often it flags. Under a null model, the flagging proportion estimates the false-positive rate of the rule. Under an alternative, it estimates that rule’s detection probability. Repeating this across realistic sample sizes reveals whether a diagnostic that works with hundreds of cells remains useful with a small sample.

1
Write down the data-generating model.
Specify count distributions, parameters, cell exposures, and any dependence or cell-level variation relevant to the diagnostic.
2
Generate replicate datasets.
For each replicate, generate the same number of cells and the same design as in the intended analysis.
3
Repeat the diagnostic procedure.
Re-estimate parameters in each replicate if the real procedure estimates them, and calculate the same statistic or flag.
4
Summarize performance and uncertainty.
Report the fraction flagged, compare null and alternative results, and check whether Monte Carlo error is small enough for the conclusion.

For \(B\) simulated datasets, if \(b\) are flagged, the estimated flagging rate is \(b/B\). Its approximate Monte Carlo standard error is \(\sqrt{(b/B)(1-b/B)/B}\). This quantifies simulation noise, not uncertainty from having observed only one biological dataset. Using more replicates reduces Monte Carlo error but cannot repair an unrealistic simulation model.

Worked Example: A Small Sample Can Blur the Distinction

Consider \(n=30\) independent cells and a hypothetical observation of 10 zero counts. The observed zero fraction is \(10/30=0.3333\). Under a Poisson model with known \(\lambda=2\), \(p_0=0.135335\), so the expected number of zeros is \(30(0.135335)=4.060\). The binomial standard error of the zero fraction is

$$ SE(\hat p_0)=\sqrt{\frac{p_0(1-p_0)}{n}} =\sqrt{\frac{0.135335(0.864665)}{30}} \approx 0.06246. $$

The observed fraction is about \((0.3333-0.135335)/0.06246=3.17\) standard errors above the model expectation. This is a substantial discrepancy for this simple benchmark. But under the negative-binomial alternative with \(\mu=2\) and \(r=2\), \(p_0=0.25\), and \(SE(\hat p_0)=\sqrt{0.25(0.75)/30}=0.07906\). The observed fraction is only \((0.3333-0.25)/0.07906=1.05\) standard errors above that alternative’s expectation.

The calculation indicates why simulation should compare the diagnostic under both mechanisms at the actual sample size. It is not a complete model assessment: the standard-error comparison treats the probabilities as known and focuses only on zeros. If the mean is estimated, exposures vary, or cells are dependent, simulate those steps and features rather than treating this calculation as a calibrated test.

Check More Than One Feature, Without Overreading It

A zero fraction alone discards information in the nonzero counts. Useful prespecified summaries can include the sample mean, variance, upper-tail counts, and zero fraction. A simulation-based check can compare the observed summaries with the distributions generated under each candidate model. If the observed data are unusual on several distinct summaries, the model may be inadequate in more than one way; if only the zero count is unusual, the discrepancy is specifically about zeros relative to that model.

Small samples can make even visibly different mechanisms hard to separate. The next example uses a hypothetical count vector and calculates a simple zero-count comparison. It also illustrates why simulation distributions, rather than a single summary or a visual impression, are useful.

Worked Example: Interpreting an Eight-Cell Dataset

Suppose eight hypothetical cell counts are \(0,1,2,0,4,1,3,0\). There are \(K=3\) zeros, for a zero fraction of \(3/8=0.375\). The sum is 11, so the sample mean is \(11/8=1.375\). The sum of squared counts is 31, giving sample variance

$$ s^2=\frac{\sum_i Y_i^2-(\sum_i Y_i)^2/n}{n-1} =\frac{31-121/8}{7} =\frac{15.875}{7} \approx 2.268. $$

As a simple zero-count benchmark, suppose the Poisson mean is known to be 2. Then \(K\sim\operatorname{Binomial}(8,e^{-2})\), and \(P(K\ge3)\) is about 0.08199. Under a negative-binomial model with \(\mu=2,r=2\), \(K\sim\operatorname{Binomial}(8,0.25)\), and \(P(K\ge3)\) is about 0.321. Three observed zeros are therefore not especially unusual under either reference at this small \(n\). The same data’s mean and variance can be compared with simulated distributions as additional, prespecified checks; the observed summaries alone do not establish which mechanism generated the counts.

A practical simulation would repeatedly draw eight counts under each model, calculate all chosen summaries, and see how often simulated datasets are at least as discrepant as the observed one. With only eight cells, such comparisons will often be uncertain. Increasing the number of simulation replicates improves the estimated reference distribution, but it does not increase the information in the eight observed cells.

Calibration, Diagnostics, and Alternatives

When a simulation check disagrees with observed data, inspect the assumptions before assigning a mechanism. A diagnostic can appear miscalibrated because the simulation used a fixed mean when the analysis estimates it, assumed equal cell exposures when they vary, or ignored shared sources of dependence. Conversely, matching a zero fraction does not guarantee that the distribution of positive counts or the upper tail is adequate. When exposures differ, for example, even a Poisson model can yield different zero probabilities across cells; the resulting zero count is not generally binomial with one common \(p_0\).

Parameter fitting matters especially for model-based reference distributions. If the real analysis estimates \(\lambda\) from the observed counts, generating every replicate at a fixed value and comparing the result with that reference does not reproduce the full procedure. A parametric bootstrap-style check instead fits the model, simulates a dataset from the fitted model, refits it, and recalculates the diagnostic. Whether parameters should be refit depends on the diagnostic and the intended question; the important point is to apply the same estimation steps in simulations as in the analysis being evaluated.

ApproachWhat it can assessMain caution
Analytic binomial benchmarkZero-count behavior under independent cells with a common known zero probabilityMay not apply with estimated parameters, unequal probabilities, or dependence
Model-based simulation checkFinite-sample behavior under a specified data-generating and analysis procedureConclusions depend on whether the simulated mechanism and fitting steps are realistic
Residual or summary plotsPatterns in means, spread, zeros, or tails that a single threshold may missVisual agreement is not a calibrated error rate and can be subjective

One limitation is that a simulation check answers a conditional question: what would this diagnostic do if the specified model generated the data? A model that passes a finite set of checks is not proven correct, and different mechanisms can produce similar observed summaries at realistic sample sizes. Excess zeros, for instance, do not by themselves prove structural zero inflation. Biological gene-off states, capture variation, and sampling noise may all contribute to non-detection. The simulation should therefore clarify what patterns a model can reproduce, not turn one diagnostic into a causal claim about why zeros occurred.

Common pitfall: “The observed zero fraction is high, so the gene must have structural zeros.” This treats one summary as a unique mechanism identifier. Instead, compare the observed statistic with simulations under well-specified alternatives, include other relevant count features, and report how reliably the diagnostic distinguishes those alternatives at the available sample size.

Check Your Understanding

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

  1. For independent Poisson counts with known \(\lambda=1.5\) in \(n=50\) cells, what is the zero probability, and what distribution does the number of zeros follow?
  2. Why might a binomial reference for the number of zeros be inappropriate when cell exposures differ?
  3. A simulation under a null model flags 48 of 1,000 datasets. Calculate the estimated flagging rate and its approximate Monte Carlo standard error.
  4. Why should parameter estimation be repeated in simulated datasets when it is part of the real analysis?
  5. What conclusion is justified if a zero-count diagnostic flags frequently under a Poisson simulation but rarely under a negative-binomial simulation at the observed sample size?

Solutions. (1) The zero probability is \(e^{-1.5}\approx0.2231\). With independent cells and the same probability in each, \(K\sim\operatorname{Binomial}(50,0.2231)\). (2) Different exposures generally imply different Poisson means and thus different zero probabilities, so the zero indicators are not identically distributed; their sum is not the stated common-probability binomial count. (3) The estimated rate is \(48/1000=0.048\), or 4.8%. Its approximate Monte Carlo standard error is \(\sqrt{0.048(0.952)/1000}\approx0.00676\), or 0.676 percentage points. (4) Refitting reproduces the actual analysis procedure and its finite-sample effects; omitting fitting can give a reference distribution that is not calibrated for the procedure used on the observed data. (5) The diagnostic distinguishes those two specified models reasonably often under the simulated conditions. It does not prove which mechanism generated the observed data, nor does it establish that other mechanisms are absent.

Clinical Biostats