Laplace Approximation in Dynamic Prediction

Why Approximate the Random-Effects Integral?

A dynamic prediction for a patient uses both the patient’s observed history and uncertainty about the patient-specific random effects. Conditional on a particular random-effects value, survival predictions may be straightforward to calculate. The difficulty is that a patient’s random effects are not known: the prediction must average over their distribution given the observed longitudinal measurements and survival to the landmark.

Likelihood-based integration expresses that average as an integral over the random-effects vector. Gauss-Hermite quadrature, covered in the earlier material on Gauss-Hermite quadrature for random-effects integration, evaluates such integrals with weighted nodes. The Laplace approximation instead represents the integrand near its dominant mode by a Gaussian-shaped function. It can be computationally efficient, but its accuracy depends on how well that local shape represents the full integral.

Prediction Integral and Laplace Formula

Let $b_i\in\mathbb{R}^q$ denote patient $i$’s random-effects vector, with density $p(b_i)$. Let $y_i(\leq s)$ denote the patient’s longitudinal measurements observed by landmark time $s$, and let $S_i(t\mid b_i)$ be the survival probability through time $t$, conditional on $b_i$ and the model parameters. For this tutorial, the model parameters are treated as fixed at their fitted values. The patient must have remained event-free through $s$.

The likelihood contribution for the observed history, conditional on $b_i$, is proportional to $p\{y_i(\leq s)\mid b_i\}S_i(s\mid b_i)p(b_i)$. Consequently, the conditional survival probability from $s$ to $s+u$ is the ratio below. The survival factor in the numerator includes survival through the entire interval from the origin to $s+u$; dividing by the corresponding integral through $s$ conditions on being event-free at the landmark.

P(T_i>s+u\mid T_i>s,y_i(\leq s)) = \frac{A_i(s+u)}{A_i(s)}, \qquad A_i(t)=\int_{\mathbb{R}^q} p\{y_i(\leq s)\mid b_i\}S_i(t\mid b_i)p(b_i)\,db_i.

Define the log-integrand $\ell_i(t;b_i)$ as the logarithm of the integrand in $A_i(t)$, and let $\hat b_i(t)$ be an interior mode of that log-integrand. Write $J_i(t)=-\nabla^2\ell_i\{t;\hat b_i(t)\}$, where $\nabla^2$ is the Hessian matrix of second derivatives with respect to $b_i$.

Definition: If $\ell_i(t;b_i)$ is twice differentiable near an interior mode $\hat b_i(t)$, and $J_i(t)$ is positive definite, the single-mode Laplace approximation is $$ A_i(t)\approx \exp\!\left[\ell_i\{t;\hat b_i(t)\}\right] (2\pi)^{q/2}|J_i(t)|^{-1/2}. $$ Here $|J_i(t)|$ is the determinant. Positive definiteness means the mode is nondegenerate: the log-integrand curves downward in every direction at the mode.

Why the Approximation Has This Form

The approximation follows from replacing the log-integrand near its mode by its second-order Taylor expansion. At an interior differentiable mode, the gradient is zero. Writing $d=b_i-\hat b_i(t)$, the expansion is

\ell_i(t;b_i) \approx \ell_i\{t;\hat b_i(t)\} -\frac{1}{2}d^\mathsf{T}J_i(t)d.

Exponentiating gives a constant at the mode times a Gaussian-shaped kernel. Integrating that kernel over $\mathbb{R}^q$ yields the determinant factor in the formula:

\begin{aligned} A_i(t) &\approx \exp[\ell_i\{t;\hat b_i(t)\}] \int_{\mathbb{R}^q} \exp\!\left(-\frac12 d^\mathsf{T}J_i(t)d\right)\,dd\\ &= \exp[\ell_i\{t;\hat b_i(t)\}] (2\pi)^{q/2}|J_i(t)|^{-1/2}. \end{aligned}

The last equality follows by transforming $d$ to $z=J_i(t)^{1/2}d$. The Jacobian contributes $|J_i(t)|^{-1/2}$, and the remaining integral is the $q$-dimensional standard Gaussian integral, $(2\pi)^{q/2}$. This argument justifies the formula locally; it does not establish that the local approximation captures all important mass in the original integral.

A twice-differentiable local maximum has a negative-semidefinite Hessian, but that alone does not ensure the Laplace formula is finite: the negative Hessian may be singular. Positive definiteness of $J_i(t)$ is an additional nondegeneracy condition. In the usual asymptotic justification, the log-integrand becomes increasingly concentrated around a unique interior mode as information grows, is sufficiently smooth there, and has negligible integral contribution away from that mode. For an individual patient with limited data, these conditions may be poor approximations rather than guarantees.

Worked Conditional-Risk Calculation

Worked Example: A Quadratic Log-Integrand

Consider a hypothetical patient with a one-dimensional random effect $b$. Suppose the two log-integrands, including the observed-history likelihood, survival factor, and random-effects density, are

\begin{aligned} \ell_i(s;b)&=-b^2+K_s, & K_s&=\log\!\left(\frac{0.951}{\sqrt{\pi}}\right),\\ \ell_i(s+u;b)&=-1.4b^2+K_{s+u}, & K_{s+u}&=\log\!\left(\frac{0.5999}{\sqrt{\pi/1.4}}\right). \end{aligned}

Both modes are at $b=0$. The negative second derivatives are $J_i(s)=2$ and $J_i(s+u)=2.8$, respectively, so both are positive. The first integral is

A_i(s)\approx e^{K_s}\sqrt{\frac{2\pi}{2}} = \frac{0.951}{\sqrt{\pi}}\sqrt{\pi} =0.951.

At the later time, the calculation is

A_i(s+u)\approx e^{K_{s+u}}\sqrt{\frac{2\pi}{2.8}} = \frac{0.5999}{\sqrt{\pi/1.4}}\sqrt{\frac{\pi}{1.4}} =0.5999.

Because these log-integrands are exactly quadratic, the Gaussian integral used by the Laplace approximation is exact. The estimated conditional survival probability is $0.5999/0.951\approx0.6308$. Thus the conditional risk over the interval is $1-0.6308=0.3692$, or about 36.92%. The risk refers to the probability of an event in $(s,s+u]$, conditional on survival to $s$ and the stated longitudinal history.

This construction is also consistent with a conditional survival factor no greater than one at each $b$: the ratio of the two integrands is $\exp(K_{s+u}-K_s-0.4b^2)$. Its maximum occurs at $b=0$ and is less than one, since $\exp(K_{s+u}-K_s)\approx0.7464$.

When Local Curvature Is Not Enough

The formula uses curvature at a mode to estimate the volume of nearby integrand mass. A mode can be a poor summary when the peak is flat, the integrand is strongly skewed, or substantial mass lies near another mode. These problems are especially relevant in individual prediction: a patient may have few measurements, and the likelihood contribution for that patient may not sharply identify the random effects.

Worked Example: A Strict Maximum With Zero Curvature

Let the one-dimensional integrand be $g(b)=e^{-b^4}$, so that $\ell(b)=-b^4$. The point $b=0$ is a strict global maximum, but $\ell''(0)=0$. Therefore $J=-\ell''(0)=0$, and the nondegenerate Laplace formula would require division by $\sqrt{0}$; it is not applicable.

The integral itself is finite. By symmetry and the substitution $x=b^4$ on the positive half-line,

\int_{-\infty}^{\infty}e^{-b^4}\,db = 2\int_0^\infty e^{-b^4}\,db = \frac12\Gamma\!\left(\frac14\right) \approx1.813.

This example shows why a strict local maximum does not, by itself, imply positive curvature. The integrand has a flatter-than-quadratic peak, so its local second-order expansion cannot supply a valid Gaussian approximation.

Worked Example: Two Modes and Missed Integral Mass

For a second hypothetical integrand, take

g(b)=e^{-(b-3)^2/2}+e^{-(b+3)^2/2}.

Each term has integral $\sqrt{2\pi}$, so the exact integral is $2\sqrt{2\pi}\approx5.0133$. The two peaks are well separated and have nearly equal height. A local Laplace calculation around just the positive peak uses height approximately 1 and curvature approximately 1, giving $\sqrt{2\pi}\approx2.5066$. That is about half the total integral: the other mode contributes comparable mass.

If a numerical optimizer returns only one mode, the approximation can therefore be seriously low even though the curvature at that mode is positive. When distinct modes are identified, a local approximation at each can be computed and the contributions summed. This example is a warning to check for multiple modes, not a claim that a single-mode formula automatically detects them.

Practical Checks, Alternatives, and Limitations

Before using a Laplace-based prediction, examine the optimization and the approximation, not just the final risk. Useful checks include whether optimization from different starting values finds the same mode, whether the negative Hessian is positive definite and reasonably conditioned, and whether the fitted log-integrand appears strongly asymmetric or has other peaks. Large changes in the result under modest changes in numerical tolerances or starting values are reasons to investigate further.

Curvature: A zero or negative eigenvalue of $J_i(t)$ means the stated nondegenerate approximation is unavailable. A very small positive eigenvalue can make the determinant factor unstable and signal a broad direction.
Shape: Compare the log-integrand to its quadratic expansion over the region contributing appreciable mass. Strong skewness or heavy tails may make a local Gaussian shape misleading.
Modes: Use multiple starting points and inspect the objective values. A single optimizer output does not establish that the mode is unique.
Prediction ratio: Check that estimated survival is between zero and one, so conditional risk is also between zero and one. A violation signals numerical or approximation trouble; satisfying the bounds does not prove accuracy.

Gauss-Hermite quadrature evaluates the integrand at a set of nodes with assigned weights, whereas Laplace approximation summarizes the integrand through its mode and local curvature. Quadrature can be a useful comparison when dimension and computation permit, though its accuracy also depends on the integration setup and numerical choices. Direct plug-in prediction at an estimated random-effects value is faster but does not integrate over the patient-specific random-effects uncertainty; it is not a substitute for validating an integral approximation.

Laplace approximation is most attractive when the integrand is smooth, dominated by one well-curved interior mode, and expensive to integrate directly. It can be unreliable with sparse patient histories, weak identification, near-boundary modes, substantial skewness, flat directions, or multiple modes. It also approximates the random-effects integral at fixed model parameters; it does not by itself incorporate uncertainty in those parameters.

Key takeaway: The Laplace formula replaces a well-behaved, locally dominant peak by a Gaussian integral. Its validity requires more than finding a maximum: the negative Hessian must be positive definite, and the region around that mode must account for nearly all relevant integral mass.

Check Your Understanding

Answer each question, then compare with the solutions below.

For a two-dimensional log-integrand with mode value $-2$ and $J=\operatorname{diag}(4,9)$, write its Laplace approximation.
Why is a strict local maximum insufficient to guarantee that the nondegenerate Laplace formula applies?
In the quadratic worked example, what conditional risk follows from the displayed integrals?
What diagnostic issue is demonstrated by the two-Gaussian integrand, and why can positive curvature at one mode fail to reveal it?

Solutions. (1) The determinant is $4\cdot9=36$, so the approximation is $e^{-2}(2\pi)/6\approx0.1417$. (2) A twice-differentiable local maximum only guarantees a negative-semidefinite Hessian; it may be singular, as for $\ell(b)=-b^4$ at zero. The formula requires a positive-definite negative Hessian. (3) The conditional survival estimate is $0.5999/0.951\approx0.6308$, so the conditional risk is $1-0.6308\approx0.3692$, or 36.92%. (4) The integrand has two separated modes with comparable mass. Curvature measures the local shape at one mode and does not establish that other modes are absent.