Deriving The GLM Score Equation

Why the GLM Score Depends on the Link

A generalized linear model connects a response distribution to predictors through a linear predictor and a link function. To estimate its regression coefficients by likelihood, we differentiate the log likelihood with respect to those coefficients. The resulting score vector measures how the log likelihood changes locally as each coefficient changes. Setting it equal to zero gives the likelihood equations for an interior maximum-likelihood estimate.

The calculation has a useful general form. The response distribution determines a variance function, while the link determines how the fitted mean changes with the linear predictor. For a canonical link, those pieces cancel in the score. For a noncanonical link, they do not, so using the simpler canonical-link expression would give the wrong result.

Definition: Suppose independent responses \(Y_1,\ldots,Y_n\), conditional on predictor vectors \(\mathbf{x}_i\), have an exponential-family density or probability mass function \(f(y_i;\theta_i,\phi)=\exp\{[y_i\theta_i-b(\theta_i)]/a(\phi)+c(y_i,\phi)\}\). Here \(\theta_i\) is the natural parameter, \(\phi\) is a common dispersion parameter, and \(a(\phi)\) is the dispersion factor. The conditional mean is \(\mu_i=b'(\theta_i)\), and the variance is \(a(\phi)b''(\theta_i)\). Write \(V(\mu_i)=b''(\theta_i)\), so the variance is \(a(\phi)V(\mu_i)\). A GLM specifies \(g(\mu_i)=\eta_i=\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta}\), where \(g\) is the link, \(\eta_i\) is the linear predictor, and \(\boldsymbol{\beta}\) is the coefficient vector.

The score is the gradient of the log likelihood with respect to \(\boldsymbol{\beta}\). It is a vector with one component per coefficient. A score evaluated at a particular coefficient vector need not be zero; the score equation is the equation that sets the vector to zero.

Deriving the General Score Equation

Independence means the joint likelihood is the product of the individual likelihoods, so its log is their sum. Terms \(c(y_i,\phi)\) do not depend on the regression coefficients when \(\phi\) is fixed. Differentiating the remaining terms first with respect to \(\theta_i\) gives \((y_i-\mu_i)/a(\phi)\). The chain rule then accounts for how \(\theta_i\) changes with the linear predictor and how the linear predictor changes with the coefficients.

Write the log likelihood.
For independent observations, sum their log densities or log mass functions. The part depending on the natural parameters is \(\sum_i [y_i\theta_i-b(\theta_i)]/a(\phi)\).

Differentiate with respect to the natural parameter.
Since \(b'(\theta_i)=\mu_i\), the derivative of \(y_i\theta_i-b(\theta_i)\) is \(y_i-\mu_i\).

Apply the chain rule.
Use \(d\theta_i/d\mu_i=1/V(\mu_i)\) and \(d\mu_i/d\eta_i=1/g'(\mu_i)\). Also, \(d\eta_i/d\boldsymbol{\beta}=\mathbf{x}_i\).

\begin{aligned} \ell(\boldsymbol{\beta})&=\sum_{i=1}^n \left\{\frac{y_i\theta_i-b(\theta_i)}{a(\phi)}+c(y_i,\phi)\right\},\\ \mathbf{U}(\boldsymbol{\beta}) =\frac{\partial \ell}{\partial \boldsymbol{\beta}} &=\sum_{i=1}^n \mathbf{x}_i\, \frac{y_i-\mu_i}{a(\phi)V(\mu_i)g'(\mu_i)}. \end{aligned}

This formula requires differentiable \(b\) and \(g\), and a valid model at the parameter values under consideration. In particular, the inverse link must give admissible means, and \(g'(\mu_i)\) must be defined and nonzero for the displayed derivative to apply. The formula uses a common dispersion factor; it does not cover likelihoods with observation-specific dispersion factors.

Canonical-link result: For a canonical link, \(\eta_i=\theta_i\), so \(g'(\mu_i)=d\theta_i/d\mu_i=1/V(\mu_i)\). Therefore \(\mathbf{U}(\boldsymbol{\beta})=\sum_i \mathbf{x}_i(y_i-\mu_i)/a(\phi)\). This simplification does not apply to an arbitrary link.

The cancellation in the canonical case follows directly from the chain rule: \(V(\mu_i)g'(\mu_i)=V(\mu_i)/V(\mu_i)=1\). The score is therefore built from residuals \(y_i-\mu_i\), multiplied by the predictor vectors. With a noncanonical link, the factor \(1/[V(\mu_i)g'(\mu_i)]\) remains and adjusts each observation's contribution.

Worked Score Calculations

Worked Example: Poisson Regression With a Log Link

A hypothetical insurer models the number of service requests \(Y_i\) from three policy groups using a Poisson GLM. The predictor values are \(x=(0,1,2)\), and the observed counts are \(y=(1,4,2)\). Use the log link, \(\log(\mu_i)=\beta_0+\beta_1x_i\), and evaluate the score at \((\beta_0,\beta_1)=(\log 2,\log 2)\). The Poisson distribution has \(a(\phi)=1\), \(\theta_i=\log\mu_i\), and \(V(\mu_i)=\mu_i\). Find both score components.

Method. The log link is canonical for the Poisson model, so the score components are \(\sum_i(y_i-\mu_i)\) and \(\sum_i x_i(y_i-\mu_i)\).

Step 1: Calculate fitted means. The linear predictors are \(\eta_1=\log 2\), \(\eta_2=2\log 2=\log 4\), and \(\eta_3=3\log 2=\log 8\). Exponentiating gives \(\mu=(2,4,8)\).

Step 2: Calculate residuals. The residuals are \(y_i-\mu_i=(-1,0,-6)\).

Step 3: Calculate the score. The intercept score is \(-1+0-6=-7\). The slope score is \(0(-1)+1(0)+2(-6)=-12\). Thus \(\mathbf{U}=(-7,-12)^\mathsf{T}\).

Check and trap. Direct differentiation of \(\sum_i[y_i\eta_i-\exp(\eta_i)-\log(y_i!)]\) gives the same sums. The score is not zero at the proposed coefficients, so this point does not satisfy the likelihood equations. Do not report the evaluated score as though it were an estimate of the coefficients.

Worked Example: Bernoulli Regression With a Logit Link

A hypothetical study records whether each of three applicants completes a training module. Let \(Y_i=1\) for completion and \(Y_i=0\) otherwise. The predictor values are \(x=(0,1,2)\), and the outcomes are \(y=(1,0,1)\). Use \(\log[p_i/(1-p_i)]=\beta_0+\beta_1x_i\) and evaluate the score at \((0,\log 3)\). Find the score vector.

Method. A Bernoulli response has \(a(\phi)=1\), \(V(p)=p(1-p)\), and the logit is its canonical link. Hence the score is the predictor-weighted residual sum.

Step 1: Calculate fitted probabilities. The linear predictors are \(0,\log 3,2\log 3\), giving odds \(1,3,9\), respectively. Therefore \(p=(1/2,3/4,9/10)\).

Step 2: Calculate residuals. The residuals \(y_i-p_i\) are \(1/2,-3/4,1/10\).

Step 3: Calculate both components. The intercept score is \(1/2-3/4+1/10=(10-15+2)/20=-3/20=-0.15\). The slope score is \(0(1/2)+1(-3/4)+2(1/10)=-3/4+1/5=(-15+4)/20=-11/20=-0.55\). Thus \(\mathbf{U}=(-0.15,-0.55)^\mathsf{T}\).

Check and trap. The second component weights each residual by its predictor value; the observation with \(x=0\) contributes nothing to the slope score. The probability residuals here are \(y_i-p_i\), not residuals formed by treating the response as a count other than zero or one.

Worked Example: Gamma Regression With a Log Link

A hypothetical actuary models positive claim-processing times \(Y_i\) with a Gamma GLM. For this parameterization, the conditional mean is \(\mu_i\), the shape is \(1/\phi\), and the scale is \(\phi\mu_i\). Let \(\phi=0.5\), \(x=(0,1)\), and \(y=(3,2)\). Use \(\log(\mu_i)=\beta_0+\beta_1x_i\) and evaluate the score at \((\log 2,\log 2)\).

Method. For the Gamma exponential-family model, \(V(\mu)=\mu^2\). The log link has \(g'(\mu)=1/\mu\), so the general score contribution is \(\mathbf{x}_i(y_i-\mu_i)/[\phi\mu_i]\). This is a noncanonical link for the Gamma family.

Step 1: Calculate means and residuals. The linear predictors are \(\log 2\) and \(2\log 2=\log 4\), so \(\mu=(2,4)\). The residuals are \(3-2=1\) and \(2-4=-2\).

Step 2: Calculate the adjusted contributions. For the first observation, \((y_1-\mu_1)/(\phi\mu_1)=1/(0.5\cdot2)=1\). For the second, \((y_2-\mu_2)/(\phi\mu_2)=-2/(0.5\cdot4)=-1\).

Step 3: Sum for each coefficient. The intercept score is \(1+(-1)=0\). The slope score is \(0(1)+1(-1)=-1\). Therefore \(\mathbf{U}=(0,-1)^\mathsf{T}\).

Check and trap. Substituting the Gamma variance function and log-link derivative into the general formula gives the same adjusted contributions. Using the canonical-link expression would instead give the unadjusted sums \(1+(-2)=-1\) and \(0(1)+1(-2)=-2\). That is the wrong score because the log link is not canonical for this Gamma model.

Using the Score Equation Correctly

For a coefficient vector that is an interior maximum-likelihood estimate, the score must satisfy \(\mathbf{U}(\boldsymbol{\beta})=\mathbf{0}\). This gives one equation for each coefficient. In a canonical-link model, those equations set predictor-weighted residual sums to zero (up to the common dispersion factor). They do not generally say that every residual is zero. Nor do they say that an unweighted residual sum is zero for every model: that particular conclusion requires an intercept represented by a constant predictor and the canonical score form.

For a noncanonical link, retain the variance and link-derivative factor for every observation. The score equation is then generally nonlinear in \(\boldsymbol{\beta}\), because the fitted means in that factor depend on the coefficients. A numerical method may be needed to solve it. A zero score is a necessary condition for an interior maximum under regularity conditions, but a stationary point should not be assumed to be a maximum without checking the likelihood and the model's conditions.

State the distribution and link first. They determine \(V(\mu)\), \(g'(\mu)\), and whether the link is canonical.
Keep the predictor vector attached to each residual. Each score component weights the observation's contribution by the corresponding predictor value.
Separate general and canonical formulas. The simple residual form is justified only when the link is canonical (or when an equivalent cancellation has been established).
Distinguish evaluating from solving. Computing \(\mathbf{U}\) at proposed coefficients is not the same as finding coefficients that make \(\mathbf{U}=\mathbf{0}\).

Key takeaway: The general GLM score contribution is the residual \(y_i-\mu_i\), multiplied by the predictor vector and adjusted by \(1/[a(\phi)V(\mu_i)g'(\mu_i)]\). For a canonical link, the variance and link derivative cancel; for a noncanonical link, keep the adjustment.

Check Your Understanding

For each problem, calculate the requested score using the stated distribution and link. Score vectors are ordered as intercept, then slope where applicable.

A Poisson log-link model has \(x=(0,1)\), \(y=(2,1)\), and proposed coefficients \((0,\log 2)\). Find the score vector. (Answer: \((0,-1)^\mathsf{T}\))
A Bernoulli logit model has \(x=(0,1,2)\), \(y=(0,1,0)\), and proposed coefficients \((0,0)\). Find the score vector. (Answer: \((-0.5,-0.5)^\mathsf{T}\))
A Gamma model with log link has \(\phi=0.25\), \(x=(0,1)\), \(y=(1,6)\), and proposed coefficients \((\log 2,\log 2)\). Find the score vector. (Answer: \((0,2)^\mathsf{T}\))
A Poisson log-link model has \(x=(0,1,2)\), \(y=(0,3,1)\), and proposed coefficients \((0,\log 2)\). Find the score vector. (Answer: \((-3,-5)^\mathsf{T}\))
A Gamma model with identity link has \(\phi=0.5\), observations \((x,y)=(0,3),(1,4)\), and proposed coefficients \((2,1)\), so the fitted means are positive. Find the score vector. (Answer: \((13/18,2/9)^\mathsf{T}\))

What You'll Learn