Observed Hessians For Noncanonical GLMs

Why the Observed Hessian Has Two Terms

The regression score tells us the direction in which the log-likelihood changes as the coefficient vector changes. Its derivative, the observed Hessian, describes how that direction itself changes. For a noncanonical link, differentiating the score produces two distinct contributions: one determined by the fitted mean and variance, and another multiplied by the observed residual. The second contribution is easy to miss, but it can materially affect the curvature at the data actually observed.

We derive the Hessian for independent observations in a one-parameter exponential-dispersion model. The dispersion parameter is treated as fixed. The design vectors are fixed, and the log-likelihood is assumed to be twice differentiable in the regression coefficients on the parameter region under consideration.

Definition: Let $Y_1,\ldots,Y_n$ be independent observations with $E[Y_i]=\mu_i$ and $\operatorname{Var}(Y_i)=\phi V(\mu_i)$, where $\phi>0$ is fixed and $V(\mu)>0$. Let $\eta_i=x_i^T\beta=g(\mu_i)$, with inverse link $m=g^{-1}$, so $\mu_i=m(\eta_i)$. The observed Hessian is the matrix $H(\beta)$ of second derivatives of the log-likelihood: $H(\beta)=\partial^2\ell(\beta)/(\partial\beta\,\partial\beta^T)$.

Assume that the mean-dependent part of each observation's log-likelihood has derivative $(y_i-\mu_i)/\{\phi V(\mu_i)\}$ with respect to $\mu_i$, as it does in the usual one-parameter exponential-dispersion likelihood. Let $m$ be twice differentiable, let $g$ be twice differentiable with $g'(\mu)\ne 0$ on the relevant mean range, and let $V$ be differentiable there. Write $r_i=y_i-\mu_i$ for the observed residual. The general regression score can then be written using the link-scale weight

q_i=\frac{m'(\eta_i)}{\phi V(\mu_i)} =\frac{1}{\phi g'(\mu_i)V(\mu_i)}.

The score is $U(\beta)=\partial\ell(\beta)/\partial\beta=\sum_{i=1}^n x_iq_ir_i$. Here $q_i$ depends on the fitted mean, so it changes with $\beta$; the residual $r_i$ changes with $\beta$ as well. Both changes must be differentiated.

Deriving the Observed Hessian

Theorem: Under the stated differentiability conditions, define $W_i=m'(\eta_i)^2/\{\phi V(\mu_i)\}$, and let $\dot q_i=dq(\eta_i)/d\eta_i$. Then $$ H(\beta)=\sum_{i=1}^n x_ix_i^T\left\{r_i\dot q_i-W_i\right\} =X^T\operatorname{diag}(r_i\dot q_i-W_i)X. $$ The term $-W_i$ is the fitted-weight contribution, while $r_i\dot q_i$ is the residual-dependent curvature contribution.

To prove the result, regard $q_i=q(\eta_i)$ and $r_i=y_i-m(\eta_i)$. Since $\eta_i=x_i^T\beta$, the chain rule gives the row derivatives

\frac{\partial q_i}{\partial\beta^T}=\dot q_i x_i^T, \qquad \frac{\partial r_i}{\partial\beta^T}=-m'(\eta_i)x_i^T.

Differentiate the score contribution $x_iq_ir_i$. The design vector $x_i$ is fixed, so its derivative is zero. Applying the product rule to the other two factors gives

\frac{\partial(x_iq_ir_i)}{\partial\beta^T} =x_i\left(r_i\dot q_i x_i^T-q_i m'(\eta_i)x_i^T\right) =x_ix_i^T\left(r_i\dot q_i-q_i m'(\eta_i)\right).

By definition, $q_i=m'(\eta_i)/\{\phi V(\mu_i)\}$. Therefore $q_i m'(\eta_i)=W_i$. Summing the differentiated score contributions over $i$ proves the stated formula. Stacking the rows $x_i^T$ into the design matrix $X$ gives the equivalent matrix expression. $\square$

For the noncanonical-link derivative in this formula, the inverse-function rule gives $m'(\eta)=1/g'(\mu)$, where $\mu=m(\eta)$. Differentiating once more with respect to $\eta$ gives $m''(\eta)=-g''(\mu)/[g'(\mu)]^3$. The cubic power in this expression matters when computing $\dot q$.

Derivative identity: If $g$ is twice differentiable, $g'(\mu)\ne0$, and $V$ is positive and differentiable, then $$ \dot q=\frac{dq}{d\eta} =-\frac{1}{\phi}\left\{ \frac{g''(\mu)}{[g'(\mu)]^3V(\mu)} +\frac{V'(\mu)}{[g'(\mu)]^2[V(\mu)]^2} \right\}. $$ Here $V'(\mu)$ and $g'(\mu)$ are derivatives with respect to $\mu$, while $\dot q$ is a derivative with respect to $\eta$.

For a direct derivation, start from $q=m'/\{\phi V(m)\}$. Differentiating with respect to $\eta$ gives

\dot q =\frac{m''(\eta)}{\phi V(\mu)} -\frac{[m'(\eta)]^2V'(\mu)}{\phi[V(\mu)]^2}.

Substitute $m'=1/g'(\mu)$ and $m''=-g''(\mu)/[g'(\mu)]^3$. The first term becomes $-g''(\mu)/\{\phi[g'(\mu)]^3V(\mu)\}$, and the second becomes $-V'(\mu)/\{\phi[g'(\mu)]^2[V(\mu)]^2\}$. Combining them gives the identity. $\square$

Worked Hessian Calculations

Worked Example: A Gaussian Model with a Log Link

Suppose $Y_i\sim N(\mu_i,1)$ independently, with log link $\eta_i=\log\mu_i=x_i^T\beta$. Thus $m(\eta)=e^\eta=\mu$, $\phi=1$, and $V(\mu)=1$. For one observation with $x=1$, $y=3$, and $\beta=0$, calculate the score and observed Hessian.

At $\beta=0$, $\eta=0$, so $\mu=e^0=1$ and $r=y-\mu=2$. Since $m'(\eta)=\mu$, we have $q=\mu=1$, and the score is $xqr=1\cdot1\cdot2=2$. The fitted weight is $W=(m')^2/V=1$. Also, $\dot q=d(e^\eta)/d\eta=e^\eta=1$. Hence

H=x^2(r\dot q-W)=1^2(2\cdot1-1)=1.

This can be checked directly. Apart from a constant, the log-likelihood is $\ell(\beta)=-\tfrac12(3-e^\beta)^2$. Its first derivative is $e^\beta(3-e^\beta)$, and its second derivative is $e^\beta(3-e^\beta)-e^{2\beta}$. At $\beta=0$, these are $2$ and $2-1=1$, respectively.

Worked Example: A Two-Coefficient Hessian

Use the same Gaussian model and log link, now with $x_1=(1,0)^T$, $x_2=(1,1)^T$, observed values $y_1=2$, $y_2=1$, and coefficients $\beta=(0,\log 2)^T$. Calculate the score and Hessian.

The linear predictors are $\eta_1=0$ and $\eta_2=\log 2$, so the fitted means are $\mu_1=1$ and $\mu_2=2$. The residuals are $r_1=1$ and $r_2=-1$. For this link and variance function, $q_i=\mu_i$, $\dot q_i=\mu_i$, and $W_i=\mu_i^2$.

The score contributions are $x_1q_1r_1=(1,0)^T$ and $x_2q_2r_2=(-2,-2)^T$. Therefore $U=(-1,-2)^T$. For the Hessian, the scalar coefficient multiplying $x_ix_i^T$ is $r_i\dot q_i-W_i$. It is $1\cdot1-1^2=0$ for observation 1 and $(-1)\cdot2-2^2=-6$ for observation 2. Thus

H =0\begin{pmatrix}1&0\\0&0\end{pmatrix} -6\begin{pmatrix}1&1\\1&1\end{pmatrix} =\begin{pmatrix}-6&-6\\-6&-6\end{pmatrix}.

The result is symmetric, as a twice continuously differentiable log-likelihood Hessian should be. It is singular here: both observations together do not give a full-rank Hessian at this parameter value. In particular, a computed Hessian need not be invertible merely because the score and its derivatives are well defined.

Worked Example: Checking the Link Derivative

Let the response be Gaussian with $\phi=1$ and $V(\mu)=1$, and use the noncanonical link $g(\mu)=\mu^2$ on the domain $\mu>0$. For one observation with $x=1$, $y=2$, and $\eta=\beta=1$, verify $\dot q$ and the Hessian.

Since $\eta=\mu^2$, the fitted mean at $\eta=1$ is $\mu=1$. Here $g'(\mu)=2\mu$, $g''(\mu)=2$, and $q=1/g'(\mu)=1/(2\mu)=1/2$. The derivative identity gives

\dot q=-\frac{g''(\mu)}{[g'(\mu)]^3} =-\frac{2}{2^3}=-\frac14.

This agrees with $q(\eta)=1/(2\sqrt{\eta})$, whose derivative is $-1/(4\eta^{3/2})$, equal to $-1/4$ at $\eta=1$. The residual is $r=2-1=1$, and the fitted weight is $W=[m'(\eta)]^2=1/4$. Therefore

H=r\dot q-W=1\left(-\frac14\right)-\frac14=-\frac12.

A direct check uses $\ell(\beta)=-\tfrac12(2-\sqrt{\beta})^2$, up to a constant, on $\beta>0$. Its first derivative is $(2-\sqrt{\beta})/(2\sqrt{\beta})$, and its second derivative at $\beta=1$ is $-1/2$, as required.

Interpreting the Two Contributions

The fitted-weight part $W_i=m'(\eta_i)^2/\{\phi V(\mu_i)\}$ is nonnegative when the variance is positive. By itself, it contributes the negative-semidefinite matrix $-X^T\operatorname{diag}(W_i)X$ to the Hessian. The residual-dependent part, $X^T\operatorname{diag}(r_i\dot q_i)X$, can have either sign. Consequently, the observed Hessian is not generally just the negative of a weighted cross-product, and it need not be negative semidefinite at every parameter value.

The residual term is also genuinely data-dependent: two data sets with the same fitted means and design matrix can have different observed Hessians because their residuals differ. At a score solution, the score components sum to zero, but that does not imply that each residual is zero or that the residual-dependent Hessian term disappears. The score condition is a set of weighted sums; the Hessian contains different, observation-specific weights.

Common mistake: Keeping only $-X^T\operatorname{diag}(W_i)X$ because it resembles a weighted least-squares curvature matrix. That drops the derivative of $q_i$, even though $q_i$ changes with the fitted mean under a general link. The correction is to differentiate both $q_i$ and $r_i$: the former supplies $r_i\dot q_i$, and the latter supplies $-W_i$.

A second frequent error is to differentiate $q$ with respect to $\eta$ but use the derivative $d/d\mu$ without multiplying by $d\mu/d\eta$. For example, $m''(\eta)=-g''(\mu)/[g'(\mu)]^3$, not a formula with a fourth power of $g'(\mu)$ in its denominator. Keeping the differentiation variable visible at each step prevents this chain-rule error.

Key takeaway: For a noncanonical GLM, the observed Hessian combines a fitted-weight term and a residual-dependent curvature term. Differentiate the score as a product, and track carefully whether each derivative is with respect to the mean $\mu$ or the linear predictor $\eta$.

Check Your Understanding

Use the definitions and derivation above to answer these questions.

In the general Hessian formula, which part comes from differentiating the residual, and which part comes from differentiating the score weight?
For a noncanonical link, why does $\dot q$ generally depend on derivatives of both the link and the variance function?
Does a zero regression score imply that every residual-dependent contribution to the observed Hessian is zero? Explain.
For $g(\mu)=\mu^2$ and $V(\mu)=1$, evaluate $\dot q$ at $\mu=1$ and identify the power of $g'(\mu)$ in its denominator.
Can the observed Hessian fail to be negative semidefinite even when every fitted weight $W_i$ is positive? Explain using the two contributions in the formula.

Observed Hessians For Noncanonical GLMs

What You'll Learn