Laplace Approximation of the Marginal Likelihood

Why Approximate the Random-Effects Integral?

In a joint model, each subject’s observed longitudinal measurements and event-time outcome contribute information about that subject’s unobserved random effects. The subject-level likelihood contribution used for marginal inference integrates over those effects. This integral is often not available in closed form, particularly when the event hazard depends nonlinearly on the longitudinal trajectory. The Laplace approximation replaces the integral with a calculation based on the integrand near its dominant point.

The approximation is local: it uses the height of the integrand at its mode and the curvature of its log there. It is therefore important to find the mode of the whole integrand, not merely the mode of the random-effects distribution or the conditional likelihood. The earlier material on marginal likelihoods and numerical integration provides the setting for this calculation.

Definition: For subject $i$, let $b_i$ be a $q$-dimensional random-effects vector, let $L_i(\theta,b_i)$ be the subject’s conditional joint likelihood contribution given those effects, and let $b_i$ have density $\phi_q(b_i;0,D)$. The marginal contribution is $$ f_i(\theta)=\int_{\mathbb{R}^q} L_i(\theta,b_i)\phi_q(b_i;0,D)\,db_i. $$ Define the log integrand $$ h_i(b_i;\theta)=\log L_i(\theta,b_i)+\log\phi_q(b_i;0,D). $$ Its mode is $\hat b_i(\theta)=\arg\max_{b_i}h_i(b_i;\theta)$. If $h_i$ is twice differentiable there and $$ J_i(\theta)=-\nabla^2_{b_i}h_i(\hat b_i;\theta) $$ is positive definite, the Laplace approximation is $$ f_i(\theta)\ \approx\ \widetilde f_i(\theta) =\exp\{h_i(\hat b_i;\theta)\}(2\pi)^{q/2}|J_i(\theta)|^{-1/2}. $$

The Gaussian density of the random effects is part of $h_i$, including its normalizing constant. Omitting that density, or its constant, changes the integral being approximated. When subjects are independent under the model, the sample marginal likelihood is the product of their marginal contributions; applying the approximation subject by subject gives $\widetilde L(\theta)=\prod_i\widetilde f_i(\theta)$. The corresponding approximate log likelihood is the sum of the subject-level log approximations.

Why Mode and Curvature Give an Approximation

At an interior mode, the gradient of $h_i$ is zero. A second-order Taylor expansion around $\hat b_i$, with $\delta=b_i-\hat b_i$, gives the local approximation $$ h_i(b_i;\theta)\ \approx\ h_i(\hat b_i;\theta)-\frac{1}{2}\delta^\mathsf{T}J_i(\theta)\delta. $$ Exponentiating turns this local log-quadratic expression into the mode height multiplied by a Gaussian-shaped kernel. Integrating that kernel over $q$ dimensions gives $(2\pi)^{q/2}|J_i|^{-1/2}$, which yields the stated formula.

This justification also explains the roles of the conditions. A positive-definite $J_i$ means the log integrand curves downward in every direction at the mode, so the local Gaussian integral is finite. Smoothness permits the Taylor expansion. For the approximation to represent the full integral well, the region contributing most of the integral should be adequately represented by this local shape. A valid local expansion alone does not guarantee that condition.

Form the complete log integrand.
Add the conditional joint log likelihood and the random-effects log density, retaining constants that depend on model parameters or affect the integral.

Find and check the mode.
Solve $\nabla_{b_i}h_i=0$, then check that the candidate is a maximum, not a minimum or saddle point.

Evaluate local curvature.
Calculate $J_i=-\nabla^2_{b_i}h_i$ at the mode and verify it is positive definite.

Combine height and spread.
Use the log-integrand height and curvature determinant in the Laplace formula; repeat for each subject before summing log contributions.

Worked Calculations

Worked Example: A Scalar Nonlinear Integrand

Consider a hypothetical one-dimensional subject contribution whose log integrand, up to any stated constant, is $h(b)=-b^2/2+4b-\exp(b)$. Approximate $I=\int_{-\infty}^{\infty}\exp\{h(b)\}\,db$. The term $-b^2/2$ supplies Gaussian-type tail decay, and the exponential term makes the log integrand nonlinear.

The first and second derivatives are $h'(b)=-b+4-\exp(b)$ and $h''(b)=-1-\exp(b)$. The mode solves $h'(b)=0$. Numerically, $\hat b\approx1.074$, since $-1.074+4-\exp(1.074)$ is approximately zero. Also $h''(b)<0$ for every $b$, so this stationary point is the unique maximum.

At the mode, $h(1.074)\approx-1.074^2/2+4(1.074)-\exp(1.074)\approx0.7922$, and $J=-h''(1.074)=1+\exp(1.074)\approx3.927$. The one-dimensional formula gives $$ \widetilde I=\exp(0.7922)\sqrt{\frac{2\pi}{3.927}} \approx 2.208\times1.265 \approx 2.793. $$ This is the approximate integral, not just the integrand’s value at its mode. The curvature factor accounts for the width of the region around the mode that contributes to the integral.

Worked Example: A Gaussian Case as a Check

Suppose a hypothetical subject has a scalar random intercept $b\sim N(0,1)$ and a measurement $y=1.2$ with conditional distribution $Y\mid b\sim N(b,0.25)$. Here the conditional likelihood is a proper normal density, so all normalizing constants are included. Find the marginal density of this measurement using Laplace approximation.

The log integrand is $$ h(b)=-\frac12\log(2\pi)-\frac{b^2}{2} -\frac12\log(2\pi\cdot0.25)-\frac{(1.2-b)^2}{2(0.25)}. $$ Its derivative is $-b+4(1.2-b)$, which is zero at $\hat b=0.96$. The second derivative is $-1-4=-5$, so $J=5>0$. Substituting the mode gives $h(0.96)\approx-1.7207$. Thus $$ \widetilde f(1.2)=\exp(-1.7207)\sqrt{\frac{2\pi}{5}} \approx 0.2006. $$

In this case the log integrand is exactly quadratic, so the Laplace calculation is exact apart from rounding. As a check, the marginal distribution is $Y\sim N(0,1.25)$, giving density $(2\pi\cdot1.25)^{-1/2}\exp\{-1.2^2/(2\cdot1.25)\}\approx0.2006$. The mode is the posterior mode of $b$ given this measurement, not the prior mode $0$.

Worked Example: Two Correlated Random-Effects Directions

For a hypothetical two-dimensional integral, suppose that near its unique mode the log integrand is $$ h(b_1,b_2)=0.4-2x^2-xy-1.5y^2-0.05(x^4+y^4), \qquad x=b_1-1,\quad y=b_2+0.5. $$ The quartic terms ensure the log integrand falls in the tails; the cross-product term represents local dependence between directions. Since the first derivatives vanish at $x=y=0$, the mode is $(\hat b_1,\hat b_2)=(1,-0.5)$, with height $h=0.4$.

At the mode, the negative Hessian is $$ J=\begin{pmatrix}4&1\\1&3\end{pmatrix}, \qquad |J|=4(3)-1(1)=11. $$ It is positive definite: its first leading principal minor is $4>0$ and its determinant is $11>0$. The approximation is therefore $$ \widetilde I=\exp(0.4)\frac{2\pi}{\sqrt{11}} \approx1.4918\times1.8944 \approx2.826. $$ The off-diagonal entries matter: using only the product of the diagonal curvatures would replace the determinant $11$ with $12$, producing the wrong local volume factor.

Using and Assessing the Approximation

In practice, the mode and Hessian depend on the model parameters $\theta$. To evaluate an approximate marginal likelihood at a proposed $\theta$, optimize $h_i$ over $b_i$, evaluate its height and Hessian at that solution, and assemble the subject contribution. A numerical optimizer’s success flag is not enough: a small gradient residual, a sensible mode, and positive eigenvalues of $J_i$ are useful checks. In multivariate problems, a nearly zero eigenvalue signals a very broad direction; the determinant factor can then be highly sensitive to small numerical errors.

When maximizing the approximate sample likelihood, the subject modes and curvature terms must be recomputed as $\theta$ changes. Holding them fixed generally does not maximize the stated approximation. This tutorial concerns the approximation itself; differentiating the marginal integral and constructing optimization updates require separate methods.

Laplace approximation is computationally attractive because it uses an optimization and a Hessian rather than a large grid of integration points. It differs from Gauss-Hermite, adaptive, tensor-product, and sparse-grid quadrature, which approximate the integral by weighted evaluations at multiple points. Those quadrature methods can represent more of the integrand’s shape, at greater computational cost, especially as the dimension grows. A Gaussian integrand is a favorable case for Laplace; for a skewed or heavy-tailed integrand, a local quadratic representation may miss important mass. Conversely, high-dimensional quadrature can itself become expensive or difficult to tune.

Potential warning signs include a non-converged optimization, multiple plausible modes, a Hessian that is not positive definite, pronounced skewness around the mode, or a posterior with substantial mass far from it. A single-mode Laplace calculation does not automatically account for other modes. If distinct modes contribute appreciably, a single local expansion may underestimate or otherwise misrepresent the integral. Comparing against an appropriate quadrature calculation on representative subjects is a useful diagnostic when feasible; it is not a proof that the approximation is accurate for every subject or parameter value.

Common Mistakes and Limitations

Using the prior mode instead of the integrand mode. The mode must maximize $\log L_i(\theta,b_i)+\log\phi_q(b_i;0,D)$. The data contribution can shift it substantially.
Dropping a normalizing constant. A constant in $b_i$ still multiplies the integral. Keep it if it depends on parameters or is part of the density whose value is required.
Using the Hessian with the wrong sign. The formula uses $J_i=-\nabla^2h_i$, which must be positive definite at a proper local maximum.
Ignoring dependence between random-effects directions. Use the full determinant $|J_i|$, including off-diagonal entries.
Treating the approximation as an exact identity. It is exact for a quadratic log integrand, but otherwise its accuracy depends on how well the local shape represents the integral’s important region.

The approximation also inherits limitations of the model and data. Sparse information for a subject, weak curvature, or a complicated conditional contribution can make its local expansion unstable. A numerically calculated approximation should therefore be understood as an approximation to the specified marginal likelihood, not as an independent check of model assumptions or data quality.

Key takeaway: For each subject, combine the conditional likelihood with the random-effects density, find the mode of their log product, and use the negative Hessian there to estimate the integrand’s local volume. Mode height alone is not a marginal likelihood.

Check Your Understanding

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

What function must be maximized to find the random-effects mode in a subject-level marginal likelihood?
At a scalar mode, the second derivative of the log integrand is $-4$, and its value is $h(\hat b)=0$. Write the Laplace approximation.
Why is a positive-definite negative Hessian required in the multivariate formula?
What determinant should be used for $J=\begin{pmatrix}5&2\\2&2\end{pmatrix}$?
Name one situation in which a single-mode Laplace approximation may be unreliable and one diagnostic or alternative that can help assess it.

Solutions. (1) Maximize $h_i(b_i;\theta)=\log L_i(\theta,b_i)+\log\phi_q(b_i;0,D)$, including the random-effects density. (2) Here $J=4$, so $\widetilde I=\exp(0)\sqrt{2\pi/4}=\sqrt{\pi/2}\approx1.253$. (3) It ensures downward curvature in every direction and a finite local Gaussian volume factor. (4) $|J|=5(2)-2(2)=6$. (5) Multiple important modes or strong skewness can undermine a single local expansion; inspect optimization and curvature diagnostics and, when feasible, compare with quadrature.