Calibrating Dispersion Checks With Simulation

Why Calibrate a Dispersion Check by Simulation?

A Pearson or deviance statistic can indicate that a fitted binomial or Poisson model leaves more variation unexplained than its mean–variance specification allows. In many settings, a reference distribution based on a large-sample approximation is used to judge whether the statistic is unusually large. That approximation can be doubtful when the sample is small, the data are sparse, or the fitted model has features that make the approximation unreliable.

Parametric simulation offers a direct alternative calibration. It asks: if the fitted model generated new responses under the same design, how often would the same diagnostic statistic be at least as large as the one observed? The procedure is useful only if the simulation reproduces the model being assessed and the fitting and diagnostic steps are repeated consistently.

Definition and Simulation Procedure

Suppose the observed responses are \(Y_1,\ldots,Y_n\), with fixed covariates or design, and a binomial or Poisson GLM is fitted. Write \(\hat{\mu}_i\) for the fitted mean of a Poisson response and \(\hat{p}_i\) for the fitted success probability of a binomial response with \(m_i\) trials. The usual Pearson statistics are

The corresponding residual deviances are

In these formulas, a term of the form \(0\log(0/a)\), for \(a>0\), is defined to be zero by continuity. Larger values are evidence of a discrepancy between the observations and the fitted model, but they do not by themselves identify the reason for it.

For each of \(B\) simulations, generate responses from the fitted model conditional on the observed design, refit the same model using the same fitting rules, and calculate the same statistic. Denote the resulting values by \(T_1^*,\ldots,T_B^*\), and the observed value by \(T_{\mathrm{obs}}\). A common Monte Carlo upper-tail estimate is

The added one prevents a zero estimate when no simulated value reaches the observed statistic. It does not make a fitted-model simulation an exact test for a composite model: the fitted parameter values are estimates, not known null parameters.

Worked Examples

Why Refit Every Replicate?

The statistic being calibrated is calculated after estimating model parameters from data. If simulated responses are generated from the fitted model but the original estimates are reused for every replicate, the simulated statistics omit the variability introduced by fitting. Refitting mimics the complete analysis applied to the observed data: generate responses, estimate parameters, then calculate the diagnostic.

If a replicate fit fails or reaches a boundary where the statistic is undefined, that event must not be silently discarded. The fitting procedure and a rule for handling such replicates should be stated; frequent failures can indicate that this calibration is not reliable for the problem at hand.

What the Calibration Does—and Does Not—Establish

A simple-null Monte Carlo rank calculation has a precise validity property. Suppose the null distribution is fully specified, the observed statistic and \(B\) simulated statistics are independent draws from that same distribution, and ties are handled by counting simulated values greater than or equal to the observed value. Then the plus-one upper-tail value is conservative: its probability of being at most \(\alpha\) is no greater than \(\alpha\). To see this, let \(N=B+1\) and randomly break any ties among the \(N\) exchangeable statistics. The rank of the observed statistic in this ordering is uniform on \(1,\ldots,N\). The count using “greater than or equal to” is at least that tie-broken rank, so the plus-one tail value is at least the tie-broken rank divided by \(N\). Consequently, the event that the tail value is at most \(\alpha\) is contained in the event that the uniform rank divided by \(N\) is at most \(\alpha\), whose probability is \(\lfloor \alpha N\rfloor/N\leq\alpha\). Thus the test is conservative, including when ties occur. \(\square\)

In a fitted GLM, however, the observed data and the simulated data are not generally exchangeable under a composite null: the simulation uses estimated parameters. The rank argument therefore does not establish exact finite-sample validity for the plug-in calibration. Its purpose is to approximate the diagnostic’s null behavior under the fitted model, especially when a standard asymptotic reference is doubtful.

The Monte Carlo error calculation also has a direct justification. If the model and fitted parameters are held fixed and each replicate is generated independently, define \(I_b=I(T_b^*\geq T_{\mathrm{obs}})\). Each \(I_b\) is Bernoulli with success probability \(q=P(T^*\geq T_{\mathrm{obs}})\), so \(E[I_b]=q\), \(\operatorname{Var}(I_b)=q(1-q)\), and independence gives

This describes simulation error conditional on the fitted model; it does not account for uncertainty in the fitted parameters or for whether the model itself is adequate.

Common Pitfall and Practical Interpretation

Wrong reasoning: “The observed statistic is large, so the fitted model must have overdispersion.” A large statistic relative to a simple reference can instead reflect the behavior of the fitted statistic in a small or sparse sample, or a departure from the assumed mean structure. Simulation calibrates how unusual the statistic is under the chosen fitted model; it does not diagnose the cause of a discrepancy.

Correction: State the model used to generate replicates, keep the design fixed, refit every replicate, and interpret the upper-tail estimate as a model-conditional diagnostic. A small simulated tail indicates that the observed statistic is unusual under that procedure. It is a reason to investigate the model and data, not a standalone proof of extra variation.