Which optimization algorithm is used in glm function in R?

Question

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

Answer 1

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.

Answer 2

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 

Answer 3

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 
             # adjusted response 
             # = 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.