# Which optimization algorithm is used in glm function in R?

One can perform a logit regression in R using such code:

> library(MASS)
> data(menarche)
> glm.out = glm(cbind(Menarche, Total-Menarche) ~ Age,
+                                              family=binomial(logit), data=menarche)
> coefficients(glm.out)
(Intercept)         Age
-21.226395    1.631968

It looks like the optimization algorithm has converged - there is information about steps number of the fisher scoring algorithm:

Call:
glm(formula = cbind(Menarche, Total - Menarche) ~ Age, family = binomial(logit),
data = menarche)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-2.0363  -0.9953  -0.4900   0.7780   1.3675

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -21.22639    0.77068  -27.54   <2e-16 ***
Age           1.63197    0.05895   27.68   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 3693.884  on 24  degrees of freedom
Residual deviance:   26.703  on 23  degrees of freedom
AIC: 114.76

Number of Fisher Scoring iterations: 4

I am curious about what optim algorithm it is? Is it Newton-Raphson algorithm (second order gradient descent)? Can I set some parameters to use Cauchy algorithm (first order gradient descent)?

• Do you mind citing where a Newton-Raphson algorithm is referred to as "gradient descent level II"? Commented Oct 24, 2015 at 18:00
• FIsher scoring itself is related to, but different from Newton-Raphson, in effect replacing the Hessian with its expected value under the model. Commented Oct 25, 2015 at 0:10
• @CliffAB sory. I ment that Newton's method is a second order gradient descent method. Commented Oct 25, 2015 at 15:21
• @hxd1011, you should not edit to change someone else's question. It is one thing to edit when you think you know what they mean, but their question is unclear (perhaps English isn't their native language, eg), but you should not make their question different (eg, more general) than they had wanted. Instead, ask a new question with what you want. I am rolling the edit back. Commented Jul 8, 2016 at 21:07
• @MarcinKosiński Newton's method is Newton-Raphson, Raphson merely built on Newton's ideas for a more general case. Commented Jul 8, 2016 at 22:27

You will be interested to know that the documentation for glm, accessed via ?glm provides many useful insights: under method we find that iteratively reweighted least squares is the default method for glm.fit, which is the workhorse function for glm. Additionally, the documentation mentions that user-defined functions may be provided here, instead of the default.

• You can also just type the function name glm or fit.glm at the R prompt to study the source code. Commented Oct 24, 2015 at 19:03
• @MatthewDrury I think you mean the workhorse glm.fit which will not be entirely reproducible since it relies on C code C_Cdqrls. Commented Jul 8, 2016 at 17:39
• Yah, you're right, i mix up the order in R a lot. What do you mean not reproducible though? Commented Jul 8, 2016 at 17:45

The method used is mentioned in the output itself: it is Fisher Scoring. This is equivalent to Newton-Raphson in most cases. The exception being situations where you are using non-natural parameterizations. Relative risk regression is an example of such a scenario. There, the expected and observed information are different. In general, Newton Raphson and Fisher Scoring give nearly identical results.

Another formulation of Fisher Scoring is that of Iteratively Reweighted Least Squares. To give some intuition, non-uniform error models have the inverse variance weighted least squares model as an "optimal" model according to the Gauss Markov theorem. With GLMs, there is a known mean-variance relationship. An example is logistic regression where the mean is $p$ and the variance is $p(1-p)$. So an algorithm is constructed by estimating the mean in a naive model, creating weights from the predicted mean, then re-estimating the mean using finer precision until there is convergence. This, it turns out, is Fisher Scoring. Additionally, it gives some nice intuition to the EM algorithm which is a more general framework for estimating complicated likelihoods.

The default general optimizer in R uses numerical methods to estimate a second moment, basically based on a linearization (be wary of curse of dimensionality). So if you were interested in comparing the efficiency and bias, you could implement a naive logistic maximum likelihood routine with something like

set.seed(1234)
x <- rnorm(1000)
y <- rbinom(1000, 1, exp(-2.3 + 0.1*x)/(1+exp(-2.3 + 0.1*x)))
f <- function(b) {
p <- exp(b[1] + b[2]*x)/(1+exp(b[1] + b[2]*x))
-sum(dbinom(y, 1, p, log=TRUE))
}
ans <- nlm(f, p=0:1, hessian=TRUE)

gives me

> ans$estimate [1] -2.2261225 0.1651472 > coef(glm(y~x, family=binomial)) (Intercept) x -2.2261215 0.1651474 • what's the difference comparing to a gradient decent/newton's method/ BFGS? I think you explained, but I am not quite follow. Commented Jul 8, 2016 at 17:40 • @hxd1011 see "Numerical Methods for Unconstrained Optimization and Nonlinear Equations" Dennis, J. E. and Schnabel, R. B. Commented Jul 8, 2016 at 17:43 • @hxd1011 as far as I can tell, Newton Raphson does not require or estimate a Hessian in the steps. The Newton-Type method in nlm estimates the gradient numerically then applies Newton Raphson. In BFGS I think the gradient is required as with Newton Raphson, but successive steps are evaluated using a second order approximation which requires an estimate of the Hessian. BFGS is good for highly nonlinear optimizations. But for GLMs, they are usually very well behaved. Commented Jul 8, 2016 at 17:58 • The existing answer mentioned "iteratively reweighted least squares", how does that enter the picture, compared to the algorithms you mentioned here? Commented Jul 8, 2016 at 23:58 • @AdamO could you address amoeba's comments? Thanks Commented Jul 11, 2016 at 15:07 For the record, a simple pure R implementation of R's glm algorithm, based on Fisher scoring (iteratively reweighted least squares), as explained in the other answer is given by: 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
# = linearised version of log-likelihood function ℓ around η
W        = weights * as.vector(gprime^2 / var)
# = working weights
betaold  = beta
beta     = solve(crossprod(X,W*X), crossprod(X,W*z))
# coefficient update based on quadratic approximation of log likelihood
# using weighted least square regression = as.matrix(coef(lm.wfit(x=X, y=z, w=W)), ncol=1)
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)) } Example: ## Dobson (1990) Page 93: Randomized Controlled Trial : y <- counts <- c(18,17,15,20,10,20,25,13,12) outcome <- gl(3,1,9) treatment <- gl(3,3) X <- model.matrix(counts ~ outcome + treatment) coef(glm.fit(x=X, y=y, family = poisson(log))) (Intercept) outcome2 outcome3 treatment2 treatment3 3.044522e+00 -4.542553e-01 -2.929871e-01 -7.635479e-16 -9.532452e-16 coef(glm_irls(X=X, y=y, family=poisson(log))) (Intercept) outcome2 outcome3 treatment2 treatment3 [1,] 3.044522 -0.4542553 -0.2929871 -3.151689e-16 -8.24099e-16 A good discussion of GLM fitting algorithms, including a comparison with Newton-Raphson (which uses the observed Hessian as opposed to the expected Hessian in the IRLS algorithm) and hybrid algorithms (which start with IRLS, as these are easier to initialize, but then finish with further optimization using Newton-Raphson) can be found in the book "Generalized Linear Models and Extensions" by James W. Hardin & Joseph M. Hilbe. A pure R implementation of Newton-Raphson to fit a GLM 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))) # diagonal weight matrix

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

Newton-Raphson has the disadvantage that it calculates the full observed Hessian matrix, which is slow & that it needs coefficient starting values, which IRLS / Fisher scoring doesn't need.