Determinant and Normalization Terms in Laplace Calculations

Why the “constant” terms matter

A Laplace calculation combines the height of an integrand at its conditional mode with the local volume of the integrand around that mode. The height is recorded on a log scale, while the volume contributes both a power of $2\pi$ and a Hessian determinant. In a GLMM, the integrand also includes the density of the random effects, whose Gaussian normalization may depend on covariance parameters. These pieces are not interchangeable: dropping one changes the approximate marginal likelihood, and dropping a parameter-dependent term can also change its estimates.

The focus here is bookkeeping and interpretation of those terms. The starting point is the log integrand $h(b)$, where $b$ is a $q$-dimensional random-effects vector. Let $\hat b$ be an interior mode and let $H=-\nabla^2h(\hat b)$ be the negative Hessian there. The usual Laplace calculation requires $h$ to be sufficiently smooth near the mode and $H$ to be positive definite, so the mode has locally downward curvature in every direction.

Definition: Under these conditions, the $q$-dimensional Laplace approximation to $I=\int \exp\{h(b)\}\,db$ is $$ I\ \approx\ \exp\{h(\hat b)\}(2\pi)^{q/2}|H|^{-1/2}, \qquad \log I\ \approx\ h(\hat b)+\frac{q}{2}\log(2\pi)-\frac{1}{2}\log|H|. $$ Here $|H|$ is the determinant of the negative Hessian at the mode. The logarithm of the determinant is the curvature correction; it is not an optional adjustment.

Where the terms come from

A second-order Taylor expansion around $\hat b$ has no linear term because $\nabla h(\hat b)=0$. With $u=b-\hat b$, it gives the local approximation $h(b)\approx h(\hat b)-\tfrac12u^\mathsf{T}Hu$. Substituting this quadratic expression into the integral leaves a multivariate Gaussian integral. To evaluate it, diagonalize the positive-definite matrix $H$: its eigenvalues $\lambda_1,\ldots,\lambda_q$ are all positive. In the corresponding orthogonal coordinates the integral separates into one-dimensional Gaussian integrals, each contributing $\sqrt{2\pi/\lambda_j}$. Multiplying these terms gives the determinant and normalization factors:

\int \exp\!\left\{h(\hat b)-\frac12u^\mathsf{T}Hu\right\}\,du = \exp\{h(\hat b)\} \prod_{j=1}^{q}\sqrt{\frac{2\pi}{\lambda_j}} = \exp\{h(\hat b)\}(2\pi)^{q/2}|H|^{-1/2}.

The factor $(2\pi)^{q/2}$ comes from integrating over $q$ directions. The determinant is the product of the curvatures, $|H|=\prod_j\lambda_j$. Larger curvature makes the integrand narrower, so it reduces the approximated integral. A determinant that is zero or nonpositive is not a valid curvature correction for this interior-mode formula: it signals a singular or non-positive-definite $H$, and the calculation needs investigation rather than an absolute value used to conceal the problem.

Normalization in a GLMM integrand

For a Gaussian random-effects vector $b\sim N(0,D)$, where $D$ is positive definite, its log density is

\log\phi(b;0,D) = -\frac{q}{2}\log(2\pi) -\frac{1}{2}\log|D| -\frac{1}{2}b^\mathsf{T}D^{-1}b.

If the conditional log likelihood is $\ell(b)$, the log integrand for the marginal likelihood is $h(b)=\ell(b)+\log\phi(b;0,D)$. Thus $h(\hat b)$ includes the Gaussian density’s normalization as well as the conditional likelihood and the quadratic penalty. The Laplace factor then contributes a separate $\tfrac{q}{2}\log(2\pi)$ and $-\tfrac12\log|H|$. The two appearances of $2\pi$ have different origins; one must not cancel or omit them by inspection.

Likelihood constants also deserve attention. For example, a Poisson log likelihood includes $-\log(y!)$, and a binomial log likelihood includes its combinatorial term. Such terms may be omitted for some parameter estimates if they do not depend on the parameters being optimized. They still affect an absolute likelihood. In contrast, $-\tfrac12\log|D|$ depends on covariance parameters when $D$ does, so it generally cannot be dropped while estimating those parameters.

Worked Example: A Poisson Random Intercept and Its Normalizer

Suppose a hypothetical cluster has two Poisson observations $y_1=y_2=1$, each with unit exposure and conditional mean $\exp(b)$. Let its random intercept have distribution $b\sim N(0,4)$. The conditional log likelihood is $\ell(b)=2b-2\exp(b)$, since both factorial terms are $\log(1!)=0$.

The log integrand is $h(b)=\ell(b)+\log\phi(b;0,4)$. At $b=0$, the likelihood derivative is $2-2\exp(0)=0$, and the prior derivative is also zero, so $\hat b=0$ is a stationary point. The second derivative of the likelihood there is $-2$; that of the log prior is $-1/4$. Therefore $H=2+1/4=2.25>0$, verifying positive curvature for this one-dimensional calculation.

The log likelihood at the mode is $\ell(0)=-2$. The Gaussian log density there, including its normalizer, is

\log\phi(0;0,4) = -\frac12\log(2\pi)-\frac12\log(4)-\frac12\left(\frac{0^2}{4}\right) = -0.91894-0.69315+0 = -1.61209.

Thus $h(0)=-2-1.61209=-3.61209$. With $q=1$ and $H=2.25$, the approximate log marginal likelihood is

\log I\approx -3.61209+\frac12\log(2\pi)-\frac12\log(2.25) =-3.61209+0.91894-0.40547 =-3.09862.

The corresponding approximate integral is $\exp(-3.09862)\approx0.0451$. If the prior normalizer is omitted, the approximate log integral becomes $-2+0.91894-0.40547=-1.48653$, which is too large by $1.61209$, exactly the omitted negative log-density contribution. The mode and curvature in this example do not change when that constant is omitted, but the marginal likelihood does.

Reading the determinant in more than one dimension

In several dimensions the determinant combines curvature across directions. It is not generally adequate to look at individual diagonal entries of $H$, especially when the Hessian has off-diagonal terms. A useful computational check is a Cholesky factorization $H=RR^\mathsf{T}$, with positive diagonal entries in $R$. Then $\log|H|=2\sum_j\log R_{jj}$. This calculates the log determinant without first forming a potentially unstable determinant.

Worked Example: A Two-Dimensional Curvature Correction

Consider a hypothetical GLMM fit with two random effects. Suppose the evaluated log integrand at its reported mode is $h(\hat b)=-5.000$, and the negative Hessian is

H= \begin{pmatrix} 2&0\\ 0&8 \end{pmatrix}.

The matrix is positive definite because its two eigenvalues are $2$ and $8$, both positive. Its determinant is $|H|=2(8)=16$, and $q=2$. Therefore,

\log I\approx -5.000+\frac{2}{2}\log(2\pi)-\frac12\log(16) =-5.000+1.83788-1.38629 =-4.54842.

The determinant correction is $-1.38629$: the integrand is relatively narrow in the combined two-dimensional volume. If the determinant term is simply omitted, the result is $-5.000+\log(2\pi)=-3.16212$, too high by $1.38629$. Using only one curvature value would likewise fail to account for both directions. These supplied local quantities are enough to check the arithmetic; in an applied fit, $h(\hat b)$ and $H$ must be evaluated from the full log integrand, with its constants retained.

Changing coordinates: the Jacobian belongs in the integral

The determinant correction depends on the coordinates used to express $b$. That does not mean the approximated integral should depend on a mere linear re-expression of the same random effects. If $b=Lz$, where $L$ is invertible, then the volume elements satisfy $db=|L|\,dz$. The transformed log integrand for integration with respect to $dz$ is $h(Lz)+\log|L|$. At a stationary point, its negative Hessian is $H_z=L^\mathsf{T}HL$, so $|H_z|=|L|^2|H|$. The Jacobian and determinant consequently offset one another in the Laplace calculation:

h(\hat b)+\log|L|+\frac{q}{2}\log(2\pi)-\frac12\log|H_z| = h(\hat b)+\frac{q}{2}\log(2\pi)-\frac12\log|H|.

This identity requires the Jacobian and the transformed Hessian to be used together. It is not permission to omit the Jacobian when changing variables. The result also presumes a linear invertible transformation and an interior stationary point; more general transformations require careful treatment of their local Jacobian and curvature.

Worked Example: Rescaling a Random Intercept

Use the one-dimensional Poisson calculation above, where $\hat b=0$, $h(\hat b)=-3.61209$, and $H_b=2.25$. Define $b=2z$. The mode in $z$ is also zero, and the transformed negative Hessian is $H_z=2^2H_b=9$. The Jacobian is $|db/dz|=2$, so the transformed log integrand’s mode height is $-3.61209+\log 2$.

The transformed calculation gives

-3.61209+\log 2+\frac12\log(2\pi)-\frac12\log(9) = -3.61209+0.69315+0.91894-1.09861 = -3.09861.

This agrees with the calculation in $b$-coordinates. If the Jacobian were omitted, the result would be lower by $\log 2=0.69315$. The Hessian determinant changed under rescaling, so comparing determinant corrections alone would not establish invariance; the volume factor must be included as well.

Checks, alternatives, and limitations

A reliable calculation keeps a short ledger: the full value of $h(\hat b)$, the dimension $q$, and $\log|H|$, followed by the two Laplace terms. Check that the Hessian is the negative Hessian of the entire log integrand, including the random-effects density, and that it is positive definite. Confirm the coordinate system and integration measure used to obtain it. When $D$ depends on covariance parameters, also check that its log determinant remains in the integrand.

The determinant is sensitive to scale: rescaling random-effect coordinates changes its numerical value. A mismatch in units or a missing Jacobian can therefore produce an apparently plausible but incorrect approximation. A non-positive Hessian, an unstable factorization, or a large discrepancy between implementations using equivalent coordinates calls for investigation of mode convergence, derivative calculations, and parameterization. The formula assumes a suitable isolated interior mode with locally quadratic behavior; a local determinant cannot by itself certify that the integrand is well represented over the region contributing to the integral.

Direct numerical integration and quadrature are alternatives for evaluating the random-effects integral. They can provide useful comparisons when the integral is tractable, while Laplace is attractive when a local approximation is computationally convenient. Neither comparison removes the need to define the same integrand and measure, including all relevant constants. This tutorial’s formula also does not guarantee that the approximation is accurate merely because its determinant is calculated correctly.

Do not drop Gaussian normalizers automatically. Their dependence on covariance parameters can affect estimation, and they always affect an absolute marginal likelihood.
Do not substitute a determinant of the wrong matrix. The formula uses the negative Hessian of the complete log integrand at the mode.
Do not change coordinates without changing the measure. Include the Jacobian and recompute curvature in the new coordinates.
Do not force a non-positive determinant into the formula. Check mode finding and curvature assumptions rather than masking a failure.

Key takeaway: A GLMM Laplace log likelihood is the sum of the complete log-integrand height, the Gaussian integration normalization, and the negative half log determinant of the negative Hessian. In transformed coordinates, include the Jacobian as well.

Check Your Understanding

Answer each question, then compare your work with the solutions below.

For $q=1$, $h(\hat b)=-2.4$, and $H=4$, calculate the Laplace approximation to $\log I$.
A two-dimensional negative Hessian has eigenvalues $3$ and $5$. What is its determinant correction, $-\tfrac12\log|H|$?
If a Gaussian random-effects covariance matrix depends on an estimated variance parameter, may its log-determinant normalizer be discarded when estimating that parameter? Explain.
Under $b=3z$ in one dimension, what Jacobian term must be added to the transformed log integrand, and what is the transformed negative Hessian if $H_b=2$?

Solutions. (1) $-2.4+\tfrac12\log(2\pi)-\tfrac12\log4=-2.4+0.91894-0.69315=-2.17421$. (2) The determinant is $3(5)=15$, so the correction is $-\tfrac12\log15=-1.35403$. (3) No. The Gaussian log density includes $-\tfrac12\log|D|$, which varies with the covariance parameter and can change the likelihood being optimized. (4) Add $\log|3|=\log3=1.09861$; the transformed negative Hessian is $H_z=3^2H_b=18$.