Maximum likelihood in a GLM

Section of Introduction to generalised linear models, part of Generalised linear models

Maximum likelihood in a GLM

Type Section
Taught in Introduction to generalised linear models
Part Generalised linear models
Questions using it 1

For a single Bernoulli trial the likelihood contribution is p^{y}(1-p)^{1-y}, with y\in\{0,1\}. Independence multiplies those contributions. The log-likelihood for a sample of n is therefore \begin{align*} \ell(\boldsymbol{\beta}) &=\sum_{i=1}^n \Bigl\{y_i\log p_i+(1-y_i)\log(1-p_i)\Bigr\}, \end{align*} where p_i=\mu({\boldsymbol{X}}_i^\top\boldsymbol{\beta}) depends on \boldsymbol{\beta} through whichever inverse link we chose. The binomial coefficient, if you are in fact looking at a binomial count rather than a single trial, does not depend on \boldsymbol{\beta} and does not change the maximiser.

Differentiate with respect to \boldsymbol{\beta}. The chain rule runs through p_i and then through the linear predictor. Write \mu'(\eta) for the derivative of the mean function with respect to \eta={\boldsymbol{X}}^\top\boldsymbol{\beta}, so \partial p_i/\partial\boldsymbol{\beta}=\mu'({\boldsymbol{X}}_i^\top\boldsymbol{\beta})\,{\boldsymbol{X}}_i. Then \begin{align*} \frac{\partial\ell}{\partial\boldsymbol{\beta}} &=\sum_{i=1}^n\left\{\frac{y_i}{p_i}-\frac{1-y_i}{1-p_i}\right\}\mu'({\boldsymbol{X}}_i^\top\boldsymbol{\beta})\,{\boldsymbol{X}}_i\\ &=\sum_{i=1}^n\frac{y_i-p_i}{p_i(1-p_i)}\,\mu'({\boldsymbol{X}}_i^\top\boldsymbol{\beta})\,{\boldsymbol{X}}_i. \end{align*} Set the score equal to zero. The denominator p_i(1-p_i) is \mathrm{Var}(Y_i\mid {\boldsymbol{X}}_i). The same algebra, with the variance of the family in that slot, is the score equation for every GLM we have written down: \sum_{i=1}^n\frac{y_i-\mu_i}{\mathrm{Var}(Y_i\mid {\boldsymbol{X}}_i)}\,\mu'({\boldsymbol{X}}_i^\top\boldsymbol{\beta})\,{\boldsymbol{X}}_i={\boldsymbol{0}}. You should be able to recover Equation (11) in Introduction to generalised linear models from the Poisson log-likelihood and from the gamma log-likelihood. That is an exercise, not a new idea. For the normal identity-link model, \mu_i={\boldsymbol{X}}_i^\top\boldsymbol{\beta}, the derivative \mu' is 1, and the variance is the constant \sigma^2. Equation Equation (11) in Introduction to generalised linear models collapses to \frac{1}{\sigma^2}\sum_{i=1}^n (y_i-{\boldsymbol{X}}_i^\top\boldsymbol{\beta})\,{\boldsymbol{X}}_i={\boldsymbol{0}}, or {\color{blue}\boldsymbol{X}}^\top({\boldsymbol{y}}-{\color{blue}\boldsymbol{X}}\boldsymbol{\beta})={\boldsymbol{0}} once the rows {\boldsymbol{X}}_i^\top are stacked into the matrix {\color{blue}\boldsymbol{X}}. Those are the normal equations. The GLM score is the same residual, {\boldsymbol{y}}-\boldsymbol{\mu}, reweighted by the variance and by the slope of the mean function.

The weights depend on \boldsymbol{\beta}, and so does \boldsymbol{\mu}. Equation Equation (11) in Introduction to generalised linear models is not a linear system. We solve it by iteration.

Fisher scoring

Main article: Fisher scoring

Fix a current value \boldsymbol{\beta}^{(0)}. Evaluate the mean, the variance and the slope of the mean at that value, and write \mu_i^{(0)}, V_i^{(0)} and \mu'_{(0),i} for those three numbers at observation i. The only unknown left in Equation (11) in Introduction to generalised linear models is the next \boldsymbol{\beta}. Replace the mean by its tangent at the current linear predictor, \mu({\boldsymbol{X}}_i^\top\boldsymbol{\beta}) \approx \mu_i^{(0)}+\mu'_{(0),i}\,{\boldsymbol{X}}_i^\top\bigl(\boldsymbol{\beta}-\boldsymbol{\beta}^{(0)}\bigr). Substitute that approximation into Equation (11) in Introduction to generalised linear models. The equation becomes linear in the unknown step \boldsymbol{\beta}-\boldsymbol{\beta}^{(0)}: \sum_{i=1}^n\frac{y_i-\mu_i^{(0)}}{V_i^{(0)}}\,\mu'_{(0),i}\,{\boldsymbol{X}}_i = \sum_{i=1}^n\frac{\bigl(\mu'_{(0),i}\bigr)^2}{V_i^{(0)}}\,{\boldsymbol{X}}_i{\boldsymbol{X}}_i^\top\bigl(\boldsymbol{\beta}-\boldsymbol{\beta}^{(0)}\bigr). Name the two sides. The working score and the working information at \boldsymbol{\beta}^{(0)} are \begin{align*} {\boldsymbol{U}}^{(0)} &=\sum_{i=1}^n\frac{y_i-\mu_i^{(0)}}{V_i^{(0)}}\,\mu'_{(0),i}\,{\boldsymbol{X}}_i,\\ {\color{blue}\boldsymbol{I}}^{(0)} &=\sum_{i=1}^n\frac{\bigl(\mu'_{(0),i}\bigr)^2}{V_i^{(0)}}\,{\boldsymbol{X}}_i{\boldsymbol{X}}_i^\top. \end{align*} The update is one step of weighted least squares, \boldsymbol{\beta}^{(1)}=\boldsymbol{\beta}^{(0)}+\bigl({\color{blue}\boldsymbol{I}}^{(0)}\bigr)^{-1}{\boldsymbol{U}}^{(0)}. Iterate until \boldsymbol{\beta} stops moving. The matrix {\color{blue}\boldsymbol{I}} is the expected negative Hessian of \ell, which is the Fisher information. The algorithm is Newton–Raphson with that expected second derivative in place of the observed one, and for that reason it is called Fisher scoring. It is the same algorithm for every exponential family we put in the random component. A model outside that class does not inherit it; you would write a separate optimiser.

For the binomial this rarely takes more than five or six iterations, and the start does not matter much. The cholesterol fit below takes four. A first reading does not require the iteration by hand; working a few steps by hand is still a reasonable exercise. Either way, the update is Equation (12) in Introduction to generalised linear models, and a GLM fit is that update and not a closed form.

What the software prints

Main article: Reading glm output in R

After the score equation, the program is only evaluating Equation (12) in Introduction to generalised linear models. For a binary disease column, a sex indicator, and cholesterol,

chol_file <- "data/cholest.csv"
if (!file.exists(chol_file)) {
  stop("Could not find data/cholest.csv. Set the working directory to the folder that contains it.")
}
chol <- read.csv(chol_file)
fit <- glm(disease ~ gender + cholesterol, data = chol, family = binomial)
fit
fit$iter

Naming binomial and not naming a link selects the logit. The coefficients printed for this fit are the ones already interpreted in Equation (10) in Introduction to generalised linear models:

Coefficients:
(Intercept)     gender   cholesterol
    -9.3203    -0.1094        1.5842

Degrees of Freedom: 99 Total (i.e. Null);  97 Residual
Null Deviance:      137.6
Residual Deviance:  114     AIC: 120

One hundred patients and three coefficients leave 97 residual degrees of freedom. The deviance and the AIC are on this printout because the software always writes them. We will say what they measure when we come to goodness of fit. The line that belongs to the algorithm above is the iteration count: four Fisher-scoring steps.

Estimating \boldsymbol{\beta} is the less interesting half of a statistical analysis. The half that matters next is the precision of those three numbers, and the test of whether sex can be removed. That is the next chapter. The fit above is the point at which this chapter stops.

  1. Show that the exponential density f(y;\lambda)=\lambda e^{-\lambda y}, for y>0 and \lambda>0, is an exponential family. Identify a(y), the natural parameter, and the dispersion.

  2. Take the gamma density, with \phi fixed, written as f(y;\lambda)=\frac{y^{\phi-1}\lambda^{\phi}e^{-\lambda y}}{\Gamma(\phi)}. Here \lambda is sometimes called a scale parameter. The formula is the shape-and-rate density: \lambda is a rate, and the scale is 1/\lambda. Using the density as written, show that it is an exponential family in \lambda with \phi held fixed, and identify the natural parameter. Then match it to the mean-and-variance form used in the clotting example.

  3. With r known, show that f(y;\theta)=\binom{y+r-1}{r-1}\theta^{r}(1-\theta)^{y}, \qquad y=0,1,2,\ldots, is an exponential family in \theta. State the natural statistic and the natural parameter. (This is a negative binomial: the number of failures before r successes.)

  4. Repeat the differentiation in the proof of Note for the general form Equation (6) in Introduction to generalised linear models, in which a(y) and b(\theta) need not be the identity. The mean is no longer -c'(\theta).

  5. Show that Equation (9) in Introduction to generalised linear models maps every real linear predictor to a probability in (0,1).

  6. Starting from the Poisson log-likelihood, and then from the gamma log-likelihood in the mean-and-variance parameterisation, show that the score equation for \boldsymbol{\beta} has the form Equation (11) in Introduction to generalised linear models.

Used in

See also

Exponential family · Mean and variance of an exponential family · Generalised linear model · Canonical link · What a regression statement is hiding · Three examples the linear model does not fit · Runs in the past week · Clotting-time data

Glossary and index