I was wondering why R uses the Iterated Re-weighted Least Square in the GLM case?
When defining the problem analytically using Newton-Raphson / Fisher scoring, one comes to the following (vector notation):
$$\beta_{t+1} = \beta_t + (X^TWX)^{-1}X^TWM(y-\mu)$$
Then a "trick" is employed, to transform this to:
$$\beta_{t+1} = (X^TWX)^{-1}X^TW(X\beta_t + M(y-\mu)) := (X^TWX)^{-1}X^TWZ_t$$
And the solution is found using the Re-Weighted Least Squares algorithm.
I wonder why this extra trick/step is needed? Is it just because the algorithm used to solve RWLS is more efficient and numerically stable (uses QR decomposition)? Or is there some other underlying reason?