---
title: "Intro to multiple linear 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 = "Color blind friendlier color palettes"}

The above `setup` chunk in this qmd includes code for using a more color blind friendly palette *throughout* this activity, without having to specify the palette anew for each `ggplot()`.
When you render the doc to an html, the palette will be used.
If you want to use the palette within the qmd (not just html), run the above chunk.

:::



\
\


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

- Identify some limitations of **simple linear regression**, i.e. utilizing only *one* predictor.
- Understand the general goals behind **multiple linear regression**.
- Explore strategies for *visualizing* and *interpreting* multiple linear regression models of Y vs 2 predictors.



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

:::



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

Today is a day to discover ideas.
There are no readings or videos.

:::




\
\
\
\



# Warm-up {-}

Before starting, load a couple datasets that we'll use in this activity:

```{r}
library(tidyverse)

# Load the cps data
cps <- read_csv("https://mac-stat.github.io/data/cps_2018.csv") %>% 
  select(wage, age, industry, marital) %>% 
  filter(age >= 18, age <= 34, wage < 250000)

# Load penguins data
# Get rid of Adelie penguins
# To get rid of an outlier, keep only penguins with bill lengths under 57mm
data(penguins)
penguins <- penguins %>% 
  filter(species != "Adelie", bill_len < 57)
```




\
\
\
\


**EXAMPLE 1: Review**

The monthly U.S. Current Population Survey (CPS), run by the Census Bureau and Bureau of Labor Statistics (BLS), is the "primary source of labor force statistics" in the U.S..
It's how the government tracks dynamics like wages and unemployment:

- [Learn more about the survey here.](https://www.census.gov/programs-surveys/cps.html)
- [Find the codebook for our data here.](https://cps.ipums.org/cps/resources/codebooks/cpsmar18.pdf)

Above, we imported `cps` data from 2018.
We limited the rows and columns of this dataset to study wages among *18-34 year olds* that make *under \$250,000*.

```{r}
# Check it out
head(cps)
```



\

Let's start by exploring the relationship of `wage` with job `industry`:

```{r}
# What's the relationship of wage with job industry?
# Fix this visualization
cps %>% 
  ggplot(aes(y = wage, x = industry)) + 
  geom_point() + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1))
```

```{r}
# Calculate the average wage in each industry
cps %>% 
  ___(___) %>% 
  summarize(mean(___))
```

```{r}
# Build & summarize the linear regression model of wage by industry
cps_model_1 <- lm(wage ~ industry, cps)
coef(summary(cps_model_1))
```

The resulting estimated model formula is...

E[`wage` | `industry`] = 27410 + 10242 `industryconstruction` + 8042 `industryinstallation_production` + 21332 `industrymanagement` - 2906 `industryservice` + 2676 `industrytransportation`


**FOLLOW-UP**

a. What's the "reference level" of the `industry` predictor?

b. Interpret the `(Intercept)` coefficient.

c. Interpret the `industryservice` coefficient.

d. Which industry has the highest average wage? Predict the wage for somebody in this industry.

```{r}

```

e. This model can't be represented by a *line*. So what does the model "look" like?!?

```{r}
# NOTE: DO NOT BE DISTRACTED BY THIS CODE :)
# This is merely for demonstrating a concept
cps_model_1 %>% 
  ggplot(aes(y = wage, x = industry)) + 
  geom_boxplot() + 
  geom_point(aes(x = industry, y = .fitted), color = "red") + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1))
```


f. Does the relationship of wage with industry appear to be strong? Provide numerical and graphical evidence.


\
\
\
\

::: {.callout-note title = "ANOVA & t-tests"}

A linear regression model of Y with 1 predictor X, where X is categorical, boils down to a set of average Y values, 1 for each category of X.
This makes sense!
For example:

- If 2 people are in the same job industry, and that's all we know about them, their predicted wage should be the same.
- And it makes sense to predict that a person's wage is at the average for their industry. 

This type of analysis, essentially a comparison of means among different groups or categories, is very, very common.
If you've taken AP stats: this type of model is the basis of things like "2-sample t-tests" and "ANOVA".
 
:::
 
\
\
\
\


::: {.callout-note title = "Interpreting models with ONE categorical predictor"}

Let Y be a quantitative response X be a categorical predictor with categories (A, B, C, D).
Then the linear regression model of Y by X is:

E[Y | X] = $\beta_0$ + $\beta_1$ XB + $\beta_2$ XC + $\beta_3$ XD

where...

- A is the **reference category** (by default because it's the first alphabetically)
- intercept $\beta_0$
    - contextual: the expected / average value of Y for category A
    - technical: the expected / average value of Y when XB, XC, XD are all 0 (hence X is category A!)
- XB coefficient $\beta_1$
    - not a slope!
    - contextual: the difference of the expected / average value of Y for category B vs category A
    - technical: the change in the expected / average value of Y when we increase XB by 1 from 0 to 1 (which corresponds to going from category A to B!)
- $\beta_2$ and $\beta_3$ are interpreted similarly

:::



\
\
\
\


**EXAMPLE 2: Motivation for multiple linear regression**

Next, consider the model of `wage` by `marital` status:

```{r}
# Build & summarize the model
wage_mod <- lm(wage ~ marital, data = cps)
coef(summary(wage_mod))
```

What do you / don't you conclude from this model?
How does it highlight the limitations of a simple linear regression model?



\
\
\
\






**EXAMPLE 3: More motivation**

Let's explore Chinstrap and Gentoo penguins:

```{r}
# Check it out
head(penguins)

# We took Adelies out of the data
penguins %>% 
  count(species)
```

Our goal is to build a model that we can use to predict flipper lengths.
First, consider a simple linear regression model of `flipper_len` by `bill_len`:

```{r}
penguins %>% 
  ggplot(aes(y = flipper_len, x = bill_len)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE)
```

Thoughts?
What's going on here?
How does this highlight the limitations of a simple linear regression model?







\
\
\
\



**EXAMPLE 4: Even more motivation**

Let's build 2 separate models of `flipper_len`:

```{r}
# Model flipper_len by bill_len
mod_bill <- lm(flipper_len ~ bill_len, data = penguins)
summary(mod_bill)
```

```{r}
# Model flipper_len by species
mod_species <- lm(flipper_len ~ species, data = penguins)
summary(mod_species)
```

The model of `flipper_len` by `bill_len` has an R-squared of 0.02881.
The model of `flipper_len` by `species` has an R-squared of 0.7014.
If we *had* to pick just 1 of these predictors, we would pick `species` -- 
it explains more about the variability in `flipper_len`.

BUT what if we didn't have to pick just 1 predictor?!
Suppose we model `flipper_len` by *both* `bill_len` and `species`.
What do you think is a reasonable range of possible values for the new R-squared?





\
\
\
\



::: {.callout-note title="Reflection: Why are multiple regression models so useful?"}

We can put more than 1 predictor into a regression model!

- Descriptive viewpoint   
    Adding predictors to a model can help us better understand the isolated (causal) effect of a variable by holding constant potential *confounders*.
    
- Predictive viewpoint    
    Adding predictors to a model can help us better predict the response variable.

:::



\
\
\
\



::: {.callout-note title="Multiple linear regression model formula"}

In general, a **multiple linear regression model** of $Y$ with multiple predictors $(X_1, X_2, ..., X_p)$ is represented by the following formula:

$$E[Y \mid X_1, X_2, ..., X_p] = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + ... + \beta_p X_p$$

:::







\
\
\
\



# Exercises {-}

In the coming weeks, we'll explore how to visualize, interpret, build, and evaluate multiple linear regression models.
Today, we'll focus on some foundations using a model of penguin `flipper_len` by just 2 predictors: `bill_len` (quantitative) and `species` (categorical).





\
\
\
\



## Exercise 1: Visualizing the relationship {-}

To explore the relationship of `flipper_len` with *both* `bill_len` and `species`, let's start with a visualization.
Modify the plot below to include information about `species`.
THINK: How can you change each point / dot to indicate the corresponding penguin's species?

```{r}
penguins %>% 
  ggplot(aes(y = flipper_len, x = bill_len)) + 
  geom_point()
```


\
\
\
\





## Exercise 2: Visualizing the model {-}


### Part a {-}

We've learned that a *simple* linear regression model of `flipper_len` by `bill_len` alone can be represented by a line:

```{r}
penguins %>% 
  ggplot(aes(y = flipper_len, x = bill_len)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE)
```

Reflecting on your plot of `flipper_len` by `bill_len` and `species` in Exercise 1, how do you think a *multiple* regression model of `flipper_len` using *both* of these predictors would be represented?


### Part b {-}

Check your intuition by modifying the code below to include `species` in this plot, as you did in Exercise 1.

```{r}
penguins %>%
  ggplot(aes(y = flipper_len, x = bill_len, ___ = ___)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)
```



\
\
\
\





## Exercise 3: Reflection {-}

Your plot in Exercise 2 demonstrated that the **multiple linear regression model** of `flipper_len` by `bill_len` and `species` is represented by 2 lines. 
Let's interpret the punchlines!
For each question below, provide an answer along with evidence from the model lines that supports your answer.

### Part a {-}

What's the relationship between `flipper_len` and `bill_len`, no matter a penguin's `species`?

What's your evidence?


### Part b {-}

Does the rate of increase in `flipper_len` with `bill_len` appear to differ between the two `species`?

What's your evidence?


### Part c {-}

What's the relationship between `flipper_len` and `species`, no matter a penguin's `bill_len`? 

What's your evidence?


### Part d {-}

The formula for this linear regression model is represented by the below formula:

E[`flipper_len` | `species`, `bill_len`] = $\beta_0$ + $\beta_1$ `bill_len` + $\beta_2$ `speciesGentoo`

Based on parts a-c: 

- Will $\beta_1$ be positive or negative? Why?

- Will $\beta_2$ be positive or negative? Why?




\
\
\
\


## Exercise 4: Model formula {-}

Of course, there's a formula behind the multiple regression model.
We can obtain this using the usual `lm()` function.

```{r}
# Build & summarize the model
penguin_mod <- lm(flipper_len ~ bill_len + species, data = penguins)
summary(penguin_mod)
```

### Part a {-}

Reflect on the code.
In the `lm()` function, how did we communicate that we wanted to model `flipper_len` by both `bill_len` and `species`?
(How is the `lm()` code different than when we used just 1 predictor?)


### Part b {-}

Complete the following model formula.
(Was your answer in Exercise 3 part d correct?!)
    
E[`flipper_len` | `bill_len`, `speciesGentoo`] = ___ + ___ `bill_len` + ___ `speciesGentoo`




\
\
\
\




## Exercise 5: Sub-model formulas {-}

We now have a single *formula* for the model:

E[`flipper_len` | `bill_len`, `speciesGentoo`] = 127.75 + 1.40 * `bill_len` + 22.85 * `speciesGentoo`

And we observed earlier that this formula is represented by two lines: one describing the relationship between `flipper_len` and `bill_len` for `Chinstrap` penguins and the other for `Gentoo` penguins.
Let's bring these ideas together.
Utilize the model formula to obtain the equations of these two lines, i.e. to obtain the *sub-model formulas* for the 2 species.
Simplify these as much as possible by combining constants when relevant!
Hint: Plug `speciesGentoo` = 0 and `speciesGentoo` = 1.        
    
Chinstrap: E[`flipper_len` | `bill_len`] = ___ + ___ `bill_len`

Gentoo:    E[`flipper_len` | `bill_len`] = ___ + ___ `bill_len`


```{r}
# Put any scratch calculations here

```

\
\
\
\




## Exercise 6: coefficients -- physical interpretation {-}

Reflecting on Exercise 5, let's interpret what the model coefficients tell us about the *physical* properties of the two 2 sub-model lines.
Choose the correct option given in parentheses:        

- The intercept coefficient, 127.75, is the intercept of the line for (Chinstrap / Gentoo) penguins.

- The `bill_len` coefficient, 1.40, is the (intercept / slope) of both lines.

- The `speciesGentoo` coefficient, 22.85, indicates that the (intercept / slope) of the line for Gentoo is 22.85mm higher than the (intercept / slope) of the line for Chinstrap. Similarly, since the lines are parallel, the line for Gentoo is 22.85mm higher than the line for Chinstrap at any `bill_len`.
    
    


\
\
\
\



## Exercise 7: coefficients -- contextual interpretation {-}

Next, interpret each coefficient in *context*.
What do they tell us about penguin flipper lengths?!

### Part a {-}

Interpret 127.75, the intercept of the Chinstrap line, *in context*.



### Part b {-}

Interpret 1.40, the slope of both lines, *in context*.
For *both* Chinstrap and Gentoo penguins, the expected flipper length...



### Part c {-}

Interpret 22.85, *in context*.
At any `bill_len`, we expect...
    
    


\
\
\
\




## Exercise 8: Prediction {-}

Now that we better understand the model, let's use it to predict flipper lengths!
Recall the model and visualization:

E[`flipper_len` | `bill_len`, `speciesGentoo`] = 127.75 + 1.40 * `bill_len` + 22.85 * `speciesGentoo`

```{r}
penguins %>% 
  ggplot(aes(y = flipper_len, x = bill_len, color = species)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE)
```

### Part a {-}

Predict the flipper length of a *Chinstrap* penguin with a *50mm* long bill.

```{r}
127.75 + 1.40*___ + 22.85*___
```    

To check your work, mark your prediction on the plot using a red dot:
Make sure it's where it should be!

```{r}
penguins %>% 
  ggplot(aes(y = flipper_len, x = bill_len, color = species)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE) + 
  geom_point(aes(x = 50, y = ____), color = "red", size = 4)
```

### Part b {-}

Predict the flipper length of a *Gentoo* penguin with a *50mm* long bill.

```{r}
127.75 + 1.40*___ + 22.85*___
```    

To check your work, mark your prediction on the plot:

```{r}
penguins %>% 
  ggplot(aes(y = flipper_len, x = bill_len, color = species)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE) + 
  geom_point(aes(x = 50, y = ____), color = "red", size = 4)
```

### Part c {-}

Use the `predict()` function to confirm your predictions in Parts a and b.

```{r}
# Confirm the calculation in part a
predict(penguin_mod,
        newdata = data.frame(bill_len = ___, species = "___"))

# Confirm the calculation in part b
predict(penguin_mod,
        newdata = data.frame(bill_len = ___, species = "___"))
```




\
\
\
\




## Exercise 9: R-squared {-}

Finally, recall that improving our predictions was one motivation to build a multiple linear regression model (using 2 predictors instead of 1).
To this end, recall the R-squared values of the simple linear regression models that used just *one* predictor at a time:

```{r}
# Model of flipper_len by bill_len
# R-squared = 0.02881
summary(mod_bill)

# Model of flipper_len by species
# R-squared = 0.7014
summary(mod_species)
```


### Part a {-}

What do you *guess* is the R-squared of `penguin_mod`, our multiple regression model that uses *both* `bill_len` and `species` as predictors?
(We discussed this same question in the Warm-up, but your opinion may have changed!)


### Part b {-}

Check your intuition.
Report the R-squared value of `penguin_mod` and discuss how it compares to the 2 simple linear regression models.
Why might one be surprised about this result?

```{r}
summary(penguin_mod)
```





\
\
\
\


## Exercise 10: Challenge Part 1 {-}

In this activity, we explored the model of `flipper_len` by a quantitative predictor (`bill_len`) *and* a categorical predictor (`species`).
Challenge yourself to anticipate what visualizations / models might be like with 2 *quantitative* predictors: `bill_len` and `bill_dep`.


### Part a {-}

First, let's visualize the relationship of `flipper_len` with `bill_len` and `bill_dep`.
To do so, start with the scatterplot of `flipper_len` and `bill_len`, and *modify* it.
THINK: How might you change the points to reflect their corresponding `bill_dep`?

```{r}
# Do NOT include a geom_smooth (which wouldn't make sense here)
penguins %>% 
  ggplot(aes(y = flipper_len, x = bill_len)) + 
  geom_point()
```

### Part b {-}

In Part a, you represented a 3D point cloud in 2D.
BONUS (just for fun & for better understanding that this is a 3D relationship): Visualize this same data in something more like 3D:

```{r}
# Don't worry about this code!
# It's specialized to this exercise and we'll rarely use it
library(plotly)
penguins %>% 
  plot_ly(x = ~bill_len, y = ~bill_dep, z = ~flipper_len,
          type = "scatter3d")
```


### Part c {-}

We can think of your plot above as a 3D point cloud.
And the estimated formula for the linear regression model is

E[`flipper_len` | `bill_len`, `bill_dep`] = 179.2 + 2.38 `bill_len` - 5.15 `bill_dep`

```{r}
quant_model <- lm(flipper_len ~ bill_len + bill_dep, penguins)
coef(summary(quant_model))
```

**QUESTION:** 

The model of `flipper_len` by `bill_len` and `species` is represented by 2 parallel lines.
How would you draw / represent the model here?!


### Part d {-}

Challenge: Interpret the `(Intercept)` and `bill_len` coefficients.
What *physical* meaning do they have, i.e. where do they come into the "drawing" of the model?
What *contextual* meaning do they have?



\
\
\
\


## Exercise 11: Challenge Part 2 {-}

Next, challenge yourself to anticipate what visualizations / models might be like with 2 *categorical* predictors: `species` and `sex`.


### Part a {-}

First, let's visualize the relationship of `flipper_len` with `species` and `sex`.
To do so, start with the side-by-side boxplots of `flipper_len` by `species`, and *modify* them.

```{r}
penguins %>% 
  ggplot(aes(y = flipper_len, x = species)) + 
  geom_boxplot()
```

### Part b {-}

The estimated formula for the linear regression model is

E[`flipper_len` | `species`, `sex`] = 191.8 + 21.07 `speciesGentoo` + 8.39 `sexmale`

```{r}
cat_model <- lm(flipper_len ~ species + sex, penguins)
summary(cat_model)
```

**QUESTION:** 

The model of `flipper_len` by `species` alone would only produce 2 possible predictions, one for Chinstraps and one for Gentoos.
What about the model here?!
How many possible predictions would it make?
HINT: Checking out the plot again might help.


### Part c {-}

Challenge: Interpret the `(Intercept)` and `speciesGentoo` coefficients.
What *physical* meaning do they have, i.e. where do they come into the "drawing" of the model?
What *contextual* meaning do they have?


\
\
\
\




# Reflection {-}

You've now explored your first multiple regression model!
Thus you likely have a lot of questions about what's to come.
What are they?






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


