Integration Diagnostics from Order Refinement and Mode Checks

Why Integration Needs Diagnostics

In a generalized linear mixed model, the likelihood for the observed responses is obtained by integrating over the random effects. That integral is rarely available in closed form. Laplace approximation and adaptive quadrature replace it with numerical approximations, so a fitted model can look well behaved while its likelihood, estimates, or standard errors still depend materially on the integration calculation.

The central diagnostic is to repeat the fit with progressively more quadrature nodes and check whether the results stabilize. This check is more informative when paired with inspection of the conditional modes and the curvature of the log integrand: adaptive integration depends on locating the important region of the integrand and approximating its local shape. These checks diagnose numerical integration; they do not establish that the model is scientifically appropriate.

Definition: For cluster $j$, let $y_j$ denote its observed responses, $b_j\in\mathbb{R}^q$ its random effects, and $\theta$ all model parameters. Define the complete log joint density $$ h_j(b_j;\theta)=\log f(y_j\mid b_j;\theta)+\log f(b_j;\theta). $$ The cluster’s marginal likelihood contribution is $$ L_j(\theta)=\int_{\mathbb{R}^q}\exp\{h_j(b_j;\theta)\}\,db_j. $$ All terms depending on $\theta$, including normalizing terms in the random-effects density, must be retained. The conditional mode $\hat b_j(\theta)$ is a maximizer of $h_j(b_j;\theta)$ over $b_j$.

The definition assumes the specified conditional response distribution and random-effects distribution define a proper joint density, and that the integral exists. In the usual clustered GLMM, responses are conditionally independent given the random effects and covariates, and the random effects have their specified distribution. Conditional modes and curvature are computed for fixed $\theta$; during model fitting, $\theta$ itself changes.

What Order Refinement Checks

Adaptive quadrature places nodes around the region that contributes most to the integral, commonly using the conditional mode and local curvature. At each chosen quadrature order, the likelihood is approximated and the model parameters are optimized. Order refinement means repeating that calculation with more nodes while holding the model, data, and relevant optimization settings fixed. For a multivariate integral, the computational cost can rise quickly with order and dimension, so the refinement should be deliberate.

A useful comparison includes the maximized log likelihood, important fixed-effect estimates, variance or covariance estimates, and quantities derived from them that matter for the analysis. If the log likelihood appears stable but a variance estimate or standard error is still changing, the result is not uniformly stable. There is no universal numerical cutoff: the size of a change matters relative to the precision needed for the scientific conclusion and the numerical tolerances of the fit.

Worked Example: Comparing Fits Across Quadrature Orders

Suppose a hypothetical binary random-intercept GLMM is fitted to the same data using adaptive quadrature at orders 5, 9, 15, 25, and 35. The table reports the fitted coefficient for a treatment indicator, the random-intercept variance, and the maximized log likelihood.

Order	Treatment coefficient	Random-intercept variance	Maximized log likelihood
5	0.681	0.64	-31.842
9	0.699	0.71	-31.617
15	0.702	0.73	-31.603
25	0.702	0.73	-31.602
35	0.702	0.73	-31.602

From order 9 to 15, the treatment estimate changes by $0.702-0.699=0.003$, and the log likelihood changes by $-31.603-(-31.617)=0.014$. From order 15 to 25, the treatment estimate changes by $0.702-0.702=0$ at the displayed precision, the variance changes by $0.73-0.73=0$, and the log likelihood changes by $-31.602-(-31.603)=0.001$. From order 25 to 35, all displayed values are unchanged.

The latter comparisons are reassuring for these reported quantities, at this precision. They do not prove that the integral is exact or that no other quantity is sensitive to order. A careful analyst would also check that each optimization converged adequately and that standard errors or other target quantities are stable enough for their intended use. Because the values are hypothetical, this example illustrates the calculation rather than reporting an empirical result.

Conditional Modes and Curvature

An adaptive calculation relies on a useful local description of $h_j$. At a candidate mode, compute its gradient and Hessian with respect to $b_j$. A stationary point satisfies $\nabla_{b_j}h_j=0$. Define the negative Hessian at that point as $$ A_j=-\nabla^2_{b_j}h_j(\hat b_j;\theta). $$ If $h_j$ is twice differentiable nearby, the candidate is stationary, and $A_j$ is positive definite, then it is a strict local maximum. Positive curvature in this sense is not enough on its own: without stationarity, the point need not be a maximum. Positive definiteness also establishes only a local maximum; it does not rule out another, higher mode elsewhere.

The local quadratic expansion explains why curvature matters. For a small displacement $u$ from a stationary point, $$ h_j(\hat b_j+u;\theta)\approx h_j(\hat b_j;\theta)-\tfrac12 u^\mathsf{T}A_ju. $$ If $A_j$ is positive definite, the quadratic approximation has a concentrated Gaussian shape. Integrating that local approximation gives $$ L_j(\theta)\approx \exp\{h_j(\hat b_j;\theta)\}(2\pi)^{q/2}|A_j|^{-1/2}. $$ This is the local Laplace approximation. It also gives a diagnostic interpretation: very small eigenvalues of $A_j$ indicate a broad direction in which the integrand changes slowly, while large eigenvalues indicate sharper curvature. Neither fact alone guarantees that a particular quadrature order is adequate.

Worked Example: Checking a Candidate Mode and Its Curvature

Suppose a hypothetical fitted random-intercept-and-slope model gives a candidate conditional mode for one cluster at $\tilde b_j=(0.31,-0.18)^\mathsf{T}$. At that point, numerical differentiation gives the gradient $(0.020,-0.010)^\mathsf{T}$ and negative Hessian $$ A_j= \begin{pmatrix} 4&1\\ 1&2 \end{pmatrix}. $$ The eigenvalues of $A_j$ are $3+\sqrt{2}\approx4.414$ and $3-\sqrt{2}\approx1.586$, so $A_j$ is positive definite. But the gradient is not zero. Thus this candidate has positive curvature but has not been shown to be a local maximum. The optimizer should continue, or the calculation should be investigated if it reports convergence.

Suppose a refined optimization instead returns $\hat b_j=(0.314,-0.176)^\mathsf{T}$, with gradient $(0.0002,-0.0001)^\mathsf{T}$ and a negative Hessian with the same displayed entries. The gradient is now close to zero at the stated numerical precision, and both eigenvalues of $A_j$ are positive. Together these are reassuring local checks, subject to the software’s convergence tolerances. They support the interpretation of this point as a strict local maximum when the exact gradient is zero; they do not prove it is the only relevant mode. The small eigenvalue also signals a less sharply curved direction than the first, which can make integration more sensitive to numerical details.

In practice, inspect the mode and curvature for clusters or random-effect directions that appear numerically difficult, not just a single typical cluster. Warning signs include a large gradient at a reported mode, a non-positive eigenvalue of the negative Hessian at a purported stationary maximum, eigenvalues close to zero, or large changes in modes as $\theta$ changes slightly. These signs have different implications: a large gradient suggests failed optimization; non-positive curvature may indicate a saddle point or local minimum; near-zero curvature can indicate a broad or weakly identified direction that is difficult to integrate accurately.

Worked Example: Inspecting a Random-Intercept Mode

Consider one hypothetical cluster with 10 binary responses, 8 of which are successes. Let the fixed linear predictor be zero, use a logistic response model with random intercept $b$, and suppose $b\sim N(0,1)$. For this cluster, the complete log joint density is $$ h(b)=8\log\{\operatorname{expit}(b)\}+2\log\{1-\operatorname{expit}(b)\}-\tfrac12b^2-\tfrac12\log(2\pi), $$ where $\operatorname{expit}(b)=1/(1+\exp(-b))$. The final two terms are the log density of the random intercept, including its normalizing constant.

Writing $p=\operatorname{expit}(b)$, the derivative and negative second derivative are $$ h'(b)=8-10p-b,\qquad -h''(b)=10p(1-p)+1. $$ At $b=0.897$, $p\approx0.7103$, so $h'(b)\approx8-10(0.7103)-0.897\approx0.000$. This is an approximate stationary point. Its negative second derivative is approximately $10(0.7103)(0.2897)+1\approx3.058$, which is positive. Thus, at this stationary point, the log integrand has negative curvature and a strict local maximum. A numerical fit should verify the gradient more precisely than this rounded hand calculation.

The curvature corresponds to a local Gaussian standard deviation of approximately $1/\sqrt{3.058}\approx0.572$. That value describes the local quadratic approximation, not the exact spread or a guarantee of quadrature accuracy. Increasing the quadrature order and checking the fitted results remains necessary when integration accuracy matters.

Using the Checks Together

Fit at a baseline order.
Record the objective value, parameter estimates, standard errors or other target quantities, optimizer status, and numerical tolerances.

Refine the order.
Repeat the fit with more nodes, keeping the data and model unchanged. Use comparable optimization settings so differences primarily reflect integration rather than altered fitting conditions.

Compare results and inspect difficult modes.
Compare successive fits at relevant precision; check gradients, curvature eigenvalues, and whether nearby parameter values produce plausible changes in the conditional modes.

Investigate instability before interpreting.
If results keep changing, check optimization convergence and the mode calculations, then use a higher order or an alternative integration calculation as a comparison. Report the sensitivity rather than treating one unstable fit as definitive.

The main alternative is to compare with another integration approximation, such as Laplace approximation, or with a different numerical integration strategy already available for the model. Agreement can be reassuring, but different approximations can share biases, and disagreement does not by itself identify which answer is accurate. Simulation comparisons can assess performance under specified data-generating conditions, but they do not replace checks on the fitted data and model. Quadrature-order refinement directly tests sensitivity to the chosen order; mode and curvature checks help explain why sensitivity may occur.

Failure Modes and Limits

Declaring convergence from one order. A single successful optimizer status does not test integration accuracy. Compare multiple orders.
Checking only the fixed effects. Variance estimates, standard errors, and likelihood-based quantities may still change. Compare the quantities that support the intended conclusions.
Calling any positive-curvature point a mode. A positive-definite negative Hessian implies a strict local maximum only at a stationary point. Check both the gradient and curvature.
Treating stable rounded values as proof. Values can appear identical after rounding while differing at greater precision. Inspect adequate digits and consider the scale relevant to the analysis.
Confusing integration diagnostics with model diagnostics. Stable numerical integration does not show that the link, random-effects distribution, covariates, or dependence structure is appropriate.

These diagnostics have limits. An order sequence can appear stable before reaching the desired accuracy, a local mode search can miss another mode, and numerical derivatives can themselves be inaccurate. In higher dimensions, the cost of increasing the order can be substantial. Consequently, document the orders and tolerances examined, which quantities stabilized, and any unresolved sensitivity. Numerical stability is necessary for trustworthy inference, but it is not sufficient for valid scientific conclusions.

Check Your Understanding

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

Why must a cluster log integrand retain the normalizing term of a random-effects density when that term depends on a model parameter?
A candidate point has a positive-definite negative Hessian but a nonzero gradient. What conclusion is justified?
Between quadrature orders 15 and 25, an estimate is 0.702 at both orders to three decimals. What is the change at the displayed precision, and what does that not establish?
What does a very small positive eigenvalue of the negative Hessian suggest, and what does it not prove?
Why should a stable quadrature result not be taken as evidence that the GLMM is correctly specified?

Solutions. (1) Omitting a parameter-dependent term changes the integrand by a parameter-dependent factor, so it can change the marginal likelihood and its maximizer. (2) The point is not established as a local maximum; verify or improve stationarity before using the curvature conclusion. (3) The change is $0.702-0.702=0$ at the displayed precision. This does not establish exact equality or adequate accuracy for unreported quantities. (4) It suggests a relatively flat local direction and possible numerical sensitivity, but does not by itself prove that quadrature is inaccurate. (5) Integration checks assess numerical calculation under the fitted model; they do not test whether the model’s assumptions or scientific structure are appropriate.