A One-Dimensional Laplace Calculation

From a Random-Effect Integral to a Number

In a nonlinear mixed-effects model, a subject’s likelihood contribution is obtained by integrating over the subject-specific random effect. Even with only one random effect, this integral may not have a closed-form solution. A one-dimensional Laplace calculation replaces the log integrand near its conditional mode by a quadratic function, then integrates the corresponding Gaussian-shaped approximation. This tutorial works through the calculation itself, from model statement to numerical result.

Consider one subject with scalar random effect $\eta$ and observations $y_1,\ldots,y_m$. Suppose that, conditional on $\eta$, the observations are independent normal variables with means $f_j(\eta)$ and variances $\sigma_j^2$, and that $\eta$ has a mean-zero normal distribution with variance $\omega^2$. The parameters in $f_j$, $\sigma_j^2$, and $\omega^2$ are held fixed during this calculation.

Definition: The scalar marginal likelihood contribution is $$ L=\int_{-\infty}^{\infty} \exp\{h(\eta)\}\,d\eta, $$ where $h(\eta)$ is the log of the conditional observation density multiplied by the random-effect density: $$ h(\eta)= -\frac{1}{2}\sum_{j=1}^{m}\left[ \log(2\pi\sigma_j^2)+\frac{\{y_j-f_j(\eta)\}^2}{\sigma_j^2} \right] -\frac{1}{2}\left[ \log(2\pi\omega^2)+\frac{\eta^2}{\omega^2} \right]. $$ If $\hat{\eta}$ is an interior mode of $h$, with $h'(\hat{\eta})=0$ and $h''(\hat{\eta})<0$, define $J=-h''(\hat{\eta})>0$. The one-dimensional Laplace approximation is $$ L\approx \exp\{h(\hat{\eta})\}\sqrt{\frac{2\pi}{J}}. $$

The density normalizing constants belong in $h$ when the goal is the marginal likelihood, not just a likelihood proportional to it. Omitting a constant that does not depend on $\eta$ leaves the mode and $J$ unchanged, but changes $\exp\{h(\hat{\eta})\}$, and therefore changes the approximated integral.

Why the Curvature Formula Appears

The calculation uses a local approximation at a mode. A second-order Taylor expansion of $h$ around $\hat{\eta}$ is $$ h(\eta)\approx h(\hat{\eta})+h'(\hat{\eta})(\eta-\hat{\eta}) +\frac{1}{2}h''(\hat{\eta})(\eta-\hat{\eta})^2. $$ At an interior mode, $h'(\hat{\eta})=0$. Since $h''(\hat{\eta})=-J$, this reduces to $$ h(\eta)\approx h(\hat{\eta})-\frac{J}{2}(\eta-\hat{\eta})^2. $$ Exponentiating and integrating the resulting Gaussian-shaped function gives $$ \int_{-\infty}^{\infty} \exp\left\{h(\hat{\eta})-\frac{J}{2}(\eta-\hat{\eta})^2\right\}\,d\eta = \exp\{h(\hat{\eta})\}\sqrt{\frac{2\pi}{J}}. $$

This derivation requires a finite, interior mode and positive curvature $J$. It also assumes that the local quadratic approximation is useful over the part of the integrand that contributes substantially to the integral. The formula is an approximation unless $h$ is exactly quadratic.

A Reproducible Calculation

Worked Example: A Gaussian Benchmark With All Constants Included

Suppose a hypothetical subject has $\eta\sim N(0,0.25)$ and one observation $y\mid\eta\sim N(\eta,0.25)$. Let $y=\sqrt{0.21765}\approx0.4665$, in standardized response units. The marginal contribution is the integral of the observation density times the random-effect density.

1. Form the log integrand. The observation variance and random-effect variance are both $0.25$, so $$ h(\eta)=-\log(2\pi)-\frac{1}{2}\left[ \frac{(y-\eta)^2}{0.25}+\frac{\eta^2}{0.25} \right]. $$ The two normalizing constants together are $1/(2\pi\cdot0.25)=2/\pi$.

2. Find the mode. Differentiating the $\eta$-dependent terms gives $$ h'(\eta)=\frac{y-\eta}{0.25}-\frac{\eta}{0.25}. $$ Setting this equal to zero gives $y-2\eta=0$, so $\hat{\eta}=y/2\approx0.2333$. The second derivative is $h''(\eta)=-8$, which is negative everywhere; the stationary point is the unique maximum.

3. Evaluate the height and curvature. At the mode, both squared residual terms in the bracket equal $y^2/4$. Thus the sum inside the bracket is $2y^2$, and $$ h(\hat{\eta})=-\log(2\pi)-y^2=-\log(2\pi)-0.21765. $$ Therefore, $\exp\{h(\hat{\eta})\}=(2/\pi)\exp(-0.21765)\approx0.5121$. The negative second derivative is $J=8$.

4. Calculate the approximation. Substitution gives $$ L\approx0.5121\sqrt{\frac{2\pi}{8}} =0.5121(0.8862)\approx0.4538. $$ The marginal response density is approximately $0.4538$ per standardized response unit. As a check, the exact marginal distribution is $y\sim N(0,0.25+0.25)$; its density at $y=\sqrt{0.21765}$ is $\exp(-0.21765)/\sqrt{\pi}\approx0.4538$. This agreement is expected because the log integrand is exactly quadratic in this Gaussian benchmark.

Worked Example: Mode, Curvature, and Integral for a Given Log Integrand

For a second hypothetical scalar calculation, suppose the fully normalized log integrand, including its likelihood and random-effect density constants, is $$ h(\eta)=-0.2-\frac{1}{2}(\eta-0.3)^2. $$ Here $\eta$ is measured in standardized random-effect units. This example isolates the steps of the calculation; it does not require recovering the underlying observation model from $h$.

Find the mode. Differentiation gives $h'(\eta)=-(\eta-0.3)$. The unique stationary point is $\hat{\eta}=0.3$, and $h''(\eta)=-1$, so it is a maximum and $J=1$.

Evaluate the height. Since $h(0.3)=-0.2$, the height is $\exp\{h(\hat{\eta})\}=\exp(-0.2)\approx0.8187$.

Apply the formula. $$ L\approx\exp(-0.2)\sqrt{\frac{2\pi}{1}} =0.8187(2.5066)\approx2.052. $$ The integral is in units of $\eta$. Its value exceeding one is not an error: an integral of a likelihood times a density need not itself be a probability.

A Nonquadratic Log Integrand

Worked Example: A Nonlinear Mean With a Centered Mode

Suppose a hypothetical subject has $\eta\sim N(0,0.25)$ and one observation $y\mid\eta\sim N(1+\eta^2,0.25)$. Let $y=1$, in standardized response units. The mean depends nonlinearly on the random effect, and the complete log integrand is $$ h(\eta)=-\log(2\pi\cdot0.25)-2\eta^4-2\eta^2. $$ The constant term includes the normalizing constants of both normal densities.

Locate and check the mode. Differentiation gives $h'(\eta)=-8\eta^3-4\eta=-4\eta(2\eta^2+1)$. The only real stationary point is $\hat{\eta}=0$. For every $\eta\ne0$, $h(\eta)-h(0)=-2\eta^4-2\eta^2<0$, so this is the global maximum.

Calculate curvature and height. The second derivative is $h''(\eta)=-24\eta^2-4$, hence $J=-h''(0)=4>0$. At the mode, $\exp\{h(0)\}=1/(2\pi\cdot0.25)=2/\pi\approx0.6366$.

Calculate the Laplace approximation. $$ L\approx\frac{2}{\pi}\sqrt{\frac{2\pi}{4}} =\frac{2}{\pi}\sqrt{\frac{\pi}{2}} \approx0.7979. $$ The mode and local curvature are well-defined, but $h$ contains a quartic term and is not exactly quadratic. The displayed value is therefore the Laplace approximation, not an exact integral obtained merely from the local curvature.

Checks, Alternatives, and Limitations

A careful calculation separates the mode search from the final integration. In a numerical implementation, a root finder or optimizer may be used to locate $\hat{\eta}$, but the result should still be checked against the defining conditions. Evaluate $h'(\hat{\eta})$ to confirm that it is close to zero, verify that $h''(\hat{\eta})<0$, and inspect the log integrand over a range around the mode. If the purported mode is a minimum, a boundary point, or only one of several local maxima, the simple single-mode formula may not be appropriate as applied.

The formula also depends on using the same random-effect coordinate in the integral and the curvature calculation. If a calculation changes coordinates, the integral’s density must include the corresponding change-of-variable factor. Ignoring it changes the integral, even if an optimizer still finds a plausible-looking mode.

When the integral is analytically available, exact marginalization is preferable as a check; the Gaussian benchmark above illustrates this. Numerical quadrature is another direct alternative in one dimension: it approximates the integral from evaluations of the integrand over a range of $\eta$, rather than using only the local quadratic behavior at the mode. First-order Taylor-based approximations take a different route by linearizing predictions and approximating the resulting likelihood. These approaches are not interchangeable: the Laplace calculation specifically uses the height and curvature of the log integrand at its conditional mode.

A common mistake is to use $h''(\hat{\eta})$ itself in the denominator as if it were positive, or to omit the square root factor. The curvature in the formula is $J=-h''(\hat{\eta})>0$, and the multiplier is $\sqrt{2\pi/J}$. Another error is to report only the mode: the mode identifies the peak location but does not determine the integral without both its height and curvature. Finally, a positive $J$ verifies local concavity, not global accuracy; a narrow local peak does not by itself establish that the integrand elsewhere is negligible.

Key takeaway: For a scalar random effect, calculate the Laplace approximation in three linked pieces: the complete log-integrand height at the conditional mode, the negative second derivative $J$ at that mode, and the multiplier $\sqrt{2\pi/J}$. Check the mode and curvature before interpreting the result as an approximation to the marginal likelihood.

Check Your Understanding

Use the definitions and calculation steps above. Solutions are included so that each numerical result can be checked.

For $h(\eta)=c-\frac{1}{2}(\eta+0.4)^2$, identify the mode and $J$. What Laplace integral follows from this log integrand?
If $h'(\hat{\eta})=0$ but $h''(\hat{\eta})=3$, can the stated single-mode Laplace formula be used at $\hat{\eta}$? Explain.
In the first worked example, what changes if a normalizing constant independent of $\eta$ is dropped from $h$?
For $h(\eta)=c-\eta^2-\eta^4$, find the mode and calculate the Laplace approximation in terms of $c$.

Solutions. (1) Differentiation gives $h'(\eta)=-(\eta+0.4)$, so $\hat{\eta}=-0.4$, $h''=-1$, and $J=1$. Thus $L\approx\exp(c)\sqrt{2\pi}$. (2) No. The positive second derivative indicates local convexity, not a maximum; $J=-3<0$ cannot be used as the positive Gaussian curvature in the formula. (3) The mode and $J$ stay the same, but the height and approximate integral are multiplied by the reciprocal of the omitted constant, so the result is no longer the correctly normalized marginal likelihood. (4) Since $h'(\eta)=-2\eta-4\eta^3=-2\eta(1+2\eta^2)$, the only real stationary point is $\hat{\eta}=0$; it is the global maximum because $-\eta^2-\eta^4<0$ for $\eta\ne0$. Also $h''(0)=-2$, so $J=2$ and $L\approx\exp(c)\sqrt{\pi}$.