Fisher Scoring and IRWLS: A Unified Approach

From Newton-Raphson to Fisher Scoring

Most numerical optimization methods are based on the Newton-Raphson algorithm. You might remember it from calculus as a method for finding the roots of a function. We want to find the roots of the score function, $S(\mathbf{\beta })$.

The Newton-Raphson update step to find the maximum of a function is:
${\mathbf{\beta }}^{(t+1)}={\mathbf{\beta }}^{(t)}-[H({\mathbf{\beta }}^{(t)}){]}^{-1}S({\mathbf{\beta }}^{(t)})$

Where:

• ${\mathbf{\beta }}^{(t)}$ is the vector of parameter estimates at iteration $t$.
• $S({\mathbf{\beta }}^{(t)})$ is the score vector (first derivative of log-likelihood) evaluated at ${\mathbf{\beta }}^{(t)}$.
• $H({\mathbf{\beta }}^{(t)})$ is the Hessian matrix (second derivative of log-likelihood) evaluated at ${\mathbf{\beta }}^{(t)}$.

The Hessian matrix, $H(\mathbf{\beta })$, measures the curvature of the log-likelihood function. While Newton-Raphson is very powerful, calculating the full Hessian matrix can be computationally difficult and unstable.

The Fisher Scoring Algorithm

The Fisher Scoring algorithm is a modification of Newton-Raphson. Instead of using the observed Hessian matrix, $H(\mathbf{\beta })$, it uses its expected value, which is the negative of the Fisher Information Matrix, $I(\mathbf{\beta })$.

Key Idea: Replace $H(\mathbf{\beta })$ with $E[H(\mathbf{\beta })]=-I(\mathbf{\beta })$.

The Fisher Scoring update rule then becomes:
${\mathbf{\beta }}^{(t+1)}={\mathbf{\beta }}^{(t)}+[I({\mathbf{\beta }}^{(t)}){]}^{-1}S({\mathbf{\beta }}^{(t)})$

This small change has two big advantages:

1. Simpler Calculations: The Fisher Information matrix often has a simpler form than the Hessian.
2. Guaranteed Direction: The information matrix, $I(\mathbf{\beta })$, is always positive semi-definite, which ensures that each step moves towards the maximum of the likelihood function.

The following video provides a great visual and mathematical overview of the Newton-Raphson method and its evolution into Fisher Scoring.

As you saw in the video, a key result is that for the exponential family of distributions, the expected value of the second-order derivative term simplifies nicely. The term involving the data, $(y_i-{\mu }_i)$, has an expectation of zero, which makes the calculation of the Fisher Information matrix much easier than the observed Hessian.

From Fisher Scoring to Iteratively Reweighted Least Squares (IRWLS)

Now, here's where the magic happens. The Fisher Scoring update rule can be rearranged to look exactly like the solution to a weighted least squares (WLS) problem. This gives us a practical, step-by-step algorithm to implement.

Let's look at the components again:

• The score function (in matrix form):
  $S(\mathbf{\beta })={\mathbf{X}}^T\mathbf{W\Delta }(\mathbf{y}-\mathbf{\mu })$
• The Fisher Information matrix:
  $I(\mathbf{\beta })={\mathbf{X}}^T\mathbf{WX}$

Here, $\mathbf{W}$ is a diagonal matrix of weights, and $\mathbf{\Delta }$ is a diagonal matrix related to the derivative of the link function.

Substituting these into the Fisher Scoring update rule:

This equation can be cleverly rearranged to:

Where ${\mathbf{z}}^{(t)}$ is a new, constructed vector called the working response or adjusted dependent variable:

...and the weights are:

Does this look familiar? The equation for ${\mathbf{\beta }}^{(t+1)}$ is precisely the solution to a weighted least squares problem, where you are regressing the working response ${\mathbf{z}}^{(t)}$ on the predictors $\mathbf{X}$ using weights ${\mathbf{W}}^{(t)}$.

Because we have to re-calculate the weights $\mathbf{W}$ and the working response $\mathbf{z}$ at each step (since they depend on ${\mathbf{\beta }}^{(t)}$), the process is called Iteratively Reweighted Least Squares (IRWLS).

The following resources provide a more formal, step-by-step definition of the algorithm.

The IRWLS Algorithm in Practice

So, how does this work in a step-by-step fashion?

The IRWLS Algorithm:

1. Initialization (k=0):

   • Choose an initial estimate for the parameters, ${\mathbf{\beta }}^{(0)}$. A good starting point is often obtained by fitting a standard linear model of $g(y_i)$ on $x_i$.
   • From ${\mathbf{\beta }}^{(0)}$, calculate initial values for the linear predictor ${\mathbf{\eta }}^{(0)}$ and the mean ${\mathbf{\mu }}^{(0)}$.

2. Iteration (k):

   • Calculate the working response $z_i$: $$z_i^{(k)}={\eta }_i^{(k)}+(y_i-{\mu }_i^{(k)})g^{\prime }({\mu }_i^{(k)})$$ This can be seen as the current linear predictor adjusted by the "residual" $(y_i-{\mu }_i)$ scaled by the rate of change of the link function.
   • Calculate the weights $w_i$: $$w_i^{(k)}=\frac{1}{V({\mu }_i^{(k)})[g^{\prime }({\mu }_i^{(k)}){]}^2}$$ Notice that observations with higher variance (larger $V({\mu }_i)$) will receive less weight in the estimation.

3. Update:

   • Solve the weighted least squares problem to get the new parameter estimates ${\mathbf{\beta }}^{(k+1)}$: $${\mathbf{\beta }}^{(k+1)}=({\mathbf{X}}^T{\mathbf{W}}^{(k)}\mathbf{X}{)}^{-1}{\mathbf{X}}^T{\mathbf{W}}^{(k)}{\mathbf{z}}^{(k)}$$

4. Check for Convergence:

   • Repeat steps 2 and 3 until the change in the log-likelihood or the parameter estimates between iterations is very small (below a set tolerance).

To see how this works for specific models, review the examples for logistic and Poisson regression in the resource below.

Conclusion

In this lesson, we demystified the process of fitting a Generalized Linear Model. We saw that because the score equations are often non-linear, we need iterative methods to find the maximum likelihood estimates.

Key takeaways:

• Fisher Scoring is a common algorithm for this, which approximates the complex Hessian matrix with the more stable Fisher Information matrix.
• For GLMs, the Fisher Scoring algorithm is equivalent to Iteratively Reweighted Least Squares (IRWLS).
• The IRWLS algorithm involves creating a "working response" variable and "weights" at each step, and then solving a simple weighted least squares problem to update the parameter estimates until they converge.

Now that we understand how the model parameters are estimated, our next step is to learn how to evaluate the model's fit and perform hypothesis tests on the coefficients. We'll cover topics like deviance, and the Likelihood Ratio, Wald, and Score tests.