---
title: "Interaction (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"}

- Critically think through whether an interaction term makes sense, or should be included in a multiple linear regression model.
- Write a model formula for a multiple linear regression model with an interaction term.
- Interpret the coefficients in a multiple linear regression model with an interaction term.

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

:::



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

**Required video**

- [video](https://youtu.be/PIQk4ROJEcs)

**Optional**

- Read: Section 3.9.3 in the [STAT 155 Notes](https://mac-stat.github.io/Stat155Notes/)
- Watch: [Interaction variables](https://voicethread.com/share/15085486/)

:::






\
\
\
\

# References {-}

Let $Y$ be a response variable and $X_1$ and $X_2$ be two predictors.
Our previous multiple linear regression models have assumed that the relationship between $Y$ and $X_1$ is independent of or unrelated to $X_2$.
In contrast, if the relationship between $Y$ and $X_1$ *differs* depending upon the values / categories of $X_2$, we say that our 2 *predictors* $X_1$ and $X_2$ **interact**.
In this case we can include an interaction term in our model.
Mathematically, this is the product of the 2 predictors:

$$
E[Y | X_1, X_2] = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_1*X_2
$$


\
\
\
\


# Warm-up {-}


A multiple linear regression model includes more than 1 predictor / explanatory variable X.
Different explanatory variables can play different roles.
For example, in a model of recovery time (Y) by treatment (predictor of interest $X_1$) we should control for symptom severity (**confounder** $X_2$).
If we *don't*, it might lead to misleading conclusions (eg: that cold treatment A is "better" when really it just looked better because people with less severe symptoms took it).

```{r fig.width = 4, fig.height = 1.5, echo = FALSE, eval = TRUE}
# Don't get distracted by this code!!!
par(mar = rep(0,4))
plot(1, type = "n", xaxt = "n", yaxt = "n", bty = "n", xlab = "", ylab = "", xlim = c(0,6), ylim = c(1,4))
text(c("treatment", "severity", "recovery time"), x = c(1,3,5), y = c(1,4,1), cex = 1.1)
arrows(x0 = c(3,3), y0 = c(4,4)-0.2, x1 = c(1,5), y1 = c(1,1)+0.2, angle = 25, lwd = 4)
arrows(x0 = 1.7, y0 = 1, x1 = 4, y1 = 1, angle = 25, lwd = 4)
```

In a model of wages (Y) by education (predictor of interest $X_1$) we should be mindful that this relationship might *vary* / be *impacted* by job industry (**effect modifier** $X_2$).
Another way to say this is that education and job industry **interact**.
Unless we just want to study one job industry, we need to adjust for this effect modifier by including an **interaction term** between it and education:

```{r fig.width = 4, fig.height = 1.5, echo = FALSE, eval = TRUE}
# Don't get distracted by this code!!!
par(mar = rep(0,4))
plot(1, type = "n", xaxt = "n", yaxt = "n", bty = "n", xlab = "", ylab = "", xlim = c(0,6), ylim = c(1,4))
text(c("education", "industry", "wage"), x = c(1,3,5), y = c(1,4,1), cex = 1.1)
arrows(x0 = 3, y0 = 4-0.2, x1 = 3, y1 = 1+0.2, angle = 25, lwd = 4)
arrows(x0 = 1.7, y0 = 1, x1 = 4, y1 = 1, angle = 25, lwd = 4)
```

\

**In general**

If the relationship between $Y$ and $X_1$ *differs* depending upon the values / categories of $X_2$, then $X_1$ and $X_2$ **interact**.

\
\
\
\



**EXAMPLE 1: More than 2 predictors**

Draw a causal diagram that represents the relationship of wages (Y) with marital status ($X_1$) when *accounting for* age ($X_2$), experience ($X_3$), and hours worked per week ($X_4$).
Are any of these *confounders* or *effect modifiers*?



\
\
\
\



**EXAMPLE 2: Identify the confounder / effect modifier**

It's our job to identify important confounders and effect modifiers.
Context is important!


a. Your friend wants to explore the relationship of bikeshare ridership among casual riders with temperature.    
    - Is weekend status a potential *confounder* or *effect modifier* in this relationship?
    - Represent your ideas using a causal diagram.

```{r echo = FALSE, eval = TRUE, fig.width = 3, fig.height = 2}
library(tidyverse)
bikes <- read_csv("https://mac-stat.github.io/data/bikeshare.csv") %>% 
  filter(temp_actual < 80)

bikes %>%
  ggplot(aes(y = riders_casual, x = temp_actual, color = weekend)) + 
  geom_point()
```

b. Another friend claims that ice cream causes sun burns, noting that sun burn rates increase as ice cream sales increase. They are mistaken!   
    - What important variable did they *not* consider in their analysis?
    - Is this a potential *confounder* or *effect modifier*?
    - Represent your ideas using a causal diagram.




\
\
\
\


**EXAMPLE 3: Interaction or no interaction?**

In addition to context, we can use plots to help identify potential interactions among our predictors, i.e. potential effect modifiers.

a. In a model of penguin bill length, is there a "notable" interaction between bill depth and sex? Provide evidence and explain your answer in context.

```{r echo = FALSE, eval = TRUE}
data(penguins)
penguins %>% 
  na.omit() %>% 
  ggplot(aes(y = bill_len, x = bill_dep, color = sex)) + 
  geom_smooth(method = "lm", se = FALSE)
```

b. [Benoit and Marsh (2008)](http://www.kenbenoit.net/pdfs/ajps_348.pdf) explore the 2002 Irish Dail elections, Ireland's version of the U.S. House of Representatives. Specifically, they explore how the number of `votes` a candidate received is associated with their campaign `spending` (in 1,000s of Euros) and their `incumbent` status: `Yes` = candidate is an incumbent (they are running for RE-election), or `No`. Is there a "notable" interaction between campaign spending and incumbent status? Provide evidence and explain your answer in context.

```{r echo = FALSE, eval = TRUE}
campaigns <- read_csv("https://mac-stat.github.io/data/campaign_spending.csv") %>% 
  select(wholename, district, votes, incumbent, spending) %>% 
  mutate(spending = spending / 1000) %>% 
  filter(!is.na(spending))

campaigns %>% 
  ggplot(aes(y = votes, x = spending, color = incumbent)) +
  geom_smooth(method = "lm", se = FALSE) 
```



c. In a model of bikeshare ridership, is there a "notable" interaction between weekend and user status (whether a user is a *casual* or *registered* rider)? Provide evidence and explain what your answer means in context.

```{r echo = FALSE, eval = TRUE}
# Don't worry about this wrangling code!
new_bikes <- bikes %>%
  select(riders_casual, riders_registered, weekend, temp_actual) %>%
  pivot_longer(cols = riders_casual:riders_registered, 
               names_to = "user", names_prefix = "riders_",
               values_to = "rides") %>%
  mutate(weekend = factor(weekend))

new_bikes %>% 
  ggplot(aes(y = rides, x = user, fill = weekend)) + 
  geom_boxplot()
```


d. In a model of bikeshare ridership among registered riders, is there a "notable" interaction between weekend status and season? Provide evidence and explain your answer in context.


```{r echo = FALSE, eval = TRUE}
bikes %>% 
  ggplot(aes(y = riders_registered, x = season, fill = weekend)) + 
  geom_boxplot()
```


\
\
\
\


**EXAMPLE 4: Implementing interaction terms**

Let's get into the technical details of effect modifiers / including interaction terms in our models.
Below is a model of `wage` by `education`, `industry`, and an *interaction* between `education` and `industry`.
The industries in our analysis are: `ag`, `construction`, `service`.
Identify what part of our `lm()` code is *new*:

```{r eval = TRUE}
cps <- read_csv("https://mac-stat.github.io/data/cps_2018.csv") %>% 
  filter(industry %in% c("ag", "construction", "service"))
```

```{r eval = TRUE}
wage_model <- lm(wage ~ industry * education, cps)
coef(summary(wage_model))
```

\
\
\
\


**EXAMPLE 5: Differences by industry**

The estimated model formula, rounded, is written out and visualized below:

`E[wage | ed, industry] = 31476 - 14427 construction - 55463 service + 62 ed + 2296 ed*construction + 4385 ed*service`

```{r echo = FALSE, eval = TRUE}
cps %>% 
  ggplot(aes(y = wage, x = education, color = industry)) + 
  geom_smooth(method = "lm", se = FALSE)
```

a. Prove the related sub-model formulas for 2 of the 3 industries:    
    - `ag`: $\text{E[wage | education]} = 31476 + 62 \text{ education}$    
    
    - `service`: $\text{E[wage | education]}= -23987 + 4447 \text{ education}$   

b.  Try to answer the following using only the model formula:

    - In what industry do average wages increase the *most* per additional year of education? What is that rate of increase?
    
    - In what industry do average wages increase the *least* per additional year of education? What is that rate of increase?



\
\
\
\


**EXAMPLE 6: Interpretation**

Interpret the following model coefficients and represent them on the plot above.
Remember: *Use non-causal language, include units, and talk about averages rather than individual cases.*


a. Intercept coefficient = 31476    
    This is the expected wage for someone with 0 years of education in...

b. `industryservice` coefficient = -55463   
    The expected wage for someone with 0 years of education is -55462.76 *lower* in...

c. `education` coefficient = 62   
    For every extra year of education, the average wage increases by $61.95 in...

d. `education:industryservice` = 4385   
    The increase in average wage with every extra year of education is...


e. Note that we didn't and can't "hold constant" education when interpreting the industry coefficient, or vice versa. Why?



\
\
\
\




# Exercises {-}

**GOAL:** Explore interaction terms among various types of variables, categorical and/or quantitative!



## Exercise 1: Interaction between 1 Q & 1 C {-}

In exploring the relationship of a candidate's `votes` with their campaign `spending` (in 1,000s of Euros), we must consider their `incumbent` status (an *effect modifier*):

```{r echo = FALSE, eval = TRUE}
campaigns %>% 
  ggplot(aes(y = votes, x = spending, color = incumbent)) +
  geom_smooth(method = "lm", se = FALSE) 
```

The corresponding estimated model formula is   
    
`E[votes | spending, incumbent] = 691 + 209.7 spending + 4814 incumbentYes - 125.9 spending*incumbentYes`

```{r eval = TRUE}
vote_model <- lm(votes ~ spending * incumbent, campaigns)
coef(summary(vote_model))
```

### Part a {-}

Write out the sub-model formulas for incumbents and non-incumbents.
Simplify / combine like terms.

- For non-incumbents:  E[votes | spending] = 

- For incumbents: E[votes | spending] = 

### Part b {-}

Interpret the `(Intercept)` and `incumbentYes` coefficients, 691 and 4814, in context.



### Part c {-}

Interpret the `spending` *and* interaction (`spending:incumbentYes`) coefficients, 209.7 and -125.9, in context.
Pro tip: First focus on the *signs* of these coefficients, then interpret the details.



\
\
\
\


## Exercise 2: But is it meaningful? {-}

When interpreting coefficients, we should reflect on whether they are *scientifically meaningful*.
That is, is the *magnitude* of the coefficient meaningful or notable in the context of the analysis?

### Part a {-}

Is the `spending` coefficient of 209.66 scientifically meaningful?
That is, is the scale of the increase in expected votes with spending among non-incumbents "notable"?


### Part b {-}

Is the interaction term between incumbent status and spending scientifically meaningful?
That is, is the difference in the rate of increase in votes with spending for incumbents vs non-incumbents notable?



\
\
\
\



## Exercise 3: Interaction between 2 C {-}

Next, let's explore the relationship of daily bike ridership (Y) and whether a user is a casual or registered rider **and** whether the day falls on a weekend.
We'll use the same bike data as before, but wrangled into a format that allows us to answer this question:

```{r eval = TRUE, echo = TRUE}
head(new_bikes, 3)
```

In the warm-up, we noted that user and weekend status predictors **interact**, i.e. the relationship of ridership with weekend status varies by user group:

```{r eval = TRUE, echo = TRUE}
new_bikes %>% 
  ggplot(aes(y = rides, x = user, fill = weekend)) + 
  geom_boxplot()
```

### Part a {-}

The estimated model of this relationship is

`E[rides | weekend, user] = 625 + 777 weekendTRUE + 3300 userregistered - 1714 weekendTRUE*userregistered`

```{r eval = TRUE, echo = TRUE}
bike_model <- lm(rides ~ weekend * user, new_bikes)
coef(summary(bike_model))
```

Use this model to make 2 predictions:

- Predict casual ridership on a weekday = 

- Predict registered ridership on a weekend = 


### Part b {-}

Use the appropriate coefficient estimate(s) to answer the following.
(This is not easy! The boxplots can help you assess whether your answers make sense!)

- On weekdays, how many *more* registered vs casual riders would we expect?


- Among casual riders, how many *more* riders would we expect on a weekend vs weekday?


- Among registered riders, how many *fewer* riders would we expect on a weekend vs weekday? NOTE: The answer is NOT 1714.


- Among casual riders, what's the expected ridership on a weekday?



\
\
\
\



## Exercise 4: Interactions between 2 Q {-}

Finally, let's explore the relationship of `price` (Y) with 2 *quantitative* predictors: the `mileage` and `age` of a used car.
In doing so, it's important to keep in mind that `mileage` and `age` **interact**: the relationship of price with mileage *differs* depending upon the age of the car.

### Part a {-}

The lines below represent the relationship of price with mileage for different aged cars: 0 years old (solid), 10 years old (dashed), and 20 years old (dotted).

```{r eval = TRUE}
# Don't worry about this code
cars <- read_csv("https://mac-stat.github.io/data/used_cars.csv") %>%
  mutate(mileage = milage %>% str_replace(",","") %>% str_replace(" mi.","") %>% as.numeric(),
         price = price %>% str_replace(",","") %>% str_replace("\\$","") %>% as.numeric(),
         age = 2025 - model_year) %>% # 2025 so that yr. 2024 cars are one year old
  select(-milage)

cars %>% 
  ggplot(aes(y = price, x = mileage)) + 
  geom_abline(aes(intercept = 90965, slope = -0.66)) + 
  geom_abline(aes(intercept = 64314, slope = -0.46), linetype = 2) + 
  geom_abline(aes(intercept = 37664, slope = -0.26), linetype = 3) + 
  lims(x = c(0, 100000), y = c(0, 100000))
```

Discuss in words how the relationship between price and mileage varies by the aged of a used car.


### Part b {-}

The model of `price` by `mileage` and `age`, with an interaction term, is represented by: 

E[price | mileage, age] = $\beta_0$ + $\beta_1$ mileage + $\beta_2$ age + $\beta_3$ mileage*age

Using the above plot as a guide, what are the *signs* of the 4 $\beta$ coefficients (+/-)?




### Part c {-}

Below is the actual model summary table:

```{r eval = TRUE}
car_model <- lm(price ~ mileage * age, cars)
coef(summary(car_model))
```

Use this to:

- Check your answer to Part b.

- Specify the sub-model formula of E[price | mileage] for 0 year old (new) cars.

- Specify the sub-model formula of E[price | mileage] for 10 year old (used) cars.




\
\
\
\




## Exercise 5: Cars continued {-}

### Part a {-}

Interpret the mileage coefficient, -0.66, in context.
Think: This tells us about the rate of decrease of price with mileage for cars of what age?!

### Part b {-}

Interpret the interaction coefficient, 0.02, in context.
Think: What does it tell us about the rate of decrease of price with mileage as cars age?!


### Part c {-}

To check out plots of the raw data behind this model, visit the qmd or online manual.
Neither of the options are very helpful (it's tough to visualize 3D relationships, especially in the presence of interaction)!

```{r}
cars %>% 
  ggplot(aes(x = mileage, y = price, col = age)) +
  geom_point(alpha = 0.5) + 
  scale_color_viridis_c(option = "H") # a fun, colorblind-friendly palette!
```

```{r}
# BONUS: visualize this in "3D"
library(plotly)
cars %>% 
  plot_ly(x = ~mileage, y = ~age, z = ~price, color = ~price,
          type = "scatter3d", 
          mode = "markers",
          marker = list(size = 5, colorscale = "Viridis"))
```





\
\
\
\


## Next steps {-}

You're encouraged to type up your written notes in the provided qmd!
This is a good way to study.





## 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!


