The General GLM Regression Score

From the Likelihood to the Regression Score

A regression score measures how the log-likelihood changes when the regression coefficients change. In a generalized linear model, changing a coefficient changes each linear predictor, which changes a fitted mean, which in turn changes the likelihood contribution from that observation. The score formula keeps these steps visible: each response residual is scaled by the variance at its fitted mean, the dispersion, and the derivative of the link.

A variance function alone does not determine an exact likelihood score for every possible response distribution. The derivation below assumes independent responses from a one-parameter exponential-dispersion family. It therefore applies to likelihood-based GLMs under that model; a model specified only through its mean and variance does not automatically have this likelihood score.

Definition: For observations $Y_1,\ldots,Y_n$ with fixed design matrix $\mathbf X$ of dimension $n\times p$, write $x_i^T$ for row $i$ of $\mathbf X$. A GLM specifies $\eta_i=x_i^T\beta$ and $g(\mu_i)=\eta_i$, where $\mu_i=E[Y_i]$, $g$ is the link, and $\beta\in\mathbb R^p$. Assume independent responses with variance $\operatorname{Var}(Y_i)=\phi V(\mu_i)$, where $\phi>0$ is held fixed when differentiating with respect to $\beta$. The regression score is the column vector $U(\beta)=\partial\ell(\beta)/\partial\beta$.

For the exact likelihood calculation, use the exponential-dispersion form $f(y_i;\theta_i,\phi)=\exp\{[y_i\theta_i-b(\theta_i)]/\phi+c(y_i,\phi)\}$. Assume the natural parameter values are in the interior of their parameter space, $b$ is twice differentiable there, and differentiation is valid for the stated density or mass function. The mean and variance relations are $E[Y_i]=b'(\theta_i)=\mu_i$ and $\operatorname{Var}(Y_i)=\phi b''(\theta_i)=\phi V(\mu_i)$. Also assume the inverse link is differentiable at the fitted predictors, with $g'(\mu_i)\ne 0$, so that the chain rule below applies.

Deriving the General Score Formula

First differentiate one observation's log-likelihood with respect to its natural parameter. Since $c(y_i,\phi)$ does not depend on $\theta_i$, and $b'(\theta_i)=\mu_i$,

\frac{\partial \ell_i}{\partial \theta_i} =\frac{y_i-b'(\theta_i)}{\phi} =\frac{y_i-\mu_i}{\phi}.

The derivative of the mean with respect to the natural parameter is $d\mu_i/d\theta_i=b''(\theta_i)=V(\mu_i)$. Because $V(\mu_i)>0$, the inverse derivative is $d\theta_i/d\mu_i=1/V(\mu_i)$. Applying the chain rule gives the derivative with respect to the mean:

\frac{\partial \ell_i}{\partial \mu_i} = \frac{\partial \ell_i}{\partial \theta_i} \frac{d\theta_i}{d\mu_i} = \frac{y_i-\mu_i}{\phi V(\mu_i)}.

Next, $g(\mu_i)=\eta_i$, so $d\mu_i/d\eta_i=1/g'(\mu_i)$. The linear predictor has derivative $\partial\eta_i/\partial\beta=x_i$. Thus the contribution from observation $i$ to the coefficient score is

\frac{\partial \ell_i}{\partial\beta} = \frac{y_i-\mu_i}{\phi V(\mu_i)g'(\mu_i)}x_i.

Theorem: Under the independent exponential-dispersion GLM assumptions above, with fixed $\phi$, the regression score is $$ U(\beta) = \sum_{i=1}^n x_i\frac{y_i-\mu_i}{\phi V(\mu_i)g'(\mu_i)} = \mathbf X^T \operatorname{diag}\!\left( \frac{1}{\phi V(\mu_i)g'(\mu_i)} \right)_{i=1}^n (y-\mu), $$ where $y=(y_1,\ldots,y_n)^T$ and $\mu=(\mu_1,\ldots,\mu_n)^T$. The result holds at parameter values where the stated differentiability and link conditions hold.

Proof. Independence makes the full log-likelihood the sum $\ell(\beta)=\sum_i\ell_i$. The chain rule calculation gives each vector contribution $\partial\ell_i/\partial\beta$ above. Summing these contributions gives the first expression. In the matrix expression, the $i$th diagonal entry multiplies $y_i-\mu_i$, and multiplication by $\mathbf X^T$ then adds that scaled residual times $x_i$ for each observation. This is exactly the same sum. The result is a $p$-vector, as required for derivatives with respect to $p$ coefficients. $\square$

The formula uses $g'(\mu_i)$, the derivative of the link with respect to the mean. It does not use the derivative of the inverse link in that position: the inverse-link derivative is $1/g'(\mu_i)$. Also, if $\phi$ is unknown, this is still the partial score for $\beta$ with $\phi$ held fixed; it is not by itself a joint score for $(\beta,\phi)$.

Worked Applications

Worked Example: Inverse Gaussian Regression With a Log Link

Suppose the response has an inverse Gaussian distribution with mean $\mu$ and variance $\phi\mu^3$, so $V(\mu)=\mu^3$. Use the log link $g(\mu)=\log\mu$, known dispersion $\phi=0.5$, design rows $x_1^T=(1,0)$, $x_2^T=(1,1)$, and coefficient candidate $\beta=(\log 2,\log 2)^T$. Let the observed responses be $y_1=1$, $y_2=6$.

The linear predictors are $\eta_1=\log 2$ and $\eta_2=2\log 2$, so the fitted means are $\mu_1=2$ and $\mu_2=4$. Since $g'(\mu)=1/\mu$, the score multiplier for observation $i$ is $1/\{\phi\mu_i^3(1/\mu_i)\}=1/(\phi\mu_i^2)$.

For observation 1, the multiplier is $1/(0.5\cdot 2^2)=0.5$, so its scalar residual contribution is $(1-2)(0.5)=-0.5$. For observation 2, the multiplier is $1/(0.5\cdot 4^2)=0.125$, giving $(6-4)(0.125)=0.25$. Multiplying by the design rows and adding,

U(\beta) = \begin{pmatrix}1\\0\end{pmatrix}(-0.5) + \begin{pmatrix}1\\1\end{pmatrix}(0.25) = \begin{pmatrix}-0.25\\0.25\end{pmatrix}.

These are the two coefficient derivatives at the candidate value. In particular, the intercept derivative includes both observations, while the second coefficient derivative includes only the second row's contribution.

Worked Example: Poisson Regression With a Square-Root Link

Let $Y_i$ be Poisson with mean $\mu_i$, so $\phi=1$ and $V(\mu)=\mu$. Use $g(\mu)=\sqrt{\mu}$, with design rows $x_1^T=(1,0)$, $x_2^T=(1,1)$, coefficient candidate $\beta=(1,1)^T$, and responses $y_1=0$, $y_2=5$. Both fitted predictors are positive, as required for this link's mean range.

The predictors are $\eta_1=1$ and $\eta_2=2$, and the fitted means are $\mu_1=1$ and $\mu_2=4$. The derivative is $g'(\mu)=1/(2\sqrt{\mu})$. For observation 1, $V(\mu_1)g'(\mu_1)=1/2$, so the scaled residual is $(0-1)/(1/2)=-2$. For observation 2, $V(\mu_2)g'(\mu_2)=4(1/4)=1$, so the scaled residual is $(5-4)/1=1$. Therefore,

U(\beta) = \begin{pmatrix}1\\0\end{pmatrix}(-2) + \begin{pmatrix}1\\1\end{pmatrix}(1) = \begin{pmatrix}-1\\1\end{pmatrix}.

The example shows why the link derivative matters: the residuals are not simply added to form the coefficient score; each is scaled according to the mean and the link.

Worked Example: Gamma Regression With an Inverse Link

For a Gamma response with mean $\mu$ and variance $\phi\mu^2$, take $V(\mu)=\mu^2$. Use the inverse link $g(\mu)=1/\mu$, known dispersion $\phi=0.5$, and a single observation with $x_1=1$, $y_1=1$, and candidate coefficient $\beta=0.5$. The predictor is $\eta_1=0.5$, giving $\mu_1=2$.

Here $g'(\mu)=-1/\mu^2$, so $g'(2)=-1/4$. The denominator in the score multiplier is $0.5\cdot V(2)\cdot g'(2)=0.5\cdot4\cdot(-1/4)=-0.5$. The response residual is $y_1-\mu_1=1-2=-1$. Hence

U(\beta)=1\cdot\frac{-1}{-0.5}=2.

The score is positive at this coefficient value. The negative link derivative is not an error: it records that increasing the inverse-link predictor decreases the fitted mean.

What a Zero Score Does—and Does Not—Tell You

A likelihood maximizer that lies in the interior of an open coefficient parameter space must have zero score. This is a necessary condition, not a guarantee that a solution is a maximum. A zero score can also occur at a minimum or a saddle point, and a maximum at the boundary need not have zero score.

Result: If $\ell(\beta)$ is differentiable on an open subset of $\mathbb R^p$ and $\beta^\ast$ is an interior local maximum, then $U(\beta^\ast)=0$.

Proof. Let $a\in\mathbb R^p$ be any direction. For all sufficiently small positive and negative $t$, openness ensures that $\beta^\ast+ta$ remains in the parameter space. The one-variable function $h(t)=\ell(\beta^\ast+ta)$ has a local maximum at $t=0$. Differentiability therefore gives $h'(0)=0$. By the chain rule, $h'(0)=a^TU(\beta^\ast)$. Since this holds for every direction $a$, choose each coordinate unit vector in turn to conclude that every component of $U(\beta^\ast)$ is zero. $\square$

In practice, the score is used to identify candidate coefficient estimates and to assess how a likelihood changes locally. The formula also provides a useful diagnostic for a derivation or implementation: every observation contributes its design row multiplied by a scaled residual. A missing factor of $\phi$, $V(\mu_i)$, or $g'(\mu_i)$ changes the score and therefore changes the estimating equation.

Common Mistake

Dropping the link derivative. It is tempting to differentiate the log-likelihood with respect to the mean and then treat that derivative as the coefficient score. But coefficients change the mean through $\eta_i=x_i^T\beta=g(\mu_i)$. The chain rule contributes $d\mu_i/d\eta_i=1/g'(\mu_i)$. The correct multiplier is $1/\{\phi V(\mu_i)g'(\mu_i)\}$, not $1/\{\phi V(\mu_i)\}$.
Assuming variance and mean alone guarantee this likelihood score. The derivation used the exponential-dispersion likelihood. If only the first two moments are specified, an exact likelihood score cannot be inferred from those moments alone.
Treating a zero score as proof of a maximum. It is a necessary condition at an interior differentiable maximum. Additional analysis is needed to determine whether a stationary point is actually a maximum.

Key takeaway: Under the exponential-dispersion GLM likelihood, each response contributes its design row times the residual $y_i-\mu_i$, scaled by $1/\{\phi V(\mu_i)g'(\mu_i)\}$. This is the general regression score with dispersion held fixed.

Check Your Understanding

Try these questions using the assumptions and score formula above.

Which distributional assumption makes the displayed expression an exact likelihood score rather than a formula based only on the mean and variance?
For a link $g$, what derivative connects a change in the linear predictor to a change in the mean?
Why does the score have $p$ components when the design matrix has $p$ columns?
In the Gamma inverse-link example, why is the score multiplier's denominator negative?
Why is a zero regression score necessary at an interior local maximum, but not sufficient to establish one?