Approximation Error in Fixed-Effect and Variance Estimates

Why Integration Accuracy Affects Estimates

In a generalized linear mixed model, the likelihood for the observed responses is obtained by integrating over the random effects. Except in special cases, this integral is not available in closed form. A numerical integration rule therefore approximates the likelihood, and the optimizer estimates parameters by maximizing that approximation. The resulting numerical error can shift both fixed effects and random-effects variance parameters; it does not have to affect them equally.

The distinction matters in practice. A coefficient can appear stable while a variance estimate changes appreciably, or the reverse. Moreover, a small difference between fits is evidence about numerical sensitivity, not proof that either fit is statistically unbiased or that the model is correct.

Likelihood Approximation and Parameter Estimates

Suppose there are \(m\) independent clusters. Let \(Y_i\) be the response data in cluster \(i\), \(b_i\) its random effect, and \(\beta\) the vector of fixed effects. Let \(D\) denote the random-effects covariance matrix. The marginal likelihood integrates the conditional response likelihood against the random-effects distribution:

L(\beta,D)=\prod_{i=1}^{m}\int p(Y_i\mid b_i;\beta)\,\phi(b_i;D)\,db_i, \qquad \ell(\beta,D)=\log L(\beta,D).

Here \(\phi(b_i;D)\) is the Gaussian random-effects density. A numerical integration rule with accuracy setting \(K\) replaces each integral by an approximation, giving an approximate log likelihood \(\ell_K(\theta)\), where \(\theta\) collects the model parameters. The estimate \(\hat\theta_K\) maximizes \(\ell_K\). Changing \(K\) can therefore change the objective function and its maximizer.

Definition: Approximation error in an estimate is the difference between an estimate obtained by maximizing a numerically approximated likelihood and the estimate that would be obtained by maximizing the exact marginal likelihood, under the same data and model. Since the exact maximizer is usually unknown, changes across integration accuracies are used to assess sensitivity; they are not automatically the exact approximation error.

This comparison assumes the same observed data, conditional model, random-effects distribution, parameterization, and optimization target at every accuracy setting. The optimizer must have converged at each setting. Otherwise, a change may reflect an optimization failure or a model change rather than integration error. When possible, use the same scaling, starting-value strategy, and convergence tolerances. For adaptive rules, update the integration center and scale consistently with the parameter values being optimized.

How a Small Likelihood Error Moves the Maximizer

A local calculation explains why a seemingly small likelihood error can move an estimate. Write the approximate log likelihood as \(\ell_K(\theta)=\ell(\theta)+e_K(\theta)\), where \(e_K\) is its numerical error. Suppose the exact likelihood has an interior maximizer \(\hat\theta\), so its score is zero there. Let \(H\) be the Hessian matrix of second derivatives of \(\ell\) at \(\hat\theta\), and let \(g\) be the gradient of \(e_K\) there. A first-order expansion of the approximate score around \(\hat\theta\) gives:

0\approx \nabla\ell_K(\hat\theta+\Delta) \approx H\Delta+g, \qquad \Delta\approx-H^{-1}g.

This is a local approximation, not a guarantee that the shift is small. It shows that the shift depends on the derivative of the likelihood error, not just on its overall size, and on the curvature of the exact likelihood. A direction with weak curvature can be sensitive to a modest score error. For variance parameters near a boundary, this interior quadratic argument may not apply.

Worked Example: A Local Shift in a Fixed Effect and a Variance

Consider a hypothetical two-parameter fit with \(\theta=(\beta,\eta)\), where \(\beta\) is a log-odds coefficient and \(\eta=\log(\sigma)\) for a random-intercept standard deviation \(\sigma\). Suppose the exact maximizer is \((0.80,\log(0.60))\), the local Hessian is \(H=\begin{pmatrix}-100&0\\0&-25\end{pmatrix}\), and the approximate likelihood error has gradient \(g=(2,-1)^\mathsf{T}\) at that point. These stipulated values are for illustration.

Using \(\Delta\approx-H^{-1}g\), the fixed-effect change is \(-(-1/100)(2)=0.020\), and the change in \(\eta\) is \(-(-1/25)(-1)=-0.040\). Thus the locally approximated parameter values are \(\beta_K\approx0.80+0.020=0.820\) and \(\eta_K\approx\log(0.60)-0.040\).

Converting back to the standard-deviation scale gives \(\sigma_K\approx0.60\exp(-0.040)=0.5765\). In this illustration, the fixed effect moves upward by 0.020 log-odds units, while the random-intercept standard deviation moves downward by about 0.0235. The variance-scale change is not obtained by simply subtracting 0.040: the perturbation was specified on the log-standard-deviation scale.

A Controlled Comparison Across Quadrature Accuracies

A useful numerical comparison changes integration accuracy while holding everything else fixed. For illustration, consider a hypothetical binary random-intercept model:

Y_{ij}\mid b_i\sim\operatorname{Bernoulli}(p_{ij}),\qquad \operatorname{logit}(p_{ij})=\beta_0+\beta_1x_{ij}+b_i,\qquad b_i\sim N(0,\sigma^2).

Suppose the hypothetical data consist of 40 clusters with 8 observations per cluster; each cluster has four observations at \(x=0\) and four at \(x=1\). The total number of responses equal to 1 is 62 of 160 at \(x=0\) and 88 of 160 at \(x=1\). The cluster-level arrangement also affects the fitted random-intercept model, so these totals alone do not determine its fit. The following table is a stipulated, illustrative set of fit outputs, not a claim that the totals uniquely reproduce them. \(K\) denotes the number of quadrature nodes used per cluster for this one-dimensional random effect.

Nodes per cluster, \(K\)	\(\hat\beta_0\)	\(\hat\beta_1\)	\(\hat\sigma\)	Maximized approximate log likelihood
3	-1.02	0.72	0.91	-218.463
7	-0.98	0.76	0.83	-217.912
15	-0.97	0.77	0.80	-217.884
31	-0.97	0.77	0.80	-217.883

Worked Example: Quantifying Changes Between Quadrature Orders

Compare the 7-node and 31-node fits in the table. For the fixed slope, \(\hat\beta_{1,7}-\hat\beta_{1,31}=0.76-0.77=-0.010\) log-odds units. For the random-intercept standard deviation, \(\hat\sigma_7-\hat\sigma_{31}=0.83-0.80=0.03\). Relative to the 31-node value, that is \(0.03/0.80=0.0375\), or 3.75%. The maximized log likelihood changes by \(-217.912-(-217.883)=-0.029\).

The 15-node and 31-node estimates agree to the displayed precision, while the 7-node standard deviation differs more noticeably. In this constructed comparison, the variance estimate is more sensitive to the lower integration accuracy than the fixed slope. The pattern is informative about numerical stability for these fits; agreement at 15 and 31 nodes is not, by itself, proof of exact integration.

Reading the Size of a Numerical Shift

Parameter changes need context. A fixed-effect difference can be compared with its standard error to assess whether the numerical change is small relative to the estimated sampling uncertainty. That comparison does not turn numerical error into sampling uncertainty; it provides a scale for judging practical impact. For variance parameters, standard errors can be especially difficult to interpret near zero, so also report the scale and boundary behavior of the estimate.

Worked Example: Comparing a Shift With Standard Errors

Continue with the hypothetical quadrature fits above. Suppose the 31-node fit has \(SE(\hat\beta_1)=0.20\) log-odds units and an estimated standard error of 0.18 for \(\hat\sigma\). The 7-to-31-node fixed-slope shift is \(0.010/0.20=0.050\) of its standard error. The standard-deviation shift is \(0.03/0.18=0.167\) of its reported standard error.

These ratios suggest that, in this illustration, both numerical changes are small relative to the reported standard errors, with the variance shift larger on that scale. They do not establish that the standard errors include quadrature error, nor that the 31-node estimates equal the exact-likelihood estimates. If the estimate is close to zero, a relative percentage can also become unstable or misleading; report the absolute change and the boundary context.

Practical Checks, Alternatives, and Limitations

For an estimation comparison to be interpretable, fit the same model at multiple integration accuracies, confirm optimizer convergence for each fit, and record the fixed effects, variance parameters, and maximized approximate log likelihood. Check whether conclusions important to the analysis change, rather than judging only from rounded agreement in one coefficient. If estimates fail to stabilize, do not treat the lowest-order result as definitive. Increased integration accuracy may help, but computational cost also rises, particularly for higher-dimensional random effects.

Laplace approximation and adaptive Gaussian quadrature are alternatives for evaluating the marginal likelihood, with different accuracy and computational trade-offs. Sparse-grid rules can be useful in moderate dimensions. A controlled comparison should use methods and settings appropriate to the random-effects dimension; estimates from different integration methods are not automatically comparable if other parts of the fitting procedure also differ. None of these numerical methods removes uncertainty from finite data, and convergence across numerical settings does not diagnose misspecification of the conditional model or random-effects distribution.

Common pitfall: “The high-order quadrature fit gives a larger variance, so low-order quadrature always underestimates random-effects variance.” This conclusion is not justified. The direction and size of numerical shifts depend on the data, model, parameter values, and approximation; in another fit, the change could go the other way. Compare controlled fits, report the observed direction, and avoid treating one example as a general rule.

Check Your Understanding

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

If an approximate log likelihood differs from the exact log likelihood by a small constant that does not depend on the parameters, does that constant shift the maximizer?
In the local shift formula, what two features determine the first-order estimate change?
A random-intercept standard deviation is 0.50 at one accuracy and 0.46 at a higher accuracy. Compute the absolute and relative change, taking 0.46 as the reference.
Why is agreement between two high quadrature orders not proof that the model is correct?

Solutions. (1) No. A parameter-independent constant changes the log-likelihood value but not its maximizer, because it has zero gradient. (2) The gradient of the likelihood approximation error and the curvature matrix of the exact log likelihood, through \(\Delta\approx-H^{-1}g\). (3) The absolute difference is \(0.50-0.46=0.04\); relative to 0.46 it is \(0.04/0.46\approx0.0870\), or 8.70%. (4) The comparison assesses numerical sensitivity under the fitted model. It does not test whether the conditional response model, random-effects distribution, or other modeling assumptions are appropriate.