Tutorials › Joint Modeling of Longitudinal and Time-to-Event Data › Deriving the Laplace Approximation

Home  /  Tutorials  /  Maximum Likelihood and EM Estimation

Maximum Likelihood and EM Estimation · Tutorial 262 of 697

Deriving the Laplace Approximation

Learn to derive the Laplace formula for a joint-model marginal likelihood and interpret the determinant of the negative Hessian at the integrand’s mode.

Advanced 10 min read

What You'll Learn

  • Express a subject’s marginal likelihood as an integral of an exponentiated log integrand.
  • Derive the local quadratic approximation around an interior mode.
  • Explain why the negative Hessian must be positive definite for the nondegenerate formula.
  • Obtain the determinant factor by transforming the quadratic form to a standard Gaussian integral.
  • Work through one-dimensional and multivariate hypothetical examples.
  • Recognize when curvature, multiple modes, or tails may make the approximation unreliable.

Why Curvature Appears in a Marginal Likelihood

In a joint model, a subject’s longitudinal and event-time observations are linked through latent random effects. The subject-specific marginal likelihood is obtained by integrating those effects out. Even when the conditional likelihood and random-effects density are straightforward to evaluate, this integral may not have a closed form. The Laplace approximation replaces the log integrand near its dominant mode by a quadratic function, whose exponential has a Gaussian integral.

The determinant in the resulting formula is not an optional adjustment. It measures the local volume of the integrand around its mode: a sharply curved peak contributes less integrated mass than a flatter peak with the same height. The derivation below shows exactly where this term comes from and what conditions it requires.

Definition: For subject \(i\), let \(b_i\in\mathbb{R}^q\) denote the \(q\)-dimensional random-effects vector, and write the marginal likelihood contribution as \(L_i(\theta)=\int_{\mathbb{R}^q}\exp\{h_i(b_i;\theta)\}\,db_i\). Here \(h_i\) is the log of the full integrand, including both the conditional contribution of the observed data and the random-effects density. Let \(\hat b_i\) be an interior mode of \(h_i\), and define \(H_i=-\nabla^2 h_i(\hat b_i;\theta)\), the negative Hessian evaluated at that mode.

The integrand must be defined with respect to the same coordinates and measure used in the integral. In particular, constants in the conditional data contribution or random-effects density belong in \(h_i\) if they depend on \(b_i\), and constants independent of \(b_i\) must still be retained when computing the marginal likelihood’s numerical value. A change of coordinates also changes the integration measure and can contribute a Jacobian.

Deriving the Formula

At an interior differentiable mode, the gradient is zero: \(\nabla h_i(\hat b_i;\theta)=0\). A second-order Taylor expansion around \(\hat b_i\), with \(\delta=b_i-\hat b_i\), gives

$$ h_i(\hat b_i+\delta;\theta) = h_i(\hat b_i;\theta) +\frac{1}{2}\delta^\mathsf{T}\nabla^2h_i(\hat b_i;\theta)\delta +o(\|\delta\|^2). $$

For the nondegenerate Laplace formula, assume that the Hessian at the mode is negative definite. Then \(H_i=-\nabla^2h_i(\hat b_i;\theta)\) is positive definite, and the local expansion becomes

$$ h_i(\hat b_i+\delta;\theta) \approx h_i(\hat b_i;\theta)-\frac{1}{2}\delta^\mathsf{T}H_i\delta. $$

Exponentiating this approximation separates the height of the integrand at its mode from its local shape. If the integral is dominated by this mode, replacing the log integrand by its quadratic expansion yields

$$ L_i(\theta) \approx \exp\{h_i(\hat b_i;\theta)\} \int_{\mathbb{R}^q} \exp\left\{-\frac{1}{2}\delta^\mathsf{T}H_i\delta\right\}\,d\delta. $$

To evaluate the remaining integral, use the Cholesky factorization \(H_i=R^\mathsf{T}R\), where \(R\) is nonsingular. Set \(u=R\delta\). Then \(\delta=R^{-1}u\), and the absolute Jacobian determinant is \(|\det(R^{-1})|=|H_i|^{-1/2}\), because \(|H_i|=|R|^2\). Also \(\delta^\mathsf{T}H_i\delta=u^\mathsf{T}u\). Therefore

$$ \int_{\mathbb{R}^q}\exp\left\{-\frac{1}{2}\delta^\mathsf{T}H_i\delta\right\}\,d\delta = |H_i|^{-1/2} \int_{\mathbb{R}^q}\exp\left\{-\frac{1}{2}u^\mathsf{T}u\right\}\,du = \frac{(2\pi)^{q/2}}{\sqrt{|H_i|}}. $$

The last equality follows by multiplying the \(q\) one-dimensional standard Gaussian integrals, each equal to \(\sqrt{2\pi}\). Combining the pieces gives the Laplace approximation:

Result: If the integral is dominated by an interior mode \(\hat b_i\), and the negative Hessian \(H_i\) there is positive definite, then the nondegenerate Laplace approximation is \(L_i(\theta)\approx \exp\{h_i(\hat b_i;\theta)\}(2\pi)^{q/2}/\sqrt{|H_i|}\). The corresponding approximate log contribution is \(h_i(\hat b_i;\theta)+(q/2)\log(2\pi)-(1/2)\log|H_i|\).

An interior local maximum has a negative semidefinite Hessian, but that alone is not enough for this formula. A strict maximum can still have zero curvature in a direction; for example, \(h(b)=-b^4\) has a strict maximum at zero and a zero second derivative there. Positive definiteness of \(H_i\) is the additional nondegeneracy condition that makes the quadratic Gaussian integral finite and supplies the stated determinant factor.

Worked Examples

Worked Example: A One-Dimensional Quadratic Log Integrand

Hypothetical data. Suppose the full log integrand for one subject’s random intercept, denoted \(b\), is \(h(b)=-2-(b-1)^2\). Find the Laplace approximation to the marginal likelihood contribution.

The derivative is \(h'(b)=-2(b-1)\), so the mode is \(\hat b=1\). The second derivative is \(h''(b)=-2\), hence \(H=-h''(1)=2>0\). The one-dimensional formula gives

$$ L\approx \exp\{h(1)\}\frac{\sqrt{2\pi}}{\sqrt{H}} =e^{-2}\frac{\sqrt{2\pi}}{\sqrt{2}} =e^{-2}\sqrt{\pi}. $$

Numerically, \(e^{-2}\approx0.135335\) and \(\sqrt{\pi}\approx1.77245\), so \(L\approx0.135335(1.77245)\approx0.23988\). Because this log integrand is exactly quadratic, the Gaussian calculation is exact here; for a nonquadratic log integrand, it is generally an approximation.

Worked Example: A Nonquadratic Log Integrand

Hypothetical data. Suppose a subject’s log integrand for a scalar random effect is \(h(b)=-0.7-\frac{1}{2}(b-1)^2-\frac{1}{4}(b-1)^4\). Apply the Laplace approximation and identify what feature it does not capture.

Let \(x=b-1\). Then \(h'(b)=-x-x^3=-x(1+x^2)\), which is zero only at \(x=0\), or \(b=1\). Since \(h(b)\) decreases without bound as \(|b|\) increases and the quartic term is nonpositive, this is the global mode. At that mode \(h(1)=-0.7\), and \(h''(b)=-1-3(b-1)^2\), so \(H=-h''(1)=1>0\).

$$ L\approx e^{-0.7}\frac{\sqrt{2\pi}}{\sqrt{1}} \approx 0.496585(2.50663) \approx 1.2448. $$

This calculation satisfies the interior-mode and positive-curvature conditions. Nevertheless, it replaces the log integrand locally by \(-0.7-\frac{1}{2}(b-1)^2\) and omits the quartic term. The approximation therefore uses the peak height and curvature but does not fully describe the integrand’s shape away from the mode; its accuracy is not established by this calculation alone.

Worked Example: The Hessian Determinant in Two Dimensions

Hypothetical data. Suppose the log integrand for a subject’s two random effects \(b=(b_1,b_2)^\mathsf{T}\) is \(h(b)=-3-\frac{1}{2}(b-m)^\mathsf{T}A(b-m)\), where \(m=(0.2,-0.1)^\mathsf{T}\) and \(A=\begin{pmatrix}4&1\\1&2\end{pmatrix}\). Calculate the Laplace contribution.

The mode is \(\hat b=m\), with \(h(m)=-3\), and \(H=A\). This Hessian is positive definite: its leading diagonal entry is \(4>0\), and its determinant is \(4(2)-1(1)=7>0\). Thus the required curvature condition holds. Since \(q=2\), the approximation is

$$ L\approx e^{-3}\frac{(2\pi)^{2/2}}{\sqrt{7}} =e^{-3}\frac{2\pi}{\sqrt{7}} \approx 0.0497871(2.37482) \approx 0.11824. $$

The determinant term is \(1/\sqrt{7}\), not the reciprocal of either diagonal entry separately. It accounts for the joint curvature in both directions, including the off-diagonal relationship between the random effects. Since this log integrand is exactly quadratic, the approximation is again exact.

Conditions, Diagnostics, and Alternatives

The derivation uses a smooth log integrand near the mode, an interior stationary point, and a positive-definite negative Hessian there. It also relies on the selected mode representing the dominant part of the integral. The Taylor expansion justifies a local approximation; it does not by itself show that regions far from the mode contribute negligibly. A common justification is asymptotic concentration: if \(h(b)=n f(b)\), with a unique dominant mode and a positive-definite negative Hessian for \(f\), the peak becomes increasingly concentrated as \(n\) grows. The curvature of \(h\) is then \(n\) times the curvature of \(f\), and the determinant contributes a factor proportional to \(n^{-q/2}\).

In practice, check that the numerical optimizer has found a stationary point, that the Hessian is symmetric to numerical precision, and that its eigenvalues are all positive after changing sign to form \(H_i\). Inspect the log integrand around the mode for skewness, flat directions, heavy tails, or another mode of similar height. These features may signal that the local quadratic shape is a poor summary of the full integral. If an eigenvalue of \(H_i\) is close to zero, the determinant is sensitive to small numerical errors and the approximation may be unstable.

Direct numerical integration, including quadrature, does not replace the integrand by a local quadratic expansion and can be useful for checking a Laplace result when the random-effects dimension is manageable. Its computational cost can grow quickly with dimension, while Laplace is often attractive for repeated subject-level integrals. The trade-off is that a computationally convenient approximation can be inaccurate when the integrand is skewed, weakly identified, or multimodal. Neither method removes the need to define the full integrand correctly.

Common pitfall: A strict local maximum does not necessarily imply a negative-definite Hessian. At an interior local maximum the Hessian is negative semidefinite; the nondegenerate Laplace formula additionally requires it to be negative definite. Also, using the Hessian with the wrong sign makes a maximum look like negative curvature in the Gaussian exponent. Define \(H_i=-\nabla^2h_i(\hat b_i;\theta)\), verify it is positive definite, and use \(\sqrt{|H_i|}\) in the denominator.

Check Your Understanding

Try each question, then compare your reasoning with the solutions below.

  1. For \(h(b)=c-\frac{a}{2}(b-m)^2\), where \(a>0\), identify the mode and negative Hessian, and write the Laplace integral.
  2. Why does a two-dimensional approximation involve \(\sqrt{|H|}\), rather than the product of the two diagonal entries of \(H\)?
  3. At an interior strict maximum, what additional condition is needed for the nondegenerate Laplace formula?
  4. Give two features of a log integrand that can make a local quadratic approximation unreliable, even if the mode and Hessian have been computed accurately.

Solutions. (1) The mode is \(m\), \(H=a\), and \(L\approx e^c\sqrt{2\pi/a}\). (2) The full determinant accounts for curvature in all directions and their dependence; off-diagonal entries affect the volume of the quadratic peak. (3) The Hessian at the mode must be negative definite, so \(H\) is positive definite. Strictness alone is insufficient. (4) Examples include substantial skewness, a nearly flat direction, heavy tails, or a second mode with comparable height; each can make the integral’s important mass poorly represented by the quadratic shape at one mode.

Previous

Next

Back to Tutorials
Clinical Biostats