r_tutorials/logistic_regression_interpretation.Rmd

   1 ---
   2 title: "Interpreting Logistic Regression Coefficients with Examples in R"
   3 author: "[Benjamin Mako Hill](https://mako.cc/) ([mako@atdot.cc](mailto:mako@atdot.cc)) and [Aaron Shaw](http://aaronshaw.org) ([aaronshaw@northwestern.edu](mailto:aaronshaw@northwestern.edu))"
   4 output:
   5   pdf_document:
   6     toc: yes
   7     toc_depth: '3'
   8   html_document:
   9     toc: yes
  10     toc_depth: 3
  11     toc_float:
  12       collapsed: true
  13       smooth_scroll: true
  14     theme: readable
  15   header-includes:
  16   - \newcommand{\lt}{<}
  17   - \newcommand{\gt}{>}
  18 ---
  19
  20 ```{r setup, include=FALSE}
  21 knitr::opts_chunk$set(echo = TRUE)
  22 ```
  23
  24 Because many people in this course wind up conducting and interpreting logistic regressions, I wanted to provide a quick overview of how to do that. I strongly recommend [this page at UCLA](http://stats.idre.ucla.edu/r/dae/logit-regression/) that covers how to fit and interpret logistic regression as well as how to create model-predicted probabilities with R. In this document, I'll show how to do it both manually and with R in a bit less detail and complexity than the UCLA statistical consulting folks do.
  25
  26 ## Fitting a logistic model
  27
  28 First, let's use the built-in `mtcars` dataset to create a simple logistic regression on whether or not a car will have greater-than-average gas mileage. This means we're creating a new, dichotomous version of the `mpg` variable that just indicates whether each value is greater or less than the mean of the original `mpg` distribution.
  29
  30 ```{r}
  31 mako.cars <- mtcars
  32 mako.cars$mpg.gtavg <- mako.cars$mpg > mean(mako.cars$mpg)
  33 ```
  34
  35 Now we can fit a logistic regression in which we regress our dichotomous "better than average gas mileage" variable on horsepower and number of gears.
  36 ```{r}
  37 m <- glm(mpg.gtavg ~ hp + gear, data=mako.cars, family=binomial("logit"))
  38 summary(m)
  39 ```
  40
  41 As you can see, the standard errors are quite large. That said, we can still intrepret the coefficients for pedagogical purposes.
  42
  43 ## Intrepreting Coefficients
  44
  45 Interpret coefficients in logistic regression is different from an ordinary least squares model, but still relatively straightforward. If we look at the parameter estimate of the variable `hp`, we can see that the coefficient is -0.3387. Coefficients in logistic regression are *logged odds ratios*. The easiest/usual way to intrepret these is to exponentiate them using the `exp()` function to turn the coefficient into a normal odds ratio.
  46
  47 For example, if we look at the estimate for `hp` we would get:
  48
  49 ```{r}
  50 exp(-0.3378)
  51 ```
  52
  53 *Interpretation:* In this case, we could say that our model estimates that a one unit increase in horsepower is associated with odds of being above average in gas mileage that are 0.71 times as large (i.e., 71% of the odds) as without the increased horsepower. That's a pretty substantial change!
  54
  55 ## Predicted Probabilities By Hand
  56
  57 Odds ratios are way easier to interpret than log-odds ratios, but can still be difficult (unless you come from a gambling family?). I (and others) always suggest that folks interpret logistical regression in terms of specific probabilities instead of just in terms of odds ratios.
  58
  59 Converting to probabilities from logistic regression is a bit complicated. The idea is basically that you will plug numbers into your fitted logistic regression and report the probabilities that your model predicts for what we might think of as hypothetical — or prototypical — individauls, who you can represent as combinations of values for the model predictors. Let's walk through how to do that and then we'll come back to the conceptual side of it.
  60
  61 The standard logistic function is:
  62
  63 $$\frac{1}{1+e^{-k}}$$
  64
  65 In other words, to turn the results from our logistic regression above back into a probability, we plug our model in for $k$:
  66
  67 $$\frac{1}{1+e^{-1(\beta_0 + \beta_1\mathrm{hp} + \beta_2\mathrm{gear})}}$$
  68
  69 In our fitted model above, we know that $\beta_0 = 26.7692$, $\beta_1 = -0.3387$, $\beta_2 = 3.2290$. In order to convert these into probabilities, we use the [logistic function](https://en.wikipedia.org/wiki/Logistic_function). There's good information on how this works in practice in the Wikipedia article on [logistic regression](https://en.wikipedia.org/wiki/Logistic_regression).
  70
  71 For this example, let's say your interested in the predicted probabilities for two "prototypical" cars: both with three gears and one with 100 horsepower and one with 120 horse power. First, lets plug in the numbers:
  72
  73 ```{r}
  74 26.7692 + (-0.3387 * 100) + (3.2290 * 3) # a car with 100 hp and 3 gears
  75 26.7692 + (-0.3387 * 120) + (3.2290 * 3) # a car with 120 hp and 3 gears
  76 ```
  77
  78 These numbers are equivalent to $k$ in the first equation. We can now plug each value for $k$ into that equation (and remember we have to use *negative* $k$):
  79
  80 ```{r}
  81 1/(1+exp(-1*2.5862))
  82 1/(1+exp(-1*-4.1878))
  83 ```
  84
  85 **Interpretation:** In other words, our model predicts that a car with three gears and 100 horse power will have above average mileage 93% of the time and a car with 120 horsepower will have above average mileage 1.5% of the time.
  86
  87 ## Created predicted probabilities with `predict()` in R
  88
  89 You can do the same thing in R using the `predict()` function. First, we make a dataset that includes the predictor values for the prototypical individuals whose outcomes would like to predict. For example, if I wanted to represent the two kinds of cars described above, I could do:
  90
  91 ```{r}
  92 prototypical.cars <- data.frame(gear=3, hp=c(100, 120))
  93 prototypical.cars
  94 ```
  95
  96 If I had more variables, I would need columns for each of them. In general, it's a good idea to hold any control variables at the median values observed in the sample for all of my "prototypical" individuals. This allows me to focus on the predicted change in the outcome associated with a change in just one (or a few) of my key predictors.
  97
  98 I can now use the `predict()` function to create a new variable. We use the option `type="response"` to have it give us predicted probabilities:
  99
 100 ```{r}
 101 predict(m, prototypical.cars, type="response")
 102 ```
 103
 104 These numbers look extremely similar. They're not exactly the same because of rounding.