Why Compare Integration Methods by Simulation?
When a generalized linear mixed model is fitted, numerical integration is part of the likelihood calculation. Different integration methods can therefore produce different parameter estimates from the same data. A simulation study can reveal when these differences matter, but only if the comparison distinguishes numerical approximation from sampling variability: even a perfectly evaluated likelihood will not make every finite-sample estimate equal to the parameter that generated the data.
The central design is to generate many data sets from one specified model, fit each data set using Laplace approximation, fixed Gauss–Hermite quadrature, adaptive Gauss–Hermite quadrature, and a higher-accuracy integration benchmark, and compare the results replicate by replicate. The benchmark is a numerical reference, not the true parameter and not an exact likelihood. The known generating parameter supplies a separate reference for assessing statistical performance.
Simulation Model and Comparison Targets
For a concrete design, consider a binary random-intercept logistic model. Let \(i=1,\ldots,m\) index clusters and \(j=1,\ldots,n_i\) index observations within cluster \(i\). Generate
The random intercepts are independent across clusters, and responses are conditionally independent given the random effects and covariates. For the simulation, set \(\beta_0=-1\), \(\beta_1=0.8\), and \(\sigma=0.9\), with \(m=100\) clusters and \(n_i=4\) observations per cluster. Let \(x_{ij}\) alternate between 0 and 1 within each cluster. These choices, including the covariate pattern and cluster sizes, are part of the design, not details to leave implicit.
Each fitting method approximates the cluster-integrated likelihood. For example, the exact contribution of cluster \(i\) is the integral of its conditional likelihood against the random-intercept density. Laplace approximation, fixed-order quadrature, and adaptive quadrature approximate this integral differently. Use the same model specification, data, parameter constraints, and optimization tolerances for every method. For the fixed-quadrature fit, specify the number of nodes in advance; for the adaptive fit, specify its node count as well as its adaptation convention.
The benchmark should evaluate the same marginal likelihood with substantially higher numerical accuracy than the methods being compared. In this one-dimensional random-intercept setting, a high-order adaptive Gauss–Hermite calculation is a practical choice. Its accuracy should be checked independently—for example, by confirming that still-higher-order evaluation changes the maximized log likelihood and parameter estimates by less than a prespecified tolerance. The benchmark remains an approximation, so report its construction and accuracy target rather than describing it as exact.
Worked Example: Generating One Hypothetical Cluster
Consider one cluster with covariates \(x_{ij}=(0,1,0,1)\). Suppose its simulated random intercept is \(b_i=0.45\). With the specified parameters, the linear predictors are
Applying the inverse logit gives \(p_{ij}=1/(1+\exp(-\eta_{ij}))\), so the probabilities are approximately \((0.366,0.562,0.366,0.562)\). To make the response generation explicit, take hypothetical uniform draws \(u=(0.70,0.40,0.80,0.30)\) and set \(Y_{ij}=1\) when \(u_j\le p_{ij}\), and 0 otherwise. The resulting responses are \((0,1,0,1)\). Repeating the same generation process with fresh random effects and uniform draws gives one complete simulated data set.
What to Measure
For each method \(M\), let \(\hat\theta_{M,r}\) be its estimate on replicate \(r\), and let \(\hat\theta_{B,r}\) be the benchmark estimate on that same replicate. Here \(\theta\) can be a fixed effect or the random-effect standard deviation. The paired numerical discrepancy \(\hat\theta_{M,r}-\hat\theta_{B,r}\) assesses how much the method differs from the benchmark on identical data. Summarize it with its mean and root mean squared discrepancy.
These bias and RMSE measures compare estimates with the known generating value \(\theta_0\), not with the benchmark. Also report interval coverage: the proportion of simulated confidence intervals that contain \(\theta_0\). State how those intervals were constructed, and give the denominator if some fits failed. A method can have small numerical discrepancies but still show finite-sample bias or poorly calibrated intervals; these are different findings.
Convergence is part of the comparison, not a reason to silently discard inconvenient replicates. Report the number of successful fits and the failure rate for each method, together with how failures were defined. If summaries use only converged fits, say so. When failure rates differ, compare both the conditional-on-convergence accuracy and the unconditional computational reliability.
A Paired Comparison and Its Rationale
Worked Example: Summarizing Five Paired Fits
Suppose five hypothetical data sets produce the following estimates of \(\beta_1\). The values are illustrative; each row uses the same data set for every method.
| Replicate | Benchmark | Laplace | Fixed quadrature | Adaptive quadrature |
|---|---|---|---|---|
| 1 | 0.74 | 0.70 | 0.73 | 0.74 |
| 2 | 0.91 | 0.88 | 0.90 | 0.91 |
| 3 | 0.65 | 0.62 | 0.64 | 0.65 |
| 4 | 0.82 | 0.79 | 0.81 | 0.82 |
| 5 | 0.77 | 0.75 | 0.76 | 0.77 |
For Laplace, the five paired differences from the benchmark are \(-0.04,-0.03,-0.03,-0.03,-0.02\). Its mean difference is \((-0.04-0.03-0.03-0.03-0.02)/5=-0.030\), and its RMSE relative to the benchmark is
The fixed-quadrature differences are all \(-0.01\), giving mean difference \(-0.010\) and RMSE \(0.010\). Adaptive quadrature matches the displayed benchmark values in all five rows, giving RMSE 0 at the table’s two-decimal reporting precision; this does not establish exact equality at greater precision. These summaries measure agreement with the benchmark, not bias relative to the generating value.
Pairing matters because method differences are computed on the same simulated data. If \(D_r=\hat\theta_{M,r}-\hat\theta_{B,r}\), then
Both estimates usually respond similarly to how informative or unusual a particular data set is, so their covariance is often positive. Pairing therefore removes much of the between-data-set sampling variation from the comparison and makes numerical differences easier to detect. Fitting the methods on unrelated simulated data sets would mix integration differences with random differences between data sets.
Worked Example: Separating Estimation Error From Numerical Error
Use the five benchmark estimates above and the known generating value \(\beta_1=0.8\). The benchmark errors are \(0.74-0.8=-0.06\), \(0.91-0.8=0.11\), \(0.65-0.8=-0.15\), \(0.82-0.8=0.02\), and \(0.77-0.8=-0.03\). Thus the benchmark’s mean estimation error is
and its RMSE against the generating value is
For Laplace, the mean estimate is \(0.748\), so its mean error against \(0.8\) is \(-0.052\). Its method-to-benchmark mean difference is \(-0.030\), as calculated above. These quantities answer different questions: the first combines numerical and sampling effects, while the second isolates average disagreement with the benchmark in these five paired fits. With only five replicates, neither is a stable estimate of long-run performance.
Interpreting a Larger Simulation
For the stated design, imagine \(R=500\) replications, with each method fitted to every generated data set. Suppose the following hypothetical results are obtained. Bias and RMSE are against the generating parameter; the numerical-discrepancy RMSE is against the benchmark. Coverage is the proportion of reported 95% Wald intervals containing the generating value among converged fits.
| Method | Bias for \(\beta_1\) | RMSE for \(\beta_1\) | RMSE numerical difference for \(\beta_1\) | Bias for \(\sigma\) | RMSE for \(\sigma\) | 95% coverage for \(\sigma\) | Failed fits |
|---|---|---|---|---|---|---|---|
| Laplace | -0.006 | 0.151 | 0.028 | -0.104 | 0.205 | 88.4% | 2/500 |
| Fixed quadrature, 5 nodes | -0.004 | 0.149 | 0.006 | -0.018 | 0.171 | 93.2% | 3/500 |
| Adaptive quadrature, 5 nodes | -0.003 | 0.148 | 0.003 | -0.008 | 0.166 | 94.0% | 3/500 |
| High-accuracy benchmark | 0.001 | 0.147 | Reference | -0.004 | 0.164 | 94.2% | 4/500 |
In this hypothetical result, the methods have similar fixed-effect RMSE, while Laplace has greater numerical discrepancy and more downward bias for \(\sigma\). The adaptive method is closer to the benchmark than the five-node fixed method on average, though neither is identical to the benchmark. The coverage values are empirical proportions: for instance, 88.4% means that 88.4% of the converged Laplace intervals contained the generating \(\sigma=0.9\) in these simulated data sets. It does not mean that a particular interval has an 88.4% probability of containing \(\sigma\).
Assumptions, Failure Modes, and Alternatives
The comparison is interpretable only if the data-generating process and fitted model are correctly specified for the question being asked, the random effects are generated from the stated distribution, and the benchmark is appreciably more accurate than the candidate methods. The numerical integration convention and optimization settings must be documented. Random-number generation should be reproducible, and the number of replications should be large enough that reported Monte Carlo summaries are not dominated by simulation noise.
A frequent mistake is to conclude that one method is superior because its average estimate is closest to the generating value in a small simulation. Sampling variation can produce that ranking even when its integration is less accurate. Compare paired numerical discrepancies to the benchmark, and separately compare bias, RMSE, coverage, and failure rates against the generating model. Another mistake is to use a low-order method as its own benchmark; that can make two inaccurate methods appear to agree.
- Too few replications: estimated bias, coverage, and failure rates can fluctuate substantially. Increase \(R\), and report Monte Carlo uncertainty when the ranking is close.
- Unpaired data sets: unrelated sampling noise obscures method differences. Apply every method to the same replicate and summarize paired differences.
- Unreported failures: dropping failed fits without a count can make a method appear more reliable or accurate than it is. Report failures and define the analysis denominator.
- Overgeneralizing one design: results for balanced clusters of size four do not establish performance for sparse outcomes, unequal cluster sizes, or random slopes. Vary design features when those conditions matter to the intended use.
The principal alternatives are not interchangeable. Laplace approximation is computationally economical but can show appreciable approximation error in some designs. Fixed Gauss–Hermite quadrature uses a common integration grid, whereas adaptive Gauss–Hermite quadrature locates and scales the grid using the integrand. A high-accuracy benchmark is used as a reference for numerical comparison, not as a routine replacement for every fitting method. Simulation can characterize these methods under chosen conditions, but it cannot establish a universal ranking or guarantee performance outside those conditions.
Check Your Understanding
Answer each question, then check the solutions below.
- Why should each integration method be fitted to the same simulated data set?
- A method has a small RMSE relative to the benchmark but poor coverage of the generating value. Are these results contradictory? Explain.
- In a simulation with 200 replications, a method fails to converge on 8 data sets. What denominator should be reported for a coverage proportion computed only on successful fits, and what additional information should accompany it?
- For five paired replicates, a method’s differences from the benchmark are \(0.02,-0.01,0.03,-0.02,0.01\). Calculate the mean difference.
Solutions. (1) Pairing holds the data set fixed, reducing sampling variation in the method-to-method contrast. (2) No. Numerical agreement with the benchmark concerns integration error; coverage concerns how often intervals contain the generating parameter and also depends on sampling behavior and interval calibration. (3) The coverage denominator is 192 successful fits. Report the 8 failures, the failure definition, and that coverage is conditional on convergence. (4) The mean is \((0.02-0.01+0.03-0.02+0.01)/5=0.03/5=0.006\).