# Are GLMs just glorified WLS regressions?

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]$$.

• Not, as it uses IRLS algorithm that can be used with different link-functions Commented Mar 20 at 13:02
• I'm curious how you arrive at calling expression 3 in your question a "weighted least squares" gradient. Surely, you're referencing page 25 of Turner's notes which point out the algorithmic approach to solving GLMs. But the general method of WLS doesn't require any expression for the weight - so $V(\mu)$ is endemic just to GLMs whereas WLS uses $W$ as a matrix to be practically anything (symmetric positive def). It may also be of interest to know that a "one step" estimator completes just one NR iteration and approximates the GLM solution but is in fact a WLS estimator. Commented Apr 17 at 16:45
• @AdamO This was 5 years ago, so I do not really remember what I was thinking at that time. But looking at it again now, for the choices $r_i ≝ y_i- μ_i = y_i-g^{-1}(x_i^𝖳β)$ and $w_i ≝ \frac{1}{σ_i^2} = \frac{1}{\operatorname{Var[y_i]}} = \frac{1}{\operatorname{V[μ_i]}}$ both expressions $∑_i \frac{y_i-μ_i}{V[μ_i]} \frac{x_i}{g'(μ_i)}$ and $∑_i \frac{r_i}{σ_i^2}∇_β r_i$ are identical. So, GLMs are engineered in a way so that the gradient of the LL can always be interpreted as a WLS Gradient, which is not the case for general probabilistic models. Commented Apr 17 at 20:28
• $$\frac{\partial\ell}{\partial\beta} = \text{an expression in which “\beta ” does not appear?}$$ Where is $\beta$ on the right side of this equality? Commented Jun 22 at 20:44
• Might the difference be that in GLMs the weights depend on the data, whereas in WLS they don't? Commented Jun 22 at 20:46

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.

• I do not see how the method of optimization is of any relevance here. I was simply observing that the gradient of the log-likelihood of a GLM can be interpreted as the gradient of a WLS scheme, indicating that both are connected. Commented Apr 29, 2019 at 11:07
• 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 GLM's conceptually is the same as WLS. Commented Apr 29, 2019 at 11:10
• @Hyperplane For regular GLMs like logistic and poisson regression, the method of IRLS is a kind of Newton Raphson procedure for the log likelihood, solve the root of the log-likelihood, obtain an extremum of the likelihood. Regular GLMs without boundary solutions have concave unimodal likelihoods. It breaks down as you consider other edge cases. Commented Apr 17 at 16:39
• @Hyperplane the problem is that which WLS scheme the GLM is equivalent to depends on where the gradient is being evaluated. If the GLM's gradient agreed with the gradient of a fixed WLS everywhere (i.e. for all $\beta$), then they would be equivalent. Commented Apr 17 at 16:50

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))
}


I think, I see your misinterpretation: "GLM is being calculated using IRLS algorithm with WLS-gradient". In general, WLS is a generalization of OLS for dataset with heteroscedasticity, but is specialization of GLM, when all the off-diagonal entries of the covariance matrix of the errors, are null, as wiki explains good

N.B.: best description of GLM, depending on link-function it can be either OLS, or WLS, or any other estimator, because can be used with other Maximum-Likelihood-algorithms (e.g. for Logistic Regression) , OLS & WLS are just special cases of MLE