Simulation Recovery as a Convergence Check

Why Check a GLMM by Simulation?

A single successful fit does not establish that an optimization procedure works reliably for the design and model at hand. The optimizer may converge to estimates that are systematically far from the generating values; it may fail for a substantial fraction of plausible data sets; or its failures may occur especially often when a particular parameter is difficult to estimate. Simulation recovery makes these issues observable: specify a clustered-data model with known parameters, repeatedly generate data under it, refit the intended GLMM, and compare the fits with the known values.

This is a computational check under specified conditions, not a guarantee about every future data set. Its conclusions depend on the chosen number of groups, cluster sizes, covariate pattern, outcome model, true parameters, and fitting procedure. A useful report therefore describes the simulation design as well as the summaries.

Define the Data-Generating Model First

For example, a random-intercept logistic GLMM for a binary outcome \(Y_{ij}\), where \(i\) indexes observations and \(j\) indexes groups, can be specified as

The generator must state the number of groups, observations per group, covariate values or their sampling rule, and the values of \(\beta_0\), \(\beta_1\), and \(\sigma_b\). It should also state the conditional distribution and link, as above. In this model, conditional on the group effect and covariates, outcomes are independent Bernoulli variables. Group effects are independent draws from a mean-zero normal distribution and are independent of the covariates. The fitted model should use the same intended structure if the goal is to assess recovery under correct specification.

“Converged” also needs an operational definition set before examining results. For example, count a fit as converged only when the fitting procedure reports convergence and returns finite estimates for the parameters being assessed. A singular fit or an estimate at a boundary should be recorded separately; it is not automatically the same thing as optimizer failure. The chosen rule should be applied consistently across all simulated data sets.

Summarize Convergence and Recovery Separately

Let \(R\) be the number of attempted simulations, and let \(C_r=1\) if replicate \(r\) meets the prespecified convergence criterion and \(C_r=0\) otherwise. The estimated convergence rate is the fraction of attempts that succeed:

If the simulation replicates are independent and the success probability is \(q\), then \(\widehat{q}\) is a sample proportion with Monte Carlo standard error \(\sqrt{\widehat{q}(1-\widehat{q})/R}\). This quantifies simulation noise in the estimated rate; it is not a standard error for a GLMM coefficient.

For a parameter \(\theta\), summarize estimates from converged replicates, while stating clearly that the summary is conditional on convergence. Useful measures include the conditional mean estimate, its bias relative to the generating value, and its root mean squared error (RMSE):

where \(R_c=\sum_{r=1}^R C_r\). Bias measures average signed error, whereas RMSE combines systematic error and variability. Also inspect the spread and shape of the estimates: an average can conceal skewness, extreme values, or a cluster of boundary estimates. If interval estimation is part of the intended procedure, assess coverage as the fraction of eligible intervals that contain the generating value, and state whether that fraction is conditional on successful fitting.

A Justification for Monte Carlo Uncertainty

The convergence indicator \(C_r\) is a Bernoulli variable: it equals 1 with probability \(q\), and 0 otherwise. Thus \(E(C_r)=q\) and \(\operatorname{Var}(C_r)=q(1-q)\). For independent replicates, the variance of their mean is the individual variance divided by \(R\):

Replacing the unknown \(q\) with \(\widehat{q}\) and taking the square root gives the Monte Carlo standard error stated above. This explains why a rate near 0.9 from 20 replicates is much less precise than the same rate from 1,000 replicates. Monte Carlo uncertainty describes how much the simulation summary might vary if the repeated experiment were run again; it does not account for whether the simulated model is appropriate for the scientific setting.

Inspect the Pattern of Difficult Fits

A recovery check is most informative when it separates several questions rather than collapsing everything into one success percentage. Record the convergence rate, the number and nature of failures, and recovery summaries for each parameter. Check whether failures or poor recovery are concentrated in replicates with particular simulated outcomes, covariate patterns, or random-effect draws. Such patterns can indicate that the computational procedure struggles in a specific region of the data it is expected to encounter.

For variance parameters, the estimated value may be on or near a boundary even when the optimizer has converged. If the generating variance is positive, repeated estimates at zero or systematic underestimation deserve attention. But a boundary estimate is not, by itself, proof of a coding error or failed optimization: finite samples can provide weak information about variation among groups. Report its frequency separately and interpret it alongside the design and the other recovery summaries.

Simulation results are also sensitive to how the generator and fitter are implemented. If the generator accidentally uses a different link, covariance structure, or covariate coding from the intended model, apparent non-recovery may reflect a mismatch rather than an optimizer weakness. Conversely, simulating only easy data sets can hide weaknesses that matter under the planned design. Confirm that the same intended parameterization is used, and retain the true values and design details with each set of simulation summaries.

What This Check Can and Cannot Establish

Small-data-perturbation checks assess whether a particular fit changes sharply under nearby changes to the observed data. Simulation recovery instead asks how a fitting procedure behaves across data sets generated under a specified model with known parameters. The two checks answer different questions: a stable fit may still be biased across repeated samples, while good average recovery in simulation does not guarantee stability for one difficult observed data set.

Compared with relying only on reported optimizer success, recovery adds a direct comparison with known values. Compared with an asymptotic standard error alone, it can reveal frequent failures, finite-sample bias, and poor interval coverage under the selected design. Its limitations are equally important: it assesses only the scenarios simulated, and good recovery under a correctly specified generator does not validate that generator for real data. Changing the number of groups, cluster sizes, covariate balance, or true variance components may change the results substantially.

A practical simulation report should therefore identify the data-generating model, design, parameter values, number of replicates, convergence rule, and software fitting procedure. It should present convergence and recovery separately, include Monte Carlo uncertainty for rates where useful, and distinguish boundary estimates from failures. The goal is not to demand perfect recovery in every finite sample; it is to learn whether the intended procedure behaves adequately under conditions relevant to its use.