Score Equations for GLM Coefficients

Why GLM Score Equations Matter

A generalized linear model (GLM) estimates coefficients by choosing values that make the observed data most plausible under the specified model. The log-likelihood summarizes that plausibility; differentiating it with respect to each coefficient gives a score equation. A maximum that lies inside the parameter space must have zero score, so solving the system of score equations is a central step in maximum-likelihood estimation.

The score equations also show how the model connects observed responses to fitted means. For canonical links, they take a particularly direct form: predictor-weighted residuals sum to zero. For other links, a derivative factor enters, so the equations generally are not just ordinary residual balances.

Definition: For a log-likelihood \(\ell(\boldsymbol{\beta})\), the score for coefficient \(\beta_j\) is \(U_j(\boldsymbol{\beta})=\partial\ell/\partial\beta_j\). The score vector is the collection of these derivatives. A score root is a coefficient vector \(\widehat{\boldsymbol{\beta}}\) satisfying \(U_j(\widehat{\boldsymbol{\beta}})=0\) for every coefficient.

From the Likelihood to the Score

Assume \(n\) responses \(Y_1,\ldots,Y_n\) are independent and each has an exponential-family density or probability 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)\) does not depend on the regression coefficients. The mean is \(\mu_i=b'(\theta_i)\). A GLM specifies \(g(\mu_i)=\eta_i=\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta}\), where \(\mathbf{x}_i\) is the vector of predictor values, including a 1 for the intercept.

The log-likelihood is the sum of the individual log-likelihood contributions because the responses are independent. Terms \(c(y_i,\phi)\) do not depend on \(\boldsymbol{\beta}\), so they vanish when differentiated with respect to a coefficient. Applying the chain rule to the remaining terms gives the score.

\ell(\boldsymbol{\beta}) =\sum_{i=1}^n\left\{ \frac{y_i\theta_i-b(\theta_i)}{a(\phi)}+c(y_i,\phi) \right\}

\frac{\partial\ell}{\partial\beta_j} =\sum_{i=1}^n \frac{y_i-b'(\theta_i)}{a(\phi)} \frac{\partial\theta_i}{\partial\eta_i} \frac{\partial\eta_i}{\partial\beta_j} =\sum_{i=1}^n \frac{y_i-\mu_i}{a(\phi)} \frac{d\theta_i}{d\eta_i}x_{ij}.

This expression is the general coefficient score under the stated model. It assumes a differentiable link-to-parameter relationship and an interior coefficient value where the likelihood can be differentiated. The maximum-likelihood equations set this derivative equal to zero for each \(j\). In vector form, the same result is \(\sum_i \mathbf{x}_i(y_i-\mu_i)(d\theta_i/d\eta_i)/a(\phi)=\mathbf{0}\).

Canonical-link result: With a canonical link, \(\eta_i=\theta_i\), so \(d\theta_i/d\eta_i=1\). If \(a(\phi)\) is common to all observations, the score equations reduce to \(\sum_i x_{ij}(y_i-\mu_i)=0\) for every coefficient \(j\). Equivalently, \(\mathbf{X}^{\mathsf T}(\mathbf{y}-\boldsymbol{\mu})=\mathbf{0}\). This simplification does not generally hold for a noncanonical link.

Why the Score Must Be Zero at an Interior Maximum

The score equations follow from a basic condition for a differentiable function to have a local maximum at an interior point. If a coefficient can be increased or decreased slightly, then at the maximum neither small change can improve the log-likelihood. The derivative in each coordinate direction must therefore be zero. This explains why a finite interior maximum-likelihood estimate is a root of the score vector.

The condition is necessary, not automatically sufficient. A score root could be a minimum or another stationary point, and some models have no finite root even though the likelihood can be made progressively larger toward a boundary. After solving the equations, check that the solution is in the parameter space and is a maximum. In simple examples, the likelihood’s curvature or direct likelihood comparison can provide that check.

Write the model and likelihood.
State the response distribution, mean structure, link, and any fixed parameters before differentiating.

Differentiate by coefficient.
Use the chain rule through \(\theta_i\), \(\mu_i\), and \(\eta_i\); retain the link derivative unless the link is canonical.

Set each score component to zero.
Solve the resulting system for the coefficient estimates, while checking that the solution is finite and allowed by the model.

Verify the candidate.
Substitute the estimates into the score equations and confirm that the candidate corresponds to a maximum, not merely a stationary point.

Worked Score-Equation Applications

Worked Example: Poisson Intercept-Only Frequency Model

A hypothetical portfolio contains three equal-exposure policy records with claim counts \(y_1=0\), \(y_2=1\), and \(y_3=3\). Assume the counts are independent Poisson random variables with common mean \(\lambda=e^\beta\). Find the maximum-likelihood estimate of \(\lambda\) and \(\beta\), using the score equation, and verify the result.

Model and likelihood. For a Poisson response, \(P(Y_i=y_i)=e^{-\lambda}\lambda^{y_i}/y_i!\). Since \(\log\lambda=\beta\), the log-likelihood is

\ell(\beta)=\sum_{i=1}^3(y_i\beta-e^\beta-\log(y_i!)) =4\beta-3e^\beta-\log 6.

Here \(\sum_i y_i=0+1+3=4\), and \(0!1!3!=6\). The Poisson log link is canonical, so the coefficient score is the sum of observed minus fitted counts:

U(\beta)=\frac{d\ell}{d\beta}=4-3e^\beta.

Set the score to zero. Thus \(4-3e^\beta=0\), giving \(e^{\widehat{\beta}}=4/3\). Therefore \(\widehat{\lambda}=4/3\) and \(\widehat{\beta}=\log(4/3)\approx0.2877\).

Check. The fitted mean total is \(3(4/3)=4\), equal to the observed count total. Substitution gives \(U(\widehat{\beta})=4-3(4/3)=0\). Also, \(\ell''(\beta)=-3e^\beta<0\) for every finite \(\beta\), so the log-likelihood is strictly concave and this root is its maximum.

Method and mistake avoided. The method is direct differentiation of the Poisson log-likelihood. A common mistake is to set \(\sum_i(y_i-\beta)=0\); the fitted mean is \(e^\beta\), not \(\beta\), so the correct residual is \(y_i-e^\beta\).

Worked Example: Binomial Logistic Model With a Binary Predictor

A hypothetical insurer groups eight independent risks into two classes. In class 0, one of four risks has a claim; in class 1, three of four have a claim. Let each individual response be Bernoulli with probability \(p_i\), and use the logistic model \(\log\{p_i/(1-p_i)\}=\beta_0+\beta_1x_i\), where \(x_i=0\) or \(1\) identifies the class. Find both coefficients and verify both score equations.

Likelihood and scores. For a Bernoulli response, the log-likelihood, omitting no coefficient-dependent terms, is \(\ell=\sum_i[y_i\log p_i+(1-y_i)\log(1-p_i)]\). The logistic link is canonical for the binomial model, so the score equations are

U_0=\sum_i(y_i-p_i)=0,\qquad U_1=\sum_i x_i(y_i-p_i)=0.

The slope equation includes only class 1. Its observed success total is 3, so it requires the fitted success total in that class to equal 3. With four risks, \(4p_1=3\), hence \(p_1=3/4\). The intercept equation then requires the fitted success total across both classes to equal the observed total \(1+3=4\). Since class 1 contributes 3, class 0 must contribute 1, giving \(4p_0=1\) and \(p_0=1/4\).

Convert fitted probabilities to coefficients. The intercept is the class 0 log-odds: \(\widehat{\beta}_0=\log\{(1/4)/(3/4)\}=-\log 3\approx-1.0986\). The class 1 log-odds are \(\log\{(3/4)/(1/4)\}=\log 3\), so \(\widehat{\beta}_1=\log 3-(-\log 3)=2\log 3\approx2.1972\).

Check. The fitted probabilities are \(p_0=1/(1+e^{\log 3})=1/4\) and \(p_1=1/(1+e^{-\log 3})=3/4\). Thus \(U_0=(1-1)+(3-3)=0\) when written as observed minus fitted class totals, and \(U_1=3-3=0\). The fitted values match the empirical proportions in each class.

Method and mistake avoided. The method uses the two canonical-link score equations, first for the class-specific total and then for the overall total. Do not interpret the intercept as the overall claim probability: with this coding, it is the log-odds for class 0.

Worked Example: Normal Identity-Link Regression

For a hypothetical severity response, suppose \(Y_i\) are independent normal random variables with \(Y_i\sim N(\mu_i,1)\), where the second parameter is the variance, and \(\mu_i=\beta_0+\beta_1x_i\). The three observed values are \((x_i,y_i)=(0,2),(1,4),(2,5)\). Find the coefficient estimates from the score equations and verify the residual balances.

Likelihood and equations. The variance is known to be \(1\), and the normal identity link is canonical. Ignoring constants independent of the coefficients, the log-likelihood is

\ell(\beta_0,\beta_1) =-\frac{1}{2}\sum_{i=1}^3(y_i-\beta_0-\beta_1x_i)^2.

Differentiating with respect to \(\beta_0\) and \(\beta_1\) gives

\sum_i(y_i-\beta_0-\beta_1x_i)=0,\qquad \sum_i x_i(y_i-\beta_0-\beta_1x_i)=0.

The first equation becomes \(11-3\beta_0-3\beta_1=0\), because \(\sum y_i=11\), \(\sum 1=3\), and \(\sum x_i=3\). The second becomes \(14-3\beta_0-5\beta_1=0\), because \(\sum x_iy_i=0(2)+1(4)+2(5)=14\) and \(\sum x_i^2=0+1+4=5\). Subtracting the first equation from the second gives \(3-2\beta_1=0\), so \(\widehat{\beta}_1=3/2\). Substitution into the first gives \(3\beta_0=11-9/2=13/2\), hence \(\widehat{\beta}_0=13/6\).

Check. The fitted means at \(x=0,1,2\) are \(13/6,11/3,31/6\). The residuals are \(-1/6,1/3,-1/6\). Their sum is \((-1+2-1)/6=0\); their predictor-weighted sum is \(0(-1/6)+1(1/3)+2(-1/6)=0\). Both score equations hold. The log-likelihood is a negative sum of squared residuals plus a constant, so this solution minimizes that sum and maximizes the likelihood.

Method and mistake avoided. The method solves the intercept and slope score equations together. A common error is to check only that residuals sum to zero; the slope score also requires their predictor-weighted sum to be zero.

Noncanonical Links and Limits of the Score Rule

For a noncanonical link, the score still follows from the same likelihood, but the factor \(d\theta_i/d\eta_i\) remains. It can change the equations’ weighting of observations. Therefore, it is not valid to replace every GLM score by a simple sum of predictor times raw residual. The canonical-link simplification is a consequence of \(\theta_i=\eta_i\), not a general property of maximum likelihood.

A score root also need not exist at finite coefficient values. For example, in logistic regression, if a predictor perfectly separates all successes from all failures, coefficient values can move toward infinity while fitted probabilities approach the observed outcomes. There is then no finite coefficient vector at which all score components equal zero. More generally, a rank-deficient design can prevent unique coefficient estimates. These cases are reasons to check existence and identifiability rather than assuming every formal score system has a unique finite solution.

Key takeaway: Derive the score from the specified likelihood, retain the chain-rule factor for a noncanonical link, and solve one equation per coefficient. A finite interior maximum-likelihood estimate must be a score root, but a root must still be checked as a valid maximum.

Check Your Understanding

For each problem, identify the relevant score equation and calculate or characterize its solution.

Five independent Poisson counts are \(1,0,2,3,4\), all with mean \(e^\beta\). Find the score root and fitted mean. (Answer: \(\widehat{\lambda}=2\), \(\widehat{\beta}=\log 2\).)
Twelve independent Bernoulli observations share a logistic intercept only, and nine are successes. Find the fitted probability and intercept. (Answer: \(\widehat p=3/4\), \(\widehat{\beta}_0=\log 3\).)
Independent Poisson observations have log means \(\beta_0+\beta_1x\). At \(x=0\), the observed counts are 0 and 2; at \(x=1\), they are 3 and 5. Find the two score roots. (Answer: \(\widehat{\beta}_0=0\), \(\widehat{\beta}_1=\log 4\).)
For independent \(Y_i\sim N(\mu_i,1)\) with \(\mu_i=\beta_0+\beta_1x_i\), the observations are \((0,1),(1,2),(2,4)\). Solve the intercept and slope score equations. (Answer: \(\widehat{\beta}_0=5/6\), \(\widehat{\beta}_1=3/2\).)
In a logistic model with a binary predictor, suppose every observation with \(x=0\) is a failure and every observation with \(x=1\) is a success. Does a finite maximum-likelihood coefficient vector satisfy the score equations? (Answer: No; the data are perfectly separated, and the likelihood approaches its supremum as the intercept tends to negative infinity and the slope tends to positive infinity.)

What You'll Learn