Hessian of the Individual Log Posterior

Why the Hessian at the Individual Mode Matters

For an individual, the random effects are not observed; their contribution to the marginal likelihood is integrated out. The earlier material on marginal likelihood as an integral and conditional random-effect modes provides the setting: the individual integrand combines the conditional data likelihood with the random-effect density. Near a well-behaved mode, the second derivatives of its log determine how sharply that integrand falls away from the mode. This tutorial focuses on calculating and interpreting those derivatives, not on evaluating the subsequent marginal-likelihood approximation.

The key matrix is the negative Hessian of the individual log integrand at its mode. Its eigenvectors identify local directions in random-effect space, while its eigenvalues quantify curvature along those directions. This is local information: it describes the shape near the mode, not necessarily the shape of the full individual distribution.

Definition: For individual $i$, let $y_i$ denote the observations, $\eta_i$ the $q$-dimensional random-effect vector, and $\theta$ the fixed-effect and residual-model parameters. Suppose the random effects have density $p(\eta_i\mid\Omega)$, and the conditional data density is $p(y_i\mid\eta_i,\theta)$. Define the individual log integrand, up to terms independent of $\eta_i$, by $$ h_i(\eta_i)=\log p(y_i\mid\eta_i,\theta)+\log p(\eta_i\mid\Omega). $$ An interior mode $\hat\eta_i$ satisfies $\nabla_{\eta_i}h_i(\hat\eta_i)=0$. The negative Hessian at that mode is $$ H_i=-\left.\nabla_{\eta_i}^2h_i(\eta_i)\right|_{\eta_i=\hat\eta_i}. $$

Calculating the Negative Hessian

Assume a Gaussian random-effect distribution $\eta_i\sim N(0,\Omega)$, where $\Omega$ is positive definite. Its log density contributes a quadratic penalty, so its negative Hessian is $\Omega^{-1}$. Consequently, for a differentiable conditional log likelihood $\ell_i(\eta_i)=\log p(y_i\mid\eta_i,\theta)$,

H_i=\Omega^{-1}-\left.\nabla_{\eta_i}^2\ell_i(\eta_i)\right|_{\eta_i=\hat\eta_i}.

This formula separates prior curvature from data curvature. The prior contributes positive precision in every direction when $\Omega$ is positive definite. The conditional likelihood contribution need not be positive in every direction; the complete matrix must be checked at the mode.

A useful explicit form arises when the observations have a Gaussian residual model with fixed covariance $R_i$, and the conditional mean vector is $f_i(\eta_i)$. Let $J_i$ be the matrix of first derivatives of $f_i$ with respect to $\eta_i$, and let $r_i=y_i-f_i(\eta_i)$. Then

\ell_i(\eta_i)=-\frac{1}{2}r_i^{\mathsf T}R_i^{-1}r_i+\text{constant},

and differentiation gives

H_i=\Omega^{-1}+J_i^{\mathsf T}R_i^{-1}J_i -\sum_{j,k}(R_i^{-1}r_i)_j\,\nabla_{\eta_i}^2 f_{i,jk}(\eta_i) \quad\text{at }\eta_i=\hat\eta_i.

Here $f_{i,jk}$ is the $k$th component of the predicted mean for observation $j$. The middle term is the familiar nonnegative Gauss–Newton contribution. The final term accounts for curvature in the nonlinear predictions when the residuals are nonzero. Dropping that term generally gives an approximation, not the exact observed negative Hessian.

Condition for a regular local maximum: If $h_i$ is twice continuously differentiable near an interior mode and $H_i$ is positive definite, then the mode is a strict local maximum. If $H_i$ is singular or has a negative eigenvalue, the quadratic curvature does not describe a regular, strictly peaked local maximum in every direction.

From Curvature to Directional Precision

The interpretation follows from a second-order Taylor expansion around the mode. Write $\delta=\eta_i-\hat\eta_i$. Since the gradient is zero at the mode,

h_i(\hat\eta_i+\delta) \approx h_i(\hat\eta_i)-\frac{1}{2}\delta^{\mathsf T}H_i\delta.

If $H_i$ is symmetric positive definite, it has an orthogonal eigendecomposition $H_i=Q\Lambda Q^{\mathsf T}$, where $\Lambda$ contains positive eigenvalues $\lambda_1,\ldots,\lambda_q$. Let $z=Q^{\mathsf T}\delta$. The quadratic decrease in log integrand is then $\frac12\sum_{j=1}^q\lambda_j z_j^2$. Thus a large $\lambda_j$ means sharper local curvature along its eigenvector, and a small positive $\lambda_j$ means a flatter direction. In this local quadratic description, the covariance is $H_i^{-1}$; its directional variances are $1/\lambda_j$.

This is a local Gaussian interpretation, not a claim that the full conditional distribution is Gaussian. A small eigenvalue can indicate weak local information, but it does not by itself establish that the full distribution is broad, unimodal, or symmetric. Eigenvalues also have units determined by the random-effect coordinates: if an effect is measured in different units, the numerical curvature changes accordingly.

Worked Calculations

Worked Example: Scalar Nonlinear Prediction at Its Mode

Suppose a hypothetical individual has one observation $y=10$, conditional mean $f(\eta)=10e^\eta$, residual standard deviation $\sigma=2$, and random effect $\eta\sim N(0,0.25)$. Consider $\eta=0$. The residual is $r=10-10e^0=0$, and the score is zero: the likelihood score is $f'(0)r/\sigma^2=10(0)/4=0$, while the prior score is $-0/0.25=0$. Thus $\eta=0$ is a stationary point.

For a scalar prediction, the negative Hessian is $\Omega^{-1}+(f')^2/\sigma^2-rf''/\sigma^2$. At zero, $f'(0)=10$, $f''(0)=10$, and $r=0$, so

H= \frac{1}{0.25}+\frac{10^2}{2^2}-\frac{0(10)}{2^2} =4+25=29.

The curvature is positive, so this stationary point is a strict local maximum. The local variance implied by the quadratic approximation is $1/29=0.03448$ in squared $\eta$-units, and its local standard deviation is approximately $0.1857$ $\eta$-units. The prior alone has variance $0.25$; the data have increased local precision at this mode.

Worked Example: Two Random Effects and Eigenvalue Interpretation

For a hypothetical individual, let the linear conditional mean be $f(\eta)=X\eta$, with $X=\begin{pmatrix}1&1\\1&-1\end{pmatrix}$, residual covariance $R=I$, and $\Omega=\begin{pmatrix}1&0\\0&4\end{pmatrix}$. Take the observed vector to be $y=(2,0)^{\mathsf T}$. The negative Hessian is constant because the mean is linear and the residual covariance is fixed:

H=\Omega^{-1}+X^{\mathsf T}X =\begin{pmatrix}1&0\\0&1/4\end{pmatrix} +\begin{pmatrix}2&0\\0&2\end{pmatrix} =\begin{pmatrix}3&0\\0&2.25\end{pmatrix}.

The mode solves $H\hat\eta=X^{\mathsf T}y$. Since $X^{\mathsf T}y=(2,2)^{\mathsf T}$, the result is $\hat\eta=(2/3,\,8/9)^{\mathsf T}$. The eigenvalues of $H$ are $3$ and $2.25$, with eigenvectors along the two coordinate axes. The local directional variances are $1/3=0.3333$ and $1/2.25=0.4444$ squared random-effect units, respectively. Both directions have positive curvature; the first is locally more precise.

Worked Example: Nonzero Residual Changes Nonlinear Curvature

Consider one hypothetical observation with $f(\eta)=\eta+\eta^2/2$, $y=2$, residual variance $1$, and prior $\eta\sim N(0,1)$. The score is $f'(\eta)\{y-f(\eta)\}-\eta$. At $\eta=1$, $f(1)=1.5$, $f'(1)=2$, and the residual is $0.5$; therefore the score is $2(0.5)-1=0$. The negative Hessian there is

$$ H=1+[f'(1)]^2-rf''(1) =1+2^2-(0.5)(1)=4.5. $$

The positive value confirms strict local maximum curvature at this stationary point. If the residual correction were omitted, the Gauss–Newton value would be $1+4=5$, which is not the exact observed curvature. In this example the nonzero residual reduces the curvature by $0.5$, so omitting it would overstate local precision.

Derivation, Diagnostics, and Alternatives

For the fixed-covariance Gaussian model, the residual is $r=y-f(\eta)$, so $\partial r/\partial\eta=-J$. Differentiating the log likelihood gives $\nabla\ell=J^{\mathsf T}R^{-1}r$. Differentiating again produces a negative term $-J^{\mathsf T}R^{-1}J$ and a residual-weighted second-derivative term:

\nabla^2\ell=-J^{\mathsf T}R^{-1}J+ \sum_{j,k}(R^{-1}r)_j\nabla^2 f_{jk}.

Adding the prior's negative Hessian $\Omega^{-1}$ and changing the sign yields the stated exact formula. The residual term vanishes at a zero-residual fit, but not merely because the mode has zero score: at a mode, likelihood and prior scores can cancel while residuals remain nonzero.

In practice, compute the Hessian at the converged individual mode using analytic derivatives, automatic differentiation, or carefully checked numerical differentiation. Confirm that the mode's score is small, the computed matrix is approximately symmetric, and its eigenvalues are stable to reasonable changes in derivative settings. A substantially negative eigenvalue indicates that the point is not a strict local maximum; a near-zero eigenvalue signals weak or numerically unresolved curvature. Possible causes include an unconverged mode, a flat direction, an inaccurate Hessian, or a nonregular local shape. Symmetrizing a noisy matrix may address tiny numerical asymmetry, but it does not repair incorrect derivatives or a genuinely non-positive curvature matrix.

The main alternative used in practice is the Gauss–Newton curvature, which keeps $\Omega^{-1}+J^{\mathsf T}R^{-1}J$ and omits the residual-weighted second derivatives of $f$. It is often simpler and is exact for a linear mean, but can differ materially for nonlinear predictions with appreciable residuals. Expected information is another alternative: under a correctly specified fixed-covariance Gaussian model, averaging the residual-weighted term over the conditional data distribution gives zero, leaving the prior precision plus the expected $J^{\mathsf T}R^{-1}J$. Expected information describes average curvature under repeated data; it is not generally the observed Hessian for the individual data at hand.

Curvature depends on parameterization. For a linear change of coordinates $\eta=A z$, with nonsingular $A$, the Hessian at the mapped mode is $H_z=A^{\mathsf T}H_\eta A$. Rescaling can therefore change eigenvalues. An orthogonal rotation changes entries and eigenvectors but preserves eigenvalues, because it gives an orthogonal similarity transformation. More general nonorthogonal transformations need not preserve eigenvalues. Accordingly, compare eigenvalues only with attention to the random-effect scales and coordinates in which they were calculated.

Common pitfall: Treating every positive eigenvalue as proof that the full individual random-effect distribution is well approximated by a Gaussian. Positive eigenvalues establish positive quadratic curvature at the mode, not global symmetry, unimodality, or negligible higher-order effects. The correction is to interpret the Hessian as local curvature and use the surrounding model behavior and diagnostics before making broader distributional claims.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

With $\eta\sim N(0,\Omega)$, what term does the prior add to the negative Hessian?
Why can the exact observed curvature differ from the Gauss–Newton curvature when predictions are nonlinear?
If a positive Hessian eigenvalue doubles, what happens to the corresponding local quadratic directional variance?
Does an orthogonal rotation of the random-effect coordinates change the Hessian eigenvalues? What about rescaling?

Solutions. (1) It adds $\Omega^{-1}$. (2) The exact Hessian includes residual-weighted second derivatives of the predictions, which Gauss–Newton omits. (3) The local directional variance, the reciprocal eigenvalue, is halved. (4) An orthogonal rotation preserves eigenvalues, although entries and eigenvectors can change; rescaling can change eigenvalues.