---
title: "Multiple logistic regression (Notes)"
subtitle: "STAT 155"
author: "Your Name"
format:
  html:
    toc: true
    toc-depth: 2
    embed-resources: true
---


```{r setup}
#| include: false
knitr::opts_chunk$set(
  collapse = TRUE, 
  warning = FALSE,
  message = FALSE,
  error = TRUE,
  fig.height = 2.75, 
  fig.width = 4.25,
  fig.env = 'figure',
  fig.pos = 'h',
  fig.align = 'center')

# Use a color blind friendly color palette throughout doc
library(tidyverse)
cb_palette <- c("black", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
scale_colour_discrete <- function(...) scale_colour_manual(values = cb_palette, ...)
scale_fill_discrete   <- function(...) scale_fill_manual(values = cb_palette, ...)
theme_set(theme_bw())
```


::: {.callout-note title = "Learning goals"}

- 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.

![](https://ajohns24.github.io/155_Fall_2025/images/concept_map_7.png)

:::

::: {.callout-note title="Additional resources"}

**Required videos**

- [Interpreting logistic regression coefficients](https://youtu.be/FY5bsbtSbpo)
- [Evaluating logistic regression models](https://www.youtube.com/watch?v=jisiKHHDz0U)



**Optional**

- Read Section 4.4 in the [STAT 155 Notes](https://mac-stat.github.io/Stat155Notes/)

:::


\
\
\
\


# References {-}

**Logistic Regression**    


Let $Y$ be a binary categorical response variable:    

$$Y = 
\begin{cases}
1 & \; \text{ if event happens} \\
0 & \; \text{ if event doesn't happen} \\
\end{cases}$$    

Further define

$$\begin{split}
p           & = \text{ probability event happens} \\
1-p         & = \text{ probability event doesn't happen} \\
\text{odds} & = \text{ odds event happens} = \frac{p}{1-p} \\
\end{split}$$    


Then a logistic regression model of $Y$ by $Y$ is

$$\log(\text{odds}) = \beta_0 + \beta_1 X$$

We can transform this to get (curved) models of odds and probability:

$$\begin{split}
\text{odds} & = e^{\beta_0 + \beta_1 X} \\
p           & = \frac{\text{odds}}{\text{odds}+1} = \frac{e^{\beta_0 + \beta_1 X}}{e^{\beta_0 + \beta_1 X}+1} \\
\end{split}$$


\
\


**Coefficient interpretation for quantitative X**    

$$\begin{split}    
\beta_0     & = \text{ LOG(ODDS) when } X=0 \\
e^{\beta_0} & = \text{ ODDS when } X=0 \\
\beta_1     & = \text{ change in LOG(ODDS) per 1 unit increase in } X \\
e^{\beta_1} & = \text{ multiplicative change in ODDS per 1 unit increase in } X \text{ (ie. the odds ratio)} \\
\end{split}$$


**Coefficient interpretation for categorical X**    

$$\begin{split}    
\beta_0     & = \text{ LOG(ODDS) at the reference category } \\
e^{\beta_0} & = \text{ ODDS at the reference category } \\
\beta_1     & = \text{ unit change in LOG(ODDS) relative to the reference} \\
e^{\beta_1} & = \text{ multiplicative change in ODDS relative to the reference } \text{ (ie. the odds ratio)} \\
\end{split}$$



\
\
\
\



# Warm-up {-}


We've learned two foundational modeling tools: "Normal" linear regression (when Y is quantitative) and logistic regression (when Y is binary).

- Normal linear regression examples   
    - model a person's reaction time (Y) by how long they slept (X)
    - model a person's reaction time (Y) by whether or not they took a nap (X)

- logistic regression examples
    - model whether or not a person believes in climate change (Y) based on their age (X)
    - model whether or not a person believes in climate change (Y) based on their location (X)


\
\
\
\


**EXAMPLE 1: Probability, odds, odds ratio review**

In this and the previous activity, we'll be modeling the *success* of a mountain climber.
Consider the following example calculations:

\

event                    prob         odds
------------------------ ----------- ------
climber A is successful  0.50          1
climber B is successful  0.66 (2/3)    2
climber C is successful  0.75          ?


\


Calculate & interpret the conditional *odds* of success for climber C:   
    
Odds(success | climber C) = 
    
    
\

Calculate & interpret the *odds ratio* of success, comparing climber C to climber B:   

Odds(success | Climber C) / Odds(success | Climber B) = 

\

Calculate & interpret the *odds ratio* of success, comparing climber A to climber B:   

Odds(success | Climber A) / Odds(success | Climber B) = 





\
\
\
\


**EXAMPLE 2: Success over time**

Load the climbers data:

```{r}
# Load data and packages
library(tidyverse)
climbers <- read.csv("https://mac-stat.github.io/data/climbers_sub.csv") %>% 
  select(peak_name, success, age, oxygen_used, year, season)
```

To better understand how climbing success rates have changed over time, let's explore success by year:

```{r}
# Bad plot!!
# Shows how year (X) varies by success (Y), not vice versa
climbers %>% 
  ggplot(aes(x = year, y = success)) +
  geom_boxplot() 
```

```{r}
# Bad plot!! How to interpret this?!
climbers %>% 
  ggplot(aes(x = year, y = success)) +
  geom_point() 
```

```{r}
# Better plot!
climbers %>% 
  group_by(year) %>% 
  summarize(success_rate = mean(success)) %>% 
  ggplot(aes(x = year, y = success_rate)) +
  geom_point() + 
  theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1, size = 8))
```

```{r}
climbers %>% 
  mutate(year_bracket = cut(year, breaks = 10)) %>%
  group_by(year_bracket) %>%
  summarize(success_rate = mean(success)) %>% 
  ggplot(aes(x = year_bracket, y = success_rate)) +
  geom_point() + 
  theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1, size = 8))

```

a. Let's now *model* that relationship...

```{r}
# Build the model
success_by_year <- ___(success ~ year, data = climbers, ___)

# Get the model summary table
coef(summary(success_by_year))
```


b. Write out the model formula.

___ = -119.36215 + 0.05934year


c. Interpret both model coefficients on the *odds* scale.


```{r}
# Exponentiate the coefficients

```




\
\
\
\


**EXAMPLE 3: Success by season**

Next, let's explore success rates by season:

```{r}
# Visualize the relationship
# Bad plot! Shows dependence of season (X) on success (Y), not vice versa
climbers %>% 
  ggplot(aes(fill = season, x = success)) +
  geom_bar(position = "fill")
```

```{r}
# Better plot!
# x = X, fill = Y
climbers %>% 
  ggplot(aes(x = season, fill = success)) +
  geom_bar(position = "fill")
```

```{r}
# Model the relationship
success_by_season <- glm(success ~ season, data = climbers, family = "binomial")
coef(summary(success_by_season))
```

Thus the estimated model formula is:

log(odds of `success` | `season`) = -0.6492953 + 0.3801896 `seasonSpring` + 14.2153618 `seasonSummer` - 1.0453004 `seasonWinter`


a. Interpret the intercept coefficient on the *odds* scale.   

```{r}
exp(-0.6492953)
```

- The odds of success are 0.52, i.e. success is roughly half as likely as failure.
- The odds of success are 0.52 in Autumn, i.e. success is roughly half as likely as failure in Autumn.
- The odds of success are roughly half (52%) as big in Autumn than in Summer.



b. Interpret the `seasonWinter` coefficient on the *odds* scale.

```{r}
exp(-1.0453004)
```

- The odds of success are 0.35 in Winter.
- The odds of success are 0.35 higher in Winter than in Autumn.
- The odds of success in Winter are 35% higher than the odds of success in Autumn.
- The odds of success in Winter are 35% as high as (i.e. 65% lower than) the odds of success in Autumn.



    


\
\
\
\





**EXAMPLE 4: Predictions**

Suppose your friend plans an expedition in Winter.
Their *probability* of success is roughly 0.16:

```{r}
# By hand
log_odds <- -0.6492953 - 1.0453004*1
log_odds

odds <- exp(log_odds)
odds

prob <- odds / (odds + 1)
prob
```

```{r}
# Check our work
# Note the use of type = "response"! If we left this off we'd get a prediction on the log(odds) scale
predict(success_by_season, newdata = data.frame(season = "Winter"), type = "response")
```


a. Suppose we had to translate the probability calculation into a yes / no prediction. Yes or no: do you predict that your friend will be successful?




b. In logistic regression, we don't measure the accuracy of a prediction by a residual. Our prediction is either right or wrong! For example, suppose that your friend *does* successfully summit their peak. Was our prediction right or wrong?





\
\
\
\




# Exercises {-}

In these exercises, you'll explore how to:

- interpret *multiple* logistic regression models (with > 1 predictor)
- evaluate the quality of logistic regression predictions

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

```{r}
climb_model_1 <- glm(success ~ age, climbers, family = "binomial")
climb_model_2 <- glm(success ~ age + oxygen_used, climbers, family = "binomial")
```


## Exercise 1: `climb_model_2`  {-}


### Part a {-}

`climb_model_1` uses only `age` as a predictor:

```{r}
# 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?


### Part b {-}

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

```{r}
# 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")
```

```{r}
# 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))
```




### Part c {-}

The estimated model formula is:

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

```{r}
# Get the model summary table
coef(summary(climb_model_2))
```


```{r}
# Exponentiated coefficients
exp(-0.520)
exp(-0.022)
exp(2.896)
```

Interpret the `age` coefficient on the odds scale.
Don't forget to "control for..."!


### Part d {-}

Interpret the `oxygen_usedTRUE` coefficient on the odds scale.
Don't forget to "control for..."!





\
\
\
\


## Exercise 2: 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:

```{r}
climbers %>% 
  mutate(.fitted = predict(climb_model_1, newdata = ., type = "response")) %>% 
  head()
```

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:

```{r}
# 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`:

```{r}
# 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()
```


### Part a {-}

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.

### Part b {-}

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:

```{r}
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 3: 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:

```{r}
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:

```{r}
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()
```

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:

```{r}
# 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"))
```

### Part a {-}

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.    

```{r}

```


### Part b {-}

`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:

```{r}

```


### Part c {-}

`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:

```{r}

```








    



\
\
\
\


## Exercise 4: 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.    

```{r}
# Confusion matrix

# Overall accuracy

# Sensitivity

# Specificity

```
    

\
\
\
\



## Exercise 5: 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


### Part a {-}

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.   


### Part b {-}

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?



\
\
\
\





## Exercise 6: Reflection {-}

What are some similarities and differences between how we interpret and evaluate linear and logistic regression models?


\
\
\
\


## Exercise 7: Choose your own adventure {-}

As a group, decide what to do next:

- Work on PP 5.
- Explore the impact of using a *higher* probability threshold of 0.5, instead of 0.25, to make binary predictions using our 2 models. Calculate the overall accuracy, sensitivity, and specificity for the 2 models.




## Done!

- Finalize your notes: (1) Render your notes to an HTML file; (2) Inspect this HTML in your Viewer -- check that your work translated correctly; and (3) Outside RStudio, navigate to your 'Activities' subfolder within your 'STAT155' folder and locate the HTML file. Open and inspect this in a web browser.
- Clean up your RStudio session: End the rendering process by clicking the 'Stop' button in the 'Background Jobs' tab in the lower left pane.
- Check the wrap-up announcements and solutions in the online course manual, at the bottom of the corresponding chapter.
- Work on homework!


