Curvature Derivatives and the Laplace Score

Why Curvature Appears in the Laplace Score

A Laplace approximation to a GLMM likelihood depends on more than the value of the integrand at its conditional mode. It also depends on the local curvature there: a sharply peaked integrand and a broadly peaked integrand have different approximate integrals even if their peak heights agree. When a model parameter changes, this curvature changes too, and the derivative of that adjustment is part of the Laplace score.

The focus here is the derivative of the log-determinant term. The earlier material on differentiating the Laplace approximation and the envelope theorem supplies useful context: the derivative of the log integrand evaluated at its mode has a simplified form, but that simplification does not remove the derivative of the curvature adjustment.

The notation \(h_b\) and \(h_{bb}\) denotes the gradient and Hessian with respect to \(b\). A subscript \(\theta\) denotes a partial derivative holding \(b\) fixed. A derivative of \(J(\theta)\), by contrast, is a total derivative: it accounts for both the explicit dependence on \(\theta\) and the fact that \(\hat b\) moves as \(\theta\) changes.

For this calculation, assume \(h\) has continuous derivatives sufficient for the displayed derivatives, the mode is interior, and \(h_b(\hat b;\theta)=0\). Assume also that \(h_{bb}(\hat b;\theta)\) is nonsingular. Since a Hessian at a local maximum is negative semidefinite, nonsingularity makes it negative definite; thus \(J\) is positive definite. A strict local maximum alone does not ensure this: a strict maximum can have zero curvature in some direction. The nonsingularity condition is what permits the implicit differentiation below and ensures the ordinary Laplace determinant is well-defined and positive.

The Matrix Identity Behind the Derivative

For a differentiable, nonsingular matrix \(J(t)\) with positive determinant, the derivative of its log determinant is

Here \(\operatorname{tr}(A)\) is the sum of the diagonal entries of a square matrix \(A\). To justify the identity, take a small increment \(\varepsilon\) and factor the determinant:

where \(K=dJ/dt\), and the last equality is the first-order expansion of the determinant about the identity matrix. Taking logarithms and dividing by \(\varepsilon\) gives the stated derivative. This identity applies to the total derivative of \(J(\theta)\), not merely its partial derivative at a fixed \(b\).

Differentiating Curvature at the Moving Mode

Differentiate the mode equation \(h_b(\hat b(\theta);\theta)=0\) with respect to \(\theta\). The chain rule gives \(h_{bb}\,d\hat b/d\theta+h_{b\theta}=0\), with derivatives evaluated at the mode. Since \(h_{bb}=-J\), the mode sensitivity is

This step uses the nonsingularity of \(h_{bb}\); the implicit function theorem then gives a locally differentiable mode branch. Now differentiate \(J(\theta)=-h_{bb}(\hat b(\theta);\theta)\). In component notation, with \(j,k\in\{1,\ldots,q\}\),

Thus the total derivative contains an explicit parameter term and a mode-movement term. Substituting the mode sensitivity into the derivative of the Laplace log integral yields the score.

The notation \(h_{bbb}[v]\) means contraction of the third-derivative array with vector \(v\): its \((j,k)\) entry is \(\sum_a h_{b_jb_kb_a}v_a\). The first term is the envelope contribution. The determinant term is separate and generally cannot be dropped. For a vector parameter, apply these formulas component by component.

Checking the Calculation and Understanding Its Limits

A practical implementation should check both the conditional mode and its curvature. At a proposed mode, verify that the gradient \(h_b\) is close to zero and that the eigenvalues of \(J\) are positive and not so small that \(J^{-1}\) is numerically unstable. Then compare the analytic or automatic-differentiation score with a centered finite difference of the Laplace log likelihood over several step sizes. The finite-difference check should recompute the mode and curvature at each perturbed parameter value; holding the mode fixed checks a different derivative.

There are important failure modes. A nearly singular Hessian can make the mode sensitivity and trace term very large, while an exactly singular Hessian invalidates the displayed inverse and ordinary nondegenerate Laplace formula. If competing modes exchange which one is selected as the global mode, the selected approximation may fail to be differentiable at the switch. Numerical instability in the mode or Hessian can also create an apparent score discrepancy. These are reasons to inspect gradients, eigenvalues, mode solutions, and sensitivity to numerical tolerances rather than trusting a score value in isolation.

The score derived here is the derivative of the Laplace-approximated log integral. It is not generally the exact score of the integrated likelihood, because the Laplace approximation itself can have error. Exact integration or adaptive quadrature provides an alternative when feasible; numerical finite differences can serve as a check, but can be sensitive to step size and mode-solving error. Omitting the determinant derivative is not a valid alternative score: it differentiates only the mode-height part of the approximation.

A frequent mistake is to argue that any strict local maximum guarantees positive-definite \(J\). That inference is false: for \(h(b)=-b^4\), \(b=0\) is a strict local maximum but \(h_{bb}(0)=0\). The correction is to require a nonsingular Hessian at the local maximum (or explicitly assume a nondegenerate maximum); then negative semidefiniteness plus nonsingularity implies negative definiteness and \(J\) is positive definite. This distinction matters both for the validity of the determinant and for differentiability of the mode.