Multiple Logistic Regression

Notes and in-class exercises

You can download the .qmd file for this activity here and open in R-studio. The rendered version is posted in the course website (Activities tab). I often experiment with the class activities (and see it in live!) and make updates, but I always post the final version before class starts. To be sure you have the most up-to-date copy, please download it once you’ve settled in before class begins.

Notes

Learning goals

By the end of this lesson, you should be able to:

Construct multiple logistic regression models in R
Interpret coefficients in multiple logistic regression models
Use multiple logistic regression models to make predictions
Evaluate the quality of logistic regression models by using predicted probability boxplots and by computing and interpreting accuracy, sensitivity, specificity, false positive rate, and false negative rate

Readings and videos

Please go through the following reading or videos before class.

Reading: Section 4.4 in the STAT 155 Notes
Videos:
- Part 1: Concepts (script)
- Part 2: R Code (script)

File organization: Save this file in the “Activities” subfolder of your “STAT155” folder.

Exercises

Context: In this activity, we’ll look at data from an experiment conducted in 2001-2002 that investigated the influence of race and gender on job applications. The researchers created realistic-looking resumes and then randomly assigned a name to the resume that “would communicate the applicant’s gender and race” (e.g., they assumed the name Emily would generally be interpreted as a white woman, whereas the name Jamal would generally be interpreted as a black man). They then submitted these resumes to job postings in Boston and Chicago and waited to see if the applicant got a call back from the job posting.

You can find a full description of the variables in this dataset here. Today, we’ll focus on the following variables:

received_callback: indicator that the resume got a call back from the job posting
gender: inferred binary gender associated with the first name on the resume
race: inferred race associated with the first name on the resume

Our research question is: does an applicant’s inferred gender and race have an effect on the chance that they receive a callback after submitting their resume for an open job posting?

library(tidyverse)
library(ggplot2)

resume <- read_csv("https://mac-stat.github.io/data/resume.csv")
head(resume)
## # A tibble: 6 × 30
##   job_ad_id job_city job_industry               job_type   job_fed_contractor
##       <dbl> <chr>    <chr>                      <chr>                   <dbl>
## 1       384 Chicago  manufacturing              supervisor                 NA
## 2       384 Chicago  manufacturing              supervisor                 NA
## 3       384 Chicago  manufacturing              supervisor                 NA
## 4       384 Chicago  manufacturing              supervisor                 NA
## 5       385 Chicago  other_service              secretary                   0
## 6       386 Chicago  wholesale_and_retail_trade sales_rep                   0
## # ℹ 25 more variables: job_equal_opp_employer <dbl>, job_ownership <chr>,
## #   job_req_any <dbl>, job_req_communication <dbl>, job_req_education <dbl>,
## #   job_req_min_experience <chr>, job_req_computer <dbl>,
## #   job_req_organization <dbl>, job_req_school <chr>, received_callback <dbl>,
## #   firstname <chr>, race <chr>, gender <chr>, years_college <dbl>,
## #   college_degree <dbl>, honors <dbl>, worked_during_school <dbl>,
## #   years_experience <dbl>, computer_skills <dbl>, special_skills <dbl>, …

# Preparing data for visualizations
resume2 <- resume %>%
  mutate(
    gender = factor(gender),
    race = factor(race),
    received_callback = factor(received_callback)
  )

Exercise 1: Graphical and numerical summaries

Our research question involves three categorical variables: received_callback (1 = yes, 0 = no), gender (f = female, m = male), and race (Black, White). Let’s start by creating a mosaic plot to visually compare inferred binary gender and callbacks:



ggplot(resume2, aes(x = gender, fill = received_callback)) +
  geom_bar(position = "fill") +
  scale_fill_manual(
    name = "Received Callback?\n(1 = yes, 0 = no)",
    values = c("0" = "lightblue", "1" = "steelblue")
  ) +
  labs(
    x = "Inferred Binary Gender (f = female, m = male)",
    y = "Proportion"
  )

In this activity, we’re also interested in looking at the relationship between inferred race and callbacks. One way we can add a third variable to a plot is to use the facet_grid function, particularly when that third variable is categorical. Let’s try that now:

ggplot(resume2, aes(x = gender, fill = received_callback)) +
  geom_bar(position = "fill") +
  facet_grid(. ~ race) +
  scale_fill_manual(
    "Received Callback? \n(1 = yes, 0 = no)",
    values = c("0" = "lightblue", "1" = "steelblue")
  ) +
  labs(
    x = "Inferred Binary Gender (f = female, m = male)",
    y = "Proportion"
  )

Here’s another way of looking at the relationship between these three variables, switching the placement of gender and race in the mosaic plot:

ggplot(resume2, aes(x = race, fill = received_callback)) +
  geom_bar(position = "fill") +
  facet_grid(. ~ gender) +
  scale_fill_manual(
    "Received Callback? \n(1 = yes, 0 = no)",
    values = c("0" = "lightblue", "1" = "steelblue")
  ) +
  labs(
    x = "Inferred Race",
    y = "Proportion"
  )

When we are comparing three categorical variables, a useful numerical summary is to calculate relative frequencies/proportions of cases falling into each category of the outcome variable, conditional on which categories of the explanatory variables they fall into. Run this code chunk to calculate the conditional proportion of resumes that did nor did not receive a callback, given the inferred gender and race of the applicant:

# corresponding numerical summaries
resume %>%
    group_by(race, gender) %>%
    count(received_callback) %>%
    group_by(race, gender) %>%
    mutate(condprop = n/sum(n))
## # A tibble: 8 × 5
## # Groups:   race, gender [4]
##   race  gender received_callback     n condprop
##   <chr> <chr>              <dbl> <int>    <dbl>
## 1 black f                      0  1761   0.934 
## 2 black f                      1   125   0.0663
## 3 black m                      0   517   0.942 
## 4 black m                      1    32   0.0583
## 5 white f                      0  1676   0.901 
## 6 white f                      1   184   0.0989
## 7 white m                      0   524   0.911 
## 8 white m                      1    51   0.0887

Write a short description that summarizes the information you gain from these visualizations and numerical summaries. Write this summary using good sentences that tell a story and do not resemble a checklist. Don’t forget to consider the context of the data, and make sure that your summary addresses our research question: does an applicant’s inferred gender or race have an effect on the chance that they receive a callback?

Exercise 2: Logistic regression modeling

Next, we’ll fit a logistic regression model to these data, modeling the log odds of receiving a callback as a function of the applicant’s inferred gender and race:

\[\log(Odds[ReceivedCallback = 1 \mid gender, race]) = \beta_0 + \beta_1 genderm + \beta_2 racewhite\]

Fill in the blanks in the code below to fit this logistic regression model.

# fit logistic model and save it as object called "mod1"
mod1 <- glm(received_callback ~ gender + race, data = resume, family = "binomial")

Then, run the code chunk below to get the coefficient estimates and exponentiated estimates:

# Original estimates
coef(mod1)
## (Intercept)     genderm   racewhite 
##  -2.6473729  -0.1270306   0.4395841

# Exponentiated estimates
exp(coef(mod1))
## (Intercept)     genderm   racewhite 
##  0.07083706  0.88070675  1.55206159

Write an interpretation of each of the exponentiated coefficients in your logistic regression model.

exp(Intercept): We estimate the odds of getting a callback among those inferred to be black females is only 0.07, meaning that the chance of getting a callback is 0.07 times as large as the chance of not getting a callback (or, inversely, the chance of not getting a callback is 1/0.07 = 14.29) times greater than the chance of getting a callback).
exp(genderm): Comparing applicants of the same inferred race, we estimate that those inferred to be male have an odds of getting a callback that is 0.88 times as high as (or, equivalently, 12% lower than) the odds of getting a callback for those inferred to be female.
exp(racewhite): We estimate that the odds of getting a callback are 1.55 times higher for applicants whose race was inferred to be white as compared to those who were inferred to be black but the same gender.

Exercise 3: Interaction terms

Do you think it would make sense to add an interaction term (between gender and race) to our logistic regression model? Why/why not?

Including an interaction term in our model would allow us to investigate whether the effect of race on getting a callback depends on your gender or, vice versa, if the effect of gender on getting a callback depends on a race. In other words, we could ask questions like: is there more of a discrepancy in callbacks between black and white males than there is among black and white females? Is there more of a discrepancy in callbacks between male and female blacks than there is among male and female whites?

Let’s try adding an interaction between gender and race. Update the code below to fit this new interaction model.

# fit logistic model and save it as object called "mod2"
mod2 <- glm(received_callback ~ gender * race, data = resume, family = "binomial")

Then, run the code chunk below to get the coefficient estimates and exponentiated estimates for this interaction model:

# Original estimates
coef(mod2)
##       (Intercept)           genderm         racewhite genderm:racewhite 
##       -2.64532337       -0.13698360        0.43609385        0.01654707

# Exponentiated estimates
exp(coef(mod2))
##       (Intercept)           genderm         racewhite genderm:racewhite 
##         0.0709824         0.8719845         1.5466539         1.0166847

(CHALLENGE) Write out the logistic regression model formula separately for males and for females. Based on this how would we interpret the exponentiated coefficients in this model?

Exercise 4: Prediction

We can use our models to predict whether or not a resume will receive a call back based on the inferred gender and race of the applicant. Run the code below to use the predict() function to predict the probability of getting a call back for four job applicants: a person inferred to be a black female, a person inferred to be black male, a person inferred to be a white female, and a person inferred to be a white male.

# set up data frame with people we want to predict for
predict_data <- data.frame(
    gender = c("f", "m", "f", "m"),
    race = c("black", "black", "white", "white")
)
print(predict_data)
##   gender  race
## 1      f black
## 2      m black
## 3      f white
## 4      m white

# prediction based on model without interaction
mod1 %>%
    predict(newdata = predict_data, type = "response")
##          1          2          3          4 
## 0.06615111 0.05872314 0.09905323 0.08828000

# prediction based on model with interaction
mod2 %>%
    predict(newdata = predict_data, type = "response")
##          1          2          3          4 
## 0.06627784 0.05828780 0.09892473 0.08869565

Report and compare the predictions we get from predict(). Do they make sense to you based on your understanding of the data? Combine insights from visualizations and modeling to write a few sentences summarizing findings for our research question: does an applicant’s inferred gender and race have an effect on the chance that they receive a callback after submitting their resume for an open job posting?

The predicted probabilities from our logistic regression models show that we estimate those inferred to be black males have the lowest chance of receiving a callback (5.87% based on mod1 and 5.83% based on mod2), followed by black females (6.62% and 6.63%), white males (8.83% and 8.87%), and then white females (9.91% and 9.89%). This matches the trend we observed that those inferred to be white have a greater chance of getting a callback, regardless of gender, and that those who are inferred to be female have a slightly higher chance of getting a callback than those inferred to be male.

Exercise 5: Evaluating logistic models with plots

We’ll fit one more model that adds on to the interaction model to also include years of college, years of work experience, and resume quality. The code below takes our fitted models and stores the predicted probabilities in a variable called .fitted. Then we use boxplots to show the predicted probabilities of receiving a callback in those who actually did and did not receive a callback.

mod3 <- glm(received_callback ~ gender*race + years_college + years_experience + resume_quality, data = resume, family = "binomial")

# mod1 predicted probabilities
resume %>% 
  mutate(.fitted = predict(mod1, newdata = ., type = "response")) %>% 
  ggplot(aes(x = factor(received_callback), y = .fitted)) +
  geom_boxplot()

  
# mod2 predicted probabilities
resume %>% 
  mutate(.fitted = predict(mod2, newdata = ., type = "response")) %>% 
  ggplot(aes(x = factor(received_callback), y = .fitted)) +
  geom_boxplot()

 
# mod3 predicted probabilities
resume %>% 
  mutate(.fitted = predict(mod3, newdata = ., type = "response")) %>% 
  ggplot(aes(x = factor(received_callback), y = .fitted)) +
  geom_boxplot()

Y-axis (.fitted): The predicted probability of getting a callback (from logistic model). So, each box represents the distribution of the model’s predicted probabilities for cases that actually had outcome 0 vs outcome 1.]

Summarize what you learn about the ability of the 3 models to differentiate those who actually did and did not receive a callback. What model seems best, and why?

All 3 models show that those who actually received a callback had higher predicted probabilities of a callback. Models 1 and 2 are very similar–although predicted probabilites of callback are high for those who did actually receive a callback, there is substantial overlap in the boxplots. There is more separation between the boxplots in the third model, perhaps model 3 is best in terms of accuracy.

If you had to draw a horizontal line across each of the boxplots that vertically separates the left and right boxplots well, where would you place them?

We would want to place the horizontal lines such that as much of the left boxplot was below the line (low predicted probabilities for those with Y = 0) and as much of the right boxplot was above the line (high predicted probabilities for those with Y = 1).

Exercise 6: Evaluating logistic models with evaluation metrics

Sometimes we may need to go beyond the predicted probabilities from our model and try to classify individuals into one of the two binary outcomes (received or did not receive a callback). How high of a predicted probability would we need from our model in order to be convinced that the person actually got a callback? This is the idea behind the horizontal lines that we drew in the previous exercise.

Let’s explore using a probability threshold of 0.08 (8%) to make a binary prediction for each case:

If a model’s predicted probability of getting a callback is greater than or equal to 8.5%, we’ll predict they got a callback.
If the predicted probability is below 8%, we’ll predict they didn’t get a callback.

We can visualize this threshold on our predicted probability boxplots:

mod1 <- glm(received_callback ~ gender + race, data = resume, family = "binomial")
mod2 <- glm(received_callback ~ gender * race, data = resume, family = "binomial")
mod3 <- glm(received_callback ~ gender*race + years_college + years_experience + resume_quality, data = resume, family = "binomial")
# Model 1 output
mod1_output <- resume %>%
  mutate(.fitted = predict(mod1, newdata = ., type = "response"))

# Model 2 output
mod2_output <- resume %>%
  mutate(.fitted = predict(mod2, newdata = ., type = "response"))

# Model 3 output
mod3_output <- resume %>%
  mutate(.fitted = predict(mod3, newdata = ., type = "response"))
ggplot(mod1_output, aes(x = factor(received_callback), y = .fitted)) +
    geom_boxplot() +
    geom_hline(yintercept = 0.08, color = "red")

ggplot(mod2_output, aes(x = factor(received_callback), y = .fitted)) +
    geom_boxplot() +
    geom_hline(yintercept = 0.08, color = "red")

ggplot(mod3_output, aes(x = factor(received_callback), y = .fitted)) +
    geom_boxplot() +
    geom_hline(yintercept = 0.08, color = "red")

Next, we can use our threshold to classify each person in our dataset based on their predicted probability of getting a callback: we’ll predict that everyone with a predicted probability higher than our threshold got a callback, and otherwise they did not. Then, we’ll compare our model’s prediction to the true outcome (whether or not they actually did get a callback).

# get binary predictions for mod1 and compare to truth
threshold <- 0.08
mod1_output %>%
    mutate(predictCallback = .fitted >= threshold) %>% ## predict callback if probability greater than or equal to threshold
    count(received_callback, predictCallback) ## compare actual and predicted callbacks
## # A tibble: 4 × 3
##   received_callback predictCallback     n
##               <dbl> <lgl>           <int>
## 1                 0 FALSE            2278
## 2                 0 TRUE             2200
## 3                 1 FALSE             157
## 4                 1 TRUE              235

mod2_output %>%
    mutate(predictCallback = .fitted >= threshold) %>%
    count(received_callback, predictCallback)
## # A tibble: 4 × 3
##   received_callback predictCallback     n
##               <dbl> <lgl>           <int>
## 1                 0 FALSE            2278
## 2                 0 TRUE             2200
## 3                 1 FALSE             157
## 4                 1 TRUE              235

mod3_output %>%
    mutate(predictCallback = .fitted >= threshold) %>%
    count(received_callback, predictCallback)
## # A tibble: 4 × 3
##   received_callback predictCallback     n
##               <dbl> <lgl>           <int>
## 1                 0 FALSE            2465
## 2                 0 TRUE             2013
## 3                 1 FALSE             159
## 4                 1 TRUE              233

We can use the count() output to fill create contingency tables of the results. (These tables are also called confusion matrices.)

Fill in the confusion matrix for Model 3.

Models 1 and 2: (Both models result in the same confusion matrix.)

	Predict callback	Predict no callback	Total
Actually got callback	235	157	392
Actually did not	2200	2278	4478
Total	2435	2435	4870

Model 3:

	Predict callback	Predict no callback	Total
Actually got callback	____	____	____
Actually did not	____	____	____
Total	____	____	____

In-class Notes from the last class!

Now compute the following evaluation metrics for the models:

Models 1 and 2:

Accuracy: P(Predict Y Correctly)
Sensitivity: P(Predict Y = 1 | Actual Y = 1)
Specificity: P(Predict Y = 0 | Actual Y = 0)
False negative rate: P(Predict Y = 0 | Actual Y = 1)
False positive rate: P(Predict Y = 1 | Actual Y = 0)

Model 3:

Accuracy: P(Predict Y Correctly)
Sensitivity: P(Predict Y = 1 | Actual Y = 1)
Specificity: P(Predict Y = 0 | Actual Y = 0)
False negative rate: P(Predict Y = 0 | Actual Y = 1)
False positive rate: P(Predict Y = 1 | Actual Y = 0)

Imagine that we are a career center on a college campus and we want to use this model to help students that are looking for jobs. Consider the consequences of incorrectly predicting whether or not an individual will get a callback. What are the consequences of a false negative? What about a false positive? Which one is worse?

False Negatives (predicting no callback, but actually got callback): this would be a lost opportunity if a student decided not to submit their resume, thinking they wouldn’t get a callback, when actually they would have. False Positives (predicting callback, but actually didn’t get callback): this would be a disappointment for the student, thinking they were going to get a callback but they ended up not getting one.

Additional exercises

Build 2 models that you’ll use below, the first of which you explored in the previous activity:

# Load data
climbers <- read.csv("https://mac-stat.github.io/data/climbers_sub.csv") %>% 
  select(peak_name, success, age, oxygen_used, year, season)

# Build 2 models
climb_model_1 <- glm(success ~ age, climbers, family = "binomial")
climb_model_2 <- glm(success ~ age + oxygen_used, climbers, family = "binomial")

Exercise 7: `climb_model_2`

climb_model_1 uses only age as a predictor:

# Plot climb_model_1
climbers %>% 
  ggplot(aes(y = as.numeric(success), x = age)) + 
  geom_smooth(method = "glm", se = FALSE, method.args = list(family = "binomial")) +
  labs(y = "probability of success")

In contrast, the multiple logistic regression model climb_model_2 includes 2 predictors of success: age and oxygen_used. In comparison to the above plot of climb_model_1, what do you anticipate climb_model_2 looking like?

Check your intuition! Summarize, in words, what you learn about the relationship of success with age and oxygen use from the below plots.

# Plot climb_model_2
climbers %>% 
  ggplot(aes(y = as.numeric(success), x = age, color = oxygen_used)) + 
  geom_smooth(method = "glm", se = FALSE, method.args = list(family = "binomial")) +
  labs(y = "probability of success")

# Just for fun zoom out!
climbers %>% 
  ggplot(aes(y = as.numeric(success), x = age, color = oxygen_used)) + 
  geom_smooth(method = "glm", se = FALSE, method.args = list(family = "binomial"),
              fullrange = TRUE) +
  labs(y = "probability of success") + 
  lims(x = c(-300, 400))

The estimated model formula is:

log(odds of success | age, oxygen_used) = -0.520 - 0.022 age + 2.896 oxygen_usedTRUE

# Get the model summary table
coef(summary(climb_model_2))
##                    Estimate  Std. Error   z value      Pr(>|z|)
## (Intercept)     -0.51957247 0.207331212 -2.506002  1.221049e-02
## age             -0.02197407 0.005459778 -4.024718  5.704361e-05
## oxygen_usedTRUE  2.89559690 0.126370801 22.913496 3.408591e-116

# Exponentiated coefficients
exp(coef(climb_model_2))
##     (Intercept)             age oxygen_usedTRUE 
##       0.5947748       0.9782656      18.0942987

Interpret the age coefficient on the odds scale. Don’t forget to “control for…”!

Interpret the oxygen_usedTRUE coefficient on the odds scale. Don’t forget to “control for…”!

Exercise 8: Model evaluation using plots

Recall our 2 models of climbing success:

model	predictors
`climb_model_1`	age
`climb_model_2`	age, oxygen_used

Let’s compare the quality of these models. To begin, the code below calculates the probability of success for every climber in our dataset using climb_model_1, and stores it as the .fitted variable:

climbers %>% 
  mutate(.fitted = predict(climb_model_1, newdata = ., type = "response")) %>% 
  head()
##    peak_name success age oxygen_used year season   .fitted
## 1 Ama Dablam    TRUE  28       FALSE 1981 Spring 0.4389240
## 2 Ama Dablam    TRUE  27       FALSE 1981 Spring 0.4448359
## 3 Ama Dablam    TRUE  35       FALSE 1981 Spring 0.3981129
## 4 Ama Dablam    TRUE  37       FALSE 1981 Spring 0.3866827
## 5 Ama Dablam    TRUE  43       FALSE 1981 Spring 0.3531755
## 6 Ama Dablam   FALSE  38       FALSE 1981 Spring 0.3810132

A visual way to evaluate the quality of our model is to compare the probability of success it assigns to climbers that were and were not actually successful:

# Compare the climb_model_1 probability of success calculations for those that were & weren't successful
climbers %>% 
  mutate(.fitted = predict(climb_model_1, newdata = ., type = "response")) %>% 
  ggplot(aes(x = success, y = .fitted)) +
  geom_boxplot()

Similarly, check out the probability calculations using climb_model_2:

# Compare the climb_model_2 probability of success calculations for those that were & weren't successful
climbers %>% 
  mutate(.fitted = predict(climb_model_2, newdata = ., type = "response")) %>% 
  ggplot(aes(x = success, y = .fitted)) +
  geom_boxplot()

Summarize what you learn from the plots about the ability of climb_model_1 and climb_model_2 to differentiate between those who were and were not actually successful. Comment on which model seems better, and why.
Focus on the plot for climb_model_2. Sometimes we may need to go beyond the probability calculations and make a binary prediction about whether a climber will succeed or fail. Suppose we used a 0.5 probability threshold: If the probability of success exceeds 0.5, then predict success. Otherwise, predict failure. In pictures:

climbers %>% 
  mutate(.fitted = predict(climb_model_2, newdata = ., type = "response")) %>% 
  ggplot(aes(x = success, y = .fitted)) +
  geom_boxplot() + 
  geom_hline(yintercept = 0.5, color = "red")

Discuss:

What do you think of the 0.5 threshold? Would it lead to good predictions? Would it do better at predicting when climbers will succeed or when they’ll fail?
If we used the 0.5 threshold, what would the specificity be: less than 25%, between 25% & 50%, between 50% & 75%, or above 75%? Answer this using only the plot.
If we used the 0.5 threshold, what would the sensitivity be: less than 25%, between 25% & 50%, between 50% & 75%, or above 75%? Answer this using only the plot.
If you don’t like the 0.5 threshold, what do you think would be a better option? Why?

Exercise 9: Model evaluation using evaluation metrics

Let’s explore using a probability threshold of 0.25 (25%) to make a binary prediction of success / failure for each climber in our dataset:

If a model’s probability of success is greater than or equal to 25%, we’ll predict that they’ll succeed.
If the probability is below 25%, we’ll predict that they’ll fail.

We can visualize this threshold on our predicted probability boxplots for our 2 models:

climbers %>% 
  mutate(.fitted = predict(climb_model_1, newdata = ., type = "response")) %>% 
  ggplot(aes(x = success, y = .fitted)) +
  geom_boxplot() + 
  geom_hline(yintercept = 0.25, color = "red")

climbers %>% 
  mutate(.fitted = predict(climb_model_2, newdata = ., type = "response")) %>% 
  ggplot(aes(x = success, y = .fitted)) +
  geom_boxplot() + 
  geom_hline(yintercept = 0.25, color = "red")

To evaluate the quality of the resulting predictions, we essentially want to count up how many climbers fell on the “correct” side of the threshold, both overall and within both groups (successful and unsuccessful climbers). Let’s start with climb_model_1. Using the 0.25 threshold, this model predicted success for the first 6 climbers – 5 of these predictions were correct (5 of the 6 were actually successful) and 1 was wrong:

threshold <- 0.25
climbers %>% 
  mutate(.fitted = predict(climb_model_1, newdata = ., type = "response")) %>% 
  mutate(predictSuccess = .fitted >= threshold) %>% 
  select(success, age, oxygen_used, .fitted, predictSuccess) %>% 
  head()
##   success age oxygen_used   .fitted predictSuccess
## 1    TRUE  28       FALSE 0.4389240           TRUE
## 2    TRUE  27       FALSE 0.4448359           TRUE
## 3    TRUE  35       FALSE 0.3981129           TRUE
## 4    TRUE  37       FALSE 0.3866827           TRUE
## 5    TRUE  43       FALSE 0.3531755           TRUE
## 6   FALSE  38       FALSE 0.3810132           TRUE

But that’s just the first 6 climbers! We can build a tabyl() to compare the observed success to the predictSuccess for all climbers in our dataset:

# Rows = actual / observed success (FALSE or TRUE)
# Columns = predicted success (FALSE or TRUE)
library(janitor)
climbers %>% 
  mutate(.fitted = predict(climb_model_1, newdata = ., type = "response")) %>% 
  mutate(predictSuccess = .fitted >= threshold) %>% 
  tabyl(success, predictSuccess) %>%
  adorn_totals(c("row", "col"))
##  success FALSE TRUE Total
##    FALSE    25 1244  1269
##     TRUE     7  800   807
##    Total    32 2044  2076

Use the confusion matrix to prove that climb_model_1 has an overall accuracy rate of 40%. That is:

P(Predict Y Correctly) = 0.40

HINT: Divide the total number of correct predictions by the total number of climbers.

climb_model_1 has a sensitivity (aka true positive rate) of 99%, hence a false negative rate of 1%

Sensitivity: P(Predict Y = success | Actual Y = success) = 0.99
Among (conditioned on) successful climbers, the probability of the model correctly predicting success is 99%.
False negative rate: P(Predict Y = failure | Actual Y = success) = 0.01 Among (conditioned on) successful climbers, the probability of the model INcorrectly predicting failure is 1%.

Prove these calculations using the confusion matrix:

climb_model_1 has a specificity (aka true negative rate) of only 2%, hence a false positive rate of 98%:

Specificity: P(Predict Y = failure | Actual Y = failure) = 0.02
Among (conditioned on) unsuccessful climbers, the probability of the model correctly predicting failure is only 2%.
False positive rate: P(Predict Y = success | Actual Y = failure) = 0.98 Among (conditioned on) unsuccessful climbers, the probability of the model INcorrectly predicting success is 98%.

Prove these calculations using the confusion matrix:

Exercise 10: More confusion matrices

Obtain the confusion matrix, overall accuracy, sensitivity (true positive rate), and specificity (true negative rate) for climb_model_2 (with a 0.25 probability threshold). NOTE: The correct metrics are outlined in the next exercise if you want to check your work.

# Confusion matrix

# Overall accuracy

# Sensitivity

# Specificity

Exercise 11: Model comparison

We’ve now considered 2 models of climbing success, and evaluated their predictive performance using a 0.25 probability threshold:

model	predictors	overall	sensitivity	specificity
1	age	0.40	0.99	0.02
2	age, oxygen_used	0.75	0.69	0.79

Fill in the blanks with the model number.

Model ___ had the higher overall accuracy.
Model ___ was better at predicting when a climber would succeed.
Model ___ was better at predicting when a climber would not succeed.

Imagine that you are planning an expedition. Consider the consequences of incorrectly predicting whether or not you will be successful. What are the consequences of a false negative? What about a false positive? Which one is worse? With that in mind, which climbing model do you prefer?

Solutions

Exercise 1: Graphical and numerical summaries

ggplot(resume) + 
    geom_mosaic(aes(x = product(gender), fill = received_callback)) +
    scale_fill_manual("Received Callback? \n(1 = yes, 0 = no)", values = c("lightblue", "steelblue")) + 
    labs(x = "Inferred Binary Gender (f = female, m = male)", y = "Received Callback? (1 = yes, 0 = no)")
## Error in make_title(..., self = self): unused arguments (list(), "Inferred Binary Gender (f = female, m = male)")

ggplot(resume) + 
    geom_mosaic(aes(x = product(gender), fill = received_callback)) +
    facet_grid(. ~ race) +
    scale_fill_manual("Received Callback? \n(1 = yes, 0 = no)", values = c("lightblue", "steelblue")) + 
    labs(x = "Inferred Binary Gender (f = female, m = male)", y = "Received Callback? (1 = yes, 0 = no)")
## Error in make_title(..., self = self): unused arguments (list(), "Inferred Binary Gender (f = female, m = male)")

ggplot(resume) + 
    geom_mosaic(aes(x = product(received_callback, race), fill = received_callback)) +
    facet_grid(. ~ gender) +
    scale_fill_manual("Received Callback? \n(1 = yes, 0 = no)", values = c("lightblue", "steelblue")) + 
    labs(x = "Inferred Race", y = "Received Callback? (1 = yes, 0 = no)")
## Error in make_title(..., self = self): unused arguments (list(), "Inferred Race")

resume %>%
    group_by(race, gender) %>%
    count(received_callback) %>%
    group_by(race, gender) %>%
    mutate(condprop = n/sum(n))
## # A tibble: 8 × 5
## # Groups:   race, gender [4]
##   race  gender received_callback     n condprop
##   <chr> <chr>              <dbl> <int>    <dbl>
## 1 black f                      0  1761   0.934 
## 2 black f                      1   125   0.0663
## 3 black m                      0   517   0.942 
## 4 black m                      1    32   0.0583
## 5 white f                      0  1676   0.901 
## 6 white f                      1   184   0.0989
## 7 white m                      0   524   0.911 
## 8 white m                      1    51   0.0887

Overall, a small proportion of applicants received a callback, with those who were inferred to be black males being least likely to get a callback (5.8%) and those inferred to be white females being most likely to get a callback (9.9%). In general, job applicants whose race was inferred to be white were more likely to receive a callback than those whose race was inferred to be black, regardless of their inferred gender. On the other hand, inferred gender does not seem to have as much of an effect on the chance of receiving a callback, with perhaps just a slight advantage for females.

Exercise 2: Logistic regression modeling

mod1 <- glm(received_callback ~ gender + race, data = resume, family = "binomial")

coef(mod1)
## (Intercept)     genderm   racewhite 
##  -2.6473729  -0.1270306   0.4395841

exp(coef(mod1))
## (Intercept)     genderm   racewhite 
##  0.07083706  0.88070675  1.55206159

exp(Intercept): We estimate the odds of getting a callback among those inferred to be black females is only 0.07, meaning that the chance of getting a callback is 0.07 times as large as the chance of not getting a callback (or, inversely, the chance of not getting a callback is 1/0.07 = 14.29) times greater than the chance of getting a callback).
exp(genderm): Comparing applicants of the same inferred race, we estimate that those inferred to be male have an odds of getting a callback that is 0.88 times as high as (or, equivalently, 12% lower than) the odds of getting a callback for those inferred to be female.
exp(racewhite): We estimate that the odds of getting a callback are 1.55 times higher for applicants whose race was inferred to be white as compared to those who were inferred to be black but the same gender.

Exercise 3: Interaction terms

Including an interaction term in our model would allow us to investigate whether the effect of race on getting a callback depends on your gender or, vice versa, if the effect of gender on getting a callback depends on a race. In other words, we could ask questions like: is there more of a discrepancy in callbacks between black and white males than there is among black and white females? Is there more of a discrepancy in callbacks between male and female blacks than there is among male and female whites?
Let’s try adding an interaction between gender and race. Update the code below to fit this new interaction model.

mod2 <- glm(received_callback ~ gender * race, data = resume, family = "binomial")

coef(mod2)
##       (Intercept)           genderm         racewhite genderm:racewhite 
##       -2.64532337       -0.13698360        0.43609385        0.01654707

exp(coef(mod2))
##       (Intercept)           genderm         racewhite genderm:racewhite 
##         0.0709824         0.8719845         1.5466539         1.0166847

Overall logistic regression model formula in terms of beta coefficients:

\[\log(Odds[ReceivedCallback = 1 \mid gender, race]) = \beta_0 + \beta_1 genderm + \beta_2 racewhite + \beta_3 genderm \times racewhite\]

The model formula for males:

\[\log(Odds[ReceivedCallback = 1 \mid gender=m, race]) = (\beta_0 + \beta_1) + (\beta_2 + \beta_3) racewhite\]

The model formula for females:

\[\log(Odds[ReceivedCallback = 1 \mid gender=f, race]) = \beta_0 + \beta_2 racewhite\] Focusing first on the female model formula, we can see that this is a simple logistic regression model.

exp(beta0): Odds of callback for black females
exp(beta2): This is the odds ratio for race among females. That is, white females have exp(beta2) times the odds of callback than black females.

Then focusing on the male model formula, we can see that this is also a simple logistic regression model.

exp(beta0+beta1): Odds of callback for black males
exp(beta2+beta3): This is the odds ratio for race among males. White males have exp(beta2+beta3) times the odds of callback than black males.

Comparing the male to the female model formula, we have:

exp(beta1): Black males have exp(beta1) times the odds of a callback than black females
exp(beta3): This tells us how many times higher the odds ratio for race is in males as compared to females.

Exercise 4: Prediction

# set up data frame with people we want to predict for
predict_data <- data.frame(
    gender = c("f", "m", "f", "m"),
    race = c("black", "black", "white", "white")
)
print(predict_data)
##   gender  race
## 1      f black
## 2      m black
## 3      f white
## 4      m white

# prediction based on model without interaction
mod1 %>%
    predict(newdata = predict_data, type = "response")
##          1          2          3          4 
## 0.06615111 0.05872314 0.09905323 0.08828000

# prediction based on model with interaction
mod2 %>%
    predict(newdata = predict_data, type = "response")
##          1          2          3          4 
## 0.06627784 0.05828780 0.09892473 0.08869565

The predicted probabilities from our logistic regression models show that we estimate those inferred to be black males have the lowest chance of receiving a callback (5.87% based on mod1 and 5.83% based on mod2), followed by black females (6.62% and 6.63%), white males (8.83% and 8.87%), and then white females (9.91% and 9.89%). This matches the trend we observed that those inferred to be white have a greater chance of getting a callback, regardless of gender, and that those who are inferred to be female have a slightly higher chance of getting a callback than those inferred to be male.

Exercise 5: Evaluating logistic models with plots

mod3 <- glm(received_callback ~ gender*race + years_college + years_experience + resume_quality, data = resume, family = "binomial")

resume %>% 
  mutate(.fitted = predict(mod1, newdata = ., type = "response")) %>% 
  ggplot(aes(x = factor(received_callback), y = .fitted)) +
  geom_boxplot()

  
resume %>% 
  mutate(.fitted = predict(mod2, newdata = ., type = "response")) %>% 
  ggplot(aes(x = factor(received_callback), y = .fitted)) +
  geom_boxplot()


resume %>% 
  mutate(.fitted = predict(mod3, newdata = ., type = "response")) %>% 
  ggplot(aes(x = factor(received_callback), y = .fitted)) +
  geom_boxplot()

All 3 models show that those who actually received a callback had higher predicted probabilities of a callback. Models 1 and 2 are very similar–although predicted probabilites of callback are high for those who did actually receive a callback, there is substantial overlap in the boxplots. There is more separation between the boxplots in the third model, perhaps model 3 is best in terms of accuracy.
We would want to place the vertical lines such that as much of the left boxplot was below the line (low predicted probabilities for those with Y = 0) and as much of the right boxplot was above the line (high predicted probabilities for those with Y = 1).

Exercise 6: Evaluating logistic models with evaluation metrics

ggplot(mod1_output, aes(x = factor(received_callback), y = .fitted)) +
    geom_boxplot() +
    geom_hline(yintercept = 0.08, color = "red")

ggplot(mod2_output, aes(x = factor(received_callback), y = .fitted)) +
    geom_boxplot() +
    geom_hline(yintercept = 0.08, color = "red")

ggplot(mod3_output, aes(x = factor(received_callback), y = .fitted)) +
    geom_boxplot() +
    geom_hline(yintercept = 0.08, color = "red")

threshold <- 0.08
mod1_output %>%
    mutate(predictCallback = .fitted >= threshold) %>% 
    count(received_callback, predictCallback) 
## # A tibble: 4 × 3
##   received_callback predictCallback     n
##               <dbl> <lgl>           <int>
## 1                 0 FALSE            2278
## 2                 0 TRUE             2200
## 3                 1 FALSE             157
## 4                 1 TRUE              235

mod2_output %>%
    mutate(predictCallback = .fitted >= threshold) %>%
    count(received_callback, predictCallback)
## # A tibble: 4 × 3
##   received_callback predictCallback     n
##               <dbl> <lgl>           <int>
## 1                 0 FALSE            2278
## 2                 0 TRUE             2200
## 3                 1 FALSE             157
## 4                 1 TRUE              235

mod3_output %>%
    mutate(predictCallback = .fitted >= threshold) %>%
    count(received_callback, predictCallback)
## # A tibble: 4 × 3
##   received_callback predictCallback     n
##               <dbl> <lgl>           <int>
## 1                 0 FALSE            2465
## 2                 0 TRUE             2013
## 3                 1 FALSE             159
## 4                 1 TRUE              233

Models 1 and 2: (Both models result in the same confusion matrix.)

	Predict callback	Predict no callback	Total
Actually got callback	235	157	392
Actually did not	2200	2278	4478
Total	2435	2435	4870

Model 3:

	Predict callback	Predict no callback	Total
Actually got callback	233	159	392
Actually did not	2013	2465	4478
Total	2246	2624	4870

Now compute the following evaluation metrics for the models:

Models 1 and 2:

Accuracy: P(Predict Y Correctly) = (235 + 2278)/(235 + 157 + 2200 + 2278) = 0.5160164

Sensitivity: P(Predict Y = 1 | Actual Y = 1) = 235/(235 + 157) = 0.5994898

Specificity: P(Predict Y = 0 | Actual Y = 0) = 2278/(2200 + 2278) = 0.5087092

False negative rate: P(Predict Y = 0 | Actual Y = 1) = 157/(235 + 157) = 0.4005102 (notice that this is equal to 1 - Sensitivity)

False positive rate: P(Predict Y = 1 | Actual Y = 0) = 2200/(2200 + 2278) = 0.4912908 (notice that this is equal to 1 - Specificity)

Model 3:

Accuracy: P(Predict Y Correctly) = (233 + 2465)/(233 + 159 + 2013 + 2465) = 0.5540041

Sensitivity: P(Predict Y = 1 | Actual Y = 1) = 233/(233 + 159) = 0.5943878

Specificity: P(Predict Y = 0 | Actual Y = 0) = 2465/(2013 + 2465) = 0.550469

False negative rate: P(Predict Y = 0 | Actual Y = 1) = 159/(233 + 159) = 0.4056122 (notice that this is equal to 1 - Sensitivity)

False positive rate: P(Predict Y = 1 | Actual Y = 0) = 2013/(2013 + 2465) = 0.449531 (notice that this is equal to 1 - Specificity)

Imagine that we are a career center on a college campus and we want to use this model to help students that are looking for jobs. Consider the consequences of incorrectly predicting whether or not an individual will get a callback. What are the consequences of a false negative? What about a false positive? Which one is worse?
- False Negatives (predicting no callback, but actually got callback): this would be a lost opportunity if a student decided not to submit their resume, thinking they wouldn’t get a callback, when actually they would have.
- False Positives (predicting callback, but actually didn’t get callback): this would be a disappointment for the student, thinking they were going to get a callback but they ended up not getting one.

Exercise 7: `climb_model_2`

# Load data
climbers <- read.csv("https://mac-stat.github.io/data/climbers_sub.csv") %>% 
  select(peak_name, success, age, oxygen_used, year, season)

# Build 2 models
climb_model_1 <- glm(success ~ age, climbers, family = "binomial")
climb_model_2 <- glm(success ~ age + oxygen_used, climbers, family = "binomial")

We should expect this to turn into 2 curves, 1 for each category of oxygen_used.

# Plot climb_model_1
climbers %>% 
  ggplot(aes(y = as.numeric(success), x = age)) + 
  geom_smooth(method = "glm", se = FALSE, method.args = list(family = "binomial")) +
  labs(y = "probability of success")

The probability of success decreases with age, both among climbers that use oxygen and among those that don’t. Further, at any age, the probability of success is higher among people that use oxygen.

# Plot climb_model_2
climbers %>% 
  ggplot(aes(y = as.numeric(success), x = age, color = oxygen_used)) + 
  geom_smooth(method = "glm", se = FALSE, method.args = list(family = "binomial")) +
  labs(y = "probability of success")


# Just for fun zoom out!
climbers %>% 
  ggplot(aes(y = as.numeric(success), x = age, color = oxygen_used)) + 
  geom_smooth(method = "glm", se = FALSE, method.args = list(family = "binomial"),
              fullrange = TRUE) +
  labs(y = "probability of success") + 
  lims(x = c(-300, 400))

When controlling for oxygen use, the odds of success are 97.8% as high (2.2% lower) for every additional year in age, on average.

# Get the model summary table
coef(summary(climb_model_2))
##                    Estimate  Std. Error   z value      Pr(>|z|)
## (Intercept)     -0.51957247 0.207331212 -2.506002  1.221049e-02
## age             -0.02197407 0.005459778 -4.024718  5.704361e-05
## oxygen_usedTRUE  2.89559690 0.126370801 22.913496 3.408591e-116

# Exponentiated coefficients
exp(coef(climb_model_2))
##     (Intercept)             age oxygen_usedTRUE 
##       0.5947748       0.9782656      18.0942987

When controlling for age, the odds of success are roughly 18 times higher for people that use oxygen.

Exercise 8: Model evaluation using plots

climbers %>% 
  mutate(.fitted = predict(climb_model_1, newdata = ., type = "response")) %>% 
  head()
##    peak_name success age oxygen_used year season   .fitted
## 1 Ama Dablam    TRUE  28       FALSE 1981 Spring 0.4389240
## 2 Ama Dablam    TRUE  27       FALSE 1981 Spring 0.4448359
## 3 Ama Dablam    TRUE  35       FALSE 1981 Spring 0.3981129
## 4 Ama Dablam    TRUE  37       FALSE 1981 Spring 0.3866827
## 5 Ama Dablam    TRUE  43       FALSE 1981 Spring 0.3531755
## 6 Ama Dablam   FALSE  38       FALSE 1981 Spring 0.3810132

# Compare the climb_model_1 probability of success calculations for those that were & weren't successful
climbers %>% 
  mutate(.fitted = predict(climb_model_1, newdata = ., type = "response")) %>% 
  ggplot(aes(x = success, y = .fitted)) +
  geom_boxplot()


# Compare the climb_model_2 probability of success calculations for those that were & weren't successful
climbers %>% 
  mutate(.fitted = predict(climb_model_2, newdata = ., type = "response")) %>% 
  ggplot(aes(x = success, y = .fitted)) +
  geom_boxplot()

climb_model_2 does a better job differentiating between those who were and were not actually successful – it tends to predict higher probabilities of success for those that were actually successful. In contrast, climb_model_1 assigns pretty similar probabilities of success for both people that were actually successful and those that weren’t.

climbers %>% 
  mutate(.fitted = predict(climb_model_2, newdata = ., type = "response")) %>% 
  ggplot(aes(x = success, y = .fitted)) +
  geom_boxplot() + 
  geom_hline(yintercept = 0.5, color = "red")

The 0.5 threshold would be bad. It would do a good job of predicting when climbers will fail, but will often predict failure for people that eventually succeed (i.e. it would have a high false negative rate).
The specificity be above 75% – the entire box is below the threshold, meaning at least 75% of failed climbers had predicted probability of success below 0.5.
The sensitivity would be between 50% & 75% – more than half of successful climbers but less than 75% had predicted probability of success above 0.5.
0.25ish better separates the 2 groups: Most climbers that failed had a probability prediction below 0.25, and most climbers that succeeded had a probability prediction above 0.25.

Exercise 9: Model evaluation using evaluation metrics

climbers %>% 
  mutate(.fitted = predict(climb_model_1, newdata = ., type = "response")) %>% 
  ggplot(aes(x = success, y = .fitted)) +
  geom_boxplot() + 
  geom_hline(yintercept = 0.25, color = "red")

climbers %>% 
  mutate(.fitted = predict(climb_model_2, newdata = ., type = "response")) %>% 
  ggplot(aes(x = success, y = .fitted)) +
  geom_boxplot() + 
  geom_hline(yintercept = 0.25, color = "red")


threshold <- 0.25
climbers %>% 
  mutate(.fitted = predict(climb_model_1, newdata = ., type = "response")) %>% 
  mutate(predictSuccess = .fitted >= threshold) %>% 
  select(success, age, oxygen_used, .fitted, predictSuccess) %>% 
  head()
##   success age oxygen_used   .fitted predictSuccess
## 1    TRUE  28       FALSE 0.4389240           TRUE
## 2    TRUE  27       FALSE 0.4448359           TRUE
## 3    TRUE  35       FALSE 0.3981129           TRUE
## 4    TRUE  37       FALSE 0.3866827           TRUE
## 5    TRUE  43       FALSE 0.3531755           TRUE
## 6   FALSE  38       FALSE 0.3810132           TRUE

# Rows = actual / observed success (FALSE or TRUE)
# Columns = predicted success (FALSE or TRUE)
library(janitor)
climbers %>% 
  mutate(.fitted = predict(climb_model_1, newdata = ., type = "response")) %>% 
  mutate(predictSuccess = .fitted >= threshold) %>% 
  tabyl(success, predictSuccess) %>%
  adorn_totals(c("row", "col"))
##  success FALSE TRUE Total
##    FALSE    25 1244  1269
##     TRUE     7  800   807
##    Total    32 2044  2076

# a
# Total number of correct predictions / total number of climbers
(25 + 800) / 2076
## [1] 0.3973988

# b
# Number of successful climbers that were predicted to be successful / total number of successful climbers
800 / 807
## [1] 0.9913259

# c
# Number of unsuccessful climbers that were predicted to be unsuccessful / total number of unsuccessful climbers
25 / 1269
## [1] 0.01970055

Exercise 10: More confusion matrices

# Confusion matrix
climbers %>% 
  mutate(.fitted = predict(climb_model_2, newdata = ., type = "response")) %>% 
  mutate(predictSuccess = .fitted >= threshold) %>% 
  tabyl(success, predictSuccess) %>%
  adorn_totals(c("row", "col"))
##  success FALSE TRUE Total
##    FALSE  1008  261  1269
##     TRUE   249  558   807
##    Total  1257  819  2076

# Overall accuracy
(1008 + 558) / 2076
## [1] 0.7543353

# Sensitivity
558 / 807
## [1] 0.6914498

# Specificity
1008 / 1269
## [1] 0.7943262

Exercise 11: Model comparison

model	predictors	overall	sensitivity	specificity
1	age	0.40	0.99	0.02
2	age, oxygen_used	0.75	0.69	0.79

Model 2 had the higher overall accuracy.
Model 1 was better at predicting when a climber would succeed. (it had higher sensitivity)
Model 2 was better at predicting when a climber would not succeed. (it had higher specificity)

In my opinion, specificity is more important. Climbing can be dangerous and so maybe it’s better to anticipate when you might not succeed.

Notes

Learning goals

Readings and videos

Exercises

Exercise 1: Graphical and numerical summaries

Exercise 2: Logistic regression modeling

Exercise 3: Interaction terms

Exercise 4: Prediction

Exercise 5: Evaluating logistic models with plots

Exercise 6: Evaluating logistic models with evaluation metrics

Additional exercises

Exercise 7: climb_model_2

Exercise 8: Model evaluation using plots

Exercise 9: Model evaluation using evaluation metrics

Exercise 10: More confusion matrices

Exercise 11: Model comparison

Solutions

Exercise 1: Graphical and numerical summaries

Exercise 2: Logistic regression modeling

Exercise 3: Interaction terms

Exercise 4: Prediction

Exercise 5: Evaluating logistic models with plots

Exercise 6: Evaluating logistic models with evaluation metrics

Exercise 7: climb_model_2

Exercise 8: Model evaluation using plots

Exercise 9: Model evaluation using evaluation metrics

Exercise 10: More confusion matrices

Exercise 11: Model comparison

Exercise 7: `climb_model_2`

Exercise 7: `climb_model_2`