The Laplace Approximation to the Marginal Likelihood

From a Conditional Model to a Marginal Likelihood

In a nonlinear mixed-effects model, a subject’s observations depend on both population parameters and that subject’s unobserved random effects. The likelihood contribution for the subject is therefore not the conditional likelihood at one selected random-effect value. It is the conditional likelihood averaged over the random-effects distribution. This average is an integral, and the Laplace approximation replaces it with a calculation based on the integrand’s value and curvature near its mode.

The key is to keep track of what is being approximated. The integrand contains both the conditional likelihood and the normalized random-effects density. Its mode determines the height term; the negative Hessian at that mode determines the local volume term. For Gaussian random effects, the density’s normalization contributes a determinant term that must not be omitted.

The Gaussian density is $$ \phi_q(\eta_i;0,\Omega)=(2\pi)^{-q/2}|\Omega|^{-1/2} \exp\left\{-\frac{1}{2}\eta_i^\mathsf{T}\Omega^{-1}\eta_i\right\}. $$ Thus, apart from constants, the log integrand is the individual log posterior kernel: the conditional log likelihood plus the Gaussian log-prior penalty. The normalized log integrand, including all constants, is $$ f_i(\eta_i;\theta)=\ell_i(\eta_i;\theta)-\frac{1}{2}\eta_i^\mathsf{T}\Omega^{-1}\eta_i-\frac{q}{2}\log(2\pi)-\frac{1}{2}\log|\Omega|. $$

Let \(\hat{\eta}_i\) be an interior mode of \(f_i\), equivalently of the individual log posterior. Define the local negative Hessian (curvature matrix) at that mode by $$ J_i=-\left.\frac{\partial^2 f_i(\eta_i;\theta)}{\partial\eta_i\,\partial\eta_i^\mathsf{T}}\right|_{\eta_i=\hat{\eta}_i}. $$ Because the Gaussian terms are quadratic, this can also be written $$ J_i=\Omega^{-1}-\left.\frac{\partial^2\ell_i(\eta_i;\theta)}{\partial\eta_i\,\partial\eta_i^\mathsf{T}}\right|_{\eta_i=\hat{\eta}_i}. $$ A valid ordinary Laplace calculation requires \(J_i\) to be positive definite, so the selected mode has locally downward curvature in every direction.

Deriving the Laplace Formula

Near an interior mode, the first derivative of \(f_i\) is zero. A second-order Taylor expansion therefore gives the local quadratic approximation $$ f_i(\eta_i;\theta)\approx f_i(\hat{\eta}_i;\theta)-\frac{1}{2}(\eta_i-\hat{\eta}_i)^\mathsf{T}J_i(\eta_i-\hat{\eta}_i). $$ The Laplace approximation substitutes this quadratic form for the log integrand in the integral. The resulting integral is a Gaussian integral: for positive-definite \(J_i\), its value over \(\mathbb{R}^q\) is \((2\pi)^{q/2}|J_i|^{-1/2}\). Hence $$ L_i(\theta)\approx \exp\{f_i(\hat{\eta}_i;\theta)\}(2\pi)^{q/2}|J_i|^{-1/2}. $$

Substituting the full expression for \(f_i\) shows exactly how the Gaussian normalization enters: $$ L_i(\theta)\approx \exp\left\{\ell_i(\hat{\eta}_i;\theta)-\frac{1}{2}\hat{\eta}_i^\mathsf{T}\Omega^{-1}\hat{\eta}_i\right\} |\Omega|^{-1/2}|J_i|^{-1/2}. $$ The factor \((2\pi)^{-q/2}\) in the normalized random-effects density cancels the \((2\pi)^{q/2}\) from integrating the local Gaussian approximation. The covariance determinant \(|\Omega|^{-1/2}\), however, remains. The log form, often more stable for calculation, is $$ \log L_i(\theta)\approx \ell_i(\hat{\eta}_i;\theta) -\frac{1}{2}\hat{\eta}_i^\mathsf{T}\Omega^{-1}\hat{\eta}_i -\frac{1}{2}\log|\Omega| -\frac{1}{2}\log|J_i|. $$

This derivation uses a local approximation to the integrand; it does not assert that the full log integrand is quadratic. The approximation is most persuasive when the integral is dominated by one smooth interior mode and the log integrand is well represented by its quadratic expansion over the region contributing appreciably to the integral. Positive definiteness is necessary for the displayed Gaussian integral, but by itself does not guarantee that the local approximation is accurate.

Normalization and Subject-Level Assembly

A frequent source of mistakes is mixing two conventions: using a normalized log integrand and then adding its normalization a second time, or using only the unnormalized kernel and forgetting the Gaussian density’s normalization. Either approach can be used, but the factors must be handled consistently. With the normalized integrand \(f_i\), use \(\exp\{f_i(\hat{\eta}_i)\}(2\pi)^{q/2}|J_i|^{-1/2}\). With the likelihood-kernel form, use the formula containing both \(|\Omega|^{-1/2}\) and \(|J_i|^{-1/2}\) above.

If subjects are conditionally independent given the population parameters and their random effects, and their random-effects vectors are mutually independent across subjects, so their joint density factors as \(\prod_{i=1}^{N}\phi_q(\eta_i;0,\Omega)\), integrating each subject’s random effects gives the marginal likelihood $$ L(\theta)=\prod_{i=1}^{N}L_i(\theta). $$ Using the Laplace approximation separately for each subject yields $$ \log L(\theta)\approx\sum_{i=1}^{N}\left[ \ell_i(\hat{\eta}_i;\theta) -\frac{1}{2}\hat{\eta}_i^\mathsf{T}\Omega^{-1}\hat{\eta}_i -\frac{1}{2}\log|\Omega| -\frac{1}{2}\log|J_i| \right]. $$ Each subject has its own mode and curvature matrix, while the same population covariance \(\Omega\) may apply across subjects. The terms must be evaluated at the parameter values under consideration; when \(\theta\) changes, the modes and Hessians generally change as well.

Conditions, Diagnostics, and Alternatives

Before interpreting or optimizing the approximation, check its mathematical and computational prerequisites. The random-effects covariance must be positive definite, the selected mode should be an interior stationary point, and \(J_i\) should be symmetric positive definite. A zero or negative eigenvalue means that the local quadratic form does not define a finite Gaussian integral in the stated way. A very small positive eigenvalue indicates weak curvature in a direction and can make the determinant and numerical calculation sensitive.

Useful checks include verifying that the gradient of \(f_i\) is close to zero at \(\hat{\eta}_i\), that the computed Hessian is symmetric to numerical tolerance, and that a Cholesky factorization of \(J_i\) succeeds. For a positive-definite matrix, the determinant can be evaluated from its Cholesky factor rather than by a potentially unstable direct determinant calculation. Repeating the mode optimization from different starting values can reveal optimization failures or different candidate modes. A mode that changes markedly with small changes in starting values, or a Hessian that is nearly singular, deserves further investigation.

The approximation is not an exact identity. It can be poor if important integrand mass lies far from the mode, if the local shape is strongly nonquadratic over the relevant region, or if there are multiple important modes. These are limitations of using a single local expansion, not errors that can be repaired merely by including the determinant term. When feasible, comparison with adaptive Gaussian quadrature or Monte Carlo integration can help assess the numerical approximation. FOCE is another commonly used approximation in NLME estimation, but it constructs an approximation through conditional linearization rather than through the local Hessian-volume calculation derived here; the two methods need not give identical objective values.

A common mistake is to use the Hessian of the conditional log likelihood alone as \(J_i\). The correct curvature is that of the full log integrand, including the Gaussian random-effects penalty: $$ J_i=\Omega^{-1}-\nabla_{\eta_i}^2\ell_i(\hat{\eta}_i;\theta). $$ Another mistake is to interpret a large \(|J_i|\) as a larger likelihood contribution by itself. In the formula, the curvature factor is \(|J_i|^{-1/2}\): greater local curvature means a narrower region around the mode and a smaller local volume, all else equal. The height, prior penalty, covariance normalization, and curvature must be considered together.