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