Why a Finite Number of Shuffles Creates Uncertainty
A permutation test compares an observed statistic with the statistics that would arise under a valid permutation scheme. If every allowable arrangement is enumerated, the resulting permutation p-value is determined by that finite reference set. In genomic applications, however, the number of allowable arrangements can be so large that analysts draw only a random subset. The resulting p-value then has Monte Carlo uncertainty: a different set of random shuffles could produce a different count of statistics as extreme as the observed one.
This uncertainty matters most when the estimated p-value is close to a decision threshold. A reported value just below \(\alpha\) may be based on very few extreme shuffles, and a small change in the simulated count can move the estimate above the threshold. The correction and uncertainty calculations in this tutorial concern the accuracy of the permutation calculation, not uncertainty in the observed data or validity of the permutation assumptions themselves.
Definitions and Sampling Assumptions
Let \(T_{\mathrm{obs}}\) be the statistic calculated from the observed data. Before shuffling, specify what counts as at least as extreme as \(T_{\mathrm{obs}}\). For an upper-tail test this might mean \(T^* \geq T_{\mathrm{obs}}\); for a two-sided test it is commonly defined using a prespecified measure such as \(|T^*| \geq |T_{\mathrm{obs}}|\). The same rule, including its treatment of ties, must be used for every shuffle.
The standard binomial calculation assumes that the permutation scheme is valid for the null hypothesis being tested and that each of the \(B\) shuffles is drawn independently from the same permutation distribution. The permutations must respect any design restrictions, such as exchangeability or blocking, established for the test. These calculations do not repair an invalid permutation scheme. They also assume that the extremeness rule was set in advance rather than chosen after looking at the simulated results.
With these assumptions, each simulated shuffle is an indicator event: it either reaches the tail or it does not. Each event has probability \(q\), so
Thus \(\hat q=K/B\) estimates the conditional permutation-tail probability. Its Monte Carlo standard error, when \(q\) is known, is \(\sqrt{q(1-q)/B}\); in practice, an estimated standard error substitutes \(\hat q\) for \(q\). This error describes variation from repeating the random simulation with the observed data held fixed.
Why Use the Plus-One Correction?
Reporting \(K/B\) alone can give a zero p-value when no sampled shuffle reaches the observed statistic. That is not a defensible conclusion that the permutation-tail probability is zero: the finite simulation may simply have missed rare tail outcomes. The plus-one value includes the observed arrangement in the comparison and prevents a zero estimate. It also accounts for the rank of the observed statistic among the observed arrangement and the simulated arrangements.
For a randomization or permutation test whose null makes the observed arrangement exchangeable with the sampled arrangements, the rank-based argument explains the correction. With no ties, the rank of the observed statistic among \(B+1\) exchangeable statistics is uniform. Counting the observed statistic itself as at least as extreme gives a tail probability no smaller than \(1/(B+1)\); ties, handled by counting equal statistics as extreme, make the test conservative rather than anti-conservative. The plus-one construction is therefore a valid finite-simulation p-value under the stated exchangeability conditions. It is not a license to use permutations that violate the design or null model.
The corrected value is not exactly the same quantity as \(\hat q\). In particular,
The correction is nonnegative and at most \(1/(B+1)\). When \(K=0\), \(p_{\mathrm{MC}}=1/(B+1)\), not zero. When the number of shuffles is large, this adjustment is usually small, but the smallest attainable corrected p-value is still determined by \(B\).
Worked Example: A Small Tail Count With 999 Shuffles
Suppose a hypothetical genomic test uses a valid, prespecified two-sided permutation statistic. The observed value is \(T_{\mathrm{obs}}=2.8\). Of \(B=999\) independent valid shuffles, \(K=18\) produce a statistic with absolute value at least \(2.8\). The significance threshold is \(\alpha=0.02\).
The uncorrected tail estimate is \(\hat q=K/B=18/999=0.01802\). The plus-one p-value is \[ p_{\mathrm{MC}}=\frac{18+1}{999+1}=\frac{19}{1000}=0.01900. \] It is below \(0.02\), but only slightly. This comparison alone does not indicate whether Monte Carlo error could plausibly affect the decision.
A Wilson interval gives a useful approximate confidence interval for \(q\). For \(z=1.96\), define \(\hat q=K/B\) and \(d=1+z^2/B\). The Wilson centre and half-width are \[ c=\frac{\hat q+z^2/(2B)}{d}, \qquad h=\frac{z\sqrt{\hat q(1-\hat q)/B+z^2/(4B^2)}}{d}. \] Here \(d=1+3.8416/999=1.003845\). The square-root term is \[ \sqrt{\frac{0.018018(1-0.018018)}{999} +\frac{1.96^2}{4\cdot999^2}} \approx 0.004321. \] The margin numerator is \(1.96(0.004321)\approx0.00847\). Substitution gives \(c\approx0.01986\), \(h\approx0.00844\), and an interval of approximately \((0.0114,\,0.0283)\). Because this interval includes \(0.02\), the simulated tail probability is not estimated precisely enough to make its position relative to that threshold clear. The interval concerns \(q\), the permutation-tail probability, not a probability that the null hypothesis is true.
Quantifying Monte Carlo Uncertainty
For a fixed \(B\), the estimated standard error is \[ \widehat{SE}_{\mathrm{MC}} =\sqrt{\frac{\hat q(1-\hat q)}{B}}. \] It is a quick measure of simulation variability, but a normal interval based on this standard error can perform poorly when \(K\) is small or zero. A Wilson interval is often preferable for a binomial proportion, particularly near zero or one. With \(z=1.96\), the interval is \((c-h,c+h)\) using the Wilson formulas in the worked example. It is an approximate 95% interval for the conditional tail probability \(q\), under the binomial sampling assumptions.
The interval does not turn a Monte Carlo p-value into a confidence interval for a scientific effect. Nor does overlap with \(\alpha\) prove that a test decision is wrong. It is a diagnostic: if the interval straddles the decision threshold, the random simulation has not clearly resolved whether the underlying permutation-tail probability is above or below that threshold. More shuffles can reduce simulation error, though they cannot change the validity of the permutation model.
Worked Example: Five Extreme Shuffles in 199 Draws
In a hypothetical analysis, \(B=199\) independent valid shuffles yield \(K=5\) statistics at least as extreme as the observed statistic. The uncorrected estimate and corrected p-value are \[ \hat q=\frac{5}{199}=0.02513, \qquad p_{\mathrm{MC}}=\frac{5+1}{199+1}=\frac{6}{200}=0.03000. \]
For a 95% Wilson interval, \(d=1+3.8416/199\approx1.019305\). The centre is \[ c=\frac{0.025126+3.8416/(2\cdot199)}{1.019305} \approx 0.03412. \] The square-root term is \[ \sqrt{\frac{0.025126(1-0.025126)}{199} +\frac{3.8416}{4\cdot199^2}} \approx 0.01214, \] so \(h=1.96(0.01214)/1.019305\approx0.02334\). The interval is approximately \((0.0108,\,0.0575)\). At a threshold of \(0.05\), this interval spans both sides of the threshold. The corrected estimate is below \(0.05\), but 199 shuffles provide limited Monte Carlo precision; a larger prespecified simulation would give a more informative estimate.
Decisions Near a Threshold and Simulation Diagnostics
When a result is close to \(\alpha\), report enough information to reveal its simulation resolution: \(K\), \(B\), the corrected p-value, and an uncertainty summary such as the Wilson interval. A p-value based on just a few exceedances should not be presented with spurious decimal precision. If the interval is well separated from \(\alpha\), Monte Carlo error is less likely to alter the threshold comparison; if it overlaps \(\alpha\), consider using more shuffles and reporting that the initial result was simulation-limited.
The scale of the problem can be anticipated. Since \(q(1-q)\leq1/4\), the Monte Carlo standard error is at most \(1/(2\sqrt{B})\). To make this worst-case standard error no larger than \(0.0025\), for example, one needs \(B\geq 1/(4(0.0025)^2)=40{,}000\). Near a specific threshold, the needed number depends on how far the underlying \(q\) is from that threshold; no finite \(B\) can guarantee a clear distinction when \(q\) is arbitrarily close to \(\alpha\).
Worked Example: An Estimate Exactly at the Threshold
Suppose a hypothetical test uses \(B=999\) independent valid shuffles and observes \(K=49\) exceedances. The corrected value is \[ p_{\mathrm{MC}}=\frac{49+1}{999+1}=0.05000. \] The raw proportion is \(\hat q=49/999\approx0.04905\), and its estimated Monte Carlo standard error is \[ \sqrt{\frac{0.04905(1-0.04905)}{999}}\approx0.00683. \] For context, the 95% Wilson interval is approximately \((0.0373,\,0.0642)\). It is much wider than the difference between the estimated tail probability and the threshold \(0.05\). The appropriate conclusion is that this simulation does not resolve which side of \(0.05\) the underlying permutation-tail probability lies on. The point estimate being exactly \(0.05\) does not make the Monte Carlo uncertainty disappear.
A useful process is to set the decision threshold and extremeness rule first, choose a simulation size suited to the required precision, and then examine the tail count and uncertainty. If a fixed simulation is inconclusive, use a larger, independently generated set of shuffles or a valid sequential procedure. Repeatedly checking an ordinary fixed-\(B\) interval and stopping as soon as it falls on one side of \(\alpha\) can invalidate its advertised coverage; optional stopping requires methods designed for sequential monitoring.
Alternatives, Limitations, and Common Mistakes
| Approach | What it provides | Main qualification |
|---|---|---|
| Exhaustive enumeration | The exact tail fraction over all allowable arrangements. | Can be computationally infeasible; the permutation scheme must still be valid. |
| Monte Carlo with plus-one correction | A finite-simulation p-value using \(K\) and \(B\). | Has simulation variability and a minimum value of \(1/(B+1)\). |
| Monte Carlo estimate with interval | An estimate and uncertainty assessment for the conditional tail probability \(q\). | Binomial intervals require independent draws from the specified permutation distribution; intervals do not address model misspecification. |
If permutations are sampled without replacement from a finite collection, \(K\) is not binomial: its distribution is hypergeometric, and the uncertainty calculation should reflect that sampling design. If draws are correlated, or if the allowed permutations are sampled with unequal probabilities without accounting for those probabilities, the simple binomial standard error and Wilson interval are not justified. Ties also matter. Counting ties as at least as extreme is a conventional conservative rule; changing the tie rule after seeing results can distort calibration.
- Mistake: reporting zero when \(K=0\). No sampled shuffle reaching the tail does not establish a zero tail probability. With the plus-one correction the reported value is \(1/(B+1)\), and the simulation still has uncertainty.
- Mistake: treating \(K/B\) and \((K+1)/(B+1)\) as interchangeable. They are different quantities. Report the corrected p-value when using the plus-one convention, and retain \(K/B\) as the observed proportion for binomial uncertainty calculations.
- Mistake: interpreting a confidence interval as uncertainty about whether the null is true. The interval quantifies Monte Carlo uncertainty in \(q\), conditional on the observed data and permutation scheme. It is not a posterior probability for a hypothesis.
- Mistake: assuming more shuffles fix a bad permutation scheme. Increasing \(B\) reduces Monte Carlo error, but cannot make invalid exchanges or restrictions appropriate for the null model.
Check Your Understanding
Answer each question using the stated sampling assumptions; concise calculations are sufficient.
- For \(B=499\) independent shuffles with \(K=0\) exceedances, calculate the plus-one p-value and explain why it is not zero.
- For \(B=999\) and \(K=9\), calculate \(K/B\) and the plus-one p-value. Which value is the corrected Monte Carlo p-value?
- In one sentence, interpret a 95% Wilson interval for the permutation-tail probability \(q\).
- Why is the binomial uncertainty calculation inappropriate if the 999 permutations are dependent draws?
Solutions. (1) \(p_{\mathrm{MC}}=(0+1)/(499+1)=0.00200\); zero simulated exceedances do not show that the tail probability is zero. (2) \(K/B=9/999\approx0.009009\), whereas \(p_{\mathrm{MC}}=10/1000=0.01000\); the latter is the corrected Monte Carlo p-value. (3) Under the specified permutation scheme and binomial sampling assumptions, the interval is an approximate procedure that covers the conditional tail probability \(q\) in 95% of repeated simulations. (4) Dependence violates the independent Bernoulli-draw assumption used to obtain \(K\sim\operatorname{Binomial}(B,q)\), so the binomial standard error and interval need not have their stated properties.