Are GLMs just glorified WLS regressions?

Question

When performing weighted least squares $L = \frac{1}{2} \sum_i w_i r_i^2$, Aitken showed that one ought to weight each sample by the inverse of its variance $w_i=1/\sigma_i^2$. This leads to gradients of the form

$$\nabla_\beta L = \sum_i \frac{r_i}{\sigma_i^2}\nabla_\beta r_i $$

In GLMs, the log-likelihood is $\ell(\beta) = \sum_i \frac{\theta_i y_i - b(\theta_i)}{\phi_i} + c(y_i, \phi_i)\qquad$ (cf. Turner's notes)

for which the gradients are $\frac{\partial \ell}{\partial \beta} =\sum_i \frac{y_{i}-\mu_{i}}{\operatorname{V}[\mu_i]} \cdot \frac{x_{i}}{g^{\prime}\left(\mu_{i}\right)} $.

This looks suspiciously similar to the WLS gradient. What are the precise connections between the two? Superficially it seems like all a GLM is really doing is

(1) transforming the prediction into a useful domain via the link function

(2) weighting the gradient updates 'optimally' under the assumption that $\sigma_i^2=V[\mu_i]$.

Answer 1

NO. What you have discovered is the idea behind IRLS (Iteratively Reweighted Least Squares): see Can you give a simple intuitive explanation of IRLS method to find the MLE of a GLM? for an exposition. But this is just one numerical method for computing the maximum likelihood estimates (really a version of Newton's method), and not necessarily the best one.

Generalized linear models (GLMs) provide a general framework for formulating statistical models, and have led to terminology and unified methods that permeate modern statistics. That should not be reduced to "just" a numerical optimization method.

So yes, they are connected, and that connection leads to IRLS, that the Newton method can be formulated as iteration of approximating weighted linear regressions. But that does not mean that using a GLM conceptually is the same as using WLS. For instance, GLMs as logistic regression can be used with binary outcomes, and it is not clear how WLS could be used directly in that case. GLMs as a modeling framework is much richer than WLS.

Answer 2

GLMs are most commonly fit using iteratively reweighted least squares, see here and references list there, and this post. This method is based on maximizing the maximum likelihood objective based on Fisher scoring, which is a variant of Newton-Raphson.

A minimal implementation is:

glm.irls = function(X, y, family=binomial, maxit=25, tol=1e-08, beta.start=rep(0,ncol(X))) {
    if (is.function(family)) family <- family()
    beta = beta.start
    for(j in 1:maxit)
    {
      eta    = X %*% beta
      g      = family$linkinv(eta)
      gprime = family$mu.eta(eta)
      z      = eta + (y - g) / gprime
      W      = as.vector(gprime^2 / family$variance(g))
      betaold   = beta
      beta      = as.matrix(coef(lm.wfit(x=X, y=z, w=W)),ncol=1) # regular weighted LS fit = solve(crossprod(X,W*X), crossprod(X,W*z))
      if(sqrt(crossprod(beta-betaold)) < tol) break
    }
    list(coefficients=beta, iterations=j, z=z, weights=W, X=X, y=y, wlmfit=lm(z~1+X[,-1], weights=W))
}

If you use a distribution with an identity link you can see that in the algorithm above z=y and each iteration just comes down to doing a weighted least squares regression with 1/variance weights. For Poisson e.g. one would then use initial weights = 1/(y+small epsilon) (since for Poisson the expected variance=the mean) and iterate this based on the predicted yhat, where your weights will then become 1/(yhat+small epsilon). With Gaussian errors (ie regular OLS regression) the weights would all just be equal to 1. So GLMs do reduce to the iterated fitting of weighted least squares regression, and with identity link do just use 1/variance weights. To estimate the 1/variance weights an iterative procedure has to be used though. If one would approximate the true ML objective using a single weighted least square analysis then this would only be approximately correct though. However, in practice, this can still be a useful approximation, see here for an example. It is not the case though that GLMs are just glorified weighted least squares models...

Example logistic regression:

data("Contraception",package="mlmRev")
R_GLM = glm(use ~ age + I(age^2) + urban + livch,
            family=binomial, x=T, data=Contraception)
IRLS_GLM = glm.irls(X=R_GLM$x, y=R_GLM$y, family=binomial)
print(data.frame(R_GLM=coef(R_GLM), IRLS_GLM=coef(IRLS_GLM))) # coefficients match with glm output
                   R_GLM     IRLS_GLM
(Intercept) -0.949952124 -0.949952124
age          0.004583726  0.004583726
I(age^2)    -0.004286455 -0.004286455
urbanY       0.768097459  0.768097459
livch1       0.783112821  0.783112821
livch2       0.854904050  0.854904050
livch3+      0.806025052  0.806025052