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, i.e. using iteratively reweighted least squares based on a quadratic approximation of the log likelihood.
A minimal implementation is:
glm_irls =
function(X, # design matrix
y, # response
weights = rep(1, ncol(X)), # prior observation weights
# e.g. total nr. of trials for proportions with family=binomial
start = rep(1, ncol(X)), # coefficient starting values
offset = rep(0, nrow(X)), # model offset
family = gaussian(identity), # distribution & link function
maxit = 25,
tol = 1e-08) {
beta = start
nobs = nrow(X)
nvars = ncol(X)
eval(family$initialize) # initializes n and mustart
eta = family$linkfun(mustart) # initialize η = g(µ)
mu = family$linkinv(eta) # predictions on response scale µ
for (i in 1:maxit)
{
var = family$variance(mu) # variance in function of the mean µ
gprime = family$mu.eta(eta) # derivative of link function w.r.t. η = d(g-1)/dη=dμ/dη
gradient = y - mu # gradient of log-likelihood with respect to η = ∂ℓ/∂η = deviance residual
z = eta - offset + gradient / gprime
# adjusted response
# = linearised version of log-likelihood function ℓ around η
W = weights * as.vector(gprime^2 / var)
# = working weights
betaold = beta
wlmfit = lm.wfit(x=X, y=z, w=W) # weighted LS fit of "adjusted response" z on X
# = solve(crossprod(X,W*X), crossprod(X,W*z))
beta = as.matrix(coef(wlmfit),ncol=1) # weighted LS fit
# coefficient update based on quadratic approximation of log likelihood
# using weighted least square regression
eta = offset + X %*% beta # linear predictor = "fitted adjusted responses"
mu = family$linkinv(eta) # predictions on response scale µ = g-1(η)
if(sqrt(crossprod(beta-betaold)) < tol) break
}
# # calculate explained variance on adjusted response scale
# ss_residual = sum(W * (z - eta) ^ 2)
# ss_total = sum(W * (z - weighted.mean(z, W)) ^ 2)
# r.squared = 1 - ss_residual / ss_total
return(list(coefficients=beta, iterations=i))
}
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
IRLS / Fisher scoring implicitly uses the expected Hessian.
Newton-Raphson, by contrast would explicitly form the observed Hessian, which is slower, and would need coefficient starting values.
A minimal implementation in R would be
glm_newton_raphson = function(X, # design matrix
y, # response
weights = rep(1, nrow(X)), # prior observation weights
start = rep(0, ncol(X)), # coefficient starting values
offset = rep(0, nrow(X)), # model offset
family = gaussian(identity), # distribution & link function
maxit = 25,
tol = 1e-08) {
beta = start
nobs = nrow(X)
nvars = ncol(X)
eval(family$initialize) # initializes n and mustart
eta = family$linkfun(mustart) # initialize η = g(µ)
mu = family$linkinv(eta) # predictions on response scale µ
for (i in 1:maxit) {
var = family$variance(mu) # variance as a function of the mean µ
gprime = family$mu.eta(eta) # derivative of the link function with respect to η = d(g^-1)/dη = dμ/dη
gradient = t(X) %*% ((y - mu) * gprime * weights) # gradient of the log-likelihood
W = diag(as.vector(weights*(gprime^2 / var))) # working weights
# Compute the Fisher information matrix = negative of observed Hessian matrix
XWX = t(X) %*% W %*% X
information = XWX # Fisher information matrix = negative of the observed Hessian (second derivatives of the log-likelihood)
betaold = beta
beta = beta + solve(information, gradient) # Newton-Raphson update step
eta = offset + X %*% beta # linear predictor, i.e., predictions on link scale
mu = family$linkinv(eta) # predictions on response scale µ = g^-1(η)
if (sqrt(crossprod(beta - betaold)) < tol) break
}
return(list(coefficients = beta, iterations = i))
}
