Why Simulate a Goodness-of-Fit P-Value?
A Pearson goodness-of-fit test compares observed counts with counts expected under a specified probability model. When expected counts are sufficiently large, a chi-square distribution can provide a useful approximation to the null distribution of the Pearson statistic. With sparse counts, that approximation may be unreliable. In that situation, we can approximate the p-value directly by generating count tables under the null model and seeing how often they produce a statistic at least as large as the one observed.
The simulation does not make the data less sparse or change the null model. It replaces a questionable reference distribution with a reference distribution generated from the stated model. The resulting p-value is still conditional on the model being appropriate: simulation cannot repair incorrect null probabilities, dependence between observations, or an incorrectly specified sampling design.
Here the null probabilities are treated as specified and fixed. Each simulated sample has the same total \(n\) as the observed sample. Its counts are generated from the multinomial model with probabilities \(p_1,\ldots,p_k\), and its Pearson statistic is calculated using the same expected counts \(np_i\). Because the probabilities are fixed, those expected counts are the same for every replicate.
The comparison is in the right tail: larger Pearson statistics indicate greater discrepancy from the null probabilities. The phrase “at least as large” matters. The null distribution is discrete, so simulated samples can tie the observed statistic; ties belong in the count of exceedances.
The Simulation Procedure
Let \(T_{\mathrm{obs}}\) denote the observed Pearson statistic. Generate \(B\) independent samples under \(H_0\), each with total \(n\), and calculate \(T_b\) for simulated sample \(b\). Let \(R\) be the number of simulated statistics satisfying \(T_b\geq T_{\mathrm{obs}}\). The simulated tail proportion \(R/B\) estimates the null tail probability.
A commonly recommended Monte Carlo p-value uses a plus-one adjustment:
The adjustment accounts for the observed data as one member of the comparison set alongside the \(B\) simulated samples. Under the fixed null model, the observed sample and simulated samples have the same distribution, so their ranks are exchangeable. Including the observed statistic avoids reporting a p-value of zero simply because a finite run happened to produce no simulated exceedances. For a conventional simulation report, state which estimate you use; examples below use the plus-one estimate.
Write down the fixed category probabilities and the observed total \(n\). Confirm that the categories and sampling model match the question being tested.
Find \(E_i=np_i\) for every category and calculate \(T_{\mathrm{obs}}=\sum (O_i-E_i)^2/E_i\).
Generate \(B\) independent multinomial count vectors, each with total \(n\) and probability vector \((p_1,\ldots,p_k)\).
Calculate the Pearson statistic for each simulated vector using the same expected counts, then count statistics at least as large as the unrounded observed statistic.
Give \(R\), \(B\), the p-value estimate, and enough information to show that this is a Monte Carlo approximation rather than an exact calculation.
Worked Examples
Worked Example: Simulating a Sparse-Count Test
Suppose 12 hypothetical observations fall into three categories. Under \(H_0\), the category probabilities are \((0.60,0.30,0.10)\), and the observed counts are \((11,1,0)\). These expected counts include values below 5, so we will use simulation rather than rely on chi-square calibration. The hypothetical observations are assumed independent, each assigned to one of the three categories, with the stated probabilities under \(H_0\).
The expected counts are \(12(0.60)=7.2\), \(12(0.30)=3.6\), and \(12(0.10)=1.2\). Thus:
Now suppose a hypothetical run generates \(B=4{,}999\) multinomial samples under the null probabilities, each with total 12, and \(R=402\) simulated statistics are at least as large as the observed statistic. Using the plus-one estimate gives:
The estimated p-value is about 0.0806. Under the specified null model, the simulation estimates that approximately 8.06% of samples of 12 observations would produce a Pearson statistic at least as large as the observed one. This is not the probability that \(H_0\) is true. It is also only a simulation-based estimate; a different run could give a slightly different result.
Worked Example: Counting Ties at the Exact Threshold
For a second hypothetical test, let \(n=10\), let the null probabilities be \((0.50,0.30,0.20)\), and let the observed counts be \((8,1,1)\). Assume independent categorical observations under this fixed null model. The expected counts are \((5,3,2)\), so the observed statistic is:
Suppose \(B=1{,}999\) simulated samples produce \(R=187\) statistics at least as large as the exact observed value \(109/30\). Then:
A simulated sample with counts \((8,1,1)\) has a statistic exactly equal to \(109/30\), so it must be included in \(R\). The threshold for the comparison is the unrounded statistic, even if a display shows \(3.6333\). Comparing against a rounded-up value such as \(3.6334\) could omit tied samples and would not implement the rule \(T_b\geq T_{\mathrm{obs}}\).
Worked Example: Interpreting No Simulated Exceedances
Suppose a hypothetical simulation uses \(B=999\) samples from a fixed null model and none has a Pearson statistic at least as large as \(T_{\mathrm{obs}}\). The model assumes independent observations and uses the same total and null probabilities as the observed data. Here \(R=0\), so the plus-one estimate is:
The run provides evidence that the null tail probability is small, but it does not establish that the probability is exactly zero. The estimate is limited by the number of simulations: with 999 replicates, seeing no exceedances cannot resolve very small tail probabilities precisely. If a decision depends on distinguishing, for example, a p-value below 0.001 from one near 0.001, a substantially larger simulation may be needed. A zero exceedance count is a reason to report the simulation limit, not to claim an impossible certainty.
Why the Simulation Estimates the Desired Probability
For each simulated sample, define \(I_b=1\) if its statistic is at least \(T_{\mathrm{obs}}\), and \(I_b=0\) otherwise. Under the null model, each \(I_b\) is a Bernoulli random variable whose success probability is the null tail probability \(q=P(X^2\geq T_{\mathrm{obs}}\mid H_0)\). Independent simulation runs make the indicators independent. Therefore, the average of the indicators estimates \(q\).
This can be justified directly. Each indicator has mean \(q\) and variance \(q(1-q)\). The average \(R/B\) consequently has mean \(q\) and variance \(q(1-q)/B\). As the number of simulations increases, its variability decreases, and the law of large numbers implies that \(R/B\) approaches \(q\).
Since \(q\) is not known, a rough estimate of the Monte Carlo standard error is \(\sqrt{\widehat{p}(1-\widehat{p})/B}\), using a reported estimate \(\widehat{p}\) near \(R/B\). This describes simulation variability, not uncertainty about whether the null model is correct. Increasing \(B\) reduces Monte Carlo variability, but does not correct a flawed model or sampling design.
What Simulation Does—and Does Not—Fix
When expected counts are sparse, the reference distribution of the Pearson statistic may not be well approximated by the chi-square distribution. Multinomial simulation instead reproduces the discreteness and sparsity implied by the fixed null model. It can therefore give a more direct estimate of the null tail probability than an unreliable asymptotic calibration.
The method still depends on the model. The total sample size must be fixed as described, all listed categories must be possible under the specified model, and the simulated observations must follow the same independence and category-probability assumptions as the real data. If probabilities are estimated from the data rather than fixed in advance, simply simulating from the fitted probabilities and using the fixed-model calculation may not reproduce the appropriate testing procedure. The simulation scheme must then account for that estimation.
A simulated p-value is also not the same as a list of simulated statistics being “the exact null distribution.” A finite simulation is a random approximation to that distribution. Reporting \(B\), \(R\), the comparison rule, and the resulting estimate makes the approximation transparent and reproducible.
Common Pitfalls
- Comparing rounded statistics. Rounding \(T_{\mathrm{obs}}\) before counting can drop ties. Compare each simulated statistic with the unrounded observed value, and include equality.
- Changing the null model during simulation. For a fixed-probability test, every replicate must use the specified probability vector and the same total \(n\). Simulating from different probabilities answers a different question.
- Reporting zero as the p-value. No exceedances in a finite run does not show that the null tail probability is zero. Report the simulation size and an appropriate nonzero estimate, such as the plus-one estimate.
- Treating simulation as a cure for all assumptions. Simulation avoids relying on the chi-square approximation; it does not validate the model, independence, or study design.
- Interpreting the result as a probability that \(H_0\) is true. The estimated p-value is a tail probability under \(H_0\), not a posterior probability of the hypothesis.
Check Your Understanding
Answer each question using the fixed-probability simulation procedure.
- Under a null model with four fixed category probabilities, what total count should each simulated sample have if the observed sample size is \(n\)?
- Why should simulated statistics equal to the observed statistic be included in the exceedance count?
- A run has \(B=4{,}999\) and \(R=49\). Calculate the plus-one Monte Carlo p-value.
- What does it mean if no simulated statistic reaches the observed statistic, and why is a p-value of zero not justified?
- Name one important assumption or model feature that simulation does not validate by itself.