26
$\begingroup$

I'm curious about why we treat fitting GLMS as though they were some special optimization problem. Are they? It seems to me that they're just maximum likelihood, and that we write down the likelihood and then ... we maximize it! So why do we use Fisher scoring instead of any of the myriad of optimization schemes that has been developed in the applied math literature?

$\endgroup$
1
  • 1
    $\begingroup$ As far as I understand it's got to do with the fact that the algorithm based on Fisher scoring (which uses the expected Hessian) does not need starting estimates of your coefficient vector, unlike regular Newton-Raphson (which uses the observed Hessian), which does... This makes Fisher scoring much easier to use. But some use hybrid algorithms, starting with IRLS and then switching to Newton-Raphson. See section 3.4 in the book of Hardin & Hilbe, gen.lib.rus.ec/… $\endgroup$ Commented Jan 21, 2020 at 16:03

3 Answers 3

15
$\begingroup$

Fisher's scoring is just a version of Newton's method that happens to be identified with GLMs, there's nothing particularly special about it, other than the fact that the Fisher's information matrix happens to be rather easy to find for random variables in the exponential family. It also ties in to a lot of other math-stat material that tends to come up about the same time, and gives a nice geometric intuition about what exactly Fisher information means.

There's absolutely no reason I can think of not to use some other optimizer if you prefer, other than that you might have to code it by hand rather than use a pre-existing package. I suspect that any strong emphasis on Fisher scoring is a combination of (in order of decreasing weight) pedagogy, ease-of-derivation, historical bias, and "not-invented-here" syndrome.

$\endgroup$
5
  • 2
    $\begingroup$ I don't think this is quite correct - the IRLS algorithm uses the expected Hessian, whereas Newton-Raphson uses the observed Hessian - see gen.lib.rus.ec/… for a detailed comparison of the 2 algorithms... $\endgroup$ Commented Aug 27, 2019 at 20:58
  • $\begingroup$ @TomWenseleers Could you perhaps elaborate in an answer? Does this mean the algorithmic complexity of beta regression is not the reason it is treated as a separate issue from GLMs? $\endgroup$ Commented Jan 21, 2020 at 2:29
  • 1
    $\begingroup$ @Frans Rodenburg Not so well into beta regression, but I believe that the standard IRLS method only works for single-parameter distributions from the exponential family, whereas beta regression is a 2-parameter exponential distribution... See stats.stackexchange.com/questions/304538/… Cox proportional hazard and negative binomial also each have an additional parameter though and they can be fit using a modified IRLS algo, so not sure really... $\endgroup$ Commented Jan 21, 2020 at 11:17
  • 1
    $\begingroup$ Other advantage of using Fisher scoring / IRLS with the expected Hessian btw is that the algo is much easier to initialize - see section 3.4 in Hardin & Hilbe's book. This contrasts with Newton Raphson where you need to have an initial guess of the coefficient vector, which is a little difficult... Sometimes people then use hybrid algorithms, and starts with IRLS with Fisher scoring and then after some iterations switch to regular Newton Raphson... $\endgroup$ Commented Jan 21, 2020 at 11:21
  • 3
    $\begingroup$ Another good point about Fisher scoring is that the expected Fisher information is always positive (semi-)definite, whereas the second derivative of the loglikelihood need not be. For typical GLMs this isn't a big issue, but for parametric survival models there is a real problem that the second derivative need not be positive semidefinite. $\endgroup$ Commented Jun 18, 2020 at 9:21
13
$\begingroup$

It's historical, and pragmatic; Nelder and Wedderburn reverse-engineered GLMs, as the set of models where you can find the MLE using Fisher scoring (i.e. Iteratively ReWeighted Least Squares). The algorithm came before the models, at least in the general case.

It's also worth remembering that IWLS was what they had available back in the early 70s, so GLMs were an important class of models to know about. The fact you can maximize GLM likelihoods reliably using Newton-type algorithms (they typically have unique MLEs) also meant that programs like GLIM could be used by those without skills in numerical optimization.

$\endgroup$
7
  • 1
    $\begingroup$ I don't think this is quite correct - the IRLS algorithm uses the expected Hessian, whereas Newton-Raphson uses the observed Hessian - see gen.lib.rus.ec/… for a detailed comparison of the 2 algorithms... $\endgroup$ Commented Aug 27, 2019 at 20:59
  • $\begingroup$ @TomWenseleers I've been wondering about it as well, would be cool to have something on this $\endgroup$
    – Firebug
    Commented Dec 9, 2020 at 19:45
  • 1
    $\begingroup$ @Firebug You can take a look in Hardin & Hilbe's book "Generalized Linear Models and Extensions", siteget.net/…, (Sections 3.4 & Sections 3.6 & 3.7, Listing 3.1 & 3.2 and Section 5.6, Listing 5.4 compare the IRLS algorithms & Newton-Raphson). The IRLS algo has the advantage of being much easier to initialise, $\endgroup$ Commented Dec 9, 2020 at 21:30
  • 1
    $\begingroup$ @Firebug That's because the IRLS algo you can initialize with an initial guess of the fitted values (which is easy), whereas Newton-Raphson you have to initialize with initial guesses of the coefficients to be estimated (which is very hard). Sometimes people use combined algorithms, where they start with IRLS, and then once they have good good coefficients estimates switch to Newton-Raphson... $\endgroup$ Commented Dec 9, 2020 at 21:32
  • 1
    $\begingroup$ @Firebug For a simple bare-bones IRLS implementation see bwlewis.github.io/GLM. $\endgroup$ Commented Dec 9, 2020 at 22:29
4
$\begingroup$

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 / Fisher scoring 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.

In pure mathematical terms & somewhat schematically, (from this Hardin & Hilbe book), the IRLS / Fisher scoring and Newton-Raphson algorithms are structured as follows: enter image description here

Code-wise, a simple pure R implementation of R's glm algorithm, based on IRLS / Fisher scoring 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 = NULL, # coefficient starting values, here optional
             offset = rep(0, nrow(X)), # model offset 
             family = gaussian(identity), # distribution & link function
             maxit = 25, 
             tol = 1e-08)
    {
      if (!is.null(start[[1]])) beta = start # no coefficient starting values needed
      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))
    }

A Newton-Raphson algorithm to fit a GLM, by contrast, would use something like

    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, here required
                                  offset = rep(0, nrow(X)), # model offset
                                  family = gaussian(identity), # distribution & link function
                                  maxit = 25,
                                  tol = 1e-08) {
    
      beta = start # coefficient starting values, here required
      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))
    }

Conceptually, the popularity and advantages of the IRLS algorithm are (1) that one merely has to perform a weighted least square regression in every update step (for which any range of state of the art least squares solvers for dense or sparse systems can be used, e.g. the Eigen solvers), (2) that the algorithm does not need coefficient starting values, unlike Newton-Raphson, and (3) that unlike Newton-Raphson, the IRLS algorithm does not form the full Hessian matrix (implicitly it uses the expected Hessian, but it doesn't actually form it; by contrast, Newton-Raphson forms the observed Hessian in every iteration, which is slow & inefficient). I haven't seen greatly improved GLM algorithms using more fancy state of the art solvers to be honest - it is hard to beat the best least square solvers once the algorithm is cast to merely use weighted least squares... IRLS using Eigen's least squares conjugate gradient solver for dense systems and Eigen's Cholesky LLT solver for sparse systems is very efficient e.g. In H20, they implement a couple of other solvers to fit GLMs, including LBFGS and coordinate descent, and claim LBFGS is more efficient for problems with very many covariates. But I am not sure what least square solver they use in the IRLS method (definitely not least square conjugate gradient from the description, which would be the fastest) & this might already make a huge difference.

$\endgroup$
5
  • $\begingroup$ I believe that this answer would be nicer if it was formulated in equations rather than computer code. It obscures the differences. $\endgroup$ Commented Jun 27 at 9:34
  • $\begingroup$ Dunno - I usually like to see an actual implementation - I always get annoyed when I see articles without any code, because small differences in implementation (e.g. what least squares solver is used etc) often make a huge difference in performance, and the question above is specifically asking what the advantages of that algo are over others. That one can readily use any state of the art least square solver is one advantage... But feel free to make any edits you like - e.g. adding the algorithmic details in purely mathematical format... $\endgroup$ Commented Jun 27 at 9:45
  • $\begingroup$ I normally like an implementation as well and sometimes it is easier to just write a code that shows what happens rather than writing it all down in detail with latex formulas. But the code here is relatively large and if you don't know what too look for, one is not able to see differences. A reader has to first understand what happens in the codes and there is not a direct explanation what the difference is. I see two times the use of a function solve(A,X) but the difference is not directly clear. I have to look back and forth between the two and there are many differences. $\endgroup$ Commented Jun 27 at 10:39
  • 1
    $\begingroup$ I added the pseudocode from the Hardin & Hilbe book - maybe that's clearer... Feel free to make any edits if you think some things are not clear... $\endgroup$ Commented Jun 27 at 12:14
  • 1
    $\begingroup$ that's very great. I might try to transform that image into text when I have time for it. Remind me after a few days when I didn't. $\endgroup$ Commented Jun 27 at 12:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.