Why Simulation-Based Checks Can Mislead
Simulation-based checks compare observed data summaries with summaries from data generated under a fitted GLMM. Their value depends on whether the simulations reproduce the aspect of the data being checked. A model can appear adequate because the simulation scheme conditions on too much, because the chosen summary is insensitive to the suspected failure, or because a graphical comparison is interpreted as a calibrated test.
The first question is not simply whether the observed statistic falls within a simulated range. It is: what data-generating process does this simulation represent, and what feature of that process should the comparison reveal? A simulation for responses within already observed clusters answers a different question from one that predicts responses in new clusters. Likewise, a check of the overall count total cannot by itself diagnose how counts are distributed across clusters.
Condition on the Right Quantities
In a GLMM, responses within a cluster are linked through a random effect. For example, a random-intercept logistic model has conditional response probability \(p_{ij}(b_j)=\operatorname{logit}^{-1}(x_{ij}^{\mathsf T}\beta+b_j)\), where \(b_j\) is the random intercept for cluster \(j\). Conditional on \(b_j\) and the model parameters, the responses are modeled as independent Bernoulli variables. The model also specifies a distribution for new random effects, commonly \(b_j\sim N(0,\sigma_b^2)\).
A conditional replication holds a fitted or otherwise specified cluster effect fixed and simulates responses within that cluster. It is useful for examining conditional response behavior given the represented clusters, but it does not reproduce the full variation expected when drawing new clusters from the random-effects distribution. A marginal replication draws new random effects and then simulates responses; it is more relevant to variation across clusters drawn from the assumed population. In a Bayesian posterior predictive check, parameter uncertainty may also be represented by drawing parameters from their posterior distribution. These schemes are not interchangeable: the question and the quantities held fixed must be stated.
The distinction matters especially for sparse outcomes. Given a cluster effect \(b_j\), the probability that all \(m_j\) binary responses in cluster \(j\) are zero is \(\prod_{i=1}^{m_j}[1-p_{ij}(b_j)]\). For a marginal prediction, this probability must also be averaged over the random-effects distribution. It is generally not equal to the same expression evaluated at the mean random effect, because the expression is nonlinear in \(b_j\).
This is an application of the law of total probability: first calculate the conditional chance of an all-zero cluster at a given \(b\), then average that chance over the random-effects distribution. Plugging in \(b=0\), or using a fitted cluster effect in place of a new draw, does not perform that averaging.
Worked Example: Expected All-Zero Clusters Under Conditional Replication
Suppose a hypothetical binary-outcome data set has 10 clusters, each with 10 observations. For a conditional replication, the fitted probability of a response is \(p_{\text{low}}=\operatorname{logit}^{-1}(-3.5)=0.0293122\) in five clusters. In the other five clusters, the fitted probability gives an all-zero probability of 0.00873. What is the expected number of all-zero clusters under this conditional scheme?
For each of the first five clusters, conditional independence of the 10 Bernoulli responses gives \(P(\text{all zero})=(1-p_{\text{low}})^{10}\). Substituting the probability, \((1-0.0293122)^{10}=(0.9706878)^{10}\approx0.7426697\), or about 0.7427.
Let \(I_j\) equal 1 if cluster \(j\) is all zero and 0 otherwise. The number of all-zero clusters is \(\sum_{j=1}^{10}I_j\), so its expectation is the sum of the cluster-specific probabilities: \(5(0.7426697)+5(0.00873)=3.7569985\), approximately 3.757 all-zero clusters.
This expectation describes the stated conditional simulation scheme; it is not a probability that the model is correct, nor does it alone determine whether an observed count is surprising. If the diagnostic question concerns new clusters, simulate new random effects and recalculate the reference distribution rather than reusing these fixed cluster effects.
Choose a Statistic That Targets the Concern
A discrepancy statistic is a numerical summary chosen to be sensitive to a feature of the data that could reveal model inadequacy. The statistic should be computed in the same way for observed and replicated data, using the same exposure, grouping, and covariate definitions. Its purpose is not to summarize everything about model fit. A statistic that is insensitive to a particular failure can look reassuring even when that failure is present.
For instance, a total count can be informative about the overall event burden but insensitive to how events are allocated across clusters. If the concern is unusually concentrated counts or excessive between-cluster variation, a cluster-level variance, upper quantile, or number of clusters exceeding a threshold may be more relevant. With unequal exposures or different covariate profiles, raw cluster counts may not be comparable; use a statistic that accounts for the expected differences under the fitted model.
Worked Example: A Matching Total Conceals Different Cluster Patterns
Consider hypothetical counts from four equally exposed clusters. The observed counts are \([0,0,2,8]\), while one replicated data set has counts \([1,1,3,5]\). Both totals are 10, so the total-count discrepancy is identical.
The observed mean cluster count is \(10/4=2.5\). Its sample variance is \(\{(0-2.5)^2+(0-2.5)^2+(2-2.5)^2+(8-2.5)^2\}/(4-1)\) \(=(6.25+6.25+0.25+30.25)/3=14.333\).
The replicated mean is also 2.5, but its sample variance is \(\{(1-2.5)^2+(1-2.5)^2+(3-2.5)^2+(5-2.5)^2\}/3\) \(=(2.25+2.25+0.25+6.25)/3=3.667\).
Thus, matching total counts does not imply matching cluster-to-cluster variation. If the diagnostic concern is concentration, compare an appropriate cluster-level statistic across many replications, not just the total. This raw variance is interpretable here because exposures are equal; for unequal exposures, a model-based standardized discrepancy would be more suitable.
A useful check begins by naming a possible failure and choosing a statistic that would respond to it. For example, examine the number of all-zero clusters for sparse binary outcomes, the upper tail of cluster counts for possible concentration, or within-cluster ordering for residual dependence. A broad collection of plots can reveal unexpected patterns, but selecting a statistic only after seeing an extreme feature makes the comparison harder to interpret as a planned diagnostic.
Predictive Plots Are Evidence, Not Automatic Tests
A predictive-check plot places an observed statistic or pattern beside its distribution across replicated data sets. It can show where the observed result lies relative to model-generated variation. A tail fraction, such as the proportion of replications with a statistic at least as large as observed, is a useful descriptive summary; by itself, it is not automatically a calibrated frequentist \(p\)-value. In posterior predictive checks, the observed data inform the posterior distribution used to generate replications, so observed and replicated quantities are not independent in the way required by a simple null reference test.
The same caution applies to visual intervals. If the observed point falls outside a central 95% interval of replicated values, that is a signal to investigate the targeted feature and the modeling choices. It does not mean there is a 95% probability that the GLMM is false, or that a formal test has rejected a precisely defined null hypothesis at level 0.05. Conversely, lying inside such an interval does not establish adequate fit; a statistic with little sensitivity to a particular failure can easily look unremarkable.
Worked Example: Interpreting a Predictive Tail Fraction
Suppose a hypothetical posterior predictive check uses the maximum cluster count as its discrepancy statistic. In 100 replications, 18 maxima are at least as large as the observed maximum. The empirical upper-tail fraction is \(18/100=0.18\).
A suitable report is: “The observed maximum cluster count was not unusual relative to these posterior predictive replications; 18 of 100 replicated maxima were at least as large.” This is a descriptive comparison for the stated statistic, conditioning scheme, and fitted model.
It is not correct to conclude that the null model has probability 0.18 of being true, or that a formal test has produced \(p=0.18\). The simulation also has Monte Carlo error: treating the 100 exceedance indicators as Bernoulli draws, the estimated standard error of the fraction is \(\sqrt{0.18(1-0.18)/100}=\sqrt{0.001476}\approx0.0384\). Thus, the tail fraction is itself imprecise at this simulation size. More replications can reduce Monte Carlo error, but cannot fix an inappropriate statistic or conditioning scheme.
A Practical Audit Before Interpreting a Check
When a simulation-based check appears reassuring or alarming, audit the setup before attributing the result to a model defect. The following sequence makes the assumptions visible.
Name the feature of the data under scrutiny, such as all-zero clusters, concentration of counts, or within-cluster dependence.
Record whether cluster effects are fixed at fitted values or newly drawn, and whether model parameters are fixed or drawn to represent uncertainty.
Check that the statistic would change if the suspected failure were present and that it respects the design, including unequal exposures when relevant.
Apply the same calculation to observed and replicated data, inspect the distribution and relevant patterns, and describe the result as diagnostic evidence rather than a verdict.
If conclusions change substantially between conditional and marginal replications, that difference is informative: the diagnostic depends on whether the question concerns represented clusters or new clusters. It is not necessarily evidence that one scheme is universally correct. Choose the scheme that corresponds to the scientific prediction or model feature being assessed, and state the choice.
Common Pitfalls and Limitations
| Pitfall | Why it misleads | Correction |
|---|---|---|
| Simulating only conditional on fitted cluster effects when assessing new-cluster behavior | It omits variation from drawing new random effects. | Use a marginal scheme with newly simulated random effects when the target is the distribution across new clusters. |
| Checking only a global total | Different cluster patterns can have the same total. | Add a statistic targeted to the suspected cluster-level feature. |
| Calling a predictive tail fraction a formal \(p\)-value | The calibration and conditioning may not match a formal null test. | Report it as a predictive comparison unless a specific calibrated test has been justified. |
| Trusting an in-range observed value as proof of fit | A weak or irrelevant statistic may not detect the failure of interest. | Assess the statistic’s relevance and inspect targeted patterns, not only whether a value lies inside an interval. |
Simulation checks have further limits. They are conditional on the assumed model structure, the data supplied to the simulation, and the implementation of the generator. A coding error, incorrect exposure, or mismatched grouping can produce a precise-looking but irrelevant comparison. Simulations also have finite Monte Carlo precision. Increasing their number addresses only that last limitation; it does not repair a misspecified model, an unsuitable statistic, or an inappropriate conditioning scheme.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- A researcher wants to assess variation in response rates across new clinics. Should the replication hold fitted clinic effects fixed or draw new clinic effects? Explain.
- Two four-cluster data sets have counts \([0,0,2,8]\) and \([1,1,3,5]\). What does comparing their totals show, and what does it fail to show?
- In a check with 200 replications, 30 replicated statistics are at least as large as the observed statistic. Calculate the empirical upper-tail fraction and state one conclusion that is not justified.
- Why might an observed statistic inside the central 95% of replicated values fail to reassure you about a suspected model problem?
Solutions. (1) Draw new clinic effects, because the target is variation across new clinics under the assumed random-effects distribution. (2) Both totals are 10, so the total comparison cannot distinguish the cluster patterns; a cluster-level statistic such as variability or concentration is needed. (3) The fraction is \(30/200=0.15\). It does not establish that the model is true with probability 0.85, nor is it automatically a formal \(p\)-value. (4) The statistic may be insensitive to the suspected failure, or the conditioning scheme may not represent the target prediction; being inside the interval is not proof of adequate fit.