Why GLM Fitting Becomes Weighted Least Squares
A generalized linear model (GLM) relates the mean response to predictors through a link function. Its likelihood is generally not a quadratic function of the regression coefficients, so maximizing it usually requires an iterative method. Iteratively reweighted least squares (IRLS) makes each step of that optimization look like a weighted linear regression.
The weights and adjusted responses change from one iteration to the next. They are not new observed data: they summarize the local likelihood around the current coefficient estimates. The derivation below starts with the likelihood score and expected information, then shows why a local quadratic approximation leads to the IRLS equations.
Here \(X\) is the design matrix with row \(\mathbf{x}_i^\mathsf{T}\), \(W\) is the diagonal matrix with entries \(w_i\), and \(\mathbf{z}\) collects the working responses. The weights are positive when the variance is positive, the dispersion is positive, and the link has a finite nonzero derivative at the current fitted mean. If a derivative is zero or undefined, or if \(X^\mathsf{T}WX\) is singular, this ordinary update cannot be applied as written.
From the Score to the IRLS Equations
For an exponential-family GLM, differentiating the log likelihood with respect to \(\eta_i\) gives \(a_i(y_i-\mu_i)\dot{\mu}_i/\{\phi V(\mu_i)\}\). Applying the chain rule, since \(d\eta_i/d\boldsymbol{\beta}=\mathbf{x}_i\), gives the coefficient score \(\mathbf{U}(\boldsymbol{\beta})=\sum_i\mathbf{x}_i a_i(y_i-\mu_i)\dot{\mu}_i/\{\phi V(\mu_i)\}\). The expected information is \(I(\boldsymbol{\beta})=X^\mathsf{T}WX\), with the weights defined above. This follows because the variance of the response residual is \(\phi V(\mu_i)/a_i\): the variance of its contribution to the linear-predictor score is \(a_i\dot{\mu}_i^2/\{\phi V(\mu_i)\}\).
A second-order local approximation to the log likelihood uses the score as its linear term and the expected information as its curvature. Maximizing that quadratic approximation gives the Fisher-scoring step. This is the justification for the update; it does not say that the actual log likelihood is globally quadratic. Newton–Raphson instead uses the observed Hessian, which need not equal the expected information for a noncanonical link.
Differentiating this quadratic expression with respect to \(\boldsymbol{\Delta}\) and setting the result to zero gives \(I(\boldsymbol{\beta})\boldsymbol{\Delta}=\mathbf{U}(\boldsymbol{\beta})\). Since \(\boldsymbol{\eta}=X\boldsymbol{\beta}\), the working response can be written as \(\mathbf{z}=X\boldsymbol{\beta}+\mathbf{r}\), where \(r_i=(y_i-\mu_i)/\dot{\mu}_i\). Multiplying \(\mathbf{z}\) by \(W\) and \(X^\mathsf{T}\) shows the link between the score equation and weighted least squares:
Therefore the weighted least-squares solution satisfies \(X^\mathsf{T}WX\boldsymbol{\beta}_{\mathrm{new}}= X^\mathsf{T}W\mathbf{z}\), so subtracting the current coefficient vector gives \(\boldsymbol{\beta}_{\mathrm{new}}-\boldsymbol{\beta}= (X^\mathsf{T}WX)^{-1}\mathbf{U}(\boldsymbol{\beta})\), when the inverse exists. This is exactly the Fisher-scoring increment. The algebra explains why each IRLS iteration is a weighted least-squares problem, even though the original response model may not be Gaussian.
Worked IRLS Steps
Worked Example: One Grouped-Logit Scoring Update
A hypothetical portfolio is divided into two risk groups. In each group, four policies are observed and the response is whether a policy has a specified claim indicator. The groups have predictor values \(x=0\) and \(x=1\), and observed claim counts 2 and 3, respectively. Fit a logistic model with intercept and slope, starting at \(\boldsymbol{\beta}=(0,0)^\mathsf{T}\). Calculate one IRLS update.
Method and current fitted values. For the logit link, \(\mu_i=p_i\), \(\dot{\mu}_i=p_i(1-p_i)\), and the binomial variance function for a response proportion is \(V(p_i)=p_i(1-p_i)\). The prior weight is the group size, \(a_i=4\), and \(\phi=1\). At the starting coefficients, both linear predictors are zero, so both fitted probabilities are \(p_i=0.5\), with \(\dot{\mu}_i=0.25\).
The observed proportions are \(2/4=0.5\) and \(3/4=0.75\). The working responses are \(z_1=0+(0.5-0.5)/0.25=0\) and \(z_2=0+(0.75-0.5)/0.25=1\). Each weight is \(4(0.25)^2/0.25=1\).
Weighted least squares. The design rows are \((1,0)\) and \((1,1)\), so \(X^\mathsf{T}WX=\begin{pmatrix}2&1\\1&1\end{pmatrix}\) and \(X^\mathsf{T}W\mathbf{z}=(1,1)^\mathsf{T}\). The determinant is \(2(1)-1(1)=1\), and solving gives \(\boldsymbol{\beta}_{\mathrm{new}}=(0,1)^\mathsf{T}\).
Check and trap avoided. The score at the starting values is \((1,1)^\mathsf{T}\): the first group contributes zero and the second contributes \(4(0.75-0.5)(1,1)^\mathsf{T}=(1,1)^\mathsf{T}\). Multiplying the information matrix by the update increment gives \(\begin{pmatrix}2&1\\1&1\end{pmatrix}(0,1)^\mathsf{T}=(1,1)^\mathsf{T}\), as required. A common mistake is to use the second group’s raw count, 3, in place of its response proportion, 0.75. That would give the incorrect working response \((3-0.5)/0.25=10\); the grouped-binomial calculation uses the observed proportion and the group-size prior weight.
Worked Example: A Poisson Log-Link Update
A hypothetical service team records event counts for two sites, with predictor values \(x=0\) and \(x=1\). The observed counts are 1 and 4. Fit a Poisson log-link model with an intercept and slope, starting at \(\boldsymbol{\beta}=(\log 2,0)^\mathsf{T}\). Assume unit prior weights and \(\phi=1\). Find one update.
Method and working quantities. With a log link, \(\eta_i=\log\mu_i\), so \(\dot{\mu}_i=\mu_i\). The Poisson variance function is \(V(\mu_i)=\mu_i\). Both current fitted means are 2. Thus the working responses are \(z_1=\log 2+(1-2)/2=\log 2-0.5\) and \(z_2=\log 2+(4-2)/2=\log 2+1\). The weights are \(w_i=2^2/2=2\).
Solve the update. The score is \(X^\mathsf{T}(y-\mu)=(1,2)^\mathsf{T}\): the intercept component is \((-1)+2=1\), while the slope component is the second residual, 2. The information matrix is \(X^\mathsf{T}WX=\begin{pmatrix}4&2\\2&2\end{pmatrix}\). Solving \(\begin{pmatrix}4&2\\2&2\end{pmatrix}\boldsymbol{\Delta}=(1,2)^\mathsf{T}\) gives \(\boldsymbol{\Delta}=(-0.5,1.5)^\mathsf{T}\). Therefore \(\boldsymbol{\beta}_{\mathrm{new}}=(\log 2-0.5,1.5)^\mathsf{T}\).
Check and trap avoided. The weighted least-squares equations have right-hand side \(X^\mathsf{T}W\mathbf{z}=(4\log 2+1,\,2\log 2+2)^\mathsf{T}\). Multiplying the information matrix by \((\log 2-0.5,1.5)^\mathsf{T}\) gives exactly that vector. A common error is to take the intercept score as \(-1\), forgetting that it sums both residuals. The correct intercept score is \(-1+2=1\).
Worked Example: Gaussian Identity Link as a Special Case
A hypothetical analyst observes responses \(y=(1,2,2)^\mathsf{T}\) at predictor values \(x=0,1,2\). Fit a Gaussian GLM with identity link, unit prior weights, and known dispersion \(\phi=1\). Starting from any coefficient vector, show what the IRLS weighted least-squares update produces.
Method. For the identity link, \(\eta_i=\mu_i\) and \(\dot{\mu}_i=1\). The Gaussian variance function is \(V(\mu_i)=1\), so \(w_i=1\) and \(z_i=\eta_i+(y_i-\mu_i)=y_i\). Thus this iteration is ordinary least squares on the original responses, regardless of the starting values.
Use design rows \((1,0),(1,1),(1,2)\). Then \(X^\mathsf{T}X=\begin{pmatrix}3&3\\3&5\end{pmatrix}\) and \(X^\mathsf{T}y=(5,6)^\mathsf{T}\), since \(1+2+2=5\) and \(0(1)+1(2)+2(2)=6\). The determinant is \(3(5)-3(3)=6\). Solving gives \(\hat{\beta}_0=(5(5)-3(6))/6=7/6\) and \(\hat{\beta}_1=(3(6)-3(5))/6=1/2\).
Check and trap avoided. The normal equations check: \(3(7/6)+3(1/2)=5\) and \(3(7/6)+5(1/2)=6\). It would be incorrect to conclude that all GLMs take one step to their final estimates. This one-step result holds here because the Gaussian identity-link likelihood is quadratic in the coefficients and the working responses and weights do not change with the iteration.
Assumptions, Interpretation, and Common Traps
IRLS is an optimization procedure, not a guarantee that every update is valid or that the likelihood has a well-behaved finite maximum. The update requires the relevant means and derivatives to be defined, positive variance weights, and a full-rank weighted design for the displayed inverse. In practice, the iterations begin at chosen starting coefficients and continue until changes in coefficients or the likelihood are sufficiently small. A single step is generally only a local improvement, not a final estimate.
The prior weight \(a_i\) and the IRLS weight \(w_i\) have different roles. The prior weight is specified by the data structure, such as a binomial group size. The IRLS weight also depends on the current fitted mean, the variance function, and the link derivative, so it is recalculated at every iteration. Do not confuse either with an arbitrary claim exposure or offset: an offset changes the linear predictor, while it is not automatically an IRLS weight.
- Omitting the link derivative. The working response uses \((y_i-\mu_i)/\dot{\mu}_i\), and the weight uses \(\dot{\mu}_i^2\). These factors arise from converting a change in mean to a change in linear predictor.
- Using the wrong variance function. For Poisson, \(V(\mu)=\mu\); for a Bernoulli proportion, \(V(\mu)=\mu(1-\mu)\). The family determines the variance function used in the weight.
- Confusing a response with a working response. \(z_i\) is recalculated from the current fit. It is not generally the original \(y_i\), except in special cases such as the Gaussian identity-link model.
- Assuming every GLM update is ordinary least squares. It is weighted least squares, and both the pseudo-response and weights usually change at each iteration.
Check Your Understanding
For each problem, identify the current working quantities and use them to obtain the requested update or conclusion.
- A Poisson log-link model has an intercept only, current mean \(\mu=3\), observed count \(y=5\), unit prior weight, and \(\phi=1\). Calculate the working response and weight. (Answer: \(z=\log 3+2/3\); \(w=3\).)
- A grouped-binomial logit model has one group of 6 trials, 4 successes, and current linear predictor zero. Calculate the observed proportion, working response, and IRLS weight. (Answer: Proportion \(=2/3\); \(z=2/3\); \(w=1.5\).)
- For a Gaussian identity-link model with known \(\phi=4\), one observation has current fitted mean 7 and observed response 9, with prior weight 1. Find its working response and IRLS weight. (Answer: \(z=9\); \(w=0.25\).)
- A Poisson log-link model has one intercept and two observations, with counts 2 and 6. At the current intercept, both fitted means are 4. Find the score, information, scoring increment, and updated intercept. (Answer: Score \(=0\); information \(=8\); increment \(=0\); updated intercept \(=\log 4\).)
- A logistic model with a single intercept has one Bernoulli observation \(y=1\), current \(\eta=0\), and \(\phi=a=1\). Calculate its score, information, working response, and next intercept. (Answer: Score \(=0.5\); information \(=0.25\); \(z=2\); next intercept \(=2\).)