Checking an IRLS Step Against the Score

Why Check an IRLS Step Against the Score?

An IRLS calculation converts a likelihood-based update into a weighted least-squares calculation. That conversion is useful computationally, but it also creates opportunities for small algebraic errors: a derivative can be inverted, a dispersion factor can be misplaced, or an offset can be omitted. A direct check is available: at the current coefficients, the weighted residual of the working response must reproduce the GLM score.

The comparison is made at one fixed coefficient vector. The working response and weights are computed from its fitted means; the residual is measured against the current linear predictor, including any offset. The resulting vector identity is exact, not an approximation.

Setup and Definitions

Consider independent observations \(Y_1,\ldots,Y_n\) in a GLM. Let \(x_i^T\) be row \(i\) of the design matrix \(X\), let \(\beta\) be the coefficient vector, and let \(o_i\) be a fixed offset. At the current coefficient vector \(\beta\), define \(\eta_i=x_i^T\beta+o_i\) and \(\mu_i=g^{-1}(\eta_i)\). Assume the mean is in the interior of its parameter space, the inverse link is differentiable there, and \(d_i=d\mu_i/d\eta_i\) is finite and nonzero. Let the variance function satisfy \(V(\mu_i)>0\), let the dispersion \(\phi>0\), and let any prior weight \(a_i\) be positive. The unweighted case has \(a_i=1\).

For precision, suppose the log-likelihood contribution has the exponential-dispersion form \(\ell_i=a_i\{y_i\theta_i-b(\theta_i)\}/\phi+c_i(y_i,\phi)\), where \(\mu_i=b'(\theta_i)\), \(V(\mu_i)=b''(\theta_i)\), and \(c_i\) does not depend on \(\beta\). This includes the usual likelihood-based GLM score calculations. A quasi-score with the same mean-variance relation gives the same algebra below.

Definition: At the current coefficients, the working response and IRLS weight for observation \(i\) are \(z_i=\eta_i+(y_i-\mu_i)/d_i\) and \(W_i=a_i d_i^2/\{\phi V(\mu_i)\}\). The weighted working residual relative to the current fit is \(W_i(z_i-\eta_i)\). In vector notation, \(W=\operatorname{diag}(W_1,\ldots,W_n)\), \(z=(z_1,\ldots,z_n)^T\), and \(\eta=X\beta+o\).

The residual here is \(z_i-\eta_i\), not the working response itself. Also, it is multiplied by \(W_i\) when forming the score contribution. In weighted least squares, the square-root weight multiplies a residual when forming a transformed data vector, but the normal equations use the weight \(W_i\) itself.

The Score Identity

Theorem: Under the conditions above, the score vector \(U(\beta)=\partial\ell(\beta)/\partial\beta\) satisfies \(U(\beta)=X^T W(z-\eta)\). Observation by observation, \(W_i(z_i-\eta_i)=a_i d_i(y_i-\mu_i)/\{\phi V(\mu_i)\}\).

Proof. First differentiate the log-likelihood contribution with respect to the mean parameter. Since \(d\mu_i/d\theta_i=b''(\theta_i)=V(\mu_i)\), the chain rule gives

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

Because \(d\mu_i/d\beta=d_i x_i\), another application of the chain rule yields the \(i\)th contribution to the coefficient score:

U_i(\beta) = \frac{a_i(y_i-\mu_i)}{\phi V(\mu_i)}d_i x_i.

Now substitute the definitions of \(z_i\) and \(W_i\). The offset is part of \(\eta_i\), so it cancels when forming the working residual:

z_i-\eta_i=\frac{y_i-\mu_i}{d_i}, \qquad W_i(z_i-\eta_i) = \frac{a_i d_i^2}{\phi V(\mu_i)} \frac{y_i-\mu_i}{d_i} = \frac{a_i d_i(y_i-\mu_i)}{\phi V(\mu_i)}.

Thus \(U_i(\beta)=x_i W_i(z_i-\eta_i)\). Summing over observations gives \(U(\beta)=\sum_i x_i W_i(z_i-\eta_i)=X^TW(z-\eta)\), as claimed. \(\square\)

This identity is a practical consistency check: once \(z_i\) and \(W_i\) have been computed from the current fit, their weighted residual must match the score contribution before any least-squares system is solved.

Worked Example: Bernoulli Logistic Regression

Take two Bernoulli observations with outcomes \(y_1=1\), \(y_2=0\), and design rows \(x_1^T=(1,0)\), \(x_2^T=(1,1)\). Use the logit link, no offsets, \(\phi=1\), and current coefficients \(\beta=(0,0)^T\). Then both fitted probabilities are \(p_i=1/2\), and both linear predictors are \(\eta_i=0\).

For a Bernoulli response, \(V(\mu)=\mu(1-\mu)\). Under the logit link, \(d_i=d\mu_i/d\eta_i=p_i(1-p_i)=1/4\). Therefore \(W_i=d_i^2/V(\mu_i)=(1/16)/(1/4)=1/4\) for both observations. The working residuals are \(z_1-\eta_1=(1-1/2)/(1/4)=2\) and \(z_2-\eta_2=(0-1/2)/(1/4)=-2\).

The weighted residuals are consequently \(1/2\) and \(-1/2\). Multiplying by the design rows gives \(X^TW(z-\eta)=(1/2-1/2,\; -1/2)^T=(0,-1/2)^T\).

Directly, the logistic score is \(\sum_i x_i(y_i-p_i)\). Its intercept component is \((1-1/2)+(0-1/2)=0\), and its slope component is \(0(1-1/2)+1(0-1/2)=-1/2\). The two calculations agree.

Worked Example: Poisson Log Link With an Offset

Consider two Poisson observations with a common coefficient whose design value is \(x_i=1\). Let the offsets be \(o_1=\log 2\) and \(o_2=0\), the outcomes be \(y_1=1\) and \(y_2=3\), and the current coefficient be \(\beta=0\). With the log link, the means are \(\mu_1=2\), \(\mu_2=1\), and the linear predictors are \(\eta_1=\log 2\), \(\eta_2=0\).

For Poisson data, \(V(\mu)=\mu\), and for the log link \(d\mu/d\eta=\mu\). With \(\phi=a_i=1\), the weights are \(W_i=\mu_i^2/\mu_i=\mu_i\), so \(W_1=2\) and \(W_2=1\). The working residuals are \(z_1-\eta_1=(1-2)/2=-1/2\) and \(z_2-\eta_2=(3-1)/1=2\).

Multiplying by the weights gives \(-1\) and \(2\); since each design value is 1, their sum is the coefficient score \(1\). Directly, the Poisson log-link score is \(\sum_i x_i(y_i-\mu_i)=(1-2)+(3-1)=1\).

The offset affects the fitted means and the values of \(\eta_i\), but it does not get added as a coefficient column in the score. Using \(z_i-\eta_i\) correctly includes the offset in the working-regression residual.

Worked Example: Grouped Binomial Data

Suppose two groups have trial counts \(m_1=4\), \(m_2=2\), and observed success counts \(y_1=3\), \(y_2=0\). Use design rows \((1,0)\) and \((1,1)\), a logit link, no offsets, and current coefficients \((0,0)^T\). Each fitted success probability is \(1/2\), so the fitted mean counts are \(\mu_1=2\) and \(\mu_2=1\).

For success counts, \(V(\mu_i)=m_i p_i(1-p_i)\), and \(d_i=d\mu_i/d\eta_i=m_i p_i(1-p_i)\). Thus \(W_i=m_i p_i(1-p_i)\), giving \(W_1=1\) and \(W_2=1/2\). The working residuals are \(z_1-\eta_1=(3-2)/1=1\) and \(z_2-\eta_2=(0-1)/(1/2)=-2\). Their weighted values are \(1\) and \(-1\).

Therefore \(X^TW(z-\eta)=(1-1,\; -1)^T=(0,-1)^T\). Directly, the grouped-binomial score is \(\sum_i x_i(y_i-m_i p_i)\): the intercept component is \((3-2)+(0-1)=0\), and the slope component is \(0(3-2)+1(0-1)=-1\). The grouped weights use trial counts because the modeled response is the success count, not a single Bernoulli outcome.

Worked Example: Gamma Log Link

For a short check using a different variance function, take two Gamma observations with outcomes \(y_1=4,y_2=2\), design values \(x_i=1\), offsets \(o_1=\log 2,o_2=0\), and current coefficient \(\beta=0\). Let \(\phi=1\). The log link gives means \(\mu_1=2,\mu_2=1\) and predictors \(\eta_1=\log 2,\eta_2=0\).

For the Gamma variance function \(V(\mu)=\mu^2\), the log-link derivative is \(d_i=\mu_i\), so \(W_i=\mu_i^2/\mu_i^2=1\). The working residuals are \(z_1-\eta_1=(4-2)/2=1\) and \(z_2-\eta_2=(2-1)/1=1\). Their weighted sum, and hence the coefficient score, is \(2\). Directly, the Gamma log-link score is \(\sum_i x_i(y_i-\mu_i)/(\phi\mu_i)=2/2+1/1=2\). The weight is constant here, but the working residual is not the raw response residual.

What the Identity Says About the Least-Squares Step

At the current coefficients, the IRLS working regression fits \(z-o\) on \(X\) with weights \(W_i\). Suppose \(W\) is positive definite, as it is under the assumptions above. The weighted least-squares normal equations for a candidate coefficient vector \(b\) are

$$ X^T W X b=X^T W(z-o). $$

Subtract \(X^TWX\beta\) from both sides and use \(\eta=X\beta+o\). The result is

X^T W X(b-\beta) = X^T W(z-o-X\beta) = X^T W(z-\eta) = U(\beta).

This gives a second useful result: the increment from the current coefficients to the working least-squares solution satisfies the current information-weighted linear system with the current score on the right-hand side. If \(X^TWX\) is nonsingular, the increment is uniquely determined. If it is singular, the normal equations still hold, but the coefficient solution need not be unique.

The working response and weights are held fixed while this particular weighted least-squares problem is solved. The identity does not say that \(z\) is newly observed data, or that the working residual at the solution equals the score evaluated there. It matches the current score to the residual at the current coefficients.

Common Mistakes and a Practical Check

Using a square-root weight in the score identity. The transformed weighted-least-squares residual may involve \(\sqrt{W_i}(z_i-\eta_i)\), but the score contribution is \(x_i W_i(z_i-\eta_i)\).
Dropping the offset. With an offset, the residual is \(z_i-(x_i^T\beta+o_i)\). Equivalently, the working-regression response is \(z_i-o_i\) when regressed on \(X\).
Mixing response conventions. A grouped-binomial count and a success proportion have different variance and weight formulas. Use the mean, variance, and working response for the response actually modeled.
Misplacing dispersion or prior weights. Under the stated convention, \(W_i\) contains \(a_i/\phi\). Omitting those factors changes the weighted residual and generally no longer reproduces the likelihood score under that convention.
Expecting the weighted residual alone to be the coefficient score. The score is the design-weighted sum \(X^TW(z-\eta)\). The design rows must still multiply the observation-level quantities.

Compute the current fit.
Find \(\eta_i=x_i^T\beta+o_i\), then \(\mu_i\) and \(d_i=d\mu_i/d\eta_i\).

Form each working quantity.
Calculate \(z_i-\eta_i=(y_i-\mu_i)/d_i\) and \(W_i=a_i d_i^2/\{\phi V(\mu_i)\}\).

Check against the score.
Verify that \(W_i(z_i-\eta_i)=a_i d_i(y_i-\mu_i)/\{\phi V(\mu_i)\}\) for each observation, then aggregate with \(X^T\).

Key takeaway: A correctly constructed IRLS working response and weight reproduce the GLM score at the current coefficients exactly: \(U(\beta)=X^TW(z-\eta)\). This identity checks the algebra before the weighted least-squares system is used.

Check Your Understanding

Use the score identity to reason through each question.

For a noncanonical link, which derivative appears in the working response, and which power of that derivative appears in the IRLS weight?
In a model with an offset, what predictor should be subtracted from the working response to form the current residual?
Why does the score use \(W_i(z_i-\eta_i)\), rather than \(\sqrt{W_i}(z_i-\eta_i)\)?
For a grouped-binomial model of success counts, what feature of the variance function affects the IRLS weight?
What equation does the difference between the weighted least-squares solution and the current coefficients satisfy?