# In GLMs, why do we solve score(beta)=0 instead of just minimizing the negative log-likelihood?

When we search for a numerical way to find $\hat{\beta}$ in a GLM (say, a logistic regression), we could do a numerical optimization (minimization) of the negative log-likelihood.

But instead, we go one step further, and first compute the score function, to then equvialently look for the value of $\beta$ that sets $s(\beta) = 0$.

Why do we do this extra step? They are both iterative numerical procedures. Is solving this equation system computationally easier than minimizing the negative log-likelihood?

• You seem to be hinting that there's a special (even unique) way of minimizing the negative log-likelihood ("they are both numerical procedures"). What particular optimization method did you have in mind? – Glen_b -Reinstate Monica Dec 18 '17 at 10:21
• While thinking about my answer, I realized some of these methods need to compute the derivative as well. So to keep my question valid, can I answer "the Nelder-Mead simplex method"? – Alexander Engelhardt Dec 18 '17 at 10:48
• Maybe because the variance of the score has other important roles to play in inference, see the variance covariance matrix. – Christoph Hanck Dec 18 '17 at 15:10
• @Alexx One obvious answer for that option is that the simplex method is generally a lot slower than Fisher scoring. If you do use derivatives (which may make speed more comparable), Christoph's answer would still come into play -- but with really large data sets you may want to consider other possibilities. – Glen_b -Reinstate Monica Dec 18 '17 at 23:51
• Finding the root of the score equation is a second-order optimization algorithm. Newton Raphson or Fisher Scoring both arrive at the maximum likelihood estimate and in far fewer iterations than any quasi-Newton procedure like BFGS. Exponential families are particularly nice for these kinds of estimation routines because of how well behaved they are. – AdamO Jan 29 '18 at 22:13

The solutions to the score equation are critical points of the objective function in your optimisation, so generally the fitted coefficient estimator should solve the score equation. This is not really a "step further" than using numerical techniques; it is just a way of characterising the fitted value, which is the point that numerical solutions should be moving towards.

Optimisation via IRLS: As a practical matter, numerical algorithms to fit GLMs generally use the default method of iteratively-reweighted least squares optimisation, which I think is equivalent to Fisher scoring for this class of models. This is the default fitting-method in the R function glm (see comments on the method variable of the function in the documentation). The estimation method can be derived via an argument involving linear approximation to the score equation, and linear approximation to the mean parameter in each iteration.

The GLM is characteristed by having some exponential family distribution with scale parameter $\phi$ and mean parameter $\mu = h (\eta) = h(\boldsymbol{\text{x}} \beta)$ where the function $h$ is the canonical inverse-link and $\eta$ is a linear function of the coefficients of interest. Using this distribution and taking a linear approximation to the score function yields:

$$s(\beta) = \frac{\partial \ell}{\partial \beta} (\beta) \approx \phi^{-1} \boldsymbol{\text{x}}^{\text{T}}(\boldsymbol{\text{y}} - \mu).$$

For some iteration value $\beta_{(k)}$ we can form the corresponding mean $\mu_{(k)} = h(\eta_{(k)}) = h(\boldsymbol{\text{x}} \beta_{(k)})$ and a linear approximation gives $\mu \approx \mu_{(k)} + \boldsymbol{\text{W}}_{(k)} \boldsymbol{\text{x}} (\beta - \beta_{(k)})$ with weight matrix $\boldsymbol{\text{W}}_{(k)} = \nabla_\eta h (\eta_{(k)})$. Now, if we substitute this latter approximation into the approximation for the score equation we get:

\begin{aligned} s(\beta) = \frac{\partial \ell}{\partial \beta} (\beta) \approx \phi^{-1} \boldsymbol{\text{x}}^{\text{T}}(\boldsymbol{\text{y}} - \mu) &= \phi^{-1} \boldsymbol{\text{x}}^{\text{T}}(\boldsymbol{\text{y}} - \mu_{(k)} - (\mu - \mu_{(k)})) \\[8pt] &\approx \phi^{-1} \boldsymbol{\text{x}}^{\text{T}}(\boldsymbol{\text{y}} - \mu_{(k)} - \boldsymbol{\text{W}}_{(k)} \boldsymbol{\text{x}} (\beta - \beta_{(k)})) \\[8pt] &= \phi^{-1} \boldsymbol{\text{x}}^{\text{T}} \boldsymbol{\text{W}}_{(k)} ( \boldsymbol{\text{W}}_{(k)} ^{-1} (\boldsymbol{\text{y}} - \mu_{(k)}) + \boldsymbol{\text{x}} \beta_{(k)} - \boldsymbol{\text{x}} \beta ) \\[8pt] &= \phi^{-1} \boldsymbol{\text{x}}^{\text{T}} \boldsymbol{\text{W}}_{(k)} ( \boldsymbol{\text{z}}_{(k)} - \boldsymbol{\text{x}} \beta ), \\[8pt] \end{aligned}

where $\boldsymbol{\text{z}}_{(k)} = \boldsymbol{\text{W}}_{(k)} ^{-1} (\boldsymbol{\text{y}} - \mu_{(k)}) + \boldsymbol{\text{x}} \beta_{(k)}$ is the adjusted response, which is calculable from the $k$th iteration. This approximation to the score function corresponds to a weighted linear regression. Solving the score equation we obtain the next iteration for the coefficient vector, which has the form of a weighted-least-squares estimator:

$$\hat{\beta}_{(k+1)} = (\boldsymbol{\text{x}}^\text{T} \boldsymbol{\text{W}}_{(k)} \boldsymbol{\text{x}})^{-1} (\boldsymbol{\text{x}}^\text{T} \boldsymbol{\text{W}}_{(k)} \boldsymbol{\text{z}}_{(k)}).$$

This method is a simple variation on the method of least-squares estimation used in multiple linear regression with Gaussian errors. Generally these algorithms start with a base iteration using least-squares estimation and then iterate the IRLS algorithm to within some tolerance of the true fitted value. Note that even once the fitted value is obtained via IRLS, the score function is still useful, since it is related to the variance of the coefficient vector.