Why Compare the Approximations?
A shared random intercept connects a subject’s repeated longitudinal measurements to their event-time outcome. After conditioning on that intercept, the two parts of the likelihood can be evaluated together; to obtain the subject’s marginal likelihood contribution, we must integrate over the intercept’s distribution. In general, this integral has no closed form.
Laplace approximation and Gauss–Hermite quadrature target the same integral, but use different information about its integrand. Laplace summarizes the integrand near its mode using a local quadratic approximation. Quadrature evaluates it at a weighted set of points. Adaptive Gauss–Hermite quadrature places those points around the posterior mode and scales them using local curvature. When the posterior is nearly symmetric and concentrated, the two methods can agree closely. When it is skewed, their difference can be important for likelihood comparisons and parameter estimation.
The Shared-Intercept Integral
For subject \(i\), let \(b_i\) be a shared random intercept with \(b_i\sim N(0,\tau^2)\). One possible longitudinal model is \(Y_{ij}\mid b_i\sim N(x_{ij}^{\mathsf T}\beta+b_i,\sigma^2)\). Suppose the event hazard conditional on the same intercept is \(\lambda_i(t\mid b_i)=\lambda_0(t)\exp(\gamma b_i)\). Write \(T_i\) for the observed event or censoring time and \(\delta_i\) for the event indicator, equal to 1 for an observed event and 0 for right censoring.
Assume that, conditional on \(b_i\) and the model parameters, the longitudinal observations and event-time outcome contribute their specified conditional likelihoods. The subject’s conditional joint contribution is their product. The marginal contribution integrates that product against the random-intercept density. In log form, define the kernel \(h_i(b;\theta)\) as the log of the conditional joint contribution multiplied by the density of \(b_i\), where \(\theta\) collects the model parameters. Then the marginal contribution is the integral of \(\exp\{h_i(b;\theta)\}\) over \(b\).
For example, the event-time part of \(h_i\), including the event indicator and censoring, is \(\delta_i\log\lambda_i(T_i\mid b)-\int_0^{T_i}\lambda_i(u\mid b)\,du\), with the appropriate modification if the observation process begins after delayed entry. The normal random-effects density also belongs inside \(h_i\). Omitting it would approximate a different integral and hence a different likelihood.
These approximations require a well-defined integral and a finite, sufficiently smooth log kernel near an interior mode. The usual local Laplace formula also requires negative curvature at that mode, so \(H_i>0\). A single-mode approximation can be unreliable when the kernel has multiple important modes, a boundary mode, or substantial skewness. For adaptive quadrature, the posterior kernel must be evaluated accurately at the quadrature points, including the random-effects density and all conditional likelihood terms.
How the Two Calculations Differ
Laplace approximates \(h_i(b;\theta)\) by its second-order Taylor expansion at \(\hat b_i\). Since the first derivative is zero at an interior mode, this expansion is \(h_i(\hat b_i;\theta)-H_i(b-\hat b_i)^2/2\). Integrating the exponential of this quadratic gives the Laplace expression. Thus, the approximation replaces the actual kernel with a symmetric Gaussian-shaped kernel determined only by its height and curvature at the mode.
For adaptive Gauss–Hermite quadrature, use the change of variable \(b=\hat b_i+\sqrt{2/H_i}\,x\). If \(x_k,w_k\) are nodes and weights for the convention \(\int e^{-x^2}f(x)\,dx\approx\sum_{k=1}^K w_k f(x_k)\), the resulting approximation is:
The \(x_k^2\) term compensates for the \(e^{-x^2}\) weight used in Gauss–Hermite quadrature. With only the quadratic approximation to \(h_i\), the exponential inside the sum is 1, and the approximation reduces to the corresponding weighted approximation to the Laplace integral. With the actual \(h_i\), the weighted evaluations can reflect departures from that quadratic shape, including skewness. This is the main reason adaptive quadrature can be more accurate when the posterior is asymmetric.
Worked Example: A Local Laplace Calculation
Suppose a hypothetical subject’s log kernel at a fixed parameter value has mode \(\hat b_i=0.30\), height \(h_i(\hat b_i;\theta)=-1.20\), and negative second derivative \(H_i=4.00\). These values satisfy the local curvature condition \(H_i>0\). The Laplace approximation is:
This calculation uses only the kernel’s height and curvature at the mode. It does not establish that the actual posterior is symmetric or that its tails are adequately represented by the local quadratic approximation. An adaptive quadrature calculation would evaluate the actual kernel at several points around \(0.30\), with scale \(1/\sqrt{4.00}=0.5000\), and could therefore give a different result.
Comparing Likelihoods in a Hypothetical Joint Model
Worked Example: Skewness Changes the Approximate Comparison
Consider one hypothetical subject with longitudinal measurements \(y_{i1}=0.2\) and \(y_{i2}=1.3\) at times 0 and 1, followed by an observed event at time \(1.4\). Suppose the longitudinal mean is \(\beta_0+\beta_1t+b_i\), the residual standard deviation is \(0.6\), the random-intercept standard deviation is \(0.8\), and the event hazard is \(0.2\exp(\gamma b_i)\). These specifications define the conditional joint model; the event contribution uses \(\delta_i=1\). Compare two hypothetical candidate parameter values, \(\theta_A\) and \(\theta_B\), with the other model components held fixed.
Suppose numerical evaluation of the subject kernel gives the following mode, curvature, and marginal-contribution approximations. The Laplace values use the formula above; the adaptive quadrature values use 15 Gauss–Hermite nodes. The table reports the resulting numerical evaluations, rounded to four significant figures.
| Candidate | Laplace contribution | Adaptive quadrature contribution |
|---|---|---|
| \(\theta_A\) | 0.7215 | 0.7410 |
| \(\theta_B\) | 0.7350 | 0.7291 |
For Laplace, the corresponding log contributions are \(\log(0.7215)=-0.3264\) and \(\log(0.7350)=-0.3079\). Since there is one subject, these are also the approximate sample log likelihoods. Laplace therefore ranks \(\theta_B\) higher: its approximate log likelihood exceeds that at \(\theta_A\) by \(-0.3079-(-0.3264)=0.0185\).
For adaptive quadrature, \(\log(0.7410)\approx-0.2998\) and \(\log(0.7291)\approx-0.3159\). Its approximate log likelihood is higher at \(\theta_A\), by about \(\log(0.7410)-\log(0.7291)=0.01619\). The ranking changes. In this hypothetical comparison, the discrepancy warns that the Laplace result may be sensitive to the nonquadratic shape of the subject kernel. It does not by itself prove which candidate is correct; quadrature should also be checked for stability as the number of nodes increases.
In a sample with multiple subjects, compute each marginal contribution and sum their logarithms. Do not sum the contributions on the likelihood scale and then take a shortcut based on one subject’s result. A small approximation error for each subject can accumulate in the sample log likelihood, and the errors need not be identical across candidate parameter values.
Checking Agreement and Recognizing Failure
Worked Example: Using Node-Count Stability as a Diagnostic
Suppose a hypothetical fit gives the following approximate sample log likelihoods at one fixed parameter vector. The same adaptive rule is used at each node count, and the Laplace value is included for comparison.
| Method | Approximate sample log likelihood |
|---|---|
| Laplace | -128.42 |
| Adaptive Gauss–Hermite, 5 nodes | -128.10 |
| Adaptive Gauss–Hermite, 9 nodes | -128.02 |
| Adaptive Gauss–Hermite, 15 nodes | -128.01 |
The change from 9 to 15 nodes is \(-128.01-(-128.02)=0.01\), while the 15-node result differs from Laplace by \(-128.01-(-128.42)=0.41\). The quadrature sequence appears much more stable at its last two reported orders than it is close to the Laplace value. This points to an approximation difference rather than evidence that the 15-node quadrature is still changing substantially. It is a diagnostic, not a proof of exactness; the adequacy of a node count depends on the integrand and should be assessed for the likelihood comparisons that matter.
A useful comparison holds the model, data, and parameter values fixed, then evaluates the same subject integrals by both methods. Check whether the mode is interior, whether the curvature is positive in the negative-second-derivative convention used here, and whether the adaptive quadrature result changes materially as nodes are added. Also inspect subject-level differences: a few highly skewed or weakly informed subjects can drive a larger sample-level discrepancy.
The introductory material on Gauss–Hermite quadrature and adaptive quadrature explains how the nodes and weights are selected. Here, the important comparison is what the approximations retain: Laplace uses local height and curvature; adaptive quadrature samples the actual kernel around that local reference. Ordinary, nonadaptive quadrature is another option, but it can waste nodes when the posterior is concentrated far from the random-effects distribution’s center. In higher dimensions, tensor-product quadrature can become expensive; sparse-grid methods address that computational burden, while Laplace remains computationally attractive.
Common Pitfall and Limitations
Wrong reasoning: “Laplace and adaptive quadrature differ, so the quadrature result must be exact.” Correction: Neither approximation is automatically exact. A discrepancy indicates that the local quadratic representation and the weighted evaluations of the actual kernel yield different results. Check quadrature stability and numerical implementation; do not treat the largest node count tried as a proof of correctness.
Another mistake is comparing a Laplace value at one parameter vector with a quadrature value at another, then attributing the difference to method accuracy. The parameter vector must be held fixed for a method comparison. Likewise, compare log likelihoods on the same scale and with the same set of subject contributions. Reporting a ranking alone can conceal a small or practically important difference, so retain the numerical likelihood difference and its scale.
Laplace can be accurate when each subject’s posterior is concentrated and close to Gaussian, and it is often much cheaper than high-order integration. It may be less accurate for skewed, heavy-tailed, broad, or multimodal kernels. Adaptive quadrature can account for local departures from a quadratic shape, but it still uses finitely many evaluations and can miss important features if its nodes or mode are unsuitable. Both methods also rely on the shared-intercept joint model and its conditional likelihood being correctly specified; numerical integration cannot repair a misspecified longitudinal or event-time model.
Check Your Understanding
Answer each question, then compare with the solutions.
- What quantities determine the one-dimensional Laplace approximation at a subject’s posterior mode?
- Why does the adaptive Gauss–Hermite formula include the term \(+x_k^2\) inside its exponential?
- In a comparison, Laplace gives contributions 0.60 and 0.62 at candidates A and B. Adaptive quadrature gives 0.63 and 0.61. What should be concluded?
- Name two features of a subject’s log kernel that can make a single-mode Laplace approximation unreliable.
Solutions. (1) The log-kernel height \(h_i(\hat b_i;\theta)\) and negative second derivative \(H_i>0\) at the mode determine the approximation \(\exp\{h_i(\hat b_i;\theta)\}\sqrt{2\pi/H_i}\). (2) The change of variable expresses the integral using the \(e^{-x^2}\) weight of Gauss–Hermite quadrature; \(+x_k^2\) cancels that weight when evaluating the remaining actual kernel. (3) The approximate ranking reverses: Laplace favors B, while adaptive quadrature favors A. Hold the parameter values fixed and check quadrature stability before deciding whether the discrepancy is consequential; neither result is automatically exact. (4) Examples include substantial skewness, multiple important modes, a boundary mode, or a broad kernel that is poorly represented by its local quadratic expansion.