Why Bootstrap a Fixed-Effect Test?
A likelihood-ratio test compares the maximized likelihoods of two nested models. For a fixed-effect hypothesis in a GLMM, its usual large-sample reference distribution is often a chi-square distribution, with degrees of freedom equal to the number of restrictions. That approximation can be unreliable when the sample is modest, the response is sparse, or the fitted model is complex. A parametric bootstrap offers a direct way to approximate the test statistic's null distribution for the design and fitted null model at hand.
The central discipline is to generate outcomes under the null hypothesis, not under the alternative model. Each generated data set is then analyzed just like the observed data: both nested models are refitted, and their likelihoods are compared. This refitting step allows the bootstrap distribution to reflect the estimation involved in the actual test.
Models, Statistic, and Assumptions
Let \(Y_{ij}\) denote the response for observation \(i\) in cluster \(j\), and let \(b_j\) denote that cluster's random effect or random-effects vector. Write the null and alternative models as \(M_0\) and \(M_1\). They must be nested: every distribution allowed by \(M_0\) must also be allowed by \(M_1\). The models have the same response observations, fixed design information, and random-effects structure; the alternative adds or frees the fixed-effect terms being tested.
For example, the conditional response distribution may be specified by a GLMM with linear predictor \(x_{ij}^{T}\beta+z_{ij}^{T}b_j\), where \(b_j\sim N(0,D)\). The likelihood used for model comparison is the marginal likelihood, obtained by integrating over the random effects. For a fixed-effect comparison, fit both models by maximum likelihood (ML), not by comparing restricted likelihoods.
Here \(\ell_1\) and \(\ell_0\) are the maximized marginal log-likelihoods under \(M_1\) and \(M_0\). The statistic is nonnegative when the models are properly nested and both fits attain their maxima. Under regularity conditions, the usual large-sample null reference is \(\chi^2_q\), where \(q\) is the number of independent fixed-effect restrictions. This tutorial uses simulation to calibrate that statistic, rather than assuming the chi-square approximation is adequate.
The bootstrap requires a specified null model, a valid simulation mechanism for its response distribution, and a fitting procedure that can be applied consistently to the observed and simulated data. The observed design—including covariates, cluster membership, and number of observations per cluster—is held fixed. The null model's nuisance parameters, including its random-effect covariance and any dispersion parameter, are estimated from the observed data and treated as fixed plug-in values when generating replicates. Thus, this is an approximation to the null distribution at the fitted null parameter values, not an exact test that removes all uncertainty in nuisance-parameter estimation.
Bootstrap Procedure and Justification
Fit \(M_0\) and \(M_1\) by ML to the observed responses using the same records and compatible optimization settings. Calculate \(T_{\mathrm{obs}}\).
Retain the ML estimates of the null fixed effects, random-effects distribution parameters, and other null-model nuisance parameters.
For each replicate, draw random effects from the fitted null random-effects distribution, then draw responses from the fitted conditional response distribution. Keep the observed design and cluster layout unchanged.
Fit both \(M_0\) and \(M_1\) to that replicate by ML and calculate \(T_b^*=2(\ell_{1b}^*-\ell_{0b}^*)\).
Repeat for \(B\) replicates and count how many statistics are at least as large as \(T_{\mathrm{obs}}\).
The generation step follows the null model's hierarchical data-generating process. For a Bernoulli GLMM, for instance, draw each cluster's random effect and then draw each binary response from its conditional Bernoulli distribution. Drawing responses only from fitted marginal means while omitting the random effects would generally fail to reproduce the fitted GLMM's dependence structure.
The indicator \(I(\cdot)\) equals 1 when the condition holds and 0 otherwise. The plus-one estimate avoids reporting a zero p-value just because no simulated statistic exceeded the observed one. It is a Monte Carlo estimate based on a plug-in null model; it should not be described as an exact finite-sample p-value for the original GLMM. The approximate Monte Carlo standard error of a tail-probability estimate \(\hat p\), based on \(B\) independent replicates, is \(\sqrt{\hat p(1-\hat p)/B}\). This measures simulation error, not uncertainty from estimating the null model or possible model misspecification.
Worked Examples
Worked Example: Observed Statistic for a Two-Parameter Test
Suppose a hypothetical clustered binary-response data set is analyzed using two nested logistic GLMMs. Both include the same random-intercept structure and are fitted by ML. The null model omits two treatment-related fixed-effect terms; the alternative includes both. Assume the fitted marginal log-likelihoods are \(\ell_0=-126.40\) and \(\ell_1=-121.85\).
The observed statistic is
There are \(q=2\) fixed-effect restrictions. Under the regular asymptotic reference, \(T_{\mathrm{obs}}\) is compared with a \(\chi^2_2\) distribution, giving \(P(\chi^2_2\geq9.10)=e^{-9.10/2}\approx0.0106\). This is the large-sample result, not a bootstrap result. To bootstrap, the null model would be used to generate outcomes on the same cluster and covariate layout, and both models would be refitted in each replicate.
Worked Example: Bootstrap Tail Estimate and Monte Carlo Error
For a hypothetical fixed-effect test, suppose \(T_{\mathrm{obs}}=7.80\). The fitted null model is used to generate \(B=1999\) data sets, and both models converge in every replicate. In 42 replicates, \(T_b^*\geq7.80\).
The plus-one estimate is
The approximate Monte Carlo standard error is
In the numerator, \(0.0215(0.9785)=0.02103775\). A rough 95% Monte Carlo range is \(0.0215\pm1.96(0.0032441)\), or approximately 0.0151 to 0.0279. This range describes simulation variability in estimating the tail probability; it is not a confidence interval for a model coefficient. At significance level \(\alpha=0.05\), this bootstrap estimate is below 0.05, so the test rejects \(H_0\) under the chosen calibration. The conclusion concerns evidence against the tested fixed-effect restrictions, not proof of a causal effect.
Worked Example: A Result Near the Decision Threshold
In a second hypothetical analysis, a one-parameter fixed-effect test has \(T_{\mathrm{obs}}=3.40\). With \(B=999\) null bootstrap replicates, 86 statistics are at least 3.40. All replicates fit successfully.
The estimated tail probability is
Its approximate Monte Carlo standard error is
The estimate exceeds 0.05, so at \(\alpha=0.05\) the test does not reject \(H_0\). This is not a decision to accept the null. Because the estimate is not far from 0.05 and its Monte Carlo error is noticeable, a larger \(B\) would make the simulated tail estimate more precise. More replicates reduce simulation error but do not resolve misspecification or poor model fitting.
How to Read the Result and Diagnose Problems
A small bootstrap p-value means that, under the fitted null GLMM and the fixed observed design, simulated likelihood-ratio statistics at least as large as the observed statistic were uncommon. It does not give the probability that \(H_0\) is true. It also cannot establish that the null model is correctly specified: the bootstrap reproduces the assumed model, including its limitations.
Check that the model comparison is genuinely nested, uses the same responses and records, and differs only in the fixed-effect restrictions being tested. Confirm that both fits use ML and the same random-effects structure. Review convergence messages, optimizer behavior, and whether fitted likelihoods are sensible. A replicate in which one model fails to converge cannot simply be treated as an ordinary small test statistic or silently removed: either choice can distort the tail estimate. Investigate the failure, improve fitting settings where justified, and report how unresolved failures affect the analysis.
A bootstrap may reveal that the nominal chi-square reference is a poor approximation for the chosen design and fitted null. However, it is not a cure for data separation, weak identification, an incorrect response distribution, an unsuitable random-effects structure, or a boundary issue in nuisance parameters. If many null replicates produce singular fits or unstable estimates, the calibration itself may be unstable and the modeling problem deserves attention before interpreting the p-value.
Comparison, Limitations, and a Common Pitfall
| Approach | What it uses | Main consideration |
|---|---|---|
| Asymptotic likelihood-ratio test | Observed likelihood-ratio statistic and a chi-square reference | Fast and standard under regularity and adequate large-sample approximation |
| Parametric-bootstrap likelihood-ratio test | Null-generated data sets, with both models refitted | Can better reflect finite-sample behavior under the fitted null, at substantial computational cost |
| Wald test | Estimate, estimated standard error, and an approximate reference distribution | Can be sensitive to local curvature and coefficient estimation; it is not the same statistic as the likelihood-ratio test |
| Score test | Likelihood slope and information evaluated under the null | Avoids fitting the full alternative in its usual form, but relies on its own approximation and implementation |
Common pitfall: Simulating from the fitted alternative and counting how often its statistics exceed \(T_{\mathrm{obs}}\) does not calibrate a null test. Those replicates are generated when the tested effect may be present. The correction is to simulate from the fitted null model, then refit both null and alternative models in every replicate. Likewise, generating responses under the null but calculating each replicate's statistic without refitting fails to reproduce the estimation step used to obtain the observed statistic.
The method is computationally intensive: each replicate requires two GLMM fits, and a reliable tail estimate may require many replicates. It is also conditional on the chosen null model and its plug-in parameter estimates. Report the number of replicates, number of exceedances, p-value formula, Monte Carlo uncertainty, fitting failures, and the model assumptions relevant to the simulation.
Check Your Understanding
Answer each question, then compare your reasoning with the brief solutions.
- A null and alternative GLMM differ by three fixed-effect restrictions. What is held fixed in each bootstrap replicate, and how many models are fitted to each replicate?
- For \(B=499\) replicates, 9 statistics are at least as large as \(T_{\mathrm{obs}}\). Calculate the plus-one bootstrap p-value.
- Why should the random effects be drawn during response generation rather than generating responses independently from marginal fitted means?
- A reported bootstrap estimate is 0.04, based on 999 replicates. Give its approximate Monte Carlo standard error.
- One model fails to converge in several replicates. Explain why silently dropping those replicates is not an adequate remedy.
Solutions. (1) Keep the observed covariates, cluster membership, and observation pattern fixed; generate responses under the fitted null, and fit both \(M_0\) and \(M_1\) to each replicate. Thus, two models are fitted per replicate. (2) \((9+1)/(499+1)=10/500=0.0200\). (3) Drawing random effects and then conditional responses preserves the dependence and variability implied by the fitted hierarchical model; independent marginal-mean generation generally does not. (4) \(\sqrt{0.04(0.96)/999}\approx0.00620\). (5) The failures may be related to the simulated statistic or data features, so excluding them can change the estimated tail distribution. Investigate and report the failures, use justified fitting improvements, and assess whether the calibration remains trustworthy.