Simulation Diagnostics For Post-Clustering P-Values

Why Diagnose the Whole Pipeline?

A marker p-value can be distorted when the same expression data help define clusters and then test genes for differences between those clusters. A simulation diagnostic asks a direct question: if the marker null were true, how would this complete analysis behave when repeated? The answer is informative only if each simulated data set passes through the same clustering, marker-selection, testing, and reporting steps as the analysis being evaluated.

This is not just a check of the test formula in isolation. Clustering may choose partitions that accentuate random expression differences; filtering may determine which genes are tested; and a rule for selecting a cluster contrast may favor extreme statistics. If a diagnostic simulation skips any of these steps, it does not measure the calibration of the procedure actually used. Earlier material on permutation tests that rerun the clustering pipeline motivates this end-to-end principle; here the focus is how to summarize and interpret repeated-simulation results.

Specify the Null and Repeat Every Data-Dependent Step

Begin by stating what is null. For example, a simulation may set a particular gene's population expression contrast to zero while allowing other genes to define structure. Alternatively, a global-null simulation may generate no true cluster-associated differences. These are different scenarios and answer different questions. In either case, define the target parameter and the null hypothesis before simulating; a zero contrast for a fixed biological grouping is not automatically the same target as a contrast between clusters selected from the same data.

A useful design preserves features of the intended analysis that affect calibration: cell number, library-size variation, gene-gene dependence, batch structure, sparsity, normalization, filtering, clustering choices, and the marker test. The simulated null must be credible for the question. If the marker gene is allowed to influence clustering in the real pipeline, it should generally be allowed to do so in the diagnostic as well. If the simulation removes that pathway, it may miss the very selection effect being assessed.

For each independent simulation replicate, regenerate the data and rerun all data-dependent steps, including choosing the number of clusters if that choice is part of the procedure. Record the resulting p-value for a pre-specified null test, or record a clearly defined summary of the family of tests. Do not treat genes within one simulated data set as independent simulation replicates: their p-values can be strongly dependent.

What Calibration Looks Like

For a valid p-value under a null hypothesis, the null rejection probability at threshold \(t\) should not exceed \(t\): \(P(p\leq t)\leq t\) for \(0\leq t\leq1\). This is called super-uniformity. An exactly uniform null p-value satisfies equality at every threshold; a conservative test can have a smaller probability. A histogram that is approximately flat is therefore a useful visual check for exact calibration, but a non-flat histogram is not by itself a complete diagnosis: the relevant inferential claim, test discreteness, simulation size, and p-value dependence all matter.

The rejection rate is a focused summary at a chosen significance level \(\alpha\). In replicate \(j\), define \(R_j=1\) if the null p-value is at most \(\alpha\), and \(R_j=0\) otherwise. With \(M\) independent replicates, the estimated rejection rate is \(\widehat{q}=M^{-1}\sum_{j=1}^{M}R_j\). Its expectation and variance require care. If replicate \(j\) has rejection probability \(q_j=P(R_j=1)\), then

where the variance expression assumes the \(R_j\) are independent. Only when the replicate outcomes are identically distributed Bernoulli\((q)\) does this simplify to \(E(\widehat{q})=q\) and \(\operatorname{Var}(\widehat{q})=q(1-q)/M\). This distinction matters if the simulation deliberately mixes different scenarios or uses differing procedures. In that case, report scenario-specific rates or use the general formulas, rather than silently applying a binomial standard error.

Under the identical-replicate assumption, a plug-in Monte Carlo standard error is \(\sqrt{\widehat{q}(1-\widehat{q})/M}\). It describes simulation noise in the estimated rejection rate, not uncertainty about the biological effect or the correctness of the null model. A small standard error cannot rescue a poorly specified simulation.

Use Histograms and Thresholds Together

A p-value histogram can reveal where departures occur. Excess mass near zero suggests anti-conservative behavior for some thresholds; a deficit near zero can indicate conservative behavior. The empirical cumulative proportion \(\widehat{F}(t)=M^{-1}\sum_{j=1}^{M} I(p_j\leq t)\) provides a threshold-specific view: compare it with \(t\), or with the appropriate conservative bound. Inspect thresholds that matter scientifically, such as the nominal test level, as well as the overall shape. A histogram can conceal inflation concentrated below a small threshold if its bins are too broad.

The histogram and rejection-rate calculations should be based on comparable, independently generated replicates—for example, one pre-specified null p-value per replicate. If each replicate contributes hundreds of gene-level p-values, those observations are generally dependent. Counting them all as independent can make Monte Carlo uncertainty look much smaller than it is. Options include focusing on a pre-specified gene, summarizing each replicate first, or estimating uncertainty by resampling whole independent replicates.

Coverage Diagnostics Need a Defined Target

Coverage is the proportion of simulated intervals that contain the true parameter. It is meaningful only when the parameter is clearly defined and fixed or known in each simulated replicate. For a post-clustering interval, specify whether the target is a contrast between fixed latent groups, a parameter conditional on a specified selection event, or another estimand. If the procedure reports intervals only for selected markers, distinguish coverage among reported intervals from the probability that a marker is selected. Those quantities answer different questions.

For \(K\) independent, identically distributed interval outcomes with coverage probability \(c\), the empirical coverage \(\widehat{c}\) has plug-in Monte Carlo standard error \(\sqrt{\widehat{c}(1-\widehat{c})/K}\). If selection determines whether an interval is reported, a conditional coverage estimate uses the number of covered intervals divided by the number reported. Its uncertainty calculation must reflect that conditional sampling and the simulation design; it should not be confused with an unconditional rate over all generated data sets. Most importantly, nominal 95% coverage is a repeated-sampling property under the stated model and target, not a 95% probability that a particular realized interval contains the parameter.

Derivation: Why the Rejection-Rate Standard Error Depends on the Design

Write \(R_j\) for the rejection indicator in replicate \(j\), so \(E(R_j)=q_j\). By linearity of expectation,

For independent replicates, the variance of a sum is the sum of variances, and a Bernoulli variable with success probability \(q_j\) has variance \(q_j(1-q_j)\). Therefore,

If all \(q_j=q\), then \(\widehat{q}\) is the average of \(M\) identically distributed Bernoulli\((q)\) variables, and the variance reduces to \(q(1-q)/M\). Independence alone is not enough for that reduction. This is why one should state whether replicates share one simulation scenario and one procedure before using a binomial standard error.

Diagnostics, Alternatives, and Limits

A simulation can fail in ways that look like reassuring calibration. Check that the code reruns clustering and every selection decision, that each null is correctly imposed, and that the test statistic and p-value calculation match the intended analysis. Record failures to cluster, empty groups, or untestable genes rather than silently dropping them; exclusion can change the population of runs being summarized. Repeat the diagnostic across plausible nuisance settings, since calibration under one cell count or dependence pattern does not guarantee calibration elsewhere.

Simulation is a diagnostic, not a universal replacement for inference. A method derived for a well-defined selection event may offer a targeted theoretical guarantee under its assumptions; an independent validation data set avoids reusing discovery data for the validation test; and a permutation approach can be appropriate when its exchangeability requirements hold. These alternatives answer different questions and each has assumptions. Simulation complements them by showing how a chosen procedure behaves under specified models, while remaining limited by the realism and scope of those models.

A common mistake is to report “the p-values are calibrated” because their histogram looks roughly uniform or because one estimated rejection rate is near \(\alpha\). The correction is to state the null model, estimand, complete analysis pipeline, simulation count, and uncertainty; check several relevant thresholds and, when applicable, coverage; and limit the conclusion to the simulated conditions. A second mistake is to pool dependent gene-level p-values and use their total number as \(M\). The correction is to use independent simulation replicates as the uncertainty units or to model the dependence explicitly.