Hessian-Based Standard Errors at the Optimum

Why Curvature Can Measure Uncertainty

An optimizer can return a parameter estimate even when the data determine that parameter only weakly. The estimate alone does not show how much nearby values are supported by the likelihood. At an interior optimum, the local curvature of the objective provides one approximation to that uncertainty: a sharply curved objective corresponds to less local uncertainty, while a shallow one corresponds to more.

Hessian-based standard errors are useful because they summarize local uncertainty for many parameters using derivatives evaluated at the estimate. They are not automatic evidence that the model is well identified or that uncertainty is symmetric. The result depends on the objective's scale, the quality of the optimum and likelihood evaluation, and whether a local quadratic approximation is adequate.

Definition: Let $\boldsymbol{\theta}$ be the vector of estimated parameters, and let $Q(\boldsymbol{\theta})=-2\ell(\boldsymbol{\theta})$, where $\ell$ is the log likelihood. At an estimated interior optimum $\hat{\boldsymbol{\theta}}$, the objective Hessian is the matrix of second derivatives $$ H_Q(\hat{\boldsymbol{\theta}})=\nabla^2 Q(\hat{\boldsymbol{\theta}}). $$ If this matrix is positive definite and the usual likelihood-based large-sample and regularity conditions hold, the approximate covariance matrix is $2H_Q^{-1}$. The approximate standard error for parameter $j$ is the square root of diagonal element $j$ of that covariance matrix.

From the Local Quadratic Approximation to Standard Errors

The key idea is to approximate the objective near its minimum by a quadratic. At an interior optimum, the gradient is zero. With $\boldsymbol{\delta}=\boldsymbol{\theta}-\hat{\boldsymbol{\theta}}$, a second-order Taylor expansion gives

Q(\boldsymbol{\theta}) \approx Q(\hat{\boldsymbol{\theta}}) + \frac{1}{2}\boldsymbol{\delta}^{T}H_Q(\hat{\boldsymbol{\theta}})\boldsymbol{\delta}.

Since $Q=-2\ell$, this is equivalent to

\ell(\boldsymbol{\theta}) \approx \ell(\hat{\boldsymbol{\theta}}) - \frac{1}{4}\boldsymbol{\delta}^{T}H_Q(\hat{\boldsymbol{\theta}})\boldsymbol{\delta}.

A multivariate normal density with covariance matrix $V$ has a quadratic exponent $-\frac{1}{2}\boldsymbol{\delta}^{T}V^{-1}\boldsymbol{\delta}$. Matching this exponent to the local log-likelihood expression gives $V^{-1}\approx H_Q/2$, hence $V\approx 2H_Q^{-1}$. This is a local, large-sample approximation to the sampling covariance of the estimator, not an exact description of the likelihood at all parameter values.

Scaling matters: If the minimized objective is the negative log likelihood $F=-\ell$, its Hessian is half the Hessian of $Q=-2\ell$, and the covariance approximation is $H_F^{-1}$. More generally, if the objective is $c(-\ell)$, its Hessian $H_F$ gives covariance approximation $cH_F^{-1}$. Establish the objective's definition and scaling before inverting its Hessian.

The usual interpretation also relies on more than a successful optimizer exit. The likelihood should be correctly specified for the intended inference, the estimate should be identifiable and interior to its parameter space, and the objective should be sufficiently smooth near the optimum. The likelihood calculation, including any numerical integration or approximation used to evaluate it, must be accurate enough that the computed curvature represents the intended objective. Large-sample likelihood theory is the basis for treating the estimator as approximately normal.

Worked Example: A One-Parameter Curvature Calculation

Suppose a hypothetical concentration model estimates a mean concentration $\mu$ from $n=25$ independent measurements with known residual standard deviation $\sigma=4$ mg/L. The observed sample mean is $12.4$ mg/L. For this normal model, the portion of $Q=-2\ell$ that depends on $\mu$ is $\frac{n}{\sigma^2}(\mu-\bar{x})^2$.

Differentiating twice with respect to $\mu$, the Hessian is the scalar

H_Q=\frac{2n}{\sigma^2} =\frac{2(25)}{4^2} =\frac{50}{16} =3.125\;(\mathrm{mg/L})^{-2}.

The approximate variance and standard error are therefore

\widehat{\mathrm{Var}}(\hat{\mu})=\frac{2}{H_Q} =\frac{2}{3.125}=0.6400\;(\mathrm{mg/L})^2, \qquad SE(\hat{\mu})=\sqrt{0.6400}=0.8000\;\mathrm{mg/L}.

This matches the model's standard error $\sigma/\sqrt{n}=4/\sqrt{25}=0.8000$ mg/L. The check confirms the factor of 2 for the $-2\ell$ objective. It relies here on the stated normal model and known $\sigma$; it is not a general guarantee that estimated pharmacometric parameters have exact normal uncertainty.

Several Parameters: Covariance, Standard Errors, and Correlation

For multiple parameters, the off-diagonal elements of the inverse Hessian matter. They describe local covariance between estimates; inverting only the diagonal elements of the Hessian generally does not give the correct marginal standard errors. The Hessian must also be positive definite for this interior-minimum approximation to produce a positive-definite covariance matrix.

Worked Example: Inverting a Two-Parameter Hessian

Suppose a hypothetical model estimates clearance $\theta_1$, in L/h, and volume $\theta_2$, in L. At the fitted interior optimum, the Hessian of $Q=-2\ell$, in the corresponding parameter units, is

H_Q= \begin{pmatrix} 8 & 2\\ 2 & 5 \end{pmatrix}.

Its determinant is $8(5)-2(2)=36>0$, and its first leading diagonal element is positive, so this symmetric matrix is positive definite. Its inverse and the approximate covariance matrix are

H_Q^{-1}=\frac{1}{36} \begin{pmatrix} 5 & -2\\ -2 & 8 \end{pmatrix}, \qquad V\approx 2H_Q^{-1}= \begin{pmatrix} 0.2778 & -0.1111\\ -0.1111 & 0.4444 \end{pmatrix}.

The standard errors are $\sqrt{0.2778}=0.5270$ L/h for clearance and $\sqrt{0.4444}=0.6667$ L for volume. The approximate correlation is

\frac{-0.1111}{\sqrt{(0.2778)(0.4444)}}=-0.3162.

Thus the local approximation indicates modest negative covariance between the two estimates. It does not establish that either parameter is adequately identified across the full likelihood surface; it summarizes the curvature around this particular optimum.

Parameter Scales and Nonlinear Transformations

Positive or bounded parameters are often estimated on a transformed scale. If the Hessian is calculated for transformed coordinates, its covariance is initially on that same scale. For a smooth transformation $\boldsymbol{\theta}=g(\boldsymbol{\eta})$, the first-order delta method maps covariance approximately as $V_{\theta}\approx J V_{\eta}J^T$, where $J$ is the matrix of first derivatives of $g$, evaluated at the estimate. This mapping can produce a symmetric standard error on the reported scale, even when the underlying transformed-scale uncertainty implies an asymmetric interval.

Worked Example: Log-Scale Curvature for Clearance

Consider a hypothetical clearance estimate obtained using $\eta=\log(CL)$, where $CL$ is measured in L/h. Suppose $\hat{CL}=10.0$ L/h, so $\hat{\eta}=\log(10.0)$, and the Hessian of $Q=-2\ell$ with respect to $\eta$ is 32. The transformed-scale variance and standard error are

\widehat{\mathrm{Var}}(\hat{\eta})=\frac{2}{32}=0.0625, \qquad SE(\hat{\eta})=\sqrt{0.0625}=0.2500.

Because $dCL/d\eta=CL$, the delta-method standard error on the clearance scale is approximately $10.0(0.2500)=2.500$ L/h. A symmetric 1.96-standard-error interval on the log scale, transformed back, is

\left[10.0\exp(-1.96(0.2500)),\; 10.0\exp(1.96(0.2500))\right] \approx [6.126,\;16.32]\;\mathrm{L/h}.

This interval is asymmetric around 10.0 L/h, as expected for a positive parameter analyzed on the log scale. The delta-method standard error of 2.500 L/h is a local, first-order summary; it does not make the log-scale interval symmetric on the clearance scale.

Asymmetry, Bounds, and Practical Checks

A Hessian-based standard error is most naturally paired with a symmetric Wald interval, such as estimate plus or minus 1.96 standard errors. Such an interval can be misleading if the likelihood is skewed, if the estimate is near a constraint, or if a parameter's sampling distribution is far from normal. A bound can also change the estimator's behavior: near a boundary, the unconstrained normal approximation may place substantial probability outside the allowed parameter space.

Worked Example: A Boundary-Adjacent Variance Estimate

Suppose a hypothetical variance parameter $\omega^2$, constrained to be nonnegative, is estimated as $0.00020$ (in squared response units). At that point, suppose the Hessian of $Q=-2\ell$ is $8{,}000{,}000$ in the corresponding inverse-squared units. The local standard error is

SE(\hat{\omega}^2)=\sqrt{\frac{2}{8{,}000{,}000}} =\sqrt{0.00000025} =0.0005000.

A symmetric 1.96-standard-error interval is $0.00020\pm1.96(0.0005000)$, or approximately $[-0.000780,\;0.001180]$. Its lower endpoint is outside the parameter space. That does not mean negative variances are possible; it signals that this symmetric local interval is poorly suited to describing uncertainty here. Simply truncating the endpoint to zero does not repair the normal approximation or account for the boundary's effect on the estimator.

Before interpreting a reported Hessian-based standard error, check that the fitted point is an interior optimum for the parameters being summarized, the Hessian corresponds to the stated objective scale, and its inverse produces finite, positive variances. Compare the standard errors with the parameter magnitudes and units, and consider whether the estimate is close to a bound or transformation limit. These checks can flag concerns, but none proves that the likelihood is locally quadratic over a useful range.

A common mistake is to invert the Hessian of $-2\ell$ and take the square root of its diagonal, omitting the factor of 2. The correction follows directly from the quadratic expansion: for $Q=-2\ell$, the covariance approximation is $2H_Q^{-1}$. A second mistake is to treat a finite standard error as proof of reliable uncertainty. When the likelihood is materially asymmetric or the estimate is boundary-adjacent, profile-likelihood intervals can represent the shape more directly, while bootstrap methods can reflect sampling behavior under their assumptions and computational design. These approaches answer related but not identical questions; they do not make a misspecified model or inaccurate likelihood calculation valid.

Key takeaway: Hessian-based standard errors are local likelihood-curvature approximations. Their scale, units, parameterization, optimum conditions, and proximity to boundaries all affect whether the resulting uncertainty summary is interpretable.

Check Your Understanding

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

A scalar objective is $Q=-2\ell$, and its Hessian at the optimum is 50. What is the approximate standard error?
For a two-parameter fit, why is it generally incorrect to take the reciprocal of each diagonal Hessian element as the parameter variances?
A positive parameter is estimated on its logarithmic scale. What derivative is used to map its log-scale standard error to a first-order standard error on the original scale?
A symmetric Hessian-based interval for a nonnegative parameter extends below zero. Does replacing the lower endpoint by zero make it a valid interval? Explain.

Solutions. (1) The variance approximation is $2/50=0.0400$, so the standard error is $\sqrt{0.0400}=0.2000$, in the parameter's units. (2) The covariance matrix is based on the inverse of the full Hessian, including off-diagonal curvature; those terms affect the diagonal elements of the inverse. (3) For $\theta=\exp(\eta)$, use $d\theta/d\eta=\exp(\eta)=\theta$ at the estimate, so $SE(\hat{\theta})\approx\hat{\theta}SE(\hat{\eta})$. (4) No. Clipping changes the interval without fixing the normal approximation; boundary-aware or likelihood-based uncertainty analysis may be more appropriate.