Why Diagnose Tests by Simulation?
A test with a nominal significance level of 0.05 is designed to reject a true null hypothesis about 5% of the time under its stated assumptions. In a generalized linear mixed model (GLMM), an asymptotic reference distribution may approximate that rate poorly when there are few clusters, a parameter lies on the boundary, or the likelihood is difficult to fit. Simulation makes the calibration question concrete: repeatedly generate data under a specified null model, apply the test exactly as planned, and estimate the fraction of rejections.
For a variance-component test, the issue is especially visible. A variance cannot be negative, so the null value zero lies on the boundary of the parameter space. A standard Wald test and a likelihood-ratio test using an ordinary chi-square reference may therefore behave differently from a likelihood-ratio test calibrated with an appropriate boundary-mixture reference. The goal is not to declare one reference universally correct: it is to assess how the specified procedures perform under specified GLMM conditions.
Specify the Simulation Before Fitting
Consider testing whether a logistic random-intercept variance is zero. For cluster \(i=1,\ldots,K\) and observation \(j=1,\ldots,n_i\), define the model
Conditional on the random intercepts, responses are independent. The covariates are treated as fixed, clusters are independent, and the random intercepts are independent draws from the stated normal distribution. The null hypothesis is \(H_0:\tau^2=0\); the alternative is \(H_a:\tau^2>0\). The fixed-effect terms are present in both models, so the comparison targets the random-intercept variance rather than a change in the fixed-effect specification.
For the simple one-component boundary case, the asymptotic likelihood-ratio reference under suitable regularity conditions is the equal-weight mixture of a point mass at zero and a chi-square distribution with one degree of freedom. The earlier material on Finite-Sample Calibration of Boundary Tests explains why finite-sample behavior can depart from this limit. The simulation must state its conditions rather than assume that the mixture is exact at every cluster count.
Use a fixed design with \(n_i=8\) observations per cluster, four with \(x_{ij}=0\) and four with \(x_{ij}=1\). One null scenario might use \(K=60\), \(\beta_0=-1\), \(\beta_1=0.3\), and \(\tau^2=0\). A corresponding alternative scenario changes only the variance to, for example, \(\tau^2=0.25\). This separation matters: a null simulation estimates type I error, while an alternative simulation estimates power. Neither rate is meaningful without the data-generating conditions attached to it.
Make the Test Implementations Explicit
A simulation can compare named procedures only after defining what each procedure actually computes. In every replicate, fit the null model with \(\tau^2=0\) and the full model with \(\tau^2\geq 0\), keeping the fixed effects the same. Let \(\ell_0\) and \(\ell_1\) be their maximized log likelihoods. With successful maximization of the nested models, the likelihood-ratio statistic is
The statistic is the same for the two likelihood-ratio procedures; their reference distributions differ. At \(\alpha=0.05\), using the ordinary \(\chi^2_1\) reference rejects when \(T>3.841\). The equal-weight boundary-mixture reference rejects when \(T>2.706\), because its upper-tail probability for a positive \(T\) is one half the \(\chi^2_1\) upper-tail probability. The mixture procedure is thus more permissive at this level.
A Wald test needs even more care because there is no parameterization-independent “Wald test for a boundary variance.” To make a comparison reproducible, define one specific naive benchmark: at an interior full-model fit, compute \(Z_W=\hat{\tau}^2/SE(\hat{\tau}^2)\) using the observed-information standard error on the variance scale, and reject when \(Z_W>1.645\), the one-sided standard-normal cutoff at 0.05. If the estimate is exactly zero, set \(Z_W=0\). If the information matrix does not yield a usable standard error, record an information failure and classify that replicate as non-rejection for the primary rejection-rate summary, while reporting the failure count separately. This operational rule is not a claim that the Wald approximation is valid at the boundary; it makes its behavior inspectable.
Use the same simulated data for all three tests in a replicate. This paired design keeps differences in rejection counts from being driven by different random data sets. Also specify the optimizer, convergence criteria, integration settings, and how failed fits are handled. A test comparison that silently drops difficult fits can give a misleadingly reassuring rejection rate.
From Rejections to Calibration
For a test, let \(I_r=1\) if replicate \(r\) rejects and \(I_r=0\) otherwise. Under independent replications generated from the same null conditions, the \(I_r\) are Bernoulli outcomes with rejection probability \(r_0\). Their mean is \(\hat{r}\), and the estimated Monte Carlo standard error is
This expression follows because a Bernoulli variable has variance \(r_0(1-r_0)\); the variance of the average of \(R\) independent indicators is \(r_0(1-r_0)/R\). Replacing the unknown \(r_0\) with its estimate gives the displayed standard-error estimate. It quantifies simulation noise, not uncertainty about whether the data-generating model is scientifically appropriate.
Worked Example: Comparing Three Tests Under the Null
Suppose a hypothetical simulation uses the null scenario above with \(K=60\), \(n_i=8\), and \(R=5{,}000\) replications. Each test is applied to the same data set within each replicate. Assume the simulation records 132 Wald rejections, 127 ordinary-\(\chi^2_1\) likelihood-ratio rejections, and 242 boundary-mixture rejections.
The empirical rejection rates are, respectively,
For the mixture procedure, the Monte Carlo standard error is
The estimated 0.0484 rejection rate is close to the nominal 0.05, and its difference from 0.05 is small compared with the Monte Carlo standard error. The two other rates are substantially lower in this hypothetical run. This pattern is consistent with the ordinary chi-square cutoff being too high when the boundary-mixture limit is relevant. It is evidence about these specified conditions and implementations, not proof that the Wald procedure or ordinary likelihood-ratio reference is always conservative.
Vary Conditions and Check Failure Modes
A single scenario cannot diagnose robustness. Vary conditions that plausibly change finite-sample behavior, especially the number of independent clusters. For example, keep the cluster size and fixed effects unchanged but repeat the null simulation with \(K=20\). A small-cluster result can reveal whether the mixture approximation itself is inaccurate, even when it is the appropriate asymptotic reference. It is also useful to vary the event frequency, covariate pattern, and nuisance fixed-effect values if they matter for the intended application.
Worked Example: A Small-Cluster Calibration Check
Suppose a second hypothetical null simulation has \(K=20\), still with eight observations per cluster, and \(R=5{,}000\). It records 94 Wald rejections, 77 ordinary-\(\chi^2_1\) likelihood-ratio rejections, and 151 mixture-reference rejections. The empirical rates are
For the mixture test, \(SE_{\mathrm{MC}}=\sqrt{0.0302(0.9698)/5000}=0.00242\). The estimated rate is about 8.2 Monte Carlo standard errors below 0.05, so ordinary simulation noise alone is an unlikely explanation for the gap. This suggests that the asymptotic mixture reference is conservative under this particular small-cluster scenario, or that some feature of the fitting and failure-handling procedure affects the result. The next step is to inspect fit diagnostics and, if justified, assess a finite-sample calibration approach—not to label the nominal level as achieved merely because the mixture is theoretically motivated.
For each scenario, report convergence failures, boundary estimates, unusable standard errors, and the number of replicates included in the denominator. Check whether likelihoods satisfy \(\ell_1\geq\ell_0\) to numerical tolerance; large violations suggest optimization problems. Inspect estimates and test statistics from failed or extreme fits. A failure rule that counts every failed fit as a non-rejection is transparent, but it can depress the reported rejection rate, so also show how many failures occurred and consider sensitivity to defensible alternative handling.
Worked Example: Separating Power from Calibration
Now suppose the \(K=60\) scenario is generated under the alternative \(\tau^2=0.25\), with all other settings and \(R=2{,}000\) unchanged. A hypothetical run records 1,024 Wald rejections, 1,198 ordinary-\(\chi^2_1\) likelihood-ratio rejections, and 1,370 mixture-reference rejections. Their estimated rejection probabilities are
These are estimated powers for this alternative, not type I error rates. The mixture procedure has the largest rejection rate here, partly because its critical value is lower than the ordinary \(\chi^2_1\) cutoff. The comparison is informative only alongside the null results: a test that rejects more often under an alternative may also reject too often under the null. For the mixture estimate, the Monte Carlo standard error is \(\sqrt{0.685(0.315)/2000}=0.0104\).
Interpreting Comparisons and Choosing Alternatives
Because all three procedures can be applied to each replicate, compare their paired rejection indicators as well as their separate rejection rates. If the mixture procedure rejects on many data sets where the ordinary chi-square procedure does not, the difference follows naturally from their cutoffs; it does not mean they produced different likelihood-ratio statistics. For a finer comparison, report the number of replicates rejected by one procedure but not another. Do not infer that a procedure is calibrated just because its rate is numerically close to 0.05 in one run; examine Monte Carlo uncertainty and repeat the assessment over conditions relevant to the analysis.
The main alternatives are larger-sample asymptotic references and parametric bootstrap calibration. The mixture reference is computationally cheap and uses the boundary asymptotics appropriate to the simple variance-component test, but can be inaccurate in finite samples. A parametric bootstrap can target a specified finite-sample null more directly, at greater computational cost and with dependence on the fitted null model and bootstrap implementation. Simulation diagnostics are complementary to either choice: they evaluate the procedure under a chosen data-generating model rather than remove the need to justify that model.
A central limitation is that the result is conditional on the simulated conditions. If the real data have unequal cluster sizes, a different event rate, overdispersion, or a covariance structure not represented in the simulation, the observed calibration need not transfer. Simulation does not establish that a model is true; it answers how a specified procedure behaves when data are generated according to a specified model.
A Practical Simulation Workflow
Specify the null or alternative, the GLMM, fixed covariates, cluster structure, parameter values, and nominal significance level.
Draw independent clusters and responses from the stated model, using the null conditions for calibration and alternative conditions for power.
Fit the same null and full models, compute each test, and follow pre-specified rules for convergence, boundary estimates, and unusable standard errors.
Report rejection counts and rates, Monte Carlo standard errors, failures, and scenario-specific results; investigate unexpected behavior before drawing conclusions.
Check Your Understanding
Answer each question, then compare with the solutions below.
- In 1,000 null simulations, a test rejects 37 times. What is the estimated rejection rate and approximate Monte Carlo standard error?
- At significance level 0.05, what likelihood-ratio cutoffs distinguish the ordinary \(\chi^2_1\) and equal-weight boundary-mixture tests?
- Why should the same simulated data sets be used to compare the three procedures?
- A test rejects 70% of the time when the variance component is positive. Does this establish that its type I error is 0.05?
Solutions. (1) The estimated rate is \(37/1000=0.037\). Its Monte Carlo standard error is \(\sqrt{0.037(0.963)/1000}=0.00597\). (2) The ordinary \(\chi^2_1\) cutoff is approximately 3.841; the equal-weight mixture cutoff is approximately 2.706, since its positive-statistic upper tail is half the \(\chi^2_1\) upper tail. (3) Using the same data pairs each procedure's result within a replicate, so differences are not confounded by random variation between separate sets of simulated data. (4) No. That is an estimated power under the stated alternative. Type I error must be assessed under the null.