Laplace Formula for a GLMM Marginal Likelihood

From a Random-Effects Integral to a Likelihood Approximation

In a generalized linear mixed model, the response distribution is specified conditionally on random effects. The marginal likelihood for the observed responses therefore requires integrating those random effects out. Even when the conditional likelihood is straightforward to evaluate, that integral is often not available in closed form. The Laplace formula replaces it with a local approximation based on the integrand's value and curvature near its mode.

For a cluster with response vector $y$, let $b$ be a $q$-dimensional random-effects vector and let $\theta$ collect the model parameters, including any covariance parameters. Write the Gaussian random-effects density as $\phi_q(b;0,D(\theta))$. The marginal likelihood contribution for this cluster is

L(\theta)=\int_{\mathbb{R}^q} f(y\mid b,\theta)\,\phi_q(b;0,D(\theta))\,db =\int_{\mathbb{R}^q}\exp\{h(b;\theta)\}\,db,

where the log integrand is

h(b;\theta)=\log f(y\mid b,\theta)+\log\phi_q(b;0,D(\theta)).

The Gaussian density, including its normalizing constant, belongs in $h$. Thus $h$ is not just the conditional log likelihood. For independent clusters, the full marginal likelihood is the product of cluster contributions; correspondingly, their marginal log likelihood contributions add.

Definition: Let $\hat b$ be an interior mode of $h(b;\theta)$, with $\nabla_b h(\hat b;\theta)=0$. Define the negative Hessian at the mode by $H=-\nabla_b^2 h(\hat b;\theta)$. If $H$ is positive definite, the $q$-dimensional Laplace approximation is $$ L_{\mathrm{Lap}}(\theta)=\exp\{h(\hat b;\theta)\}(2\pi)^{q/2}|H|^{-1/2}. $$ Equivalently, its log is $$ \log L_{\mathrm{Lap}}(\theta)=h(\hat b;\theta)+\frac{q}{2}\log(2\pi)-\frac{1}{2}\log|H|. $$

Why the Three Terms Appear

The formula combines three distinct quantities: the log integrand at the mode, the local curvature represented by the Hessian determinant, and a constant that depends on the dimension $q$. Keeping them separate makes both calculations and interpretations easier to check.

The quadratic Taylor approximation around the mode is

h(b;\theta)\approx h(\hat b;\theta)-\frac{1}{2}(b-\hat b)^\mathsf{T}H(b-\hat b).

Exponentiating and integrating this approximation gives

\begin{aligned} L(\theta) &\approx \exp\{h(\hat b;\theta)\} \int_{\mathbb{R}^q}\exp\left\{-\frac{1}{2}(b-\hat b)^\mathsf{T}H(b-\hat b)\right\}\,db\\ &=\exp\{h(\hat b;\theta)\}(2\pi)^{q/2}|H|^{-1/2}. \end{aligned}

The final equality is the multivariate Gaussian integral for positive-definite $H$. A larger determinant means stronger overall curvature around the mode and, all else equal, a smaller volume under this local Gaussian approximation. The exponent $q/2$ is set by the number of random-effect dimensions being integrated, not by the number of observations.

Worked Example: Keeping the Log-Scale Terms Separate

Suppose a numerical mode search for a hypothetical binary-response cluster at fixed $\theta$ gives $\hat b=0.4$, $h(\hat b;\theta)=-2.0$, and the one-dimensional negative Hessian $H=4$. These are stipulated outputs of the mode and curvature calculations; $H>0$, so the mode has positive local curvature in the required sense.

Here $q=1$, so

\log L_{\mathrm{Lap}} =-2.0+\frac{1}{2}\log(2\pi)-\frac{1}{2}\log(4).

Using $\tfrac12\log(2\pi)\approx0.918939$ and $\tfrac12\log(4)\approx0.693147$,

\log L_{\mathrm{Lap}}\approx-2.0+0.918939-0.693147=-1.77421.

Therefore $L_{\mathrm{Lap}}\approx\exp(-1.77421)=0.1696$. The same result follows on the likelihood scale: $e^{-2.0}\sqrt{2\pi}/\sqrt{4}\approx0.1696$. The value $-2.0$ is the log integrand at the mode, not the approximate log marginal likelihood; the curvature and dimensional terms change it.

A Binary Random-Intercept Calculation

For a binary random-intercept model, let one hypothetical cluster contain a single observed success, $y=1$, with conditional success probability $p(b)=\{1+\exp(-b)\}^{-1}$. Suppose $b\sim N(0,1)$, with model parameters held fixed. The marginal likelihood contribution is the integral of the Bernoulli probability times the standard normal density. Its log integrand is

h(b)=\log p(b)-\frac{b^2}{2}-\frac{1}{2}\log(2\pi).

The derivative is $h'(b)=1-p(b)-b$. The mode solves $b=1-p(b)$, and numerical solution gives $\hat b\approx0.4011$, with $p(\hat b)\approx0.5989$. The second derivative is $h''(b)=-p(b)\{1-p(b)\}-1$, so $H=1+p(\hat b)\{1-p(\hat b)\}\approx1.2402>0$. This verifies positive local curvature at the mode.

Worked Example: Approximating a Binary Cluster Contribution

Using the hypothetical data and model above, substitute the mode into the log integrand:

h(\hat b)\approx\log(0.5989)-\frac{(0.4011)^2}{2}-\frac{1}{2}\log(2\pi) \approx-0.5126606397-0.0804406050-0.9189385332=-1.5120397779\approx-1.5120.

There is one random effect, so the log approximation is

\log L_{\mathrm{Lap}} \approx-1.5120+\frac{1}{2}\log(2\pi)-\frac{1}{2}\log(1.2402) \approx-0.7007.

Exponentiating gives $L_{\mathrm{Lap}}\approx0.496$, a dimensionless probability contribution. In this particular example the exact marginal probability is $0.5$: symmetry of the standard normal distribution and the identity $p(b)+p(-b)=1$ imply that the average success probability is one half. Thus the approximation is close here, though this single example does not establish accuracy in other models or data sets.

More Than One Random-Effects Dimension

When $q>1$, curvature cannot generally be summarized by a single second derivative. The determinant $|H|$ combines curvature across directions, including dependencies between them. The determinant must be calculated from the full Hessian matrix, not by multiplying its diagonal entries unless the matrix is diagonal.

Worked Example: A Two-Dimensional Hessian

Suppose a hypothetical cluster has two random effects and, at the computed mode, $h(\hat b;\theta)=-4.3$ and

H= \begin{pmatrix} 5&1\\ 1&3 \end{pmatrix}.

The matrix is positive definite: its first leading principal minor is $5>0$, and its determinant is $5(3)-1(1)=14>0$. Thus $|H|=14$, $q=2$, and

\log L_{\mathrm{Lap}} =-4.3+\frac{2}{2}\log(2\pi)-\frac{1}{2}\log(14) \approx-4.3+1.83788-1.31953=-3.78165.

Consequently $L_{\mathrm{Lap}}\approx e^{-3.78165}=0.02278$. Using only the diagonal entries would give $5(3)=15$, not the correct determinant $14$, and would therefore produce a different approximation.

Checks, Alternatives, and Limitations

A Laplace calculation is most interpretable when its inputs and assumptions are checked explicitly. At fixed $\theta$, the log integrand should be sufficiently smooth near an interior mode, the mode should be located accurately, and $H$ should be positive definite. The approximation also relies on the integral being well represented by the local quadratic shape around that mode. Positive curvature at one mode does not prove that it is the only important mode or that the quadratic shape represents the integrand well over the region contributing most of the integral.

Useful checks include verifying that the gradient is close to zero at $\hat b$, that the Hessian is symmetric and positive definite, and that the computed log determinant is finite. A poor mode search can distort both $h(\hat b;\theta)$ and $H$. A nearly flat direction can make the determinant small and the curvature correction large; this can signal weak local identification or numerical instability. If the log integrand is strongly skewed, has multiple relevant modes, or has substantial non-quadratic shape, a single-mode approximation may be inaccurate even when the local Hessian is positive definite.

The principal alternative is to evaluate the random-effects integral by numerical integration, such as quadrature, when feasible. Direct integration avoids relying solely on a local quadratic approximation, but its computational cost can grow quickly with dimension. Laplace approximation is often computationally convenient, particularly when mode and Hessian calculations are already available, but it is still an approximation. It does not become exact merely because the Hessian has been calculated correctly.

A frequent mistake is to put only $\log f(y\mid b,\theta)$ in $h$, then use the formula as if the random-effects density were included. The correction is to define the full integrand first: $h=\log f(y\mid b,\theta)+\log\phi_q(b;0,D(\theta))$. Another common error is to omit $q/2$ or use the wrong determinant: the coefficient is the number of integrated random-effect dimensions divided by two, and the determinant is that of the negative Hessian of the full log integrand at the mode. On the log scale these contributions are $+(q/2)\log(2\pi)$ and $-\tfrac12\log|H|$, respectively.

Key takeaway: For a GLMM random-effects integral, evaluate the full log integrand at its mode, calculate the negative Hessian there, and combine them as $h(\hat b;\theta)+(q/2)\log(2\pi)-(1/2)\log|H|$. Keep all three terms visible and verify the mode and curvature before interpreting the result.

Check Your Understanding

Answer each question, then compare with the solutions below.

For a three-dimensional random-effects integral, what is the coefficient of $\log(2\pi)$ in the log Laplace formula?
Why must the Gaussian random-effects density be included in the log integrand?
If $h(\hat b)=-1.5$, $q=1$, and $H=2$, calculate $\log L_{\mathrm{Lap}}$ to three decimal places.
What does $H$ need to satisfy for the stated formula to apply at an interior mode?
Why is a positive-definite Hessian not, by itself, a guarantee of an accurate approximation?

Solutions. (1) The coefficient is $q/2=3/2$. (2) The marginal likelihood integrates the conditional response likelihood multiplied by the random-effects density, so omitting that density changes both the integrand and its mode and curvature. (3) $\log L_{\mathrm{Lap}}=-1.5+\tfrac12\log(2\pi)-\tfrac12\log(2)\approx-0.928$. (4) $H$, the negative Hessian of the full log integrand at the mode, must be positive definite. (5) It verifies local curvature at that mode, but does not ensure that the integrand is well approximated by a quadratic over the important integration region or that other modes are negligible.