Parametric Bootstrap for Random-Effect Likelihood-Ratio Tests

Why calibrate a variance-component test by simulation?

A likelihood-ratio test (LRT) for a random-effect variance is unusual because the null value, zero, lies on the boundary of the parameter space. For example, a random-intercept variance cannot be negative, so a null hypothesis of zero variance is not an ordinary interior-point test. A reference distribution based on a large-sample approximation can be useful, but its calibration may be imperfect for a particular finite-sample design or fitted model. A null parametric bootstrap estimates the test statistic’s null distribution by repeatedly generating data under the fitted null model and repeating the entire likelihood-ratio calculation.

The crucial phrase is “repeating the entire calculation.” The variance component must be re-estimated in every simulated data set, and nuisance parameters must be estimated again under both the null and alternative. Simulating statistics from a fixed fitted alternative, or merely simulating random effects without generating new responses, does not calibrate the null test.

Models and bootstrap procedure

Definition: Consider a GLMM with observed responses $Y$, fixed design, and random effects $b$. For a test of a single random-intercept variance, the null model fixes $\sigma_b^2=0$; the alternative estimates $\sigma_b^2\geq 0$. The likelihood-ratio statistic is $$ T_{\mathrm{obs}}=2\{\ell_1-\ell_0\}, $$ where $\ell_0$ and $\ell_1$ are the maximized log likelihoods under the null and alternative models, respectively, evaluated on the same data. Under the null parametric bootstrap, simulated data are generated from the null model fitted to the observed data. Both models are then refitted to each simulated data set, producing bootstrap statistics $T_1^*,\ldots,T_B^*$.

For example, with clusters $j=1,\ldots,J$ and observations $i=1,\ldots,n_j$, a binary random-intercept GLMM may be written

Y_{ij}\mid b_j \sim \operatorname{Bernoulli}(p_{ij}), \qquad \operatorname{logit}(p_{ij})=x_{ij}^{\mathsf T}\beta+b_j,\qquad b_j\sim N(0,\sigma_b^2).

The null sets $\sigma_b^2=0$, so the random intercept is absent; the alternative estimates a nonnegative variance. The bootstrap keeps the observed covariates, cluster sizes, offsets, and other design features fixed. For each replicate, it generates new random effects and responses according to the fitted null model. When the null variance is zero, the generated random intercepts are all zero; response generation still follows the null GLMM’s conditional response distribution.

The likelihoods used for the two models must be comparable. In the usual GLMM likelihood-ratio test, fit both by maximum likelihood using the same response and fixed-effects specification. The alternative must enforce the variance constraint $\sigma_b^2\geq 0$. Because the statistic compares maximized likelihoods, a small numerical reversal in which the alternative log likelihood is lower than the null log likelihood should trigger a fit check, not be interpreted as a negative LRT statistic.

Why the bootstrap tail estimates the test’s null probability

Under the fitted null model, each bootstrap data set represents a possible response under the null, with the observed design held fixed. Repeating the fitting procedure on that data set yields a statistic that reflects both response variability and the boundary constraint imposed on the variance estimate. The proportion of simulated statistics at least as large as the observed statistic estimates how often the fitted null model produces evidence against itself at least as extreme as the observed evidence.

This reasoning can be written as an indicator average. Let $T^*$ denote a statistic generated by the fitted null model, and let $\widehat{\theta}_0$ denote its fitted nuisance parameters. The target approximated by simulation is the conditional tail probability under that fitted model:

P_{\widehat{\theta}_0}(T^*\geq T_{\mathrm{obs}}) \ \approx\ \frac{1}{B}\sum_{r=1}^{B} I(T_r^*\geq T_{\mathrm{obs}}).

A common reported estimate uses the plus-one form

\widehat{p}_{\mathrm{boot}} = \frac{1+\sum_{r=1}^{B} I(T_r^*\geq T_{\mathrm{obs}})}{B+1}.

Including ties in the upper tail is appropriate for a test based on “at least as large.” The plus-one convention avoids reporting a zero p-value when no simulated statistic exceeds the observed statistic and accounts for the observed statistic in a Monte Carlo rank calculation. Since the null parameters are estimated from the observed data, this is a plug-in bootstrap calibration, not a guarantee of an exact finite-sample test.

Fit the observed models.
Fit the null and alternative to the observed response with the same design and compatible likelihoods. Calculate $T_{\mathrm{obs}}=2(\ell_1-\ell_0)$.

Generate under the fitted null.
For each of $B$ replicates, generate a new response from the null model using its fitted nuisance parameters and the original design.

Refit and record.
Fit both models to each generated response and calculate $T_r^*=2(\ell_{1r}^*-\ell_{0r}^*)$, enforcing the alternative’s variance constraint.

Estimate and report.
Count statistics satisfying $T_r^*\geq T_{\mathrm{obs}}$, apply the stated tail-probability convention, and assess Monte Carlo precision.

Worked calculations

Worked Example: Calculating the Observed Statistic

Hypothetical data. Suppose a binary-outcome random-intercept GLMM is fitted to 24 clusters. The fixed-effects specification is identical in both models. The null model fixes the random-intercept variance to zero and has maximized log likelihood $\ell_0=-112.60$. The alternative estimates the variance subject to $\sigma_b^2\geq 0$ and has maximized log likelihood $\ell_1=-110.90$.

The observed statistic is

T_{\mathrm{obs}}=2(\ell_1-\ell_0) =2\{-110.90-(-112.60)\} =2(1.70)=3.40.

The alternative has the larger maximized log likelihood, and the resulting statistic is nonnegative, as expected when the null model is nested within the constrained alternative and both fits are properly optimized. The value 3.40 alone does not determine the p-value: the bootstrap distribution, generated under the fitted null, supplies the calibration.

Worked Example: One Bootstrap Replicate

Hypothetical calculation. Keep the 24 cluster sizes and all observed covariates fixed. Generate a new binary response under the fitted null model. Suppose that, for this simulated response, refitting gives $\ell_0^*=-108.30$ and $\ell_1^*=-107.95$, with the alternative variance estimate on its permitted nonnegative scale.

The replicate statistic is

T^*=2\{\ell_1^*-\ell_0^*\} =2\{-107.95-(-108.30)\} =2(0.35)=0.70.

Since $0.70<3.40=T_{\mathrm{obs}}$, this replicate does not contribute to the upper-tail count. The absolute magnitudes of the two log likelihoods are not compared across data sets; only the within-replicate difference enters its statistic.

Worked Example: Estimating the Bootstrap P-Value

Hypothetical simulation result. In $B=999$ null replicates, suppose 83 statistics are at least $T_{\mathrm{obs}}=3.40$. Then

\widehat{p}_{\mathrm{boot}} =\frac{1+83}{999+1} =\frac{84}{1000} =0.0840.

At a preselected significance level $\alpha=0.05$, this estimated p-value does not provide sufficient evidence to reject the null of zero random-intercept variance. It does not establish that the variance is zero. The result is an estimated tail probability under the fitted null model, with Monte Carlo error as well as model-based uncertainty.

A rough Monte Carlo standard error for a tail estimate near 0.084 is

SE_{\mathrm{MC}} \approx \sqrt{\frac{0.084(1-0.084)}{1000}} = \sqrt{0.00007694} \approx 0.00877.

This quantifies simulation noise approximately; it does not measure uncertainty from model misspecification or from estimating the null nuisance parameters. If a decision lies close to $\alpha$, increasing $B$ can reduce the simulation component of uncertainty.

Assumptions, diagnostics, and limitations

The procedure is meaningful when the null and alternative are correctly specified as nested models, the response-generating distribution and link are appropriate, the random-effects distribution is specified adequately, and the observed design is treated as fixed for the intended inference. The likelihood calculations must use the same observations and compatible likelihood definitions. Under the null, nuisance parameters such as fixed effects and dispersion parameters are estimated and used to generate the bootstrap responses. The resulting calibration is therefore conditional on those plug-in estimates.

The simulation must reproduce the data-generation mechanism under the null, not the alternative. In a GLMM this generally means generating the random effects specified by the fitted null and then generating responses conditional on them. Retain the original pattern of covariates and cluster sizes. Do not resample individual rows as if they were independent, and do not generate responses from the fitted alternative when calibrating this null test.

Useful checks include confirming that the null and alternative have the intended nesting, examining convergence messages for both fits, checking that alternative variance estimates obey the boundary constraint, and inspecting the distribution of the simulated statistics. A large concentration at zero can be reasonable for a boundary test; a negative statistic, implausibly large values, or a substantial number of failed fits demands investigation. Check whether failures arise from the optimizer, separation, or weak identification rather than silently discarding them. If some replicates cannot be fit reliably, the reported tail count may not represent the intended procedure; resolve the fitting issue or clearly qualify the result.

The bootstrap does not repair a misspecified GLMM, guarantee reliable estimation with very few clusters, or make a poorly identified variance component informative. It can be computationally expensive because every replicate requires two model fits. Its p-value also has finite simulation precision, which should be reported or assessed when relevant. A mixture reference distribution is a faster asymptotic alternative; the introductory material on finite-sample calibration discusses why reference-law accuracy can be a concern. Simulation offers design-specific calibration under the assumed fitted null, at the cost of computation and plug-in dependence.

Common mistake

Pitfall: “No bootstrap statistic exceeded the observed statistic, so the p-value is zero.” This confuses a finite simulation with the full null distribution. With $B$ replicates, the simulation only resolves tail probabilities at a finite granularity. Report the replicate count and use a stated convention such as $(1+r)/(B+1)$, where $r$ is the number of simulated statistics at least as large as observed. A zero exceedance count is not proof that the null tail probability is zero.

Another important distinction is between the bootstrap p-value and uncertainty about the variance itself. The procedure tests a specified null using a simulated reference distribution; it is not an interval estimate for the variance component. Conclusions should remain about evidence against the null under the assumed model, not about whether the null is certainly true.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

In each null bootstrap replicate, why must both the null and alternative models be refitted?
For $B=499$ simulations, 24 statistics are at least as large as the observed statistic. Calculate the plus-one bootstrap p-value.
Why should the observed covariates and cluster sizes normally be retained when simulating the null responses?
Suppose one alternative fit returns a maximized log likelihood below its nested null fit. What should be checked?
Give one limitation of a null parametric bootstrap even when the simulation and fitting code are implemented correctly.

Solutions. (1) Refitting both models reproduces the statistic’s definition and accounts for the estimation of nuisance parameters and the constrained variance in each generated data set. (2) $\widehat{p}_{\mathrm{boot}}=(1+24)/(499+1)=25/500=0.0500$. (3) Holding the design fixed makes the generated responses reflect the observed study structure and the intended conditional calibration, rather than adding an unrelated resampling scheme. (4) Check convergence, optimizer settings, the variance constraint, and whether the fits use the same data and comparable likelihoods; a properly optimized nested alternative should not have a lower maximum. (5) The calibration relies on the fitted null model and plug-in nuisance estimates, so misspecification or weak identification can still make the inference unreliable.