Unstable Fixed-Effect Estimates and Curvature

Why Fixed-Effect Estimates Can Be Unstable

A GLMM can return finite fixed-effect estimates and still leave important uncertainty about their values. One warning is an unusually large standard error, but the more informative question is which combination of coefficients is uncertain? Two coefficients may each have large standard errors because the data estimate their sum well but their difference poorly. Looking at curvature in the likelihood makes this distinction visible.

Curvature is a local property: it describes how quickly the fitted log likelihood falls as parameters move away from their estimates. A sharply curved likelihood direction is well constrained by the data. A nearly flat direction permits substantial parameter changes with little change in fit. This tutorial uses the observed-information matrix to connect that geometry to standard errors and to diagnose which fixed-effect combinations deserve caution.

Definition: Let \(\ell(\beta,\theta)\) be the fitted GLMM log likelihood, where \(\beta\) is the vector of fixed effects and \(\theta\) denotes the random-effects covariance parameters. At an estimated parameter vector, the observed information is the negative Hessian of the log likelihood: \(I=-\nabla^2\ell\). When the likelihood is sufficiently smooth and the estimate is an interior, locally identifiable maximum, the inverse observed information approximates the covariance matrix of the estimates. Its diagonal entries give approximate variances; their square roots are approximate standard errors.

In practice, software may calculate curvature from the marginal likelihood or an approximation to it, such as a Laplace approximation. The exact implementation and any parameter transformations matter when interpreting diagnostics. The central idea is the same: small curvature means that nearby parameter values have similar likelihood, so the corresponding uncertainty is large.

From Curvature to Standard Errors

For a fixed-effect vector \(\beta\), suppose temporarily that covariance parameters are held fixed and the negative Hessian for \(\beta\) is a positive-definite matrix \(I_\beta\). A local quadratic approximation to the log likelihood is

\ell(\beta)\approx \ell(\hat\beta)-\frac{1}{2}(\beta-\hat\beta)^\mathsf{T}I_\beta(\beta-\hat\beta).

The curvature along a unit direction \(v\) is \(v^\mathsf{T}I_\beta v\). If \(v\) is an eigenvector of \(I_\beta\) with eigenvalue \(\lambda\), that curvature is \(\lambda\). Inverting the information changes the eigenvalue to \(1/\lambda\), so the variance along that direction is large when \(\lambda\) is small. This explains why a single weak direction can make multiple individual coefficient standard errors large.

The usual Hessian-based standard errors rely on a local normal approximation and regular estimation conditions: a smooth likelihood near the estimate, a locally identifiable parameterization, an interior maximum, and enough information for the approximation to be useful. In a GLMM, estimates near covariance-parameter boundaries or a likelihood that is strongly non-quadratic can make ordinary Hessian-based uncertainty unreliable. A finite inverse Hessian is not, by itself, evidence that these conditions hold.

Worked Example: Finding the Uncertain Coefficient Combination

Suppose a hypothetical GLMM fit has two fixed effects, \(\beta_1\) and \(\beta_2\), and the negative Hessian for those effects is

I_\beta= \begin{pmatrix} 100 & 99\\ 99 & 100 \end{pmatrix}.

The determinant is \(100(100)-99(99)=199\), so the approximate covariance matrix is

I_\beta^{-1}=\frac{1}{199} \begin{pmatrix} 100 & -99\\ -99 & 100 \end{pmatrix}.

Each coefficient has approximate variance \(100/199=0.5025\), and hence standard error \(\sqrt{0.5025}=0.709\). The off-diagonal covariance is \(-99/199=-0.4975\), giving an estimated correlation of \(-0.990\). Thus, the coefficient estimates move strongly in opposite directions together.

The matrix has eigenvalue 199 in direction \((1,1)/\sqrt{2}\), and eigenvalue 1 in direction \((1,-1)/\sqrt{2}\). The sum is therefore much better estimated than the difference. Indeed,

\operatorname{Var}(\hat\beta_1+\hat\beta_2) =\frac{100+100-198}{199} =\frac{2}{199}, \qquad SE(\hat\beta_1+\hat\beta_2)=0.1003,

whereas

\operatorname{Var}(\hat\beta_1-\hat\beta_2) =\frac{100+100+198}{199} =\frac{398}{199}=2, \qquad SE(\hat\beta_1-\hat\beta_2)=1.414.

For example, if the fitted values were \(\hat\beta_1=1.1\) and \(\hat\beta_2=-0.7\), the estimated sum would be \(0.4\), while the estimated difference would be \(1.8\). The difference is much less precise. Reporting only the two individual standard errors would hide this contrast between a well-estimated combination and a weak one.

Curvature When Covariance Parameters Are Estimated

The fixed effects are generally estimated together with random-effects covariance parameters. Write the joint observed-information matrix in blocks:

I= \begin{pmatrix} I_{\beta\beta} & I_{\beta\theta}\\ I_{\theta\beta} & I_{\theta\theta} \end{pmatrix}.

The off-diagonal blocks represent dependence in the likelihood curvature between fixed effects and covariance parameters. Consequently, the fixed-effect covariance matrix from joint estimation is generally the fixed-effect block of \(I^{-1}\), not simply \(I_{\beta\beta}^{-1}\). Treating \(\theta\) as fixed discards uncertainty and dependence associated with estimating it.

Result: If the joint information is positive definite and \(I_{\theta\theta}\) is invertible, the effective information for \(\beta\) after accounting for estimation of \(\theta\) is the Schur complement \(I_{\beta\beta}-I_{\beta\theta}I_{\theta\theta}^{-1}I_{\theta\beta}\). Its inverse is the fixed-effect block of \(I^{-1}\).

This result follows by solving the block linear system for the first block column of \(I^{-1}\). Let that column be \(\begin{pmatrix}A\\C\end{pmatrix}\), so that \(I\begin{pmatrix}A\\C\end{pmatrix}=\begin{pmatrix}\mathrm{Id}\\0\end{pmatrix}\), where \(\mathrm{Id}\) is the identity matrix of the fixed-effect dimension. The second block row gives \(I_{\theta\beta}A+I_{\theta\theta}C=0\), hence \(C=-I_{\theta\theta}^{-1}I_{\theta\beta}A\). Substitution into the first row gives

\left(I_{\beta\beta}-I_{\beta\theta}I_{\theta\theta}^{-1}I_{\theta\beta}\right)A=\mathrm{Id}.

Therefore \(A\), the fixed-effect block of the joint inverse, is the inverse Schur complement. The subtracted term is positive semidefinite when the joint information is positive definite. Accounting for the nuisance parameters can thus reduce effective information about fixed effects, increasing their estimated uncertainty.

Worked Example: Allowing for Covariance-Parameter Uncertainty

Consider a hypothetical fit with one fixed effect \(\beta\) and one covariance parameter \(\theta\), with joint information

I= \begin{pmatrix} 16 & 6\\ 6 & 9 \end{pmatrix}.

If \(\theta\) were incorrectly treated as known, the fixed-effect variance would be approximated by \(1/16=0.0625\), giving \(SE=0.250\). But the determinant of the joint matrix is \(16(9)-6(6)=108\), so its inverse is

I^{-1}=\frac{1}{108} \begin{pmatrix} 9 & -6\\ -6 & 16 \end{pmatrix}.

The joint-estimation variance for \(\hat\beta\) is \(9/108=0.08333\), and its standard error is \(\sqrt{0.08333}=0.2887\). The Schur-complement calculation gives the same result: \(I_{\beta\beta}-I_{\beta\theta}I_{\theta\theta}^{-1}I_{\theta\beta}=16-6(1/9)6=12\), whose inverse is \(1/12=0.08333\). Ignoring the dependence would understate the standard error in this example.

Diagnosing a Weak Direction

A useful diagnosis goes beyond asking whether a standard error is large. Examine the scale and orientation of weak curvature, then express the affected direction as a scientifically interpretable combination of coefficients. Eigenvectors of a suitably scaled information matrix can identify directions with small eigenvalues. Because eigenvectors and eigenvalues depend on parameter units, raw eigenvalues should not be compared blindly across differently scaled coefficients. Consider predictor scaling and coefficient meanings before interpreting the calculation.

Worked Example: Interpreting a Low-Curvature Direction

Suppose a hypothetical fit is expressed in two orthogonal coefficient combinations \(u=(\beta_1+\beta_2)/\sqrt{2}\) and \(v=(\beta_1-\beta_2)/\sqrt{2}\). The information in these coordinates is diagonal, with entries 80 for \(u\) and 0.5 for \(v\). The approximate standard errors are \(1/\sqrt{80}=0.1118\) for \(\hat u\) and \(1/\sqrt{0.5}=1.414\) for \(\hat v\). The difference direction is much less precisely estimated.

Under the local quadratic approximation, moving a distance \(d\) in the \(v\) direction changes the log likelihood by approximately \(-\tfrac12(0.5)d^2\). A change of \(d=2\) therefore gives an approximate log-likelihood change of \(-1\). In the \(u\) direction, a change of only \(d=\sqrt{2/80}=0.1581\) gives the same change. This is a direct numerical illustration of weak versus strong curvature.

The original coefficients are \(\beta_1=(u+v)/\sqrt{2}\) and \(\beta_2=(u-v)/\sqrt{2}\). If the two combinations are uncorrelated in this example, each original coefficient has variance \((1/80+2)/2=1.00625\), giving a standard error of \(1.003\). Their covariance is \((1/80-2)/2=-0.99375\), a correlation of approximately \(-0.988\). Thus, the unstable direction affects both coefficients, even though the sum remains comparatively precise.

Practical Checks, Alternatives, and Limitations

A large fixed-effect standard error is a signal to investigate, not a diagnosis by itself. Check whether the estimate and standard error are stable under sensible changes in parameter scaling or coding, and inspect the coefficient covariance or correlation matrix for strongly dependent estimates. Relate weak directions to the model terms and the data structure. Also check whether the fit is near a boundary or whether numerical curvature calculations are sensitive to small changes; these conditions can undermine a local quadratic approximation.

The main alternatives and complements answer related but distinct questions:

Hessian-based standard errors are efficient local summaries when the likelihood is smooth and approximately quadratic near an interior estimate. They can be misleading when curvature changes rapidly or the estimate is near a boundary.
Profile likelihood examines likelihood behavior as a parameter or contrast varies while other parameters are re-estimated. It can reveal asymmetry or a broad range that a single local standard error conceals.
Refitting and perturbation checks assess whether estimates change under reasonable numerical or data changes. These checks can reveal instability, but a stable result does not prove that the model is correct or that all uncertainty has been captured.

These approaches should not be treated as interchangeable. A weak Hessian direction describes local uncertainty under the fitted model; it does not identify the cause. For example, an unstable fixed-effect combination is not the same thing as a singular random-effects covariance matrix. Conversely, apparent computational convergence does not guarantee strong information about every fixed-effect contrast.

Common Mistake

Wrong reasoning: “Both coefficients have large standard errors, so neither their sum nor their difference can be estimated well.” This ignores covariance between the estimates. In the first worked example, each coefficient had standard error \(0.709\), but the sum had standard error \(0.1003\), while the difference had standard error \(1.414\).

Correction: Use the full covariance matrix to calculate uncertainty for the coefficient combination of interest. For a contrast \(a^\mathsf{T}\beta\), its approximate variance is \(a^\mathsf{T}\operatorname{Cov}(\hat\beta)a\). State the contrast and its units in context; individual standard errors alone do not describe all the information in the fit.

Key takeaway: Large fixed-effect standard errors often reflect weak curvature in a particular coefficient combination. Inspect the covariance structure and the corresponding direction, account for joint estimation of covariance parameters, and treat Hessian-based uncertainty cautiously when regularity or local-quadratic assumptions are doubtful.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

If an information eigenvalue is small, what does that imply about approximate variance in its eigenvector direction?
Why is the fixed-effect block of the inverse joint information generally not equal to \(I_{\beta\beta}^{-1}\)?
A fitted coefficient pair has standard errors 0.8 and a correlation of \(-0.95\). What should you inspect before concluding that every linear combination is poorly estimated?
Why should raw Hessian eigenvalues be interpreted with attention to parameter units?
Name one condition under which Hessian-based standard errors may be unreliable, and one complementary diagnostic.

Solutions. (1) The approximate variance is the reciprocal of the eigenvalue, so a small positive eigenvalue implies large variance in that direction. (2) Joint estimation allows dependence between fixed effects and covariance parameters; the relevant information is the Schur complement, not just the fixed-effect block. (3) Inspect the coefficient covariance matrix and calculate the variance of the specific sum or difference of interest; strong negative correlation can make the sum precise while the difference is uncertain. (4) Rescaling a parameter changes numerical curvature and eigenvalues even though the scientific model may be equivalent, so interpret directions on meaningful or suitably scaled parameter axes. (5) A boundary estimate or a strongly non-quadratic likelihood can invalidate the local approximation; profile likelihood is one complementary check.