Generalized Linear Models and GLMs

Generalized Linear Models

Components of the Model

A generalized linear model describes the relationship between a response variable and a set of predictors. It has three components:

• Random Component: Specifies the response variable $Y$ and its distribution, conditional on predictors $X<em>1, X</em>2, …, X_{p-1}$.
• Systematic Component: Specifies the predictors in the linear predictor $\eta = \beta<em>0 + \beta</em>1X<em>1 + \beta</em>2X<em>2 + \cdots + \beta</em>{p-1}X_{p-1}$.
• Link Function: \g(\cdot) specifies the relationship between the systematic component and the expected value of the random component.
  • Link functions are chosen based on the probability distribution of the random component.
  • The 'natural' link functions are called canonical link functions.

Canonical Link Functions

Normal error regression model uses the identity link function: $g(x) = x$.

$Y = \beta<em>0 + \beta</em>1X<em>1 + \cdots + \beta</em>{p-1}X_{p-1} + \epsilon$

$E[Y] = \beta<em>0 + \beta</em>1X<em>1 + \cdots + \beta</em>{p-1}X_{p-1}$

Examples of Generalized Linear Models

• Binomial Distribution: For a response $Y$ with a binomial distribution with parameters $n$ and $\pi = P(Y = 1)$, the generalized linear model is given by:

  $\ln\left(\frac{\pi(X<em>1, \dots, X</em>{p-1})}{1 - \pi(X<em>1, \dots, X</em>{p-1})}\right) = \beta<em>0 + \beta</em>1X<em>1 + \cdots + \beta</em>{p-1}X_{p-1}$

  where $\pi(X<em>1, \dots, X</em>{p-1}) = E[Y \mid X<em>1, \dots, X</em>{p-1}]$.

• Poisson Distribution: For a response $Y$ with a Poisson distribution with parameter $\lambda$, the generalized linear model is given by:

  $\ln(\lambda(X<em>1, \dots, X</em>{p-1})) = \beta<em>0 + \beta</em>1X<em>1 + \cdots + \beta</em>{p-1}X_{p-1}$

  where $\lambda(X<em>1, \dots, X</em>{p-1}) = E[Y \mid X<em>1, \dots, X</em>{p-1}]$.

Maximum Likelihood Estimators

• Given a random sample $Y<em>1, Y</em>2, \dots, Y<em>n$ with density functions $f(y</em>i | \theta, X)$, dependent on a parameter $\theta$ (possibly vector-valued) and predictor $X$, the likelihood function is:

  $L(\theta | y<em>1, y</em>2, \dots, y<em>n, X) = f(y</em>1 | \theta, X)f(y<em>2 | \theta, X) \cdots f(y</em>n | \theta, X)$.

• The log-likelihood function is:

  $l(\theta | y<em>1, y</em>2, \dots, y<em>n, X) = \ln(L(\theta | y</em>1, y<em>2, \dots, y</em>n, X)) = \ln f(y<em>1 | \theta, X) + \ln f(y</em>2 | \theta, X) + \cdots + \ln f(y<em>n | \theta, X) = \sum</em>{i=1}^{n} \ln f(y_i | \theta, X)$.

• The likelihood function of a random sample is a product of density functions.

• The log-likelihood function of a random sample is the sum of the logs of density functions.

• Response and predictor values are fixed (observed), and parameters are unknowns to be estimated.

Maximum Likelihood Estimate (MLE)

• The MLE of a parameter $\theta$ is the value that maximizes the likelihood of the observed data.

  $L(\text{mle } \theta | y<em>1, y</em>2, \dots, y<em>n, X) = \max</em>{\text{all possible } \theta} L(\theta | y<em>1, y</em>2, \dots, y_n, X)$.

• $\theta$ can represent multiple parameters (e.g., mean and variance).

• The MLE is usually computed using the log-likelihood function because the log function is increasing, so both maximize at the same place.

• Theoretically, you try all possible values for $\theta$ until you find the one that maximizes the (log)-likelihood function, subject to any constraints on the parameter. This is the MLE.

Example: Bernoulli Distribution

• A random sample of size three is taken from a Bernoulli distribution with unknown proportion of successes $\pi$. The observed values are one success and two failures.

  $f(y \mid \pi) = \begin{cases} \pi^y(1 - \pi)^{1-y} &amp; y = 0, 1 \ 0 &amp; 0 \text{ otherwise } \end{cases}$

  $L(\pi \mid y<em>1 = 1, y</em>2 = 0, y<em>3 = 0) = f(y</em>1 = 1 \mid \pi)f(y<em>2 = 0 \mid \pi)f(y</em>3 = 0 \mid \pi) = [\pi^1(1 - \pi)^{1-1}][\pi^0(1 - \pi)^{1-0}][\pi^0(1 - \pi)^{1-0}] = \pi(1 - \pi)^2, \quad 0 \le \pi \le 1$

  $l(\pi \mid y<em>1 = 1, y</em>2 = 0, y_3 = 0) = \ln(\pi) + 2\ln(1 - \pi)$

Likelihood Ratio Test (LRT)

• Suppose $M<em>0$ and $M</em>1$ are two models where $M<em>0$ is a restricted version of $M</em>1$.

  $M<em>0: g(\mu) = \beta</em>0 + \beta<em>1X</em>1 + \beta<em>2X</em>2 + \cdots + \beta<em>kX</em>k$ (restricted)

  $M<em>1: g(\mu) = \beta</em>0 + \beta<em>1X</em>1 + \beta<em>2X</em>2 + \cdots + \beta<em>kX</em>k + \beta<em>{k+1}X</em>{k+1} + \cdots + \beta<em>{p-1}X</em>{p-1}$ (full)

• Let $l_M$ be the maximum log-likelihood of a model $M$.

Example: Bernoulli Distribution (Continued)

• Find the maximum log-likelihood function:

  $l(\pi \mid y<em>1 = 1, y</em>2 = 0, y_3 = 0) = \ln(\pi) + 2\ln(1 - \pi)$

  $\frac{dl}{d\pi} = \frac{1}{\pi} - \frac{2}{1 - \pi} = 0$
  $\frac{\pi}{1 - \pi} = \frac{1}{2}$
  $\text{mle } \pi = \frac{1}{3}$
  $l(\text{mle } \pi \mid y<em>1 = 1, y</em>2 = 0, y_3 = 0) = \ln(\frac{1}{3}) + 2\ln(1 - \frac{1}{3})$

• The likelihood ratio test statistic is:

  $LRT = -2(l<em>{M</em>0} - l<em>{M</em>1})$.

  • If n >>> 0, the LRT statistic has an approximate $\chi^2$ distribution with $k$ degrees of freedom.
  • $k$ is the difference in the number of parameters between the two models.

• The likelihood ratio test of:

  $H<em>0: \beta</em>{k+1} = \cdots = \beta_{p-1} = 0$ (restricted)

  $H_1: \text{at least one is not equal to 0}$ (full)

• Rejects the null hypothesis if:

  1. The test statistic $LRT$ is greater than the appropriate critical value: $LRT \ge \chi^2_{\alpha, k}$.
  2. The p-value is less than $\alpha$: P(\chi^2_k \ge LRT) < \alpha.

Binary Data: Logistic Regression

• Using a linear regression model with a binary response leads to fitting a line with an infinite range to the binary response.
• Logistic regression fits a curve that estimates the probability of a response at various predictor values.

Challenger Space Shuttle Example

• On January 28, 1986, the Challenger space shuttle disaster occurred 73 seconds into the flight, killing all seven crew members.
• The investigation focused on O-rings; damage may be related to ambient temperature during the launch.

Odds of Success

• Suppose $Y<em>i$ is a Bernoulli random variable with success probability $\pi</em>i$. The odds of success are:

  $o<em>i = \frac{\pi</em>i}{1 - \pi_i}$.

Binomial Distribution

• A random variable $Y$ has a binomial distribution with parameters $n$ (number of independent trials) and $\pi$ (success probability) if its probability mass function is:

  p(y) = \binom{n}{y} \pi^y (1 - \pi)^{n-y}, \quad y = 0, 1, \dots, n, \quad 0 < \pi < 1
  dbinom(y, size = n, prob = π)

• Logistic Regression Model:

  $\ln\left(\frac{\pi(X<em>1, \dots, X</em>{p-1})}{1 - \pi(X<em>1, \dots, X</em>{p-1})}\right) = \beta<em>0 + \beta</em>1X<em>1 + \cdots + \beta</em>{p-1}X_{p-1}$

  $\pi(X<em>1, \dots, X</em>{p-1}) = \frac{1}{1 + e^{-(\beta<em>0 + \beta</em>1X<em>1 + \cdots + \beta</em>{p-1}X_{p-1})}}$

  $L(\beta<em>0, \beta</em>1, \dots, \beta<em>{p-1} | Y</em>i, X<em>{i1}, \dots, X</em>{i,p-1}, i = 1, 2, \dots, n) = \prod<em>{i=1}^{n} \left( \frac{1}{1 + e^{-(\beta</em>0 + \beta<em>1X</em>{i1} + \cdots + \beta<em>{p-1}X</em>{i,p-1})}} \right)^{Y<em>i} \left(1 - \frac{1}{1 + e^{-(\beta</em>0 + \beta<em>1X</em>{i1} + \cdots + \beta<em>{p-1}X</em>{i,p-1})}} \right)^{1-Y_i}$

Challenger Space Shuttle Example (Continued)

• Partial Dataframe:

  head(Shuttle, n = 3)

    failed temp
  1      0   66
  2      1   70
  3      0   69
  

• Fit Generalized Linear Model:

  shuttle.glm <- glm(failed ~ temp, data = Shuttle, family = binomial(link = "logit"))

• Fit Intercept Only Model:

  shuttle.intercept.only <- glm(failed ~ 1, data = Shuttle, family = binomial(link = "logit"))

• Likelihood Ratio Test:

  anova(shuttle.intercept.only, shuttle.glm, test = "LRT")

  Analysis of Deviance Table
  
  Model 1: failed ~ 1
  Model 2: failed ~ temp
    Resid. Df Resid. Dev Df Deviance Pr(>Chi)
  1        22     28.267                      
  2        21     20.315  1   7.952  0.004804 **
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  

• Summarize Results:

  sum.shuttle <- summary(shuttle.glm)

  Call:
  glm(formula = failed ~ temp, family = binomial(link = "logit"), 
      data = Shuttle)
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)
  (Intercept)   15.0429     7.3786   2.039   0.0415 *
  temp         -0.2322     0.1082  -2.145   0.0320 *
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for binomial family taken to be 1)
  Null deviance: 28.267  on 22  degrees of freedom
  Residual deviance: 20.315  on 21  degrees of freedom
  AIC: 24.315
  
  Number of Fisher Scoring iterations: 5
  

• Confidence Interval Parameter:

  confint(shuttle.glm)

                  2.5 %       97.5 %
  (Intercept)  3.3305848 34.34215133
  temp        -0.5154718 -0.06082076
  

• Estimate Odds:

  predict(shuttle.glm, newdata = data.frame(temp = c(36, 60)), type = "link")

      1        2
  6.685043 1.113137
  

• Estimate Probability:

  predict(shuttle.glm, newdata = data.frame(temp = c(36, 60)), type = "response")

          1         2
  0.9987521 0.7527135
  

Donner–Reed Party Example

• The Donner–Reed Party was a group moving to California in 1846-47 who got stuck in the Sierra Nevada mountains during the winter.
• Data on party members 15 years of age and older are in the DonnerShort dataset.

Partial Dataframe

• head(Donner, n = 2)

    age gender survived
  1  23      1        0
  2  40      0        1
  

Fit Generalized Linear Model

• donner.full <- glm(survived ~ age + gender, data = Donner, family = binomial(link = "logit"))

Fit Intercept Only Model

• donner.intercept.only <- glm(survived ~ 1, data = Donner, family = binomial(link = "logit"))

Likelihood Ratio Test

• anova(donner.intercept.only, donner.full, test = "LRT")

  Analysis of Deviance Table
  
  Model 1: survived ~ 1
  Model 2: survived ~ age + gender
    Resid. Df Resid. Dev Df Deviance Pr(>Chi)
  1        44     61.827                      
  2        42     51.256  2   10.57  0.005066 **
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  

Summarize Results

• sum.donner <- summary(donner.full)

  Call:
  glm(formula = survived ~ age + gender, family = binomial(link = "logit"), 
      data = Donner)
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)
  (Intercept)   3.23041    1.38686   2.329   0.0198 *
  age          -0.07820    0.03728  -2.097   0.0359 *
  gender       -1.59729    0.75547  -2.114   0.0345 *
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for binomial family taken to be 1)
  Null deviance: 61.827  on 44  degrees of freedom
  Residual deviance: 51.256  on 42  degrees of freedom
  AIC: 57.256
  
  Number of Fisher Scoring iterations: 4
  

Confidence Interval Parameter

• confint(donner.full, level = 0.95)

                  2.5 %       97.5 %
  (Intercept)  0.8514190 6.42669512
  age         -0.1624377 -0.01406576
  gender      -3.2286705 -0.19510198
  

Estimate Odds

• new.data <- data.frame(age = c(50, 60), gender = c(1, 0))

• predict(donner.full, newdata = new.data, type = "link")

          1          2
  -2.277083 -1.461831
  

Estimate Probability

• predict(donner.full, newdata = new.data, type = "response")

          1          2
  0.09303878 0.18818751
  

Fit Generalized Linear Model - Age Only

• donner.age.only <- glm(survived ~ age, data = Donner, family = binomial(link = "logit"))

Summarize Results

• sum.age.only <- summary(donner.age.only)

  Call:
  glm(formula = survived ~ age, family = binomial(link = "logit"), 
      data = Donner)
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)
  (Intercept)   1.81852    0.99937   1.820   0.0688 .
  age          -0.06647    0.03222  -2.063   0.0391 *
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for binomial family taken to be 1)
  Null deviance: 61.827  on 44  degrees of freedom
  Residual deviance: 56.291  on 43  degrees of freedom
  AIC: 60.291
  
  Number of Fisher Scoring iterations: 4
  

Likelihood Ratio Test - Age vs Full

• anova(donner.age.only, donner.full, test = "LRT")

  Analysis of Deviance Table
  
  Model 1: survived ~ age
  Model 2: survived ~ age + gender
    Resid. Df Resid. Dev Df Deviance Pr(>Chi)
  1        43     56.291                      
  2        42     51.256  1   5.0344  0.02485 *
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  

Likelihood Ratio Test - Age vs Intercept

• anova(donner.intercept.only, donner.age.only, test = "LRT")

  Analysis of Deviance Table
  
  Model 1: survived ~ 1
  Model 2: survived ~ age
    Resid. Df Resid. Dev Df Deviance Pr(>Chi)
  1        44     61.827                      
  2        43     56.291  1   5.5358  0.01863 *
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  

Count Data: Poisson Regression

• A random variable $Y$ has a Poisson distribution with parameter $\lambda$ if its probability mass function is:

  $p(y) = \begin{cases} \frac{e^{-\lambda}\lambda^y}{y!} &amp; y = 0, 1, 2, 3, \dots \ 0 &amp; \text{ otherwise } \end{cases}$

  dpois(y, λ)

• A Poisson distribution can be used when considering the number of occurrences of a specified event over some type of interval.

Horseshoe Crabs Example

• A study of horseshoe crabs examined factors affecting whether a female crab has male crabs (satellites) nesting near her.

• Variables measured: color, spine condition, carapace width, number of satellites, weight.

• A Poisson random variable counts the number of events occurring in a fixed interval of time, where events occur at rate $\lambda$ per interval and independently of the time since the previous event.

• Poisson Regression Model:

  $\ln(\lambda(X<em>1, \dots, X</em>{p-1})) = \beta<em>0 + \beta</em>1X<em>1 + \beta</em>2X<em>2 + \cdots + \beta</em>{p-1}X_{p-1}$

  $\lambda(X<em>1, \dots, X</em>{p-1}) = e^{\beta<em>0 + \beta</em>1X<em>1 + \cdots + \beta</em>{p-1}X_{p-1}}$

  $L(\beta<em>0, \beta</em>1, \dots, \beta<em>{p-1} | Y</em>i, X<em>{i1}, \dots, X</em>{i,p-1}, i = 1, 2, \dots, n) = \prod<em>{i=1}^{n} \frac{e^{-\lambda(X</em>{i1}, \dots, X<em>{i,p-1})} \lambda(X</em>{i1}, \dots, X<em>{i,p-1})^{Y</em>i}}{Y_i!}$

Horseshoe Crabs Example (Continued)

• head(Crabs, n = 3)

    color spine width satellite weight y
  1     3     3  28.3         8   3050 1
  2     4     3  22.5         0   1550 0
  3     2     1  26.0         9   2300 1
  

Fit Generalized Linear Model

• crab.glm <- glm(satellite ~ width, data = Crabs, family = poisson(link = "log"))

Fit Intercept Only Model

• crab.intercept <- glm(satellite ~ 1, data = Crabs, family = poisson(link = "log"))

Likelihood Ratio Test

• anova(crab.intercept, crab.glm, test = "LRT")

  Analysis of Deviance Table
  
  Model 1: satellite ~ 1
  Model 2: satellite ~ width
    Resid. Df Resid. Dev Df Deviance Pr(>Chi)
  1       172     632.79                      
  2       171     567.88  1   64.913  7.828e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  

Summarize Results

• summary(crab.glm)

  Call:
  glm(formula = satellite ~ width, family = poisson(link = "log"), 
      data = Crabs)
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)
  (Intercept) -3.30476    0.54224  -6.095 1.1e-09 ***
  width        0.16405    0.01997   8.216  < 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: 632.79  on 172  degrees of freedom
  Residual deviance: 567.88  on 171  degrees of freedom
  AIC: 927.18
  
  Number of Fisher Scoring iterations: 6
  

Estimate log lambda

• predict(crab.glm, newdata = data.frame(width = c(20, 30)), type = "link")

           1         2
  -0.0238555 1.6165954
  

Estimate lambda

• predict(crab.glm, newdata = data.frame(width = c(20, 30)), type = "response")

          1         2
  0.9764268 5.0359157