Differentiating the Laplace Approximation with Respect to Model Parameters

Why Differentiating the Approximation Requires Care

In a GLMM, the Laplace approximation replaces an integral over random effects with a quantity evaluated at the conditional mode and its local curvature. When a model parameter changes, the mode generally moves. The curvature used in the approximation can therefore change both directly with the parameter and indirectly because it is evaluated at a different mode. A derivative that holds either quantity fixed can give the wrong gradient for likelihood optimization.

Write \(u\in\mathbb{R}^q\) for the random-effects vector and \(\theta\) for one scalar model parameter, such as a fixed effect or a covariance parameter. Let \(h(u;\theta)\) be the log of the integrand, including the conditional response log likelihood and the random-effects log density. Constants in \(h\) must be retained if they depend on \(\theta\). The integral and its Laplace approximation are

where \(\hat{u}=\hat{u}(\theta)\) is the conditional mode and \(H(u;\theta)=-\partial^2 h(u;\theta)/\partial u\partial u'\) is the negative Hessian with respect to \(u\). The determinant term uses the Hessian at the mode.

Conditions and the General Derivative

For the formulas below, assume that near the parameter value of interest, \(h\) has continuous derivatives sufficient to differentiate its score and Hessian; the mode is an interior stationary point; and \(H(\hat{u};\theta)\) is positive definite. Assume also that this mode is the relevant isolated maximum for the Laplace approximation and that the integral is finite. These conditions ensure a locally smooth mode and a well-defined log determinant. A boundary mode, a singular Hessian, or competing modes can invalidate the ordinary derivative formulas.

The mode satisfies the score equation \(h_u(\hat{u};\theta)=0\), where \(h_u\) is the column vector of derivatives with respect to \(u\). Differentiating this equation gives the mode’s sensitivity to \(\theta\). In the next formula, \(h_{u\theta}\) is the vector of mixed derivatives, and all derivatives are evaluated at \((\hat{u};\theta)\).

The sign in the second equality follows from \(h_{uu}=-H\). Thus, the mode moves more sharply in directions where the curvature is weak, since those directions correspond to smaller eigenvalues of \(H\).

To differentiate the log approximation, first apply the chain rule to \(h(\hat{u}(\theta);\theta)\). The term involving the mode derivative is \(h_u' d\hat{u}/d\theta\), which is zero at the mode by its score equation. The determinant term does not have this cancellation: its Hessian is evaluated at the moving mode. For a positive definite matrix \(H\), the differential identity \(d\log|H|=\operatorname{tr}(H^{-1}dH)\) yields

Here \(H_\theta\) is the partial derivative holding \(u\) fixed, while \(H_{u_k}\) is the matrix of derivatives of \(H\) with respect to the \(k\)th coordinate of \(u\). The sum accounts for the change in the Hessian caused by movement of the mode. Equivalently, the matrix in braces is the total derivative \(dH(\hat{u}(\theta);\theta)/d\theta\). For a vector of model parameters, apply the same expression separately to each coordinate.

When the Mode’s Movement Changes Curvature

In the Gaussian example, \(H\) did not depend on \(u\), so movement of the mode did not change the curvature. For a nonquadratic log integrand, that contribution is generally present. The following smooth scalar example isolates it.

Practical Checks and Alternatives

In a GLMM, fixed effects usually enter the conditional response model, while covariance parameters enter the random-effects density and its normalizing terms. A covariance parameter may also alter the curvature in random-effects directions. The formula is unchanged: form \(h\) with all parameter-dependent terms included, differentiate the mode equation, and differentiate the Hessian at the resulting mode. If a covariance matrix is parameterized through log standard deviations or a Cholesky factor, apply the chain rule from the chosen unconstrained parameters to the covariance entries.

Useful checks include verifying that the fitted score \(h_u\) is close to zero, that the Hessian at the mode is positive definite, and that the mode sensitivity solves \(H(d\hat{u}/d\theta)=h_{u\theta}\). Compare analytic derivatives with centered finite differences of the same log Laplace approximation, re-optimizing the mode at each perturbed parameter value. A disagreement can indicate an incorrect sign, an omitted density normalizer, a failure to update the mode, or a numerical mode that has not converged. Near-zero Hessian eigenvalues make both \(H^{-1}\) and the derivative sensitive to numerical error; a nearly singular or non-positive-definite Hessian is a warning that the local approximation or mode calculation is unreliable.

Analytic differentiation is efficient and avoids finite-difference step-size choices, but it requires correct mixed and third derivatives when the Hessian depends on the mode. Finite differences are easier to implement but can be noisy and costly because each evaluation may require a new mode search. Repeated adaptive quadrature instead approximates the integrated likelihood directly and can provide a useful comparison when feasible, but requires its own numerical convergence checks. These approaches do not remove the need to check whether the mode and curvature underlying a Laplace approximation are well behaved.