Optimizing Parameters Through Repeated Adaptive Quadrature

Why Optimization Repeats the Integration

In a generalized linear mixed model, fixed effects and random-effects covariance parameters determine the marginal likelihood. That likelihood is obtained by integrating the random effects out of the conditional response model. An optimizer cannot generally assess a proposed parameter value without evaluating the marginal likelihood at that value, so integration is nested inside optimization: the optimizer proposes parameters, the integration routine evaluates the resulting likelihood, and the optimizer uses that value to choose what to try next.

This nesting matters computationally. A single fit can involve many parameter proposals, with an integral for every independent cluster at each proposal. When adaptive quadrature is used, the quadrature calculation is also tied to the current integrand. Parameters affect the conditional response probabilities and the random-effects distribution, and may affect the quadrature approximation as well.

Definition: Let $Y_i$ denote the responses in cluster $i$, let $b_i$ be its random effects, and write $\theta=(\beta,\phi)$, where $\beta$ is the fixed-effect vector and $\phi$ denotes the parameters of the random-effects covariance matrix $\Sigma(\phi)$. The cluster marginal likelihood is $$ L_i(\theta)=\int f(Y_i\mid b_i;\beta)\,f(b_i;\phi)\,db_i. $$ If clusters are independent conditional on the model parameters, the marginal log-likelihood is $$ \ell(\theta)=\sum_{i=1}^{m}\log L_i(\theta). $$ In practice, quadrature usually supplies an approximation $\widetilde L_i(\theta)$, so the optimizer works with $\widetilde\ell(\theta)=\sum_{i=1}^{m}\log\widetilde L_i(\theta)$.

The integral at a proposed value depends on that value through its integrand, and sometimes through its distributional transformation and adaptive quadrature points. Therefore, the integral must be evaluated, or obtained by a valid calculation for that proposed value. Different integrands do not, by themselves, prove that the resulting integrals differ; the computational requirement is to evaluate the objective correctly at each proposal.

The Inner and Outer Computations

For a given $\theta$, the inner computation forms each cluster’s conditional likelihood and integrates over its random-effects distribution. In multiple dimensions, whitening can express a correlated Gaussian random effect as a linear transformation of standard normal variables. A quadrature rule then approximates the expectation over those variables. For example, a one-dimensional Gaussian expectation can be approximated by a weighted sum at standard-normal quadrature points.

Adaptive quadrature concentrates its evaluation points in a region important to the current integrand, commonly by using information about the conditional integrand’s mode and curvature. The precise adaptation scheme varies. Regardless of scheme, the result at a parameter proposal is a numerical approximation to the integral, not a separate estimate of $\beta$ or $\phi$.

\widetilde L_i(\theta) \approx \sum_{k=1}^{K} w_{ik}(\theta)\, f(Y_i\mid b_{ik}(\theta);\beta)\, f(b_{ik}(\theta);\phi)\,J_{ik}(\theta)

Here the nodes $b_{ik}(\theta)$, weights $w_{ik}(\theta)$, and any transformation factor $J_{ik}(\theta)$ represent the chosen quadrature scheme. In some implementations, parts of this expression are written differently, for example by absorbing a density or transformation into the weights. The formula emphasizes the important point: the finite sum approximates the integral at the current parameter value.

The outer computation compares approximate marginal log-likelihood values across parameter proposals. A gradient-based optimizer may use derivatives of the approximate objective; another optimizer may use objective comparisons or derivative approximations. In either case, each required objective evaluation calls the integration routine for the current candidate. If derivatives are computed through a parameter-dependent quadrature representation, the implementation must account for that dependence correctly or use a valid equivalent derivative calculation.

A Change of Variables Behind Adaptive Quadrature

A short derivation shows how a mode and curvature can place quadrature points. Suppose a one-dimensional integral has the form $I=\int \exp\{h(b)\}\,db$, and $h$ has a local maximum at $b=\hat b$, with $H=-h''(\hat b)>0$. Make the change of variables $b=\hat b+\sqrt{2/H}\,x$. Then $db=\sqrt{2/H}\,dx$, and

I= \sqrt{\frac{2}{H}}\int e^{-x^2} \exp\left\{h\left(\hat b+\sqrt{\frac{2}{H}}x\right)+x^2\right\}\,dx.

An $K$-point Gauss-Hermite rule for the weight $e^{-x^2}$ therefore gives

I\approx \sqrt{\frac{2}{H}}\sum_{k=1}^{K}w_k \exp\left\{h\left(\hat b+\sqrt{\frac{2}{H}}x_k\right)+x_k^2\right\}.

The mode and curvature set the location and scale of the nodes in the original $b$-coordinate. If the parameters change, the function $h$ can change, and so can its mode, curvature, and integral. This derivation explains the structure of one common adaptive rule; it does not imply that every software package uses exactly this transformation.

Worked Example: One Cluster at Two Variance Proposals

Consider a hypothetical cluster with two ordered binary observations, one success and one failure. Use a logistic random-intercept model with fixed intercept $\beta=0$, $b\sim N(0,\tau^2)$, and conditional success probability $p(b)=1/(1+e^{-b})$. The conditional likelihood for this ordered response pair is $p(b)\{1-p(b)\}$. Thus, with $h(b)$ the log of the conditional likelihood times the normal density,

h(b)=-2\log\{2\cosh(b/2)\}-\frac{b^2}{2\tau^2} -\frac12\log(2\pi\tau^2).

The symmetry gives $\hat b=0$. At $\tau=1$, $H=-h''(0)=1+1/\tau^2=2$ would be incorrect: the likelihood curvature is $1/2$, because the second derivative of $-2\log\{2\cosh(b/2)\}$ at zero is $-1/2$. Hence $H=1/2+1=1.5$. Also $h(0)=-2\log 2-\tfrac12\log(2\pi)=-2.3052$.

For two-point Gauss-Hermite quadrature, the nodes are $x_k=\pm1/\sqrt{2}$, each with weight $\sqrt{\pi}/2$. The corresponding original-scale nodes are $b_k=\pm\sqrt{2/H}\,x_k=\pm0.8165$. At either node, $h(b_k)-h(0)=-2\log\cosh(0.4082)-0.3333\approx-0.4956$. Therefore $h(b_k)+x_k^2\approx-2.3008$, and the approximation is

\widetilde L(0,1) =\sqrt{\frac{2}{1.5}}\frac{\sqrt{\pi}}{2} \left[e^{-2.3008}+e^{-2.3008}\right] \approx 0.2050.

Now keep $\beta=0$ but propose $\tau=0.8$. The mode remains zero for this symmetric example, but $H=1/2+1/0.8^2=2.0625$, $h(0)=-2.0821$, and the nodes become $b_k=\pm1/\sqrt{2.0625}=\pm0.6963$. Substitution gives $h(b_k)+x_k^2\approx-2.0797$, so

\widetilde L(0,0.8) =\sqrt{\frac{2}{2.0625}}\frac{\sqrt{\pi}}{2} \left[e^{-2.0797}+e^{-2.0797}\right] \approx 0.2181.

These are two approximate likelihood evaluations for two different covariance proposals, not a claim that all adaptive quadrature rules must give these values. The normalizing constant, curvature scale, nodes, and integrand have all been evaluated for the relevant $\tau$. In a model with several clusters, the routine performs the corresponding calculation cluster by cluster and sums their log likelihoods.

What the Optimizer Does with These Evaluations

The outer optimizer seeks parameter values that maximize the approximate marginal log-likelihood. Its iteration is not simply “integrate once, then optimize the conditional model.” Every new parameter candidate changes the objective being evaluated, so the inner integration is repeated as needed. The following hypothetical trace illustrates a fixed-effect step while holding covariance parameters fixed.

Worked Example: Accepting a Fixed-Effect Proposal

Suppose a hypothetical data set contains 12 independent clinics, each with binary responses, and a random-intercept logistic model is being fitted. At a current fixed intercept $\beta=0$, with covariance parameters held at their current values, the quadrature routine reports $\widetilde\ell(0)=-20.000$. It reports $\widetilde\ell(-0.1)=-20.420$ and $\widetilde\ell(0.1)=-19.620$. These are illustrative outputs of the 12-cluster integration calculation, not observed study results.

A central finite-difference estimate of the slope at zero is

\frac{\widetilde\ell(0.1)-\widetilde\ell(-0.1)}{0.2} = \frac{-19.620-(-20.420)}{0.2} = \frac{0.800}{0.2}=4.000.

The positive slope suggests increasing $\beta$. If a proposed step to $\beta=0.1$ is assessed by the objective, it improves the approximate log-likelihood by $-19.620-(-20.000)=0.380$, so this step is acceptable as an ascent step. The next iteration may use new gradients or curvature information and will call the integration routine again for the new candidate. An optimizer’s eventual stopping decision also depends on its parameter-step and gradient criteria, not just on one improving step.

Propose parameters.
The outer optimizer supplies a candidate $\theta$, including fixed effects and covariance parameters.

Evaluate cluster integrals.
For each cluster, form the current conditional likelihood and approximate its random-effects integral using the selected quadrature rule.

Assemble the objective.
Take the log of each cluster contribution and sum to obtain $\widetilde\ell(\theta)$.

Choose the next action.
Use the objective and, where applicable, derivatives or a curvature approximation to select or assess another parameter proposal.

Check convergence.
Stop only when the optimizer’s stated convergence criteria are met; otherwise continue the nested evaluations.

Assumptions, Diagnostics, and Limitations

The marginal-likelihood expression assumes that the specified conditional response model and random-effects distribution are appropriate, and that the random effects have been integrated over their stated distribution. The sum of cluster log likelihoods additionally relies on the assumed independence structure across clusters. Quadrature introduces a numerical approximation on top of these statistical assumptions. Increasing optimizer precision cannot repair a poor integration approximation, and a highly accurate integral cannot repair a misspecified model.

Useful warning signs include a non-finite objective, abrupt jumps in the approximate objective for small parameter changes, sensitivity of estimates to starting values, repeated optimizer step rejection, and covariance estimates running to a boundary. These signals have several possible causes: an unsuitable parameterization, an optimization problem with weakly identified directions, a difficult integrand, or insufficient numerical accuracy. They call for investigation rather than an automatic conclusion that one particular component has failed.

The integral should be recomputed or validly recalculated at each candidate. Reusing a previous numerical value as if it were the new candidate’s integral is not a valid optimization step. Some computational work can be cached when mathematically justified, but cached quantities must still produce the correct objective for the current parameters. Likewise, convergence of the optimizer means it has met its numerical stopping conditions; it does not establish that the quadrature approximation is accurate or that the model is scientifically adequate.

Repeated adaptive quadrature is one approach, not the only one. Nonadaptive quadrature uses a fixed integration rule and may be simpler, but can place too few useful points where an integrand is concentrated. Laplace approximation replaces the integral by a local approximation around a mode and can reduce computational cost, though its accuracy depends on how well that local approximation represents the integral. Direct or simulation-based integration methods offer other trade-offs. The appropriate comparison depends on dimension, integrand shape, computing cost, and the precision needed for the inferential goal.

A practical limitation is cost: each outer iteration can require many cluster-level integrations, and derivative calculations may require additional work. In higher dimensions, the number of quadrature evaluations can grow rapidly. A fit can also be sensitive to the integration accuracy if the likelihood is flat in a covariance direction. Estimates, standard errors, and likelihood comparisons should therefore be interpreted with attention to both the optimizer’s convergence and the integration method’s numerical adequacy.

Common Mistake

Pitfall: “The integrand changed, so its integral must have changed by the same amount; I can carry the previous integral forward.” This reasoning is incorrect: different functions can have the same integral, and the relation between a change in an integrand and a change in its integral is not generally that simple. The correction is to evaluate or otherwise correctly obtain the integral for every parameter proposal required by the optimizer. Do not infer equality or inequality of integrals from pointwise integrand comparison alone.

Worked Example: Distinguishing Optimization from Integration Accuracy

Suppose an optimizer reports convergence at $\hat\theta$ using a particular quadrature order. The optimization criterion indicates that the approximate objective is locally stable under its stopping rules. It does not show that a coarser quadrature rule would give the same objective or estimates. A careful analyst can refit or reevaluate using a more accurate integration setting and compare the results. If the resulting estimates or likelihood values change materially, the original fit may have been limited by integration error even though the optimizer converged.

This check separates two questions: did the optimizer adequately maximize the objective it was given, and is that objective an adequate approximation to the marginal likelihood? Both matter, but neither answer guarantees the other.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

What are the inner and outer computations in marginal-likelihood optimization for a GLMM?
In the one-cluster example, why did the quadrature nodes change when $\tau$ changed from 1 to 0.8, even though the mode stayed at zero?
For objective values $\widetilde\ell(-0.1)=-20.420$ and $\widetilde\ell(0.1)=-19.620$, calculate the central finite-difference slope using step size 0.1.
Does optimizer convergence establish that the quadrature approximation is sufficiently accurate? Explain.
Why is it invalid to conclude that two integrals differ solely because their integrands differ?

Solutions. (1) The inner computation integrates each cluster’s random effects to approximate its marginal likelihood; the outer computation uses the summed log likelihood to choose parameter proposals. (2) The curvature changed from $H=1.5$ to $H=2.0625$, changing the scale $\sqrt{2/H}$ that maps quadrature nodes to the random-effect scale. (3) The slope is $\{-19.620-(-20.420)\}/0.2=4.000$. (4) No. Convergence concerns the optimizer’s stopping conditions for its objective; integration accuracy requires a separate numerical check. (5) Different functions can have equal integrals, so pointwise inequality alone does not establish inequality of their integrals.