Logistic Regression in R

Fitting the Model

To fit a logistic regression in R, we will use the glm function, which stands for Generalized Linear Model. Within this function, write the dependent variable, followed by ~, and then the independent variables separated by +’s. When the family is specified as binomial, R defaults to fitting a logit model. Note that we are also treating educ as a numeric variable here.

mod_log <-  tbl %>% 
  mutate(vote = as.factor(vote)) %>%
  glm(formula = vote ~ gender_factor + educ + age, family = binomial)

The broom package contains useful functions, like tidy(), for viewing the output from common models. The conf.int = T argument requests confidence intervals.

mod_log %>%
  tidy(conf.int = T) %>%
  mutate_if(is.numeric, ~ round(.,3)) %>%
  kable(align = c("l", rep("c", 4)))

Because gender_factor is a factor variable, R will automatically create dummy variables for us. educ is an ordered categorical variable, we opt here to treat its effect as linear. The last variable is age. We find that gender, age, and educ have significant effects.

We can also look at our results in terms of odds ratios by setting exponentiate = T in the tidy() function.

mod_log %>%
   tidy(exponentiate = T, conf.int = T) %>%
   mutate_if(is.numeric, ~ round(.,3)) %>%
   kable(align = c("l", rep("c", 6)))

This will produce a table where:

• estimate is the odds ratio
• std.error is the standard error of the odds ratio
• statistic is the \(z\)-statistic
• p.value is the significance
• conf.low is the lower level of the 95% confidence interval for the odds ratios
• conf.high is the upper level of the 95% confidence interval for the odds ratios

Note that odds ratios are simply the exponentiated coefficients from the logit model. For example, the coefficient for educ was -.252. The odds ratio is \(\exp(-.252) = .777\). The standard errors for the odds ratio are based on the delta method. The 95% confidence interval around the odds ratios are the exponentiated 95% confidence intervals from the logit model.

Interpretation of Odds Ratios

The coefficients returned by our logit model are difficult to interpret intuitively, and hence it is common to report odds ratios instead. An odds ratio less than one means that an increase in \(x\) leads to a decrease in the odds that \(y = 1\). An odds ratio greater than one means that an increase in \(x\) leads to an increase in the odds that \(y = 1\). In general, the percent change in the odds given a one-unit change in the predictor can be determined as

\[ \% \text{ Change in Odds} = 100(OR - 1) \]

For example, the odds of voting for Trump are \(100(.701 - 1) = -29.9\%\) lower for females compared to males. In addition, each increase on the education scale leads to a \(100(.777 - 1) = -22.3\%\) decrease in the odds of voting for Trump. Finally, each one year increase in age leads to a \(100(1.013 - 1) = 1.3\%\) increase in the odds of voting for Trump. All of these are statistically significant at \(p < .05\).

Interpretation in Terms of Predicted Probabilities

Odds ratios are commonly reported, but they are still somewhat difficult to intuit given that an odds ratio requires four separate probabilities:

\[ \text{Odds Ratio} = \left(\frac{p(y = 1 \mid x + 1)}{p(y = 0 \mid x + 1)}\right)\bigg/ \left(\frac{p(y = 1 \mid x)}{p(y = 0 \mid x)}\right) \]

It’s much easier to think directly in terms of probabilities. However, due to the nonlinearity of the model, it is not possible to talk about a one-unit change in an independent having a constant effect on the probability. Instead, predicted probabilities require us to also take into account the other variables in the model. For example, the difference in the probability of voting for Trump between males in females may be different depending on if we are talking about educated voters in their 30s or uneducated voters in their 60s.

In R we can find predicted probabilities using the augment function from the broom package, which will append predicted probabilities from our model to any data frame we provide it. Let’s say we wanted to get predicted probabilities for both genders across the range of ages 20-70, holding educ = 4 (college degree). We can specify:

at_vals <- expand.grid(age = seq(20, 75, by = 5), 
                       gender_factor = c("Male", "Female"),
                       educ = 4)

pred_prob <- augment(mod_log, type.predict = 'response', newdata = at_vals) %>%
   mutate(lower = .fitted - 1.96*.se.fit,
          upper = .fitted + 1.96*.se.fit) %>%
   mutate_if(is.numeric, ~ round(.,3)) 

The expand.grid function will create a new object that contains every combination of age, gender, and education levels that are provided as arguments. Thus, it is possible to create a new data file with the specific combination of values in which we are interested. The augment function will then append the predicted probability of voting for Trump to every one of the combinations that are specified. Note the variable names in expand.grid must exactly match the names of the independent variables specified in the previous call to glm. The resulting table is following:

pred_prob %>%
   kable(col.names = c("Age", 
                      "Gender", 
                      "Education", 
                      "Pred. Prob", 
                      "SE", 
                      "Lower",
                      "Upper"), 
         align = 'c')

In the table:

• The first three columns are the values of the independent variables for which the predicted probability is calculated.
• Pred. Prob is the predicted probability.
• SE is the standard error of the predicted probability.
• Lower and Upper give the 95% confidence interval around the predicted probability.

The probability that a 35-year-old, college-educated male votes for Trump is .432, 95% CI = [.395, .469].

This is a lot of output, so we can better represent it with a plot.

plotdata <- pred_prob %>%
  rename(
    predprob = .fitted, 
    se       = .se.fit,
    gender   = gender_factor) 

plotdata %>%
  ggplot(aes(x = age, 
             y = predprob, 
             color = gender)) + 
  geom_point() + 
  geom_line() + 
  geom_errorbar(aes(ymin = lower, 
                    ymax = upper), 
                width = .1) +
  xlab("Age") +
  ylab("Predicted Probability of Voting for Trump") + 
  labs(color = "Gender")