Parametric Bootstrap Checks for GLMs

Why Bootstrap a Goodness-of-Fit Check?

A goodness-of-fit discrepancy measures how far the observed responses are from what a fitted model would typically produce. In a generalized linear model (GLM), the fitted values depend on estimated parameters. That estimation affects the discrepancy: data are used to fit the model and are then assessed against it. A reference distribution that ignores this fitting step may therefore be a poor description of the discrepancy's variation under the model.

A parametric bootstrap addresses this by repeatedly generating responses from the fitted GLM and refitting the model to each generated data set. Each simulated discrepancy then reflects both response variation and the fitting procedure. The bootstrap provides a model-based check, not proof that the model is correct: its conclusions depend on the assumed response distribution and the way the model is refitted.

Definition and Procedure

Definition: Let \(Y\) denote the response data, \(x\) the fixed design information (including covariates, offsets, and trial sizes), and \(P_{\theta}\) a specified GLM response distribution. Let \(\hat{\theta}(Y)\) be the estimator obtained by fitting the model, and let \(T(Y,\hat{\theta}(Y))\) be a chosen goodness-of-fit discrepancy, with larger values designated as more discrepant. A parametric bootstrap goodness-of-fit check simulates \(Y^{*(b)}\) from \(P_{\hat{\theta}(Y)}\), refits the same model to each \(Y^{*(b)}\) to obtain \(\hat{\theta}^{*(b)}\), and calculates \(T^{*(b)}=T(Y^{*(b)},\hat{\theta}^{*(b)})\). The observed discrepancy is \(T_{\mathrm{obs}}=T(Y,\hat{\theta}(Y))\).

The design is held fixed because the check asks what responses the fitted model would generate for the observed covariate pattern. The model specification, estimation method, treatment of dispersion, and definition of the discrepancy should be consistent between the observed fit and the bootstrap fits. For example, if the observed statistic is the residual deviance, calculate a residual deviance after every refit, not against the original fitted means.

Fit the observed data.
Estimate the GLM parameters and calculate the observed discrepancy \(T_{\mathrm{obs}}\).

Simulate responses.
For each replicate, draw responses from the fitted response distribution using the observed covariates, offsets, and any fixed trial sizes.

Refit and recalculate.
Fit the same model to each simulated response data set and calculate its discrepancy \(T^{*(b)}\).

Compare with the simulated distribution.
Count how many bootstrap discrepancies are at least as large as \(T_{\mathrm{obs}}\), and report the count, the number of replicates, and the resulting tail estimate.

Worked Examples

Worked Example: Grouped-Binomial Logistic Model

Suppose two groups have \(n_i=10\) trials each, and the observed numbers of events are \(y_1=2\) and \(y_2=8\). Fit an intercept-only binomial GLM, so both groups have one common event probability. The maximum-likelihood estimate is

\hat p=\frac{y_1+y_2}{n_1+n_2}=\frac{10}{20}=0.5.

The binomial residual deviance compares each observed proportion with the fitted probability. For a group with \(y\) events in \(n\) trials, its contribution is

2\left\{y\log\left(\frac{y}{n\hat p}\right)+(n-y)\log\left(\frac{n-y}{n(1-\hat p)}\right)\right\}.

For the first group, the contribution is \(2\{2\log(0.4)+8\log(1.6)\}=3.8549\). The second group has the same contribution, with the event and non-event counts reversed. Thus \(T_{\mathrm{obs}}=7.7098\).

In each bootstrap replicate, generate \(Y_1^*\) and \(Y_2^*\) independently from \(\operatorname{Binomial}(10,0.5)\). Then refit the common probability as \(\hat p^*=(Y_1^*+Y_2^*)/20\), and calculate the deviance using \(\hat p^*\). For instance, if one replicate gives \((Y_1^*,Y_2^*)=(4,6)\), then \(\hat p^*=0.5\). Each group contributes \(2\{4\log(0.8)+6\log(1.2)\}=0.4027\), so \(T^*=0.8054\). This replicate is much less discrepant than the observed data. The full check compares the observed value with the distribution from many such refitted replicates.

Worked Example: Poisson Intercept-Only Model

Consider four independent Poisson counts \(y=(0,1,2,3)\) with a common mean. The fitted mean is \(\hat\mu=(0+1+2+3)/4=1.5\). The Poisson residual deviance, using \(0\log(0/\hat\mu)=0\), is

T_{\mathrm{obs}}=2\sum_{i=1}^4 y_i\log\left(\frac{y_i}{\hat\mu}\right) =2\left\{\log\left(\frac{1}{1.5}\right)+2\log\left(\frac{2}{1.5}\right)+3\log\left(\frac{3}{1.5}\right)\right\} =4.4987.

A bootstrap replicate draws four counts independently from \(\operatorname{Poisson}(1.5)\). Suppose the generated counts are \((1,1,2,2)\). Their sum is 6, so refitting gives \(\hat\mu^*=6/4=1.5\). The replicate deviance is

T^*=2\left\{2\log\left(\frac{1}{1.5}\right)+4\log\left(\frac{2}{1.5}\right)\right\} =2\{-0.81093+1.15073\}=0.6796.

This calculation illustrates the refitting rule: the simulated responses determine their own fitted mean. Here the refitted mean happens to equal the original one; in other replicates it will generally differ. Reusing \(1.5\) without refitting would not reproduce the full fitting procedure.

Worked Example: Calculating a Bootstrap Tail Estimate

Suppose the observed discrepancy is \(T_{\mathrm{obs}}=7.7098\), as in the grouped-binomial example. For illustration, a run with \(B=10\) refitted bootstrap data sets produces discrepancies \(0.42,\ 1.10,\ 2.70,\ 3.10,\ 4.20,\ 5.00,\ 6.80,\ 8.10,\ 9.40,\ 12.00\). Three of the ten values are at least as large as the observed value.

A commonly reported finite-simulation estimate is

\hat p_{\mathrm{boot}}=\frac{1+\sum_{b=1}^{B} I(T^{*(b)}\geq T_{\mathrm{obs}})}{B+1} =\frac{1+3}{10+1}=\frac{4}{11}\approx 0.364.

The added one avoids reporting a zero tail estimate when no simulated value exceeds the observation. It does not make the plug-in bootstrap an exact test for a composite model. With only ten replicates, the estimate is also coarse; a real analysis would generally use many more and report \(B\) alongside the estimate.

What the Bootstrap Comparison Establishes

The bootstrap distribution is generated under the fitted model, not under an alternative model. If \(T\) is defined so that large values indicate poor fit, an observation far into its upper tail is evidence that the fitted model does not readily reproduce a discrepancy as large as the one observed. A non-extreme value means only that this check did not find an unusually large discrepancy by this measure; it does not establish that the model is correct.

The distinction between a known parameter and an estimated parameter is important. If the model parameter were known, the observed data and independently simulated data sets would be exchangeable under the null model. In practice, the parameter is estimated from the observed data, so the simulated data are generated using a plug-in estimate. The resulting check is generally approximate, not an exact finite-sample test.

Theorem (rank calibration for a known model): Suppose \(Y^{(0)},Y^{(1)},\ldots,Y^{(B)}\) are independent data sets with the same distribution \(P_\theta\), where \(\theta\) is known, and let \(T\) be the same deterministic statistic applied to each data set. If the \(B+1\) values of \(T\) have no ties with probability one, the rank of \(T(Y^{(0)})\) among these values is uniform on \(\{1,\ldots,B+1\}\). With ties, counting ties as at least as extreme gives a conservative upper-tail rank comparison.

Proof. Independence and the common distribution make \(Y^{(0)},\ldots,Y^{(B)}\) exchangeable: permuting their labels does not change their joint distribution. Applying the same deterministic statistic to each preserves exchangeability. When there are no ties, each of the \(B+1\) values is equally likely to occupy any given rank, so the rank of the observed value is uniform. With ties, an observation tied at the boundary is counted as extreme rather than non-extreme. This can only increase the upper-tail count relative to a random ordering that breaks the tie, giving a conservative comparison. \(\square\)

For a fitted GLM, the observed statistic is calculated using \(\hat\theta(Y)\), while each simulated statistic uses \(\hat\theta(Y^{*(b)})\). But the simulated data are generated from \(P_{\hat\theta(Y)}\), not from the unknown true \(P_\theta\). The observed and simulated statistics are therefore not generally exchangeable. The rank result explains the form of the finite-simulation comparison; it does not justify treating a plug-in bootstrap result as an exact p-value.

Result (Monte Carlo approximation): Conditional on the observed data and hence on the fitted model, suppose the bootstrap replicates are independent draws from \(P_{\hat\theta(Y)}\), and the statistic is computed by the same procedure for every replicate. For any fixed \(t\), the fraction of replicates with \(T^{*(b)}\leq t\) converges almost surely, as \(B\) increases, to \(P_{\hat\theta(Y)}\{T(Y^*)\leq t\}\).

Proof. Conditional on the observed data, the indicators \(I(T^{*(b)}\leq t)\) are independent and identically distributed Bernoulli random variables with success probability \(q=P_{\hat\theta(Y)}\{T(Y^*)\leq t\}\). Their mean is the empirical fraction in the result. The strong law of large numbers gives

\frac{1}{B}\sum_{b=1}^{B} I(T^{*(b)}\leq t)\xrightarrow{a.s.}q.

Thus, increasing \(B\) reduces Monte Carlo error in approximating the distribution under the fitted model. It does not remove the approximation involved in substituting \(\hat\theta(Y)\) for the unknown parameter. \(\square\)

Assumptions, Reporting, and Common Pitfalls

A parametric bootstrap requires a fully specified probability model for generating responses. For a binomial GLM, this includes the binomial trial sizes; for a Poisson rate model, it includes how exposure enters the mean. If dispersion is estimated, the simulation and refitting procedure must specify how it is handled. A quasi-likelihood fit alone does not supply a complete response distribution for simulation; an additional distributional assumption is needed.

Fit failures are also part of the procedure, not values to silently discard. If simulated data regularly produce separation in a logistic model, boundary estimates, or non-convergence, report the frequency and the rule used to handle those replicates. Such failures may indicate that the fitting procedure is unstable for data generated under the fitted model.

Pitfall: Keep the original estimates for every replicate. This treats simulated responses as if the parameters had been known and omits the variability and adaptation caused by fitting. Refit the specified model to each simulated data set.
Pitfall: Resample observed rows instead of generating from the fitted response distribution. That is a different bootstrap procedure. A parametric GLM bootstrap simulates responses from the model while keeping the observed design information fixed.
Pitfall: Call the result an exact p-value. The plug-in estimate does not give the exact rank calibration for a known parameter. Describe it as a bootstrap tail estimate or model-based check, and state the number of replicates.

In reporting a check, state the fitted response model, what was held fixed, the discrepancy used, the number of bootstrap replicates, the observed discrepancy, and the fraction of replicates at least as extreme. This makes clear what behavior the model was asked to reproduce and how much simulation uncertainty remains.

Check Your Understanding

Use the definitions and results above to answer the following questions.

In a parametric bootstrap for a GLM with fixed covariates, which parts of the data are simulated and which are held fixed?
Why must the model be refitted to each bootstrap response data set when the observed discrepancy is based on a fitted model?
If \(B=199\) bootstrap discrepancies include 11 values at least as large as the observed discrepancy, calculate the plus-one upper-tail estimate.
What assumption makes the rank-calibration result exact, and why does fitting the parameter from the observed data generally violate that condition?
Why does increasing the number of bootstrap replicates not eliminate the plug-in approximation?