The purpose of this blog post is to review the derivation of the logit estimator and the interpretation of model estimates. Logit models are commonly used in statistics to test hypotheses related to binary outcomes, and the logistic classifier is commonly used as a pedagogic tool in machine learning courses as a jumping off point for developing more sophisticated predictive models. A secondary goal is to clarify some of the terminology related to logistic models, which – as should already be clear given the interchanging usage of “logit” and “logistic” – may be confusing. A tertiary motivation is to have a post that future posts elaborating covering causal inference and classifiers for categorical outcomes can link back to.

The following are points to keep in mind:

1. The terms “logit model”, “logistic model”, and “logistic regression model” all refer to the same thing; usage varies by discipline.
2. Logistic regression can be interpreted in many ways, but the most common are in terms of odds ratios and predicted probabilities.

• Predicted probabilities are prefered by most social scientists and the machine learning community while odds ratios are more common in biostatistics and epidemiology.
• Interpreting how much probabilities change given a change in one predictor requires setting values for all predictors.
• Interpreting how much odds change for a change in one predictor does not require taking into account other predictors, though thinking of terms of “odds” is less intuitive.

To begin, the primary reasons we prefer not to use linear regression for categorical outcomes are the following:

1. Impossible predictions, \(p(y = 1) \notin [0,1]\)
2. Heteroskedasticity (errors have non-zero variance only at \(y = 0\) and \(y = 1\))
3. Improper functional forms (assumes constant rate of change, even at extremes)

Say we are trying to determine tolerance for Justin Bieber by alcohol consumption as measured by blood/alcohol levels.

Interpreting the y-axis as the probability of having Bieber fever, there are predictions less than zero. Also, there are only observations when Bieber fever is 0 or 1, and the model assumes that each increase of .1 in blood/alcohol levels has the same effect on the probability of Bieber fever at low, medium, and high levels of consumption.

We can, however, assume that Bieber fever is a continuous concept. If it were continuous, the linear model may apply without these pathologies. As a bridge between the unobserved continuous variable and the observed dichotomization, assume that a threshold \(\tau\) exists such that, when latent Bieber fever is below it, we observe no Bieber fever, and when latent Bieber fever is above it, we observe Yes.

To generalize, denote a latent (unobserved) continuous outcome as \(y^{\ast}\). We need to map the linear model for \(y^*\) onto the range \([0, 1]\) for the probabilities. How do we get there?

We can fall back on our linear model for the latent variable \(y^{\ast}\) to determine the appropriate function that maps the linear model onto a probability.

\[ p(y=1|\mathbf{x}) = p(y^* > \tau | \mathbf{x}) \]

That is, the probability of observing a one is equal to the probability of observing a value of \(y^{\ast}\) greater than the threshold \(\tau\).

To simplify, set \(\tau = 0\). We can do this because \(y^{\ast}\) is unobserved, so its scale is arbitrary. Then

\[ \begin{aligned} p(y=1|x) &= p(y^* > \tau | x) \\ &= p(y^* > 0 | x) \end{aligned} \]

Substitute the linear model in for \(y^{\ast}\) and rearrange

\[ \begin{aligned} p(y=1|x) &= p(x\beta + \epsilon > 0 | x) \\ &= p(\epsilon > – x \beta | x ) \\ &= p(\epsilon \leq x\beta | x ) \end{aligned} \]

The last line

\[ p(y=1|x) = p(\epsilon \leq x\beta | x ) \]

says that \(p(y = 1 | x)\) is equal to the cumulative distribution function for \(\epsilon\) at a given value of \(x\).

Make an assumption that this error distribution follows a standard logistic distribution. The cdf for a standard logistic distribution is

\[ F(\epsilon) = \frac{\mbox{exp}(\epsilon)}{1 + \mbox{exp}(\epsilon)} . \]

Because

\[ p(y=1|x) = p(\epsilon \leq x\beta | x ) \]

is the cdf evaluated at \(x\beta\), we can write this as

\[ p(y=1|x) = F(x\beta) . \] The \(x\beta\) is a linear transformation. It gets stuck inside the standard logistic cdf, which maps the results onto the \([0,1]\) range.

Recall that the pmf for a Bernoulli distribution is

\[ f(y_i) = \pi_i^y(1-\pi_i)^{1-y}, \]

and that

\[ \pi_i = p(y=1|x) = p(\epsilon \leq x\beta | x ) . \]

Assuming \(N\) observations that are independent, the likelihood function is

\[ L(\beta | y, x) = \Pi_{i=1}^{N}\pi_i . \]

Combining and a little re-writing:

\[ L(\beta | y, X) = \Pi_{y=1} p(y_i = 1 | x_i)\Pi_{y=0}[1-p(y_i) = 1 | x] \]

\[ L(\beta | y, X) = \Pi_{y=1} p(y_i = 1 | x_i)\Pi_{y=0}[1-p(y_i) = 1 | x] \]

Substitute what we derived for \(p(y = 1)\)

\[ L(\beta | y, X) = \Pi_{y=1}F(x\beta)\Pi_{y=0}[1-F(x\beta)] \]

Taking logs for tractability

\[ L(\beta | y, X) = \Sigma_{y=1}\mbox{ln}F(x\beta)\Sigma_{y=0}\mbox{ln}[1-F(x\beta)] \]

Optimizing numerically gives us our estimates for \(\beta\).

An example using data from the American National Election Studies

Who voted for Trump?

Test for associations by

1. Age
2. Education
3. Gender

Start with simple associations.

Fit the logistic regression model.

trump_model <- glm(vote ~ gender + educ + age, data = anes, 
                   family = binomial(link = "logit"))

Prettying up the output:

estimates <- round(coef(trump_model), 3)

coef(summary(trump_model)) %>%
  as.data.frame %>%
  rownames_to_column("Variable") %>%
  mutate_if(is.numeric, funs(round(.,3))) %>%
  var_mapping(Variable)  %>%
  kable(align = c("l",rep("c",4)))

The var_mapping function was created for these data to map variable names to clean labels. See code at the bottom of this post.

The estimates returned by the glm() function are the coefficients for the linear part of the logit model,

\[ y^* = \alpha + \beta_1x_1 + \beta_2x_2 + \ldots + \beta_kx_k + \epsilon \]

The prediction for a 55-year-old male who finished high school but did not go to college is:

pp1 <- estimates[1]+ estimates[2]*0 + estimates[3]*1 + estimates[4]*0 + estimates[5]*0 + estimates[6]*0 + estimates[7]*55

\[ \begin{array}[r l l] \hat{y^*} &=& -1.051 + -0.374(0) + 0.655(1) + \\\ & & 0.696(0) + 0.411(0) + -0.424(0) + *0.015(55) \\ &=& 0.429 \end{array} \]

The scale of \(y^*\) is arbitray, so the meaning of this value is ambiguous.

Instead, convert to a probability.

\[ \begin{array} \mbox{Pr}(y=1|x) &= F(x\beta) \\ &= F(0.429) \\ & = \frac{\mbox{exp}(0.429)}{1 + \mbox{exp}(0.429)} \\ & = 0.61 \end{array} \]

In R,

plogis(.429)
## [1] 0.6056349

Note that the coefficient estimates are for the linear model regressing \(y^*\) on the independent variables.

The model for predicted probabilities is not linear. A change in x has a non-constant effect on the change in probability.

Let x be a vector of \(k > 1\) independent variables, and let \(\beta\) be the corresponding coefficients. The partial derivitive for a change in one independent variable \(x_k\) is

\[ \begin{eqnarray} \frac{\partial Pr(y = 1 | x_k)}{\partial x_k} &= \frac{\partial F(x\beta)}{x_k} \\ &= \frac{dF(x\beta)}{dx\beta}\frac{\partial x\beta}{\partial x_k} \\ &= f(x\beta)\beta_k . \end{eqnarray} \]

This means that the change in the predicted probability depends not just on the value of \(x_k\) but also on the values of all of the \(x\)s.

When we use the model to predict a probability, we have to assume values for all of the variables, e.g. a 55-year-old male who finished high school but did not go to college.

To get a sense of how predicted probabilities change, the base R expand.grid() function is quite helpful. We can create a grid of different combinations of independent variables:

predict_data <- expand.grid(educ    = levels(anes$educ),
                            gender  = levels(anes$gender),
                            age     = seq(20,90, by = 10))
head(predict_data, n = 15)
##                     educ gender age
## 1       HS Not Completed   Male  20
## 2           Completed HS   Male  20
## 3      College < 4 Years   Male  20
## 4  College 4 Year Degree   Male  20
## 5        Advanced Degree   Male  20
## 6       HS Not Completed Female  20
## 7           Completed HS Female  20
## 8      College < 4 Years Female  20
## 9  College 4 Year Degree Female  20
## 10       Advanced Degree Female  20
## 11      HS Not Completed   Male  30
## 12          Completed HS   Male  30
## 13     College < 4 Years   Male  30
## 14 College 4 Year Degree   Male  30
## 15       Advanced Degree   Male  30

Plug these in as the new data frame for the predict() method for glm objects.

pred_probs <- predict_data %>%
  mutate(prob_trump = predict(trump_model, newdata = predict_data, 
                              type = "response"))

Note the type argument equals response. By default, predict.glm() will return the estimate of \(y^*\).

Use these to generate plots.

ggplot(pred_probs, aes(x = age, y = prob_trump, group = educ, color = educ)) + geom_line() +
  facet_grid(~ gender) + labs(x = "Age (Years)", y = "Probability of Voting for Trump") + 
  scale_color_discrete(name="Highest\nEducation")

These are of course sample estimates but are presented without confidence intervals. How can we assess uncertainty?

1. Nonparametric bootstrap
2. Parametric bootstrap (e.g. zelig)
3. Delta method

The last of these will be covered in a subsequent blog post.

confint(trump_model) %>%
  exp
## Waiting for profiling to be done...
##                               2.5 %    97.5 %
## (Intercept)               0.2109995 0.5736639
## genderFemale              0.5820164 0.8130810
## educCompleted HS          1.2282458 3.0384294
## educCollege < 4 Years     1.3129020 3.0918562
## educCollege 4 Year Degree 0.9792916 2.3407350
## educAdvanced Degree       0.4182419 1.0295980
## age                       1.0102309 1.0201818

library(tidyverse)
library(ggplot2)
library(knitr)
library(boot)

set.seed(12345)

anes <- read_delim("data/anes_timeseries_2016_rawdata.txt", delim = "|") %>%
  select(vote = V162034a, V161270, gender = V161342, age = V161267) %>%
  mutate_all(as.numeric) %>%
  filter(vote %in% c(1,2) & gender %in% c(1,2)) %>%
  mutate(vote = factor(vote, levels = 1:2, labels = c("Clinton", "Trump")),
         educ = case_when(
           V161270 %in% 1:8 ~ 1,
           V161270 %in% 9 ~ 2,
           V161270 %in% 10:12 ~ 3,
           V161270 %in% 13 ~ 4,
           V161270 %in% 14:16 ~ 5,
           TRUE ~ -999),
         gender = factor(gender, levels = 1:2, labels = c("Male", "Female"))) %>%
  mutate(educ = factor(educ, level = 1:5, labels = c("HS Not Completed",
                                                     "Completed HS",
                                                     "College < 4 Years",
                                                     "College 4 Year Degree",
                                                     "Advanced Degree"))) %>%
  filter(!is.na(educ) & age >= 18) %>%
  dplyr::select(-V161270)


var_mapping <- function(df, var) {
  
  var <- enquo(var)
  
  df_recode <- df %>%
    mutate(!!quo_name(var) := case_when(
      (!!var) == "(Intercept)" ~ "Constant",
      (!!var) == "genderFemale" ~ "Female",
      (!!var) == "educCompleted HS" ~ "Completed HS",
      (!!var) == "educCollege < 4 Years" ~ "College < 4 Years",
      (!!var) == "educCollege 4 Year Degree" ~ "College 4 Year Degree ",
      (!!var) == "educAdvanced Degree" ~ "Advanced Degree", 
      (!!var) == "age" ~ "Age",
      TRUE ~ ""))
  
  df_recode
  
}