Chapter 5: Generalized Linear Models

“You must stick to your conviction, but be ready to abandon your assumptions.” — Dennis Waitley

Generalized Linear Models (GLM’s) are extensions of linear regression to areas where assumptions of normality and homoskedasticity do not hold. There are several versions of GLM’s, each for different types and distributions of outcomes. We are going to go through several of the most common.

This chapter is to introduce the method very briefly and demonstrate how to perform one in R. We do not delve into the details of each method much, but rather focus on showing the quirks of the coding.

We discuss:

  1. Logistic Regression
  2. Poisson Regression
  3. GLM with Gamma distribution
  4. Negative binomial
  5. Beta Regression

Logistic Regression

For binary outcomes (e.g., yes or no, correct or incorrect, sick or healthy), logistic regression is a fantastic tool that provides useful and interpretable information. Much like simple and multiple linear regression, logistic regression17 uses dummy coding and provides coefficients that tell us the relationship between the outcome and the independent variables.

Since the outcome is binary, we use a statistical transformation to make things work well. This makes it so the outcome is in “log-odds.” A simple exponentiation of the coefficients and we get very useful “odds ratios.” These are very common in many fields using binary data.

Luckily, running a logistic regression is simple in R. We first create the binary outcome variable called dep. We use a new function called mutate to create a new variable (we could do this a number of ways but this is probably the cleanest way).

## First creating binary depression variable
df <- df %>%
  mutate(dep = dpq010 + dpq020 + dpq030 + dpq040 + dpq050 +
               dpq060 + dpq070 + dpq080 + dpq090) %>%
  mutate(dep2 = ifelse(dep >= 10, 1,
                ifelse(dep < 10, 0, NA)))

Note that we added the values from the ten variables that give us an overall depression score (dep). We then use ifelse() to create a binary version of depression called dep2 with a cutoff of \(\geq 16\) meaning depressed. Because there are missing values denoted as “NA” in this variable, we use a “nested ifelse” to say:

  1. IF depression \(\geq 10\) then dep2 is 1,
  2. IF dpression \(< 10\), then dep2 is 0,
  3. ELSE dep2 is NA.

Note that these nested ifelse() statements can be as long as you want. We further need to clean up the asthma and sedentary variables.

## Fix some placeholders
df <- df %>%
  mutate(asthma = washer(mcq010, 9),
         asthma = washer(asthma, 2, value = 0)) %>%
  mutate(sed = washer(pad680, 9999, 7777))

Now let’s run the logistic regression:

l_fit <- glm(dep2 ~ asthma + sed + race + famsize,
             data = df,
             family = "binomial")
summary(l_fit)
## 
## Call:
## glm(formula = dep2 ~ asthma + sed + race + famsize, family = "binomial", 
##     data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7831  -0.4479  -0.4078  -0.3645   2.5471  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -2.6203555  0.2380770 -11.006  < 2e-16 ***
## asthma             0.5688452  0.1276326   4.457 8.32e-06 ***
## sed                0.0005638  0.0002610   2.160   0.0307 *  
## raceOtherHispanic  0.7162568  0.2328673   3.076   0.0021 ** 
## raceWhite          0.1287059  0.2116414   0.608   0.5431    
## raceBlack          0.0189205  0.2205461   0.086   0.9316    
## raceOther         -0.4901414  0.2570123  -1.907   0.0565 .  
## famsize           -0.0318309  0.0373218  -0.853   0.3937    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2706.3  on 4436  degrees of freedom
## Residual deviance: 2648.2  on 4429  degrees of freedom
##   (195 observations deleted due to missingness)
## AIC: 2664.2
## 
## Number of Fisher Scoring iterations: 5

We used glm() (stands for generalized linear model). The key to making it logistic, since you can use glm() for a linear model using maximum likelihood instead of lm() with least squares, is family = "binomial". This tells R to do a logistic regression.

Poisson Regression

As we did in logistic regression, we will use the glm() function. The difference here is we will be using an outcome that is a count variable. For example, the sedentary variable (sed) that we have in df is a count of the minutes of sedentary activity.

p_fit <- glm(sed ~ asthma + race + famsize,
             data = df,
             family = "poisson")
summary(p_fit)
## 
## Call:
## glm(formula = sed ~ asthma + race + famsize, family = "poisson", 
##     data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -27.362   -8.430   -1.477    5.823   34.507  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        5.6499871  0.0035550 1589.31   <2e-16 ***
## asthma             0.0614965  0.0021434   28.69   <2e-16 ***
## raceOtherHispanic  0.1393438  0.0040940   34.04   <2e-16 ***
## raceWhite          0.3484622  0.0033438  104.21   <2e-16 ***
## raceBlack          0.3400346  0.0034430   98.76   <2e-16 ***
## raceOther          0.3557953  0.0036273   98.09   <2e-16 ***
## famsize           -0.0188673  0.0005488  -34.38   <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: 496351  on 4436  degrees of freedom
## Residual deviance: 475428  on 4430  degrees of freedom
##   (195 observations deleted due to missingness)
## AIC: 508999
## 
## Number of Fisher Scoring iterations: 5

Sedentary may be over-dispersed (see plot) and so other methods related to poisson may be necessary. For this book, we are not going to be delving into these in depth but we will introduce some below.

Gamma

Regression with a gamma distribution are often found when analyzing costs in dollars. It is very similar to poisson but does not require integers and can handle more dispersion. However, the outcome must have values \(> 0\). Just for demonstration:

## Adjust sed
df$sed_gamma <- df$sed + .01
g_fit <- glm(sed_gamma ~ asthma + race + famsize,
             data = df,
             family = "Gamma")
summary(g_fit)
## 
## Call:
## glm(formula = sed_gamma ~ asthma + race + famsize, family = "Gamma", 
##     data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.3589  -0.4613  -0.0845   0.2926   1.6868  
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.567e-03  1.132e-04  31.515  < 2e-16 ***
## asthma            -1.604e-04  5.865e-05  -2.735  0.00626 ** 
## raceOtherHispanic -4.874e-04  1.309e-04  -3.723  0.00020 ***
## raceWhite         -1.090e-03  1.078e-04 -10.115  < 2e-16 ***
## raceBlack         -1.068e-03  1.102e-04  -9.697  < 2e-16 ***
## raceOther         -1.110e-03  1.145e-04  -9.695  < 2e-16 ***
## famsize            5.107e-05  1.552e-05   3.289  0.00101 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.2932604)
## 
##     Null deviance: 1664.8  on 4436  degrees of freedom
## Residual deviance: 1604.2  on 4430  degrees of freedom
##   (195 observations deleted due to missingness)
## AIC: 59154
## 
## Number of Fisher Scoring iterations: 5

Two-Part or Hurdle Models

We are going to use the pscl package to run a hurdle model. These models are built for situations where there is a count variable with many zeros (“zero-inflated”). The hurdle model makes slightly different assumptions regarding the zeros than the pure negative binomial that we present next. The hurdle consists of two models: one for whether the person had a zero or more (binomial) and if more than zero, how many (poisson).

To run a hurdle model, we are going to make a sedentary variable with many more zeros to illustrate and then we will run a hurdle model.

## Zero inflated sedentary (don't worry too much about the specifics)
df$sed_zero <- ifelse(sample(1:100, 
                             size = length(df$sed), 
                             replace=TRUE) %in% c(5,10,11,20:25), 0, 
                      df$sed)
## Hurdle model
library(pscl)
h_fit = hurdle(sed_zero ~ asthma + race + famsize,
               data = df)
summary(h_fit)
## 
## Call:
## hurdle(formula = sed_zero ~ asthma + race + famsize, data = df)
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5863 -1.5251 -0.2225  1.3001 10.3219 
## 
## Count model coefficients (truncated poisson with log link):
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        5.6353484  0.0037209 1514.49   <2e-16 ***
## asthma             0.0699420  0.0022275   31.40   <2e-16 ***
## raceOtherHispanic  0.1571862  0.0042678   36.83   <2e-16 ***
## raceWhite          0.3599565  0.0034881  103.20   <2e-16 ***
## raceBlack          0.3548714  0.0035945   98.73   <2e-16 ***
## raceOther          0.3643603  0.0037957   95.99   <2e-16 ***
## famsize           -0.0178646  0.0005737  -31.14   <2e-16 ***
## Zero hurdle model coefficients (binomial with logit link):
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        2.62672    0.23188  11.328   <2e-16 ***
## asthma             0.06791    0.15508   0.438    0.661    
## raceOtherHispanic -0.18188    0.26111  -0.697    0.486    
## raceWhite         -0.12325    0.21801  -0.565    0.572    
## raceBlack         -0.27244    0.22267  -1.223    0.221    
## raceOther         -0.34516    0.23388  -1.476    0.140    
## famsize           -0.02101    0.03721  -0.565    0.572    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Number of iterations in BFGS optimization: 13 
## Log-likelihood: -2.327e+05 on 14 Df

Notice that the output has two parts: “Count model coefficients (truncated poisson with log link):” and “Zero hurdle model coefficients (binomial with logit link):”. Together they tell us about the relationship between the predictors and a count variable with many zeros.

Negative Binomial

Similar to that above, negative binomial is for zero-inflated count variables. It makes slightly different assumptions than the hurdle and doesn’t use a two-part approach. In order to run a negative binomial model we’ll use the MASS package and the glm.nb() function.

library(MASS)
fit_nb <- glm.nb(sed_zero ~ asthma + race + famsize,
                 data = df)
summary(fit_nb)

Note that this model is not really appropriate because our data is somewhat contrived.

Beta Regression

For outcomes that are bound between a lower and upper bound, Beta Regression is a great method. For example, if we are looking at test scores that are bound between 0 and 100. It is a very flexible method and allows for some extra analysis regarding the variation.

For this, we are going to use the betareg package. But first, we are going to reach a little and create a ficticiously bound variable in the data set.

## Variable bound between 0 and 1
df$beta_var <- sample(seq(.05, .99, by = .01), 
                      size = length(df$asthma),
                      replace = TRUE)
library(betareg)
fit_beta <- betareg(beta_var ~ asthma + race + famsize,
                    data = df)
summary(fit_beta)
## 
## Call:
## betareg(formula = beta_var ~ asthma + race + famsize, data = df)
## 
## Standardized weighted residuals 2:
##     Min      1Q  Median      3Q     Max 
## -2.0630 -0.6913 -0.0350  0.5947  2.8995 
## 
## Coefficients (mean model with logit link):
##                    Estimate Std. Error z value Pr(>|z|)   
## (Intercept)        0.176892   0.063201   2.799  0.00513 **
## asthma             0.014604   0.043658   0.335  0.73800   
## raceOtherHispanic -0.133659   0.071931  -1.858  0.06315 . 
## raceWhite         -0.060386   0.058692  -1.029  0.30354   
## raceBlack         -0.068759   0.060930  -1.128  0.25911   
## raceOther          0.002296   0.065312   0.035  0.97196   
## famsize            0.002851   0.010809   0.264  0.79198   
## 
## Phi coefficients (precision model with identity link):
##       Estimate Std. Error z value Pr(>|z|)    
## (phi)   2.4803     0.0452   54.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 95.02 on 8 Df
## Pseudo R-squared: 0.001448
## Number of iterations: 17 (BFGS) + 1 (Fisher scoring)

The output provides coefficients and the “Phi” coefficients. Both are important parts of using beta regression but we are not going to discuss it here.

There are many resources available to learn more about beta regression and each of these GLM’s. As for now, we are going to move on to more complex modeling where there are clustering or repeated measures in the data.

Conclusions

One of the great things about R is that most modeling is very similar to the basic lm() function. In all of these GLM’s the arguments are nearly all the same: a formula, the data, and family of model. As you’ll see for Multilevel and Other Models chapters, this does not change much. Having a good start with basic models and GLM’s gets you ready for nearly every other modeling type in R.


  1. Technically, logistic regression is a linear regression model.