One Iteration of IRWLS for GLM Parameter Update

A Quick Recap: The IRWLS Loop

As a brief reminder, the IRWLS algorithm solves for the GLM parameters $\mathbf{\beta }$ by repeatedly solving a weighted least squares (WLS) problem. The update step for the parameters at iteration $(t+1)$ is given by:

${\mathbf{\beta }}^{(t+1)}=({\mathbf{X}}^T{\mathbf{W}}^{(t)}\mathbf{X}{)}^{-1}{\mathbf{X}}^T{\mathbf{W}}^{(t)}{\mathbf{z}}^{(t)}$

This looks just like the formula for ordinary least squares, but with two key differences:

1. We use a working response vector $\mathbf{z}$ instead of the actual data $\mathbf{y}$.
2. We use a weight matrix $\mathbf{W}$ to give more influence to observations with lower variance.

Since $\mathbf{z}$ and $\mathbf{W}$ depend on the current parameter estimates ${\mathbf{\beta }}^{(t)}$, we must recalculate them at each step, hence the name "Iteratively Reweighted."

For a deeper dive into the mathematical framework, the following resources provide a solid foundation.

Worked Example: One Iteration of IRWLS for Poisson Regression

Now, let's get our hands dirty. Your task is to perform one complete update of the $\mathbf{\beta }$ coefficients for a Poisson regression model.

Problem Setup:
Suppose we are modeling the number of software bugs (y) found in a code module based on its complexity score (x). We have the following data from four modules:

We propose a Poisson regression model:

• Random Component: $Y_i\sim \text{Poisson}({\mu }_i)$
• Systematic Component (Linear Predictor): ${\eta }_i={\beta }_0+{\beta }_1{x}_i$
• Link Function: $\log ({\mu }_i)={\eta }_i$ (Canonical link)

For a Poisson model, we know:

• Variance function: $V({\mu }_i)={\mu }_i$
• Derivative of the link function: $g^{\prime }({\mu }_i)=\frac{d}{d{\mu }_i}\log ({\mu }_i)=\frac{1}{{\mu }_i}$

Our initial guess for the parameters is ${\mathbf{\beta }}^{(0)}=\left(\begin{array}{c}{\beta }_0^{(0)}\\ {\beta }_1^{(0)}\end{array}\right)=\left(\begin{array}{c}0.1\\ 0.1\end{array}\right)$.

Let's execute the first iteration ($t=0$) to find ${\mathbf{\beta }}^{(1)}$.

Step 1: Calculate the Linear Predictor (${\mathbf{\eta }}^{(0)}$) and Mean Response (${\mathbf{\mu }}^{(0)}$)

First, we use our initial ${\mathbf{\beta }}^{(0)}$ to calculate the linear predictor ${\eta }_i^{(0)}$ and the expected number of bugs ${\mu }_i^{(0)}$ for each module.

• ${\eta }_i^{(0)}={\beta }_0^{(0)}+{\beta }_1^{(0)}{x}_i=0.1+0.1x_i$
• ${\mu }_i^{(0)}=\exp ({\eta }_i^{(0)})$

Step 2: Calculate the Working Response (${\mathbf{z}}^{(0)}$) and Weights (${\mathbf{w}}^{(0)}$)

Now we compute the two key components for our WLS update.

The general formula for the working response is:
$z_i={\eta }_i+(y_i-{\mu }_i)g^{\prime }({\mu }_i)$
For our Poisson model with $g^{\prime }({\mu }_i)=1/{\mu }_i$, this becomes:
$z_i^{(0)}={\eta }_i^{(0)}+\frac{y_i-{\mu }_i^{(0)}}{{\mu }_i^{(0)}}$

The general formula for the weights is:
$w_i=\frac{1}{V({\mu }_i)[g^{\prime }({\mu }_i){]}^2}$
For our Poisson model with $V({\mu }_i)={\mu }_i$ and $g^{\prime }({\mu }_i)=1/{\mu }_i$, this simplifies beautifully:
$w_i^{(0)}=\frac{1}{{\mu }_i^{(0)}(1/{\mu }_i^{(0)}{)}^2}={\mu }_i^{(0)}$

Let's calculate these values.

You can see these calculation steps outlined for a general Poisson model in the resource below.

Step 3: Perform the Weighted Least Squares Regression

We now have all the pieces to solve for ${\mathbf{\beta }}^{(1)}$ using the WLS formula: ${\mathbf{\beta }}^{(1)}=({\mathbf{X}}^T{\mathbf{W}}^{(0)}\mathbf{X}{)}^{-1}{\mathbf{X}}^T{\mathbf{W}}^{(0)}{\mathbf{z}}^{(0)}$.

Let's define our matrices:

1. Calculate ${\mathbf{X}}^T{\mathbf{W}}^{(0)}\mathbf{X}$: $${\mathbf{X}}^T{\mathbf{W}}^{(0)}\mathbf{X}=\left(\begin{array}{cc}18.312& 257.408\\ 257.408& 3781.16\end{array}\right)$$
2. Calculate the inverse $({\mathbf{X}}^T{\mathbf{W}}^{(0)}\mathbf{X}{)}^{-1}$:
   For a 2x2 matrix $\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)$, the inverse is $\frac{1}{ad-bc}\left(\begin{array}{cc}d& -b\\ -c& a\end{array}\right)$. $$({\mathbf{X}}^T{\mathbf{W}}^{(0)}\mathbf{X}{)}^{-1}=\frac{1}{2921.2}\left(\begin{array}{cc}3781.16& -257.408\\ -257.408& 18.312\end{array}\right)=\left(\begin{array}{cc}1.294& -0.088\\ -0.088& 0.00627\end{array}\right)$$
3. Calculate ${\mathbf{X}}^T{\mathbf{W}}^{(0)}{\mathbf{z}}^{(0)}$: $${\mathbf{X}}^T{\mathbf{W}}^{(0)}{\mathbf{z}}^{(0)}=\left(\begin{array}{c}25.59\\ 362.47\end{array}\right)$$
4. Finally, calculate ${\mathbf{\beta }}^{(1)}$: $${\mathbf{\beta }}^{(1)}=({\mathbf{X}}^T{\mathbf{W}}^{(0)}\mathbf{X}{)}^{-1}({\mathbf{X}}^T{\mathbf{W}}^{(0)}{\mathbf{z}}^{(0)})=\left(\begin{array}{cc}1.294& -0.088\\ -0.088& 0.00627\end{array}\right)\left(\begin{array}{c}25.59\\ 362.47\end{array}\right)=\left(\begin{array}{c}1.20\\ 0.02\end{array}\right)$$

So, after one iteration, our updated parameter estimates are ${\beta }_0^{(1)}=1.20$ and ${\beta }_1^{(1)}=0.02$.

Statistical software would now take these new $\mathbf{\beta }$ values and repeat the entire process—calculating new $\eta$, $\mu$, $z$, and $W$—until the change in $\mathbf{\beta }$ between iterations is negligible.

Conclusion

In this lesson, you walked through the core computational engine for fitting Generalized Linear Models. You have now executed a full iteration of the IRWLS algorithm, translating abstract formulas into concrete calculations.

Key Takeaways:

• The IRWLS algorithm consists of a loop with three main steps per iteration:
  1. Calculate the linear predictor $\mathbf{\eta }$ and mean response $\mathbf{\mu }$ from the current $\mathbf{\beta }$.
  2. Use these to compute the working response $\mathbf{z}$ and weights $\mathbf{W}$.
  3. Solve the weighted least squares equation ${\mathbf{\beta }}^{\text{new}}=({\mathbf{X}}^T\mathbf{WX}{)}^{-1}{\mathbf{X}}^T\mathbf{Wz}$ to get the updated parameters.
• The specific formulas for the weights and working response depend on the GLM's distribution (variance function) and link function.

You've now seen how a GLM is fit. But how do we know if the fit is any good? In our next lesson, we will address this by learning how to compute the deviance of a fitted GLM and interpret it as a measure of model fit.