Hessian Curvature at the Conditional Mode

Why Curvature at the Mode Matters

In a GLMM, the likelihood for a cluster’s observations is obtained by integrating over its random effects. Before carrying out that integration, it is useful to know how sharply the integrand peaks around its conditional mode. The negative Hessian of the log integrand records this local curvature: it describes which directions are sharply concentrated and which directions vary more gradually near the mode.

This tutorial focuses on calculating and interpreting that matrix. It does not treat curvature as a complete description of the integrand: a local calculation cannot by itself rule out skewness, distant modes, or other departures from a single-peaked shape.

Definition and Working Model

For one cluster, let $y_i$ denote response $i$, let $b$ be a $d$-dimensional random-effects vector, and let $h(b)$ be the integrand in the random-effects likelihood integral. Define the log integrand $r(b)=\log h(b)$. Terms that do not depend on $b$ may be omitted when finding derivatives. The conditional mode $\hat b$ is a value that maximizes $r(b)$, with the other model parameters held fixed.

Definition: The negative Hessian of the log integrand at the conditional mode is $$ \mathcal{H}(\hat b)=-\nabla_b^2 r(\hat b). $$ The Hessian $\nabla_b^2 r(\hat b)$ is the matrix of second partial derivatives with respect to the components of $b$. At a smooth local maximum, the Hessian is negative semidefinite, so its negative is positive semidefinite. If $\mathcal{H}(\hat b)$ is positive definite, the mode is locally strictly curved in every direction.

A common setting has conditionally independent observations and Gaussian random effects $b\sim N(0,\Sigma)$, where $\Sigma$ is positive definite. For fixed effects and covariance parameters held fixed, the log integrand is the conditional log likelihood plus the Gaussian log density:

r(b)=\sum_{i=1}^{n}\ell_i(b)-\frac{1}{2}b^\mathsf{T}\Sigma^{-1}b+C,

where $\ell_i(b)$ is observation $i$'s conditional log-likelihood contribution and $C$ does not depend on $b$. The Gaussian term contributes $-\Sigma^{-1}$ to the Hessian of $r$, or $+\Sigma^{-1}$ to its negative Hessian. It is important to include this contribution: curvature from the data alone is not the curvature of the full integrand.

Deriving the Curvature for Canonical Responses

Suppose the conditional linear predictor is $\eta_i=x_i^\mathsf{T}\beta+z_i^\mathsf{T}b$, where $z_i$ is the random-effects design vector, and the conditional response model has a canonical link and unit dispersion. Write its log-likelihood contribution, apart from terms independent of $b$, as $y_i\eta_i-a(\eta_i)$. Then the conditional mean is $\mu_i=a'(\eta_i)$, and the conditional variance is $a''(\eta_i)$.

Because $\partial\eta_i/\partial b=z_i$, the first derivative of this contribution is $(y_i-\mu_i)z_i$. Differentiating again uses $\partial\mu_i/\partial b=a''(\eta_i)z_i$, giving $-a''(\eta_i)z_i z_i^\mathsf{T}$. Adding the Gaussian-prior contribution gives the following result.

Result: For conditionally independent canonical-response observations with Gaussian random effects, $$ \nabla_b r(b)=\sum_{i=1}^{n}z_i(y_i-\mu_i)-\Sigma^{-1}b, \qquad \mathcal{H}(b)=\Sigma^{-1}+\sum_{i=1}^{n}a''(\eta_i)z_i z_i^\mathsf{T}. $$ For binomial observations with a logit link, $a''(\eta_i)=m_i p_i(1-p_i)$, where $m_i$ is the number of trials and $p_i$ is the conditional success probability. For Poisson observations with a log link, $a''(\eta_i)=\mu_i$.

Each data contribution is positive semidefinite because $a''(\eta_i)\geq 0$ and $z_i z_i^\mathsf{T}$ is positive semidefinite. The Gaussian term is positive definite when $\Sigma$ is positive definite. Thus, in this canonical setting, $\mathcal{H}(b)$ is positive definite. The result depends on the stated model and parameterization; for a noncanonical link or a different likelihood, the observed second derivatives should be calculated from that likelihood rather than substituted into this formula without justification.

Worked Curvature Calculations

Worked Example: Binary Random Intercept

Consider a hypothetical cluster with four binary responses $1,0,1,0$, a logit random-intercept model with fixed linear predictor zero, and random-intercept variance $\sigma^2=1$. The conditional probability is $p(b)=\operatorname{logit}^{-1}(b)$, and the log integrand, up to an additive constant, is $$ r(b)=2b-4\log(1+\exp(b))-\frac{b^2}{2}. $$

Its derivative is $r'(b)=2-4p(b)-b$. At $b=0$, $p(0)=0.5$, so $r'(0)=2-4(0.5)-0=0$. The second derivative is $r''(b)=-4p(b)(1-p(b))-1$, which is negative for every finite $b$. Thus $b=0$ is the unique mode, and the negative Hessian there is $$ \mathcal{H}(0)=1+4(0.5)(1-0.5)=1+1=2. $$

The value $2$ is the local curvature in the random-intercept coordinate. Its reciprocal square root, $1/\sqrt{2}\approx 0.7071$, is a local quadratic scale in the units of $b$; it is not an exact standard deviation for the full integrand.

Worked Example: Poisson Random Intercept

Suppose a hypothetical cluster has Poisson counts $0,3,1$, offsets $0,\log(2),0$, and a random intercept $b\sim N(0,0.5)$. Take the fixed-effect contribution to the linear predictor to be zero apart from the offsets. At $b=0$, the conditional means are $1,2,1$. The score for $b$ is the sum of observed counts minus fitted means, minus the prior-precision term: $$ r'(0)=(0-1)+(3-2)+(1-1)-2(0)=0. $$

The negative second derivative at a general $b$ is the prior precision plus the sum of the three conditional means. At zero it is $$ \mathcal{H}(0)=\frac{1}{0.5}+1+2+1=2+4=6. $$ The second derivative of $r$ is negative throughout this example, so the stationary point at zero is a mode. The local quadratic scale is $1/\sqrt{6}\approx 0.4082$ random-intercept units. The larger curvature than in the previous example reflects both the prior precision and the Poisson likelihood contributions at these fitted means.

Worked Example: Eigenvalues with a Random Intercept and Slope

Consider four hypothetical binary observations in a logit model. There are two observations at $x=-1$ and two at $x=1$; at each value of $x$, one response is 1 and one is 0. Let the random effects be $b=(b_0,b_1)^\mathsf{T}$, with predictor $\eta_i=b_0+x_i b_1$, and let $$ \Sigma= \begin{pmatrix} 1&0\\ 0&4 \end{pmatrix}. $$ At $b=(0,0)^\mathsf{T}$, all fitted probabilities are $0.5$. The intercept score sums to zero because there are two successes and two failures. The slope score is also zero: at each $x$, the success and failure residuals $y_i-0.5$ cancel. Thus this point is stationary.

The likelihood part of the negative Hessian is $\sum_i p_i(1-p_i)z_i z_i^\mathsf{T}$, with $z_i=(1,x_i)^\mathsf{T}$. Each weight is $0.25$. The two observations at each $x$ give $$ \sum_i 0.25 \begin{pmatrix} 1&x_i\\ x_i&x_i^2 \end{pmatrix} = \begin{pmatrix} 1&0\\ 0&1 \end{pmatrix}. $$ The prior precision is $\Sigma^{-1}=\begin{pmatrix}1&0\\0&0.25\end{pmatrix}$. Therefore $$ \mathcal{H}(0,0)= \begin{pmatrix} 2&0\\ 0&1.25 \end{pmatrix}. $$ Its eigenvalues are $2$ and $1.25$, with eigenvectors along the intercept and slope coordinates, respectively. The intercept direction has greater local curvature; the corresponding quadratic scales are $1/\sqrt{2}\approx0.7071$ and $1/\sqrt{1.25}\approx0.8944$. These scales are expressed in the chosen random-effect coordinates. Rescaling a predictor or changing coordinates changes the matrix and the numerical eigenvalues, so comparisons require attention to parameterization.

Interpreting and Checking the Matrix

For a positive-definite negative Hessian, an eigenvector identifies a direction in random-effects space and its eigenvalue measures curvature along that direction. A large eigenvalue means the log integrand falls rapidly for small movement in that direction; a small positive eigenvalue means it falls more slowly. Under a local quadratic description, movement of size $1/\sqrt{\lambda}$ along a unit eigenvector changes the log integrand by approximately $1/2$, where $\lambda$ is that eigenvalue. This is local integration-scale information, not an exact probability statement or a guarantee about the integrand far from the mode.

A practical calculation should be checked in context. Confirm that derivatives are with respect to the random effects, that all likelihood and Gaussian-density terms depending on $b$ are included, and that the Hessian is evaluated at the conditional mode rather than at zero by default. For an analytic Hessian, verify symmetry and compare selected entries with finite-difference approximations to the gradient. At a purported strict maximum, the negative Hessian should be positive definite; a Cholesky factorization or eigenvalue calculation can check this numerically.

A negative or nearly zero eigenvalue deserves investigation rather than automatic correction. Possible causes include a mode-finding algorithm that has not converged, an incorrect derivative or sign, numerical error, a nearly flat direction, or a stationary point that is not a strict maximum. A large ratio of the largest to smallest eigenvalue indicates poor conditioning: some directions are much more sharply curved than others. Poor scaling of random-effect coordinates can contribute, so interpret eigenvalues alongside the coordinate units and model design.

Common Pitfalls and Alternatives

Reversing the sign. At a maximum, the Hessian of the log integrand is usually negative definite; the negative Hessian is the positive-curvature matrix. Reporting the Hessian itself as positive curvature reverses the interpretation.
Leaving out the prior. The curvature of the conditional likelihood alone omits $\Sigma^{-1}$. The integration integrand includes the random-effects density, so its curvature includes the prior contribution.
Using expected information without saying so. The observed Hessian is evaluated at the fitted conditional means and data for the specified integrand. Fisher information is an expectation under a model. They can coincide in particular canonical GLMM calculations, but they are distinct concepts and should not be interchanged automatically.
Treating a local scale as an exact spread. The inverse square roots of eigenvalues summarize a quadratic approximation near the mode. They do not establish that the full integrand is Gaussian or symmetric.

Analytic differentiation is usually the clearest approach when the likelihood has a simple form, as in the canonical examples above. Numerical differentiation is a useful check or fallback, but step size, cancellation, and scale can distort small eigenvalues. Expected-information or Fisher-scoring matrices may be useful for mode-finding algorithms, but they answer a different calculation question from the observed negative Hessian of the log integrand at the mode. In every approach, the assumed likelihood, random-effects density, and coordinate system must be made explicit.

The matrix is most informative when the conditional mode is reliably located and the local curvature is well behaved. It has limits: a single Hessian does not detect all distant modes, describe substantial skewness, or establish the accuracy of an integration method. Those questions require examining more than the second derivatives at one point.

Check Your Understanding

Answer each question, then compare with the solutions below.

In a binomial-logit random-intercept model with four trials and $p=0.5$ at the mode, what is the likelihood contribution to the negative Hessian?
For a Gaussian random effect with variance $0.25$, what is its prior-precision contribution to the negative Hessian?
What does a small positive eigenvalue indicate locally, and what does it not establish?
If a computed negative Hessian at a claimed strict mode has a negative eigenvalue, name two checks to make before interpreting it.

Solutions. 1. The contribution is $m p(1-p)=4(0.5)(0.5)=1$. 2. The prior precision is $1/0.25=4$. 3. It indicates relatively weak curvature, and hence a broader local quadratic scale, in its eigenvector direction. It does not establish that the full integrand has that spread or is approximately Gaussian away from the mode. 4. Check that the mode-finding procedure converged to a maximum and that the Hessian sign, derivatives, and prior contribution were computed correctly; also check numerical accuracy and coordinate scaling.