Worked Score Derivations Across Three GLM Families

From Likelihood Derivatives to Fitted Residuals

A coefficient score is the derivative of the log-likelihood with respect to a regression coefficient. In a GLM, differentiating first with respect to the linear predictor and then applying the chain rule gives a direct connection between the score and the difference between each observed response and its fitted mean. This tutorial works out that connection for three common models: Gaussian identity, Bernoulli logit, and Poisson log-link.

Assume independent observations \(Y_1,\ldots,Y_n\), fixed covariate vectors \(x_i\in\mathbb{R}^p\), and linear predictors \(\eta_i=x_i^\mathsf{T}\beta\), where \(\beta\in\mathbb{R}^p\). Write \(\mu_i=E[Y_i]\) for the model mean and \(\hat\mu_i\) for its fitted value at an estimated coefficient vector \(\hat\beta\). Define the response residual \(r_i=y_i-\mu_i\) at a given coefficient vector; at the fit this becomes \(y_i-\hat\mu_i\). A score equation sets each component of the likelihood gradient equal to zero.

The calculations below use the usual maximum-likelihood convention: the Gaussian variance \(\sigma^2>0\) is common across observations, and the Bernoulli and Poisson observations have the stated distributions. The examples have finite fitted coefficient values, so the likelihood derivatives and score equations apply at their fits.

Gaussian Identity-Link Model

For independent \(Y_i\sim N(x_i^\mathsf{T}\beta,\sigma^2)\), the identity link gives \(\mu_i=\eta_i=x_i^\mathsf{T}\beta\). Treating \(\sigma^2\) as fixed, the log-likelihood, up to no omitted \(\beta\)-dependent terms, is

For coefficient \(\beta_j\), the derivative of the squared-error term is \(\partial (y_i-x_i^\mathsf{T}\beta)^2/\partial\beta_j=-2(y_i-x_i^\mathsf{T}\beta)x_{ij}\). Therefore

Thus the fitted Gaussian coefficient equations require each covariate-weighted residual sum to be zero. The common factor \(1/\sigma^2\) does not affect their roots. Estimating a common variance as well does not change this coefficient score.

Bernoulli Logit Model

For independent \(Y_i\sim\operatorname{Bernoulli}(p_i)\), the logit link specifies \(\log\{p_i/(1-p_i)\}=\eta_i=x_i^\mathsf{T}\beta\). Equivalently, \(p_i=e^{\eta_i}/(1+e^{\eta_i})\). The log-likelihood is

To differentiate, note that \(dp_i/d\eta_i=p_i(1-p_i)\). The derivative of an observation's log-likelihood with respect to \(\eta_i\) is

Since \(\partial\eta_i/\partial\beta_j=x_{ij}\), the chain rule yields \(U_j(\beta)=\sum_i x_{ij}(y_i-p_i)\). In vector form, \(U(\beta)=X^\mathsf{T}(y-p)\). Here the response residual is \(y_i-p_i\), not a difference between \(y_i\) and the linear predictor.

Poisson Log-Link Model

For independent \(Y_i\sim\operatorname{Poisson}(\mu_i)\), the log link gives \(\log\mu_i=\eta_i=x_i^\mathsf{T}\beta\), so \(\mu_i=e^{\eta_i}\). The log-likelihood is

The factorial term does not depend on \(\beta\). Differentiating the remaining terms first with respect to \(\eta_i\) gives \(y_i-e^{\eta_i}=y_i-\mu_i\); applying the chain rule gives

As in the logit model, each coefficient score is a covariate-weighted sum of response residuals, but the fitted mean is now \(e^{x_i^\mathsf{T}\beta}\).

What the Residual Checks Establish

The three models have different likelihoods and fitted mean functions, but their coefficient scores share the same residual structure. The Gaussian score has the additional common factor \(1/\sigma^2\); for Bernoulli logit and Poisson log-link, the coefficient score is exactly \(X^\mathsf{T}(y-\mu)\). The following statement makes the consequence of a zero score explicit.

Proof. In the Gaussian model the score is \(\sigma^{-2}X^\mathsf{T}(y-\mu)\), and \(\sigma^2>0\), so a zero score implies \(X^\mathsf{T}(y-\mu)=0\). In the Bernoulli logit and Poisson log-link models, the score is already \(X^\mathsf{T}(y-\mu)\), so the same implication holds directly. If the model includes an intercept, one column of \(X\) consists of ones. The component of \(X^\mathsf{T}(y-\mu)\) corresponding to that column is \(\sum_i(y_i-\mu_i)\); because every component is zero, this sum is zero. At a fitted solution, \(\mu_i=\hat\mu_i\), giving the stated fitted-residual checks. \(\square\)

These equations express a covariate-weighted balance: positive and negative residuals must cancel after weighting by each model covariate. They do not require every individual residual to be zero. Nor do they, by themselves, prove that a zero-score point is a maximum or that it is the only solution. The existence and finiteness of solutions are separate questions.

Common Mistake: Checking the Wrong Quantity

A frequent error is to set a score equal to zero using \(y_i-\eta_i\), as if the linear predictor were the fitted response. The score depends on the likelihood through the fitted mean: it uses \(y_i-\eta_i\) only when the particular mean function makes \(\mu_i=\eta_i\), as in the Gaussian identity-link model. For Bernoulli logit, the fitted mean is \(p_i=1/(1+e^{-\eta_i})\); for Poisson log-link, it is \(\mu_i=e^{\eta_i}\). Therefore the correct residuals are \(y_i-p_i\) and \(y_i-e^{\eta_i}\), respectively.

A second related mistake is to check only that the residuals sum to zero. That verifies the intercept equation when an intercept is present, but not the other coefficient equations. Each additional covariate requires its own weighted residual sum to vanish. A reliable check computes the fitted means first, forms the response residuals, and then evaluates every column of \(X^\mathsf{T}(y-\hat\mu)\).