Maximum Likelihood Estimation with Adaptive Quadrature

Why Numerical Integration Enters the Likelihood

In a shared-parameter model, an unobserved subject-level effect can link a repeated outcome process to dropout. Conditional on that effect, the two processes may be modeled separately; the observed-data likelihood for a subject still requires averaging their joint contribution over the distribution of the latent effect. That integral is generally not available in closed form, especially when the dropout model uses a nonlinear link such as the logit.

Adaptive Gauss–Hermite quadrature approximates each subject’s integral using a small set of weighted evaluations. “Adaptive” refers to centering and scaling the quadrature nodes using the shape of that subject’s log integrand. This can put more of the numerical effort where the integrand contributes most to the integral.

Definition: Let $b_i$ be a scalar shared random effect for subject $i$, with $b_i\sim N(0,\tau^2)$. Let $y_i$ denote the observed outcome history and $d_i$ the observed dropout information. If the model specifies the joint conditional contribution $p(y_i,d_i\mid b_i;\theta)$, the subject’s observed-data likelihood contribution is $$ L_i(\theta)=\int p(y_i,d_i\mid b;\theta)\,\phi(b;0,\tau^2)\,db, $$ where $\phi(b;0,\tau^2)$ is the normal density and $\theta$ collects the model parameters. For independent subjects, $L(\theta)=\prod_i L_i(\theta)$, and the maximum likelihood estimate maximizes $\ell(\theta)=\sum_i\log L_i(\theta)$.

The conditional contribution must reflect what was actually observed. For example, a dropout hazard contributes the likelihood for the observed dropout time or censoring information, not a likelihood for an unobserved future dropout time. The integration step does not repair an incorrectly specified conditional model.

Deriving the Adaptive Approximation

For a fixed parameter value, define the subject-specific log integrand $$ h_i(b;\theta)=\log p(y_i,d_i\mid b;\theta)+\log\phi(b;0,\tau^2). $$ Then $L_i(\theta)=\int \exp\{h_i(b;\theta)\}\,db$. Suppose $h_i$ has an interior mode $\hat b_i$, with negative second derivative $$ H_i=-\left.\frac{\partial^2 h_i(b;\theta)}{\partial b^2}\right|_{b=\hat b_i}>0. $$ The mode identifies where the integrand is largest, and $H_i$ describes its local curvature. In a narrow peak, $H_i$ is large; in a broad peak, it is smaller.

Let $z_k$ and $a_k$ be standard-normal Gauss–Hermite nodes and normalized weights, so that $\int g(z)\phi(z;0,1)\,dz$ is approximated by $\sum_{k=1}^K a_k g(z_k)$, with $\sum_k a_k=1$. Set $b=\hat b_i+z/\sqrt{H_i}$. The change of variable gives $db=dz/\sqrt{H_i}$. Since $dz=\sqrt{2\pi}\phi(z;0,1)\exp(z^2/2)\,dz$, the integral becomes

L_i(\theta)= \sqrt{\frac{2\pi}{H_i}}\exp\{h_i(\hat b_i;\theta)\} \int \exp\left\{ h_i\left(\hat b_i+\frac{z}{\sqrt{H_i}};\theta\right) -h_i(\hat b_i;\theta)+\frac{z^2}{2} \right\}\phi(z;0,1)\,dz.

Replacing the remaining normal expectation by its $K$-node rule yields the adaptive approximation:

\widetilde L_i(\theta)= \sqrt{\frac{2\pi}{H_i}}\exp\{h_i(\hat b_i;\theta)\} \sum_{k=1}^K a_k \exp\left\{ h_i\left(\hat b_i+\frac{z_k}{\sqrt{H_i}};\theta\right) -h_i(\hat b_i;\theta)+\frac{z_k^2}{2} \right\}.

The formula follows directly from the change of variable and the quadrature rule. As a useful check, if $h_i$ is exactly quadratic, then $h_i(\hat b_i+z/\sqrt{H_i})-h_i(\hat b_i)=-z^2/2$. The exponent inside the sum is then zero for every node, the weighted sum is one, and the approximation equals the exact Gaussian integral.

Assumptions for this calculation: The conditional model and random-effect density define the intended likelihood; the subject contribution is integrable; the log integrand has a usable interior mode with $H_i>0$; and the quadrature order is adequate for the integrand’s shape. These are numerical and model conditions, not consequences of applying the formula. If the mode is not unique, the curvature is nonpositive, or the integrand is highly skewed or heavy-tailed, the approximation needs particular scrutiny.

Worked Adaptive Quadrature Calculations

Worked Example: Transforming Nodes Around a Subject-Specific Mode

Consider a hypothetical subject contribution whose log integrand has mode $\hat b=0.675$ and curvature $H=1.447$. To illustrate the node transformation, use a two-node standard-normal rule with $z_1=-1,z_2=1$ and $a_1=a_2=0.5$. Suppose the log-integrand shape near and around the mode is represented by $h(\hat b+u)=-2-1.447u^2/2+0.08u^3-0.03u^4$. It has derivative zero and second derivative $-1.447$ at $u=0$, as required for the stated mode and curvature.

The transformed effects are $b_k=\hat b+z_k/\sqrt{H}$. Since $1/\sqrt{1.447}=0.8313$, the two effects are $b_1=0.675-0.8313=-0.1563$ and $b_2=0.675+0.8313=1.5063$. In particular, the first transformed effect is $-0.1563$, not $-0.157$.

For either node, $u^2=1/1.447=0.6911$, $u^3=0.5745$ in magnitude, and $u^4=0.4776$. Thus the log-integrand difference from the mode is $-0.5+0.08u^3-0.03u^4$. It is approximately $-0.4684$ at the positive node and $-0.5603$ at the negative node. Adding $z_k^2/2=0.5$, the correction exponents in the quadrature sum are $0.0316$ and $-0.0603$.

The weighted correction factor is $0.5\exp(0.0316)+0.5\exp(-0.0603)\approx0.9868$. The likelihood approximation is therefore $$ \widetilde L =\sqrt{\frac{2\pi}{1.447}}\exp(-2)(0.9868) \approx 2.0838(0.1353)(0.9868) \approx 0.2783. $$ The factor is close to one because the local quadratic shape is the principal contribution, but the cubic and quartic terms make the integrand nonquadratic. A two-node rule is used here to show the arithmetic, not to recommend that two nodes are adequate for an actual fit.

Worked Example: A Subject Contribution in a Joint Outcome–Dropout Model

Suppose a hypothetical joint model has a continuous outcome and a binary dropout indicator $D_i$, with conditional contribution $$ p(y_i,D_i\mid b_i;\theta)= \frac{1}{\sigma_y}\phi\left(\frac{y_i-\beta-b_i}{\sigma_y};0,1\right) \pi_i(b_i)^{D_i}\{1-\pi_i(b_i)\}^{1-D_i}, \qquad \pi_i(b)=\frac{\exp(\gamma+\alpha b)}{1+\exp(\gamma+\alpha b)}. $$ Here $\pi_i(b)$ is the conditional dropout probability, and $b_i\sim N(0,\tau^2)$ enters both processes. For a subject with observed $y_i=12$ and $D_i=1$, the dropout factor is $\pi_i(b)$, not $1-\pi_i(b)$. The random effect is integrated out only after multiplying the outcome and dropout contributions by its density.

At a trial value of $\theta$, suppose numerical maximization of $h_i$ gives $\hat b_i=0.20$, $H_i=2.00$, and $h_i(\hat b_i;\theta)=-3.00$. Assume the three-node standard-normal rule has nodes $0,\sqrt{3},-\sqrt{3}$ and normalized weights $2/3,1/6,1/6$. For this model evaluation, suppose the log-integrand differences from the mode at the transformed nodes are respectively $0,-1.30,-1.70$. These values are the numerical evaluations of the joint log integrand; the correction adds $z_k^2/2$, which is $0,1.5,1.5$.

The sum in the adaptive rule is $ (2/3)\exp(0)+(1/6)\exp(0.20)+(1/6)\exp(-0.20)$. Substituting $\exp(0.20)\approx1.2214$ and $\exp(-0.20)\approx0.8187$, the sum is approximately $0.6667+0.2036+0.1364=1.0067$. Hence $$ \widetilde L_i =\sqrt{\frac{2\pi}{2.00}}\exp(-3.00)(1.0067) \approx1.7725(0.04979)(1.0067) \approx0.08884. $$ This is one subject’s approximate observed-data likelihood contribution at the trial parameter values. Its logarithm contributes to the objective optimized in maximum likelihood estimation.

Using the Approximation for Maximum Likelihood

For each trial parameter vector, the calculation is repeated: find each subject’s mode and curvature, evaluate the integrand at transformed nodes, form $\widetilde L_i(\theta)$, and sum $\log\widetilde L_i(\theta)$ across subjects. The modes and curvatures depend on the trial parameters, so they must be updated as the optimizer changes $\theta$. Maximizing the resulting approximate log likelihood gives an approximate maximum likelihood estimate. This procedure does not treat estimated random effects as if they were observed data; it approximates their integration in the observed-data likelihood.

A log-scale implementation is usually safer than multiplying many small likelihood contributions. For positive terms $r_k$, a stable calculation of $\log\sum_k r_k$ uses the largest log term $m=\max_k\log r_k$ and evaluates $m+\log\sum_k\exp(\log r_k-m)$. This avoids numerical underflow while preserving the same likelihood.

Worked Example: Locating a Quadrature-Based Estimate

Consider a hypothetical data set of six subjects, each with a continuous outcome history and observed dropout information. Fit the joint model above, treating $\beta$ as the parameter of interest and fixing the other parameters at specified trial values for this profile calculation. Suppose a 15-node adaptive quadrature evaluation produces the following approximate log likelihoods:

Trial $\beta$ (outcome units)	Approximate log likelihood
0.20	-8.200
0.30	-8.050
0.40	-8.000
0.50	-8.050
0.60	-8.200

The largest listed value is $-8.000$ at $\beta=0.40$. The central difference in slope is $\{-8.050-(-8.050)\}/(0.50-0.30)=0$, consistent with a local maximum at $0.40$. A quadratic approximation to the curvature using the two values $0.10$ units on either side gives $$ \ell_Q''(0.40)\approx \frac{-8.050-2(-8.000)+(-8.050)}{(0.10)^2} =\frac{-0.100}{0.010}=-10.0. $$ Thus the profile is locally concave at this point, and the corresponding curvature-based standard error, if the other parameters were fixed, is approximately $1/\sqrt{10.0}=0.3162$ outcome units. The quadrature-based profile estimate is $\hat\beta\approx0.400$ outcome units. In a full fit, all parameters would be optimized jointly and uncertainty would use the relevant full observed-information matrix, not this fixed-parameter profile calculation.

Diagnostics, Alternatives, and Limitations

Quadrature is a numerical approximation, so numerical stability should be checked rather than presumed. Refit or reevaluate the likelihood with a higher quadrature order and compare the log likelihood, parameter estimates, and curvature-based uncertainty. Material changes indicate that the lower-order rule is not yet resolving the integral adequately. Also inspect the subject-specific mode searches: a nonconverged mode, a nonpositive $H_i$, a mode at an artificial search boundary, or multiple plausible modes can invalidate the local centering and scaling used by this formula.

Adaptive Gauss–Hermite quadrature is often more efficient than a nonadaptive rule when subjects’ integrands are concentrated away from zero or have different widths. A nonadaptive rule keeps nodes fixed, which can waste evaluations where the integrand is negligible. Laplace approximation uses only the mode and curvature, effectively retaining the Gaussian local approximation; it can be fast, but may be less accurate when the integrand is skewed or otherwise non-Gaussian. Adaptive quadrature samples beyond the local quadratic approximation, with accuracy controlled in part by the number of nodes. Direct Monte Carlo integration is another alternative, but its simulation error can make likelihood optimization noisy unless the numerical design is carefully controlled.

For a vector random effect of dimension $q$, the scalar curvature $H_i$ is replaced by the positive-definite negative Hessian matrix, and a matrix square-root transformation scales the nodes in all directions. A product rule with $K$ nodes per dimension requires $K^q$ evaluations per subject, so computation can grow rapidly with dimension. Moreover, quadrature addresses integration error only. It cannot establish that the shared-effect distribution, outcome model, dropout model, or conditional-independence assumptions are substantively appropriate.

Using the prior mode instead of the integrand mode. The adaptive center is the mode of the full log integrand, including both the observed-data contribution and the random-effect density.
Omitting the curvature factor. The multiplier $\sqrt{2\pi/H_i}$ is part of the change of variable. Leaving it out changes the likelihood.
Failing to update nodes during optimization. Both $\hat b_i$ and $H_i$ generally change with $\theta$; keeping them fixed does not evaluate the stated adaptive approximation at new parameter values.
Assuming more nodes fix a bad model. More nodes may improve numerical integration, but they do not correct an incorrect dropout contribution or an unsuitable latent-effect distribution.

Key takeaway: Adaptive Gauss–Hermite quadrature approximates each shared-effect likelihood integral by centering at the subject-specific log-integrand mode and scaling by its curvature. Maximum likelihood then optimizes the sum of the log approximations, with numerical adequacy checked by increasing quadrature order and examining the subject-specific calculations.

Check Your Understanding

Answer each question, then compare with the brief solutions.

In the scalar adaptive rule, what do $\hat b_i$ and $H_i$ describe?
Why does the exponent in the quadrature correction include $z_k^2/2$?
If the log integrand is exactly quadratic around its mode, what is the correction factor inside the weighted sum?
What does it suggest if estimates change materially when the quadrature order is increased?
Why must the subject-specific mode and curvature usually be recalculated as the optimizer changes model parameters?

Solutions. (1) $\hat b_i$ is the mode of the full subject log integrand, and $H_i$ is its negative second derivative there. (2) It appears when rewriting integration with respect to $dz$ as an expectation under the standard normal density. (3) The exponent is zero at every node, so the correction factor is one. (4) The lower-order rule may not adequately resolve the integral, so the likelihood and estimates are not yet numerically stable. (5) Both quantities depend on the trial parameters, and the adaptive approximation is defined using their current values.