IRLS Working Responses and Residuals

Why IRLS Creates a Working Response

In a generalized linear model, the mean is connected to the linear predictor by a link function, so the mean is generally not a linear function of the regression coefficients. Iteratively reweighted least squares (IRLS) handles this by making a local linear approximation: at the current parameter values, it constructs a working response that can be fitted by weighted least squares. The response and weights are recalculated as the fit changes.

The word “working” matters. The working response is not a new measurement of the outcome. It is an intermediate quantity built from the observed response and the current fitted mean. Its residual from the current linear predictor is determined by the local slope of the link.

Definition: For observation $i$, let $Y_i$ have mean $\mu_i$, let $\eta_i=x_i^\mathsf{T}\beta+o_i$ be its linear predictor, and let $\eta_i=g(\mu_i)$, where $g$ is a differentiable link. The IRLS working response at the current fit is $$ z_i=\eta_i+(y_i-\mu_i)g'(\mu_i). $$ If $\operatorname{Var}(Y_i)=\phi V(\mu_i)/a_i$, the corresponding IRLS weight is $$ W_i=\frac{a_i}{\phi V(\mu_i)\{g'(\mu_i)\}^2}. $$ Here $a_i$ is a prior weight, $V$ is the variance function, and $\phi$ is the dispersion parameter.

These formulas apply when the current mean is in the interior of its parameter space, $g$ is differentiable there, and $d\mu_i/d\eta_i$ is finite and nonzero. The observations are assumed independent. In a binomial model for proportions, the number of trials is commonly represented by $a_i$; equivalently, one can formulate the model for the grouped count and use its corresponding variance.

Deriving the Working Response and Weight

Write $\mu(\eta)=g^{-1}(\eta)$. At a current predictor $\eta_i$, a small change $\delta_i$ gives the first-order approximation $$ \mu(\eta_i+\delta_i)\approx \mu_i+\frac{d\mu_i}{d\eta_i}\delta_i. $$ If the observed response $y_i$ is to be matched by this local approximation, set the right side equal to $y_i$. Solving for the required change gives $\delta_i\approx (y_i-\mu_i)/(d\mu_i/d\eta_i)$. Since $d\eta_i/d\mu_i=g'(\mu_i)$, the linear predictor that would produce this local match is the working response $z_i=\eta_i+\delta_i=\eta_i+(y_i-\mu_i)g'(\mu_i)$.

Result: In the first-order approximation around the current predictor, fitting a new predictor to the response is equivalent to fitting the working response by weighted least squares, with weight equal to the inverse variance of the locally transformed response.

To see the weight, consider the response transformed locally to the predictor scale. Its change is approximately $g'(\mu_i)(Y_i-\mu_i)$, so its variance under the current fit is approximately $$ \{g'(\mu_i)\}^2\operatorname{Var}(Y_i) =\{g'(\mu_i)\}^2\frac{\phi V(\mu_i)}{a_i}. $$ The inverse of this variance is $W_i$. Also, writing $\mu'_i=d\mu_i/d\eta_i=1/g'(\mu_i)$, the local response approximation can be expressed as $$ y_i-\mu(\eta_i+\delta_i)\approx \mu'_i(z_i-\eta_i-\delta_i). $$ Thus the squared discrepancy divided by the response variance is, to first order, $$ \frac{\{y_i-\mu(\eta_i+\delta_i)\}^2}{\operatorname{Var}(Y_i)} \approx \frac{(\mu'_i)^2}{\operatorname{Var}(Y_i)} (z_i-\eta_i-\delta_i)^2 =W_i(z_i-\eta_i-\delta_i)^2. $$ This establishes the local weighted least-squares interpretation. It does not say that the original model is exactly a normal linear model: the approximation and its weights are recalculated at each iteration.

Connection to Fisher Scoring

The same working quantities give the Fisher-scoring update for a GLM. State the result for an exponential-dispersion model with independent observations, variance $\phi V(\mu_i)/a_i$, and a differentiable mean-link relationship. Assume the required likelihood derivatives exist and the expected information matrix is nonsingular. Let $X$ have rows $x_i^\mathsf{T}$, and let $W$ be the diagonal matrix with entries $W_i$.

Result: With the offset held fixed, the IRLS weighted least-squares step is the Fisher-scoring step: $$ \beta_{\mathrm{new}}=\beta+(X^\mathsf{T}WX)^{-1}U(\beta), $$ where $U(\beta)$ is the likelihood score. Equivalently, $\beta_{\mathrm{new}}$ solves the weighted least-squares normal equations for the working response after subtracting the offset.

For an exponential-dispersion model, the score for the mean parameter is $(y_i-\mu_i)/\operatorname{Var}(Y_i)$ times the derivative of the mean with respect to that parameter. Applying the chain rule to $\eta_i=x_i^\mathsf{T}\beta+o_i$ therefore gives $$ U(\beta)=\sum_i x_i\frac{a_i(y_i-\mu_i)}{\phi V(\mu_i)} \frac{d\mu_i}{d\eta_i}. $$ The expected information is $$ I(\beta)=X^\mathsf{T}WX,\qquad W_i=\frac{a_i}{\phi V(\mu_i)} \left(\frac{d\mu_i}{d\eta_i}\right)^2. $$ The latter is the same weight as above, since $d\mu_i/d\eta_i=1/g'(\mu_i)$.

For the weighted least-squares step, the proposed predictor is $X\beta_{\mathrm{new}}+o$. Its normal equations are $X^\mathsf{T}W(z-o-X\beta_{\mathrm{new}})=0$. At the current fit, $\eta=X\beta+o$, so $z-\eta$ has entries $(y_i-\mu_i)g'(\mu_i)$. Consequently, $$ X^\mathsf{T}W(z-o-X\beta) =\sum_i x_i W_i(z_i-\eta_i) =\sum_i x_i\frac{a_i(y_i-\mu_i)}{\phi V(\mu_i)} \frac{d\mu_i}{d\eta_i} =U(\beta). $$ Rearranging the normal equations gives $X^\mathsf{T}WX(\beta_{\mathrm{new}}-\beta)=U(\beta)$, which is precisely the Fisher-scoring update. This proves the equivalence. $\square$

When $\phi$ is common to every observation, omitting its factor from all the weights multiplies the whole weight matrix by a common constant. That does not change the weighted least-squares coefficient update, although it does change the numerical scale of the weights.

The Working Residual at a Fitted Observation

At a fitted observation, evaluate the working response using the fitted linear predictor $\hat\eta_i$, fitted mean $\hat\mu_i$, and observed response $y_i$. Define the working residual as the difference between that working response and the fitted linear predictor:

Definition: The working residual for observation $i$ is $$ z_i-\hat\eta_i=(y_i-\hat\mu_i)g'(\hat\mu_i). $$ It is the response residual multiplied by the derivative of the link evaluated at the fitted mean.

This identity follows directly by subtracting $\hat\eta_i$ from the fitted working response $z_i=\hat\eta_i+(y_i-\hat\mu_i)g'(\hat\mu_i)$. The factor $g'(\hat\mu_i)$ changes with the fitted mean, so the working residual is not generally equal to the response residual. It is also distinct from the weight: the residual contains one factor of the link derivative, whereas the weight contains its square in the denominator.

Worked Example: Poisson Counts with a Log Link

Suppose a Poisson observation has $y=4$, current fitted mean $\hat\mu=3$, unit prior weight, and log link $g(\mu)=\log(\mu)$. Take $\phi=1$. Compute the working response, working weight, and working residual.

The fitted linear predictor is $\hat\eta=\log(3)$, and the link derivative is $g'(\hat\mu)=1/3$. Therefore $$ z=\log(3)+(4-3)\frac{1}{3} =\log(3)+\frac{1}{3} \approx 1.4319. $$ For a Poisson variable, $V(\mu)=\mu$, so $$ W=\frac{1}{3(1/3)^2}=3. $$ The working residual is $z-\hat\eta=1/3$, which also equals $(4-3)g'(3)=1/3$. The response residual is $4-3=1$; it is not the working residual.

Worked Example: A Binomial Proportion with a Logit Link

A group has $m=10$ trials and $y=6$ successes. Model its observed proportion $y/m=0.6$, with current fitted probability $\hat p=0.5$, unit dispersion, and logit link $g(p)=\log\{p/(1-p)\}$. For a proportion, the variance function is $V(p)=p(1-p)$ and the prior weight is $a=m$.

The fitted predictor is $\hat\eta=\log(0.5/0.5)=0$. Differentiating the link gives $g'(p)=1/\{p(1-p)\}$, so $g'(0.5)=4$. Hence $$ z=0+(0.6-0.5)(4)=0.4. $$ The weight is $$ W=\frac{10}{0.5(1-0.5)(4)^2} =\frac{10}{0.25(16)} =2.5. $$ Thus the working residual is $z-\hat\eta=0.4$. The same value follows from $(0.6-0.5)g'(0.5)=0.1(4)=0.4$. This calculation is for the proportion response; the factor $m=10$ enters its IRLS weight.

Worked Example: A Gamma Mean with a Log Link

Suppose a Gamma observation has $y=8$, current fitted mean $\hat\mu=5$, unit prior weight, and log link. Use the Gamma variance function $V(\mu)=\mu^2$ and take $\phi=1$.

Here $\hat\eta=\log(5)$ and $g'(\mu)=1/\mu$, so $g'(5)=0.2$. The working response is $$ z=\log(5)+(8-5)(0.2) =\log(5)+0.6 \approx 2.2094. $$ The working weight is $$ W=\frac{1}{5^2(0.2)^2} =\frac{1}{25(0.04)} =1. $$ The working residual is $z-\hat\eta=0.6$, also obtained as $(8-5)g'(5)=3(0.2)=0.6$. Although the response residual is $3$, the local transformation to the log-predictor scale reduces it by the factor $1/5$.

Common Mistake: Reversing the Link Derivative

A frequent error is to multiply by $d\mu/d\eta$ when constructing the working response. The first-order expansion is $d\mu/d\eta$ times a change in the predictor, but solving that expansion for the required predictor change divides by $d\mu/d\eta$. The resulting multiplier is $d\eta/d\mu=g'(\mu)$, not its reciprocal. Another error is to use the response variance itself as the weight; weighted least squares uses its inverse after the local transformation, giving $W_i=a_i/\{\phi V(\mu_i)g'(\mu_i)^2\}$.

For example, with the Poisson log-link fit above, $d\mu/d\eta=3$. Multiplying the response residual $1$ by $3$ would give $3$, but the correct working residual is $1/3$, because $g'(3)=1/3$. Checking this against the definition $z-\hat\eta$ is a useful way to catch the inversion error.

Key takeaway: Linearizing the mean around the current predictor gives $z_i=\eta_i+(y_i-\mu_i)g'(\mu_i)$. The inverse variance of the locally transformed response gives $W_i=a_i/\{\phi V(\mu_i)g'(\mu_i)^2\}$. At a fitted observation, the working residual is $(y_i-\hat\mu_i)g'(\hat\mu_i)$.

Check Your Understanding

Use the current fitted mean and predictor in each calculation.

For a model with $g'(\mu)=2$, $y=5$, and $\hat\mu=4$, what is the working residual?
For a Poisson model with log link, $V(\mu)=\mu$, and prior weight $a_i$, write the IRLS weight as a function of the current mean and $a_i$, taking $\phi=1$.
Why does the working response use $g'(\mu)$, rather than $d\mu/d\eta$, as the multiplier of the response residual?
In the grouped-binomial example, what role does the number of trials play in the weight for a proportion response?
Explain why the weighted least-squares fit is an iteration within a GLM procedure, rather than an assertion that the original response follows a normal linear model.