faraway package contains the datasets in the ELMR book.
1 Introduction
Now we have learnt regression modeling of binomial, multinomial, and count responses. How about nonnegative real responses and so on? Generalized linear model (GLM) is a generic framework that encompasses normal regression, binomial regression, multinomial regression, and others. There are two essential components of the GLM framework: a distribution for response \(Y\) and a link function that relates mean of \(Y\) to covariates \(x\).
2 Exponential family distribution
2.1 Definition
In GLM, the distribution of \(Y\) is from the exponential familty of distributions of form \[
f(y \mid \theta, \phi) = \exp \left[ \frac{y \theta - b(\theta)}{a(\phi)} + c(y, \phi) \right].
\]\(\theta\) is called the canonical parameter or natural parameter and represents the location while \(\phi\) is the dispersion parameter and represents the scale. Note the canonical parameter \(\theta\) is not necessarily the mean \(\mu\).
Gamma has density \[
f(y \mid \nu, \lambda) = \frac{1}{\Gamma(\nu)} \lambda^{\nu} y^{\nu - 1} e^{-\lambda y}, \quad y > 0,
\] where \(\nu\) is the shape parameter and \(\lambda\) is the scale parameter. For the purpose of GLM, it’s convenient to reparameterize by \(\lambda = \nu / \mu\) to get \[
f(y) = \frac{1}{\Gamma(\nu)} \left( \frac{\nu}{\mu} \right)^{\nu} y^{\nu - 1} e^{-y\nu / \mu} = \exp \left\{ \frac{- y \mu^{-1} - \log \mu}{\nu^{-1}} + (\nu-1) \log y + \nu \log \nu - \log \Gamma(\nu) \right\}.
\] Now \(\mathbb{E}Y = \mu\) and \(\operatorname{Var}(Y) = \mu^2 / \nu = (\mathbb{E} Y)^2 / \nu\). So we have \[\begin{eqnarray*}
\theta &=& - \mu^{-1} \\
\phi &=& \nu^{-1} \\
a(\phi) &=& \phi \\
b(\theta) &=& - \log (- \theta) \\
c(y, \phi) &=& (\phi^{-1} - 1) \log y - \phi^{-1} \log (\phi) - \log \Gamma(\phi^{-1}).
\end{eqnarray*}\] Some books remove the minus sign in the canonical parameter/link which is fine provided we take account of this in any derivations. For the canonical link \(\eta = \mu^{-1}\), the systematic component can only be non-negative, which could cause problems. Other possible link functions are log link \(\eta = \log \mu\) and identity link \(\eta = \mu\).
Many other distributions.
2.3 Moments
Exponential family distributions have mean and variance \[\begin{eqnarray*}
\mathbb{E}Y &=& \mu = b'(\theta) \\
\operatorname{Var}Y &=& \sigma^2 = b''(\theta) a(\phi).
\end{eqnarray*}\]Show this in HW4. Thus the function \(b\) determines the moments of \(Y\).
3 Link function
Given the linear predictor or systematic component\[
\eta = \mathbf{x}^T \boldsymbol{\beta},
\] the link function, \(g\), relates the mean \(\mathbb{E} Y = \mu\) to the covariates \[
\eta = g(\mu).
\]
In principal, any monotone, continuous, and differentiable function can be a link function. But there are some convenient and common choices for the standard GLMs.
For Gaussian linear model, the identity link, \(\eta = \mu\), is the obvious choice.
For binomial model, we saw logit, probit and cloglog.
For Poisson model, a standard choice is \(\eta = \log \mu\).
The canonical link has \(g\) such that \(\eta = g(\mu) = \theta\), the canonical parameter of the exponential family distribution. This means that \(g(b'(\theta))=\theta\). If a canonical link is used, \(\mathbf{X}^T \mathbf{y}\) is sufficient for \(\boldsymbol{\beta}\). The canonical link is mathematically and computationally convenient and is often the natural choice of link.
Family
Canonical Link
Variance Function
Normal
\(\eta=\mu\)
1
Poisson
\(\eta=\log \mu\)
\(\mu\)
Bernoulli
\(\eta=\log \left( \frac{\mu}{1 - \mu} \right)\)
\(\mu (1 - \mu)\)
Gamma
\(\eta = \mu^{-1}\)
\(\mu^2\)
Inverse Gaussian
\(\eta = \mu^{-2}\)
\(\mu^3\)
4 Fisher scoring algorithm and IRWLS
GLM regreesion coefficients are estimated by MLE. Recall that the Newton-Raphson algorithm for maximizing a log-likelihood \(\ell(\boldsymbol{\beta})\) proceeds as \[
\boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} + s [- \nabla^2 \ell(\boldsymbol{\beta}^{(t)})]^{-1} \nabla \ell(\boldsymbol{\beta}^{(t)}),
\] where \(s>0\) is a step length, \(\nabla \ell\) is the score (gradient) vector, and \(-\nabla^2 \ell\) is the observed information matrix (negative Hessian). See Biostat 257 lecture notes
For GLMs with canonical links, the second term and third term cancel using the fact \[
\frac{d\mu_i}{d \eta_i} = \frac{d\mu_i}{d \theta_i} = \frac{d \, b'(\theta_i)}{d \theta_i} = b''(\theta_i) = \frac{\sigma_i^2}{a(\phi)}.
\] Therefore for canonical link the negative Hessian is positive semidefinte and Newton’s algorithm with line search is stable.
How about non-canonical link? We use the expected (Fisher) information matrix \[
\mathbb{E} [- \nabla^2 \ell(\boldsymbol{\beta})] = \sum_{i=1}^n \frac{[\mu_i'(\eta_i)]^2}{\sigma_i^2} \mathbf{x}_i \mathbf{x}_i^T = \mathbf{X}^T \mathbf{W} \mathbf{X} \succeq 0,
\] where \(\mathbf{W} = \text{diag}([\mu_i'(\eta_i)]^2/\sigma_i^2)\). This modified Newton-Raphson algorithm is called the Fisher scoring algorithm.
Take the Binomial logistic regression as an example \[\begin{eqnarray*}
\ell(\boldsymbol{\beta}) &=& \sum_{i=1}^n [y_i \log p_i + (n_i - y_i) \log (1 - p_i)] = \sum_{i=1}^n [y_i \mathbf{x}_i^T \boldsymbol{\beta} - n_i \log (1 + e^{\mathbf{x}_i^T \boldsymbol{\beta}})] \\
\nabla \ell(\boldsymbol{\beta}) &=& \sum_{i=1}^n \left( y_i - n_i \frac{\exp \mathbf{x}_i^T \boldsymbol{\beta}}{1 + \exp \mathbf{x}_i^T \boldsymbol{\beta}} \right) \mathbf{x}_i = \sum_{i=1}^n (y_i - n_i p_i) \mathbf{x}_i = \mathbf{X}^T (\mathbf{y} - \widehat{\boldsymbol{\mu}}) \\
- \nabla^2 \ell(\boldsymbol{\beta}) &=& \sum_{i=1}^n n_i p_i (1 - p_i) \mathbf{x}_i \mathbf{x}_i^T = \mathbf{X}^T \mathbf{W} \mathbf{X}, \quad \mathbf{W} = \text{diag}(w_1, \ldots, w_n), w_i = n_i p_i (1 - p_i) \\
\mathbb{E} [- \nabla^2 \ell(\boldsymbol{\beta})] &=& - \nabla^2 \ell(\boldsymbol{\beta}).
\end{eqnarray*}\] The Fisher scoring algorithmn proceeds as \[\begin{eqnarray*}
\boldsymbol{\beta}^{(t+1)} &=& \boldsymbol{\beta}^{(t)} + s(\mathbf{X}^T \mathbf{W}^{(t)} \mathbf{X})^{-1} \mathbf{X}^T (\mathbf{y} - \widehat{\boldsymbol{\mu}}^{(t)}) \\
&=& (\mathbf{X}^T \mathbf{W}^{(t)} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W}^{(t)} [\mathbf{X} \boldsymbol{\beta}^{(t)} + s (\mathbf{W}^{(t)})^{-1} (\mathbf{y} - \widehat{\boldsymbol{\mu}}^{(t)})] \\
&=& (\mathbf{X}^T \mathbf{W}^{(t)} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W}^{(t)} \mathbf{z}^{(t)},
\end{eqnarray*}\] where \[
\mathbf{z}^{(t)} = \mathbf{X} \boldsymbol{\beta}^{(t)} + s (\mathbf{W}^{(t)})^{-1} (\mathbf{y} - \widehat{\boldsymbol{\mu}}^{(t)})
\] are working responses. In this sense, the Fisher scoring algorithm for GLM is also called the IRWLS (iteratively reweighted least squares).
Recall the insecticide data bliss we used for binomial regression.
Let’s manually implement the first iteration of the Fisher scoring algorithm.
p <- bliss$dead /30eta <-logit(p)z <- etaw <-30* p * (1- p)lmod <-lm(z ~ conc, weights = w, data = bliss)coef(lmod)
(Intercept) conc
-2.302462 1.153587
The Fisher scoring algorithm converges quickly after a few iterations.
y = bliss$deadfor (iter in1:5) { eta <- lmod$fit p <-ilogit(eta) w <-30* p * (1- p) z <- eta + (y -30* p) / w lmod <-lm(z ~ conc, weights = w, data = bliss)cat(iter, coef(lmod), "\n")}
glm(cbind(dead, alive) ~ conc, family = binomial, data = bliss) %>%summary()
Call:
glm(formula = cbind(dead, alive) ~ conc, family = binomial, data = bliss)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.3238 0.4179 -5.561 2.69e-08 ***
conc 1.1619 0.1814 6.405 1.51e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 64.76327 on 4 degrees of freedom
Residual deviance: 0.37875 on 3 degrees of freedom
AIC: 20.854
Number of Fisher Scoring iterations: 4
After finding the MLE \(\widehat{\boldsymbol{\beta}}\), the variance is obtained by \[
\operatorname{Var} \widehat{\boldsymbol{\beta}} = \widehat{\phi} (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1}.
\]
When considering the choice of model for data, two extremes are the null or intercept-only model and the full or saturated model.
The null model means there’s no relation between predictors and the response. Usually it means we fit a common mean \(\mu\) for all \(y\).
The full model means data is explained exactly. Typically it means we need to use \(n\) parameters for \(n\) data points.
To assess the goodness of fit of a model, we might consider likelihood ratio test. For independent observations from exponential family with \(a_i(\phi) = \phi\), this simplifies to \[
\frac{D(\mathbf{y}, \widehat{\boldsymbol{\mu}})}{\phi} = \frac{2 \sum_i [y_i(\tilde \theta_i - \hat \theta_i) - b(\tilde \theta_i) + b(\hat \theta_i)]}{\phi},
\] where \(\tilde \theta\) are the estimates under the full model and \(\hat \theta\) are the estimates under the model of interest. \(D(\mathbf{y}, \widehat{\boldsymbol{\mu}})\) is called the deviance and \(D(\mathbf{y}, \widehat{\boldsymbol{\mu}}) / \phi\) is the scaled deviance.
An alternative measure of goodness of fit is the Pearson’s \(X^2\) statistic\[
X^2 = \sum_i \frac{(y_i - \hat \mu_i)^2}{\operatorname{Var}(\hat \mu_i)}.
\]
For goodness of fit test, we use the fact that, under certain conditions, provided the model is correct, the scaled Deviance and the Pearson’s \(X^2\) statistic are both asymptotically \(\chi^2\) with degrees of freedom equal to the difference of the numbers of identifiable parameters.
For Gaussian, \(\phi\) is unknown so this test cannot be used. For binomial and Poisson, \(\phi=1\) so the test is practical. However the accuracy of asymptotic approximation is dubious for smaller batch sizes. For binary responses, the approximation is worthless.
To compare two nested models \(\Omega\) and \(\omega\), difference of the scaled deviance \(D_\omega - D_\Omega\) is asymptotically \(\chi^2\) with degrees of freedom equal to the difference in the number of identifiable parameters in the two models. For Gaussian model and other models where the disperson \(\phi\) is unknown, we can insert an estimate of \(\phi\) and compute an \(F\) test \[
\frac{(D_\omega - D_\Omega) / (\text{df}_{\omega} - \text{df}_{\Omega})}{\hat \phi},
\] where \(\hat \phi = X^2 / (n-p)\) is a good estimate of the dispersion. For Gaussian, the F-test is exact. For other models, the F-test is approximate.
Recall the insecticide data bliss we used for binomial regression. The goodness of fit test by analysis deviance shows that this model fits data well. Comparing the fit to the null model also shows that the predictor conc is highly significant.
We mostly use the deviance residuals for diagnostics.
6.2 Leverage and influence
For GLM, the hat matrix is \[
\mathbf{H} = \mathbf{W}^{1/2} \mathbf{X} (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W}^{1/2},
\] where \(\mathbf{W}\) is the weight matrix at the fitted model. Diagonal elements of \(\mathbf{H}\) are the leverages \(h_i\). A larger value of leverage indicates that the fit may be sensitive to the response at case \(i\). Its predictor values are unusual in some way.
Leverage only measures the potential to affect the fit whereas measures of influence (change in coefficients if we omit a case) more directly access the effect of each case on the fit.
modp <-glm(Species ~ . - Endemics - Island, family = poisson, data = gala)summary(modp)
Call:
glm(formula = Species ~ . - Endemics - Island, family = poisson,
data = gala)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.155e+00 5.175e-02 60.963 < 2e-16 ***
Area -5.799e-04 2.627e-05 -22.074 < 2e-16 ***
Elevation 3.541e-03 8.741e-05 40.507 < 2e-16 ***
Nearest 8.826e-03 1.821e-03 4.846 1.26e-06 ***
Scruz -5.709e-03 6.256e-04 -9.126 < 2e-16 ***
Adjacent -6.630e-04 2.933e-05 -22.608 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 3510.73 on 29 degrees of freedom
Residual deviance: 716.85 on 24 degrees of freedom
AIC: 889.68
Number of Fisher Scoring iterations: 5
For GLM, it’s better to plot the deviance residuals against the linear predictors \(\widehat{\eta}\) rather than the predicted responses. We expect to see a constant variance of the deviance residuals.
gala %>%mutate(devres =residuals(modp, type ="deviance"),linpred =predict(modp, type ="link")) %>%ggplot() +geom_point(mapping =aes(x = linpred, y = devres)) +labs(x =expression(hat(eta)), y ="Deviance residual")
If we plot response residuals \(y_i - \widehat{\mu}_i\) against the linear predictor \(\widehat{\eta}\), then we see variance increases with \(\widehat{\mu}\) which is consistent with a Poisson model.
gala %>%mutate(resres =residuals(modp, type ="response"),linpred =predict(modp, type ="link")) %>%ggplot() +geom_point(mapping =aes(x = linpred, y = resres)) +labs(x =expression(hat(eta)), y ="Response residual")
6.4 Check linearity of predictors
Let’s plot the response Species against the predictor Area. We see majority of islands have a small area with a few exceptions.
gala %>%ggplot() +geom_point(mapping =aes(x = Area, y = Species))
Log transform of Area reveals a curvilinear relationship.
gala %>%ggplot() +geom_point(mapping =aes(x =log(Area), y = Species))
Poisson regression uses a log link, which needs to be taken into account. Plotting the linearized response (working response) \[
z = \eta + (y - \mu) \frac{d\eta}{d \mu}
\] against log(Area) shows a linear relationship.
gala %>%mutate(eta =predict(modp, type ="link"),mu =predict(modp, type ="response"),res =residuals(modp, type ="response")) %>%ggplot() +geom_point(mapping =aes(x =log(Area), y = eta + res / mu)) +labs(x ="log(Area)", y ="Linearized response")
From the same reasoning, we fit a model with log transformation of all predictors. We see a substantial reduction in deviance.
modpl <-glm(Species ~log(Area) +log(Elevation) +log(Nearest) +log(Scruz +0.5) +log(Adjacent), family = poisson, data = gala)summary(modpl)
Call:
glm(formula = Species ~ log(Area) + log(Elevation) + log(Nearest) +
log(Scruz + 0.5) + log(Adjacent), family = poisson, data = gala)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.331143 0.286899 11.611 < 2e-16 ***
log(Area) 0.350684 0.018004 19.479 < 2e-16 ***
log(Elevation) 0.032465 0.057047 0.569 0.56929
log(Nearest) -0.039966 0.014126 -2.829 0.00467 **
log(Scruz + 0.5) -0.037242 0.013783 -2.702 0.00689 **
log(Adjacent) -0.089481 0.006945 -12.885 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 3510.7 on 29 degrees of freedom
Residual deviance: 360.0 on 24 degrees of freedom
AIC: 532.83
Number of Fisher Scoring iterations: 5
6.5 Half normal plot
Q-Q plot of the residuals is the standard way to check the normality assumption on the errors. For GLM, it’s better use a half-normal plot that compares the sorted absolute residuals and the quantiles of the half-normal distribution \[
\Phi^{-1} \left( \frac{n+i}{2n + i} \right), \quad i=1,\ldots,n.
\] The residuals are not expected to be normally distributed, so we are not looking for an approximate straight line. We only seek outliers which may be identified as points off the trend. A half-normal plot is better for this purpose because in a sense the resolution of the plot is doubled by having all the points in one tail.
halfnorm(rstudent(modpl))
gali <-influence(modpl)halfnorm(gali$hat)
7 Sandwich estimation
To illustrate with the linear model, suppose we relax the assumption of the linear regression model \(\mathbf{Y} \sim \text{N}(\mathbf{X} \boldsymbol{\beta}, \sigma^2 \mathbf{I})\) to \[
\mathbf{Y} \sim \text{N}(\mathbf{X} \boldsymbol{\beta}, \boldsymbol{\Omega}),
\] where \(\boldsymbol{\Omega}\) is an unknow covariance matrix. We still use the least squares solution estimate \[
\widehat{\boldsymbol{\beta}} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}.
\] Then \[
\operatorname{Var} \widehat{\boldsymbol{\beta}} = (\mathbf{X}^T \mathbf{X})^{-1} (\mathbf{X}^T \boldsymbol{\Omega} \mathbf{X}) (\mathbf{X}^T \mathbf{X})^{-1}.
\]
For correct inference, we would like to plug in an estimate of \(\boldsymbol{\Omega}\) and use this sandwich estimator.
Generalization to GLM is similar.
library(sandwich)glm(Species ~ . - Endemics - Island, family = poisson, data = gala) %>%vcovHC() %>%diag() %>%sqrt()
