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

By the end of this lesson, you should be familiar with:

- confounding variables
- how to *control* for confounding variables in our models
- how to *represent* the role of confounding variables using causal diagrams

:::


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

**Required video**

- [Confounding (and other causal diagrams)](https://voicethread.com/myvoice/thread/15362352/96009456/96009456)   
    - Required material: Watch from 0:00 - 6:54.
    - Optional but important material: 6:55 - end.

**Optional**

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

:::


\
\
\
\


# References {-}

Types of variables in an analysis.
First:

- Y = the response variable or outcome of interest
- X = the predictor of interest

When exploring the relationship of Y with X, there are various types of other variables Z.
**We will largely focus on CONFOUNDERS though it's good to contrast these with other types!!!**


\
\




**Variables to INCLUDE in our model**

$$E[Y|X, Z] = \beta_0 + \beta_1 X + \beta_2 Z$$

- **confounders**
    - Z is a confounder if...
        - it is *causally* associated with both Y and X; or 
        - it is *causally* associated with Y, and just *associated* (not causally) with X
    - Failing to adjust for Z leads to an incorrect conclusion about the relationship between Y and X.

- **precision variable**
    - Z is a precision variable if it is *causally* associated with Y, but *not* associated* (not causally) with X.
    - We don't *have* to include Z, but it can provide more *precise* info about the relationship between Y and X, and improve predictions of Y.

- **effect modifier**   
    - Z is an effect modifier if it has an effect on the *relationship* between Y and X.
    - We can include Z as an interaction term.



\
\



**Variables to NOT include in our model**

- **mediator**    
    - Z is a mediator if it is *causally* associated with Y and *caused by* (is the result of ) X.
    - Including Z in our model can make it harder to estimate the effect of X on Y -- it can attenuate (or even eliminate) the estimated association.
    
- **collider**    
    - Z is a collider if it is *caused by* Y *and* *caused by* X.
    - Including Z in our model can introduce bias. For example, it can make it look like there's an association between Y and X even if one doesn't exist!
    
\
\



**In pictures: Causal diagrams or "directed acyclic graphs" (DAGs)**

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

    


\
\
\
\




# Warm-up {-}


**EXAMPLE 1: Review**

Let's revisit the `peaks` data:


```{r}
# Load useful packages and data
library(tidyverse)
peaks <- read_csv("https://mac-stat.github.io/data/high_peaks.csv") %>%
  mutate(ascent = ascent / 1000)

# Check it out 
head(peaks)
```

Our goal is to explore the relationship of the `time` (in hours) that it takes to complete a hike with the hike's `length` (in miles), vertical `ascent`(in 1000s of feet), and `rating` (easy, moderate, or difficult).

```{r}
# Modify this code to ATTEMPT to visualize this relationship
peaks %>% 
  ggplot(aes(y = time, x = length, ___ = ascent, ___ = rating)) + 
  geom_point()
```

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

Below is a *model* of this relationship:

```{r}
peaks_model <- lm(time ~ length + ascent + rating, data = peaks)
summary(peaks_model)
```

Thus the estimated model formula is:

E[time | length, ascent, rating] = 6.511 + 0.459 length + 0.187 ascent - 3.169 ratingeasy - 2.477 ratingmoderate

Interpret the `length` and `ratingeasy` coefficients by using the following **general strategy**: When interpreting a coefficient for variable X, compare two hikes (cases) whose values of X differ by 1 but who are identical for all other variables.





\
\
\
\


::: {.callout-note title="Interpreting multiple regression coefficients"}


**Interpreting the coefficient $\beta_Q$ for a quantitative variable Q**

Holding all other variables in the model constant, each 1-unit increase in Q is associated with $\beta_Q$ change (increase or decrease) in the expected / average value of Y.

\

**Interpreting the coefficient $\beta_C$ for an indicator variable on category C**

Holding all other variables in the model constant, the expected / average value of Y for the category C is $\beta_C$ higher/lower than that of the reference category.

:::




\
\
\
\


**EXAMPLE 2: Confounders**

In exploring the relationship of Y with our primary predictor of interest X, we may want to / need to *adjust* for potential **confounding variables** Z.
For example:

- Y = recovery time from a severe cold
- X = treatment, oral medicine or a shot / injection (received at doctor's office)

A potential **confounder** is:

- Z = severity of cold symptoms

Explain why Z might be a confounder and draw an appropriate causal diagram, i.e. directed acyclic graph (DAG).





\
\
\
\


::: {.callout-note title="2 approaches to handling a confounder"}

Let Z be a potential confounding variable in the relationship between Y and X.

1. Randomized controlled trials / experiments  
    - We can *adjust* or *control* for Z by *randomly* assigning patients to one of the 2 treatment groups.
    - In this case, we *can* draw causal conclusions (correlation implies causation).
    - This is great! But we can't always run an experiment / trial...
  
2. Observational studies   
    - Unlike experimental data, observational data is collected *without* interfering with the study subjects.
    - The saying "correlation does not imply causation" stems from unmeasured / unaccounted for confounding in observational studies.
    - We can *adjust* or *control* for Z by including it in our model of Y vs X.
    - This allows us to only assess the relationship between Y and X among groups with the same Z value, thus "removing" the issue of confounding.
    - Being able to control for confounders *without* running an experiment is one of the **superpowers of multiple linear regression**.


:::





\
\
\
\





# Exercises {-}

**Research question:** Is there a wage gap, hence possibly discrimination, by marital status among 18-34 year olds?


To explore, we'll revisit the `cps` data with employment information collected by the U.S. Current Population Survey (CPS) in 2018:

```{r}
# Import data
cps <- read_csv("https://mac-stat.github.io/data/cps_2018.csv") %>% 
  filter(age >= 18, age <= 34) %>% 
  filter(wage < 250000)

# Check it out
head(cps)
```

Recall that a **simple linear regression model** of `wage` by `marital` suggests that single workers make $17,052 *less* than married workers, on average:

```{r}
wage_model_1 <- lm(wage ~ marital, data = cps)
coef(summary(wage_model_1))
```

That's a big gap!!
BUT this model ignores important **confounding variables** that might help explain this gap.




\
\
\
\



## Exercise 1: Confounders {-}

Recall: One definition of a confounding variable is one that

- is a *cause* of the response / outcome variable (`wage`); AND
- is a *cause* of the predictor of interest (`marital` status)

We can represent this idea with a **causal diagram**:

```{r fig.width = 6, fig.height = 2, 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("marital", "confounder", "wage"), 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.5, y1 = 1, angle = 25, lwd = 4)
```

Another definition of a confounding variable is one that

- is a *cause* of the response variable (`wage`)
- is *associated* with the predictor of interest (`marital` status), but...
    - is *not* a cause of the predictor of interest
    - is *not* caused by the predictor of interest

We can represent this idea with a **causal diagram**:

```{r fig.width = 6, fig.height = 2, echo = FALSE, eval = TRUE}
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("marital", "confounder", "wage"), x = c(1,3,5), y = c(1,4,1), cex = 1.1)
arrows(x0 = 3, y0 = 4-0.2, x1 = 5, y1 = 1+0.2, angle = 25, lwd = 4)
lines(x = c(3, 1), y = c(3.8, 1.2), lwd = 4)
arrows(x0 = 1.7, y0 = 1, x1 = 4.5, y1 = 1, angle = 25, lwd = 4)
```

**QUESTION:** Name at least 2 potential confounders in exploring the relationship of `wage` with `marital` status.


\
\
\
\



## Exercise 2: How & why do confounders bias results? {-}

Unaccounted-for confounders can be a source of **bias** in our models.
When we ignore them, we often over- or under-estimate the true underlying relationship between a predictor X and response variable Y.
To explore why this is important, let's first look at how our focal predictor `marital` is associated with our response variable, `wage`:

```{r}
cps %>%
  ggplot(aes(x = marital, y = wage)) +
  geom_boxplot()
```

Next, consider `age` as a potential *confounder* (Z).
The below plot shows how `age` (Z) is associated with `marital` status (X).
The results make sense: the older a person is, the more likely they are to be married! 


```{r}
cps %>%
  ggplot(aes(x = age, y = marital)) +
  geom_boxplot()
```

Similarly, we can show how `age` (Z) is associated with `wage` (Y):
The results make sense: older people who have been in the workforce longer typically earn more!

```{r}
cps %>%
  ggplot(aes(x = age, y = wage)) +
  geom_smooth(method = "lm", se = FALSE)
```

Let's now *revisit* our initial plot of the relationship between marital status and wages:

```{r}
cps %>%
  ggplot(aes(x = marital, y = wage)) +
  geom_boxplot()
```

Since we now know that age (Z) is associated with both being married (X) *and* higher wages (Y), we know that this plot doesn't tell the full story.
People who are married could simply be earning higher wages because they tend to be older, *not necessarily because they are married!*
Age is therefore a **confounder** in the relationship between marital status and wages.

\
\
\
\


## Exercise 3: Controlling for confounders {-}

The exercise above illustrates that it is important to *control* or *adjust* for `age` when trying to understand the *actual* causal relationship between `marital` status and `wage`.
To this end, we'll have to consider that our `cps` data was *observational*, *not* collected via a *randomized trial* or *experiment*.
If we *had* collected the `cps` data via an experiment, we would have been able to control for `age` in our exploration of `wage` vs `marital` status.
BUT what might that experiment look like?!?
(Hence why are observational data so common in practice?!?)
Discuss with your group.



\
\
\
\





## Exercise 4: Age {-}

We're in luck.
Even though our `cps` data is *observational*, we can control / adjust for `age` and other confounding variables in the relationship between `wage` and `marital` status by *including them in our model*!
**That's one of the superpowers of multiple linear regression**.
Let's start by just including `age`:

```{r}
# Construct the model
wage_model_2 <- lm(wage ~ marital + age, cps)
coef(summary(wage_model_2))
```


### Part a {-}

Fill in the code below to *visualize* this model:
NOTE: The `geom_line()` layer adds in trend lines from our *actual* model whereas `geom_smooth()` would draw the model that uses "interaction terms" (which we'll learn about next class). 

```{r}
cps %>% 
  ggplot(aes(y = ___, x = ___, color = ___)) +
  geom____(size = 0.1, alpha = 0.5) +
  geom_line(aes(y = wage_model_2$fitted.values))
```

### Part b {-}

Suppose 2 workers are the *same age*, but one is married and one is single.
By how much do we expect the single worker's wage to differ from the married worker's wage? (How does this compare to the $17,052 marital gap among all workers?)


### Part c {-}

Multiple choice: How can we interpret the `maritalsingle` coefficient?

- The average wage is \$7500 less for single workers than married workers.
- When controlling for ("holding constant") age, the average wage is \$7500 less for single workers than married workers. That is, *among workers of the same age*, single workers make \$7500 less than married workers on average.


    
    
\
\
\
\




## Exercise 5: More confounders {-}

Let's control for even more potential confounders!
Model wages by marital status while controlling for `age` *and* years of `education`:   

```{r}
wage_model_3 <- lm(wage ~ marital + age + education, cps)
coef(summary(wage_model_3))
```

### Part a {-}

With so many variables, this is a tough model to visualize.
If you *had* to draw it, how would the model trend appear: 1 point, 2 points, 2 lines, 1 plane, or 2 planes?
Explain your rationale.
Hint: Pay attention to whether your predictors are quantitative or categorical.


### Part b {-}

Recall our research question: Is there a wage gap, hence possibly discrimination, by marital status?
Given this question, which coefficient is of primary interest?
Interpret this coefficient.



### Part c {-}

Interpret the two other coefficients in this model.


\
\
\
\



## Exercise 6: Even more {-}

Let's control for *another* potential confounder, the job `industry` in which one works (categorical):

```{r}
wage_model_4 <- lm(wage ~ marital + age + education + industry, cps)
coef(summary(wage_model_4))
```
    
If we *had* to draw it, this model would appear as 12 planes!
The original plane explains the relationship between wage and the 2 quantitative predictors, age and education.
Then this plane is split into 12 (2*6) individual planes, 1 for each possible combination of marital status (2 possibilities) and industry (6 possibilities).


### Part a {-}

Interpret the main coefficient of interest for our research question.


### Part b {-}

When controlling for a worker's age, marital status, and education level, which industry tends to have the highest wages?
The lowest?
Provide some numerical evidence.





\
\
\
\



## Exercise 7: Biggest model yet {-}

Build a model that helps us explore `wage` by `marital` status while controlling for: `age`, `education`, job `industry`, typical number of work `hours`, and `health` status.
Store this model as `wage_model_5` and report a model summary table.


```{r}

```





\
\
\
\




## Exercise 8: Reflection {-}

Take two workers -- one is married and the other is single.
The models above provided the following insights into the difference in expected wages for these two people:    
    
Model            Assume the two people have the same...    Wage difference
---------------- --------------------------------------- -----------------
`wage_model_1`   NA                                               -$17,052
`wage_model_2`   age                                               -$7,500
`wage_model_3`   age, education                                    -$6,478
`wage_model_4`   age, education, industry                          -$5,893
`wage_model_5`   age, education, industry, hours, health           -$4,993


### Part a {-}

Though not the case in every analysis, the `marital` coefficient got closer and closer to 0 as we controlled for more confounders.
Explain the significance of this phenomenon, in context.
What does it *mean*?


### Part b {-}

Do you still find the wage gap for single vs married people to be meaningfully "large"?
Can you think of any remaining factors that might explain part of this remaining gap?
Or do you think we've found evidence of inequitable wage practices for single vs married workers?








\
\
\
\







## Exercise 9: A new (extreme) example {-}

For a more extreme example of why it's important to control for confounding variables, let's return to the `diamonds` data:

```{r}
# Import and wrangle the data
data(diamonds)
diamonds <- diamonds %>% 
  mutate(
      cut = factor(cut, ordered = FALSE),
      color = factor(color, ordered = FALSE),
      clarity = factor(clarity, ordered = FALSE)
  ) %>% 
  select(price, clarity, cut, color, carat)
```

Our goal is to explore how the `price` of a diamond depends upon its `clarity` (a measure of quality).
Clarity is classified as follows, in order from best to worst:

clarity  description
-------- ------------------------------------------------
IF       flawless (no internal imperfections) 
VVS1     very very slightly imperfect
VVS2     " "
VS1      very slightly imperfect
VS2      " "
SI1      slightly imperfect
SI2      " "
I1       imperfect

### Part a {-}

Check out a model of `price` by `clarity`.
What clarity has the highest average price? The lowest?
(This should be surprising!)       

```{r}
diamond_model_1 <- lm(price ~ clarity, data = diamonds)

# Get a model summary
coef(summary(diamond_model_1))
```


### Part b {-}

What *confounding variable* might explain these results? What's your rationale?
    
    


\
\
\
\




## Exercise 10: Size {-}

It turns out that `carat`, the *size* of a diamond, is an important confounding variable.
Let's explore what happens when we control for this in our model:

```{r}
diamond_model_2 <- lm(price ~ clarity + carat, data = diamonds)

# Get a model summary
coef(summary(diamond_model_2))

# Plot the model
diamonds %>% 
  ggplot(aes(y = price, x = carat, color = clarity)) + 
  geom_line(aes(y = diamond_model_2$fitted.values))
```

What do you think now?
Which clarity has the highest expected price?
The lowest?
Provide numerical evidence from the *model* (not just the visualization).





\
\
\
\





## Exercise 11: Simpson's Paradox {-}

Controlling for `carat` didn't just *change* the `clarity` coefficients, hence our understanding of the relationship between `price` and `clarity`...
It flipped the *signs* of many of these coefficients.
This extreme scenario has a name: **Simpson's paradox**.
Though this paradox is not named for "The Simpson's", this graphic does summarize the paradox -- though there's a positive relationship *overall*, the relationship *within each character* is negative.
Image source: https://jollycontrarian.com/images/e/ed/Simpson_paradox.jpg

![](https://jollycontrarian.com/images/e/ed/Simpson_paradox.jpg)




CHALLENGE: Explain *why* this happened and support your argument with *graphical* evidence.
HINTS: Think about the causal diagram below.
*How* do you think `carat` influences `clarity`?
*How* do you think `carat` influences `price`?
Make 2 plots that support your answers.


```{r fig.width = 6, fig.height = 2, echo = FALSE}
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("clarity", "carat", "price"), x = c(1,3,5), y = c(1,4,1), cex = 1.1)
arrows(x0 = 3, y0 = 4-0.2, x1 = 5, y1 = 1+0.2, angle = 25, lwd = 4)
lines(x = c(3, 1), y = c(3.8, 1.2), lwd = 4)
arrows(x0 = 1.7, y0 = 1, x1 = 4.5, y1 = 1, angle = 25, lwd = 4)
```




\
\
\
\



## Exercise 12: Final conclusion {-}

What's your final conclusion about diamond prices?    
Which diamonds are more expensive: flawed ones or flawless ones?
    


\
\
\
\





## Exercise 13: Reflection {-}

Write a one-sentence warning label for what might happen if we do *not* control for confounding variables in our model.









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


