---
title: "Multiple regression principles (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"}

We're still talking about linear regression models!
But now, we have *multiple* predictors in our models.
Accordingly, the plots and models of the relationships of interest can get complicated.
Though there are no *recipes* for this process, there are some guiding *principles* that assist in long-term retention, deeper understanding, and the ability to generalize our tools in new settings.
Today we'll explore some general *principles* for... 

- incorporating additional quantitative or categorical predictors in a visualization
- how additional quantitative or categorical predictors impact the physical representation of a model
- interpreting quantitative or categorical coefficients in a multiple regression model

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

:::


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

**Required video**

- [Multivariate modeling principles](https://youtu.be/iNv4GJmvrOY) ([slides](https://drive.google.com/file/d/1QIikKI_SqXyLFDYnckvCLj56k5zEjiBI/view?usp=sharing))


**Optional**

The below videos provide extra examples for interpreting multiple linear regression models.
You're encouraged to watch them *after* today's activity.

- [Interpreting multivariate models (Part I)](https://youtu.be/DnXwu1OWgMM) and [Part II](https://youtu.be/vQ_QZcZCEMU) ([slides](https://drive.google.com/file/d/1fiJL1IbReg6RJmjxRideYzYTyQtWR3Or/view?usp=sharing))


:::




\
\
\
\




# Warm-up {-}

Let's revisit the bikeshare data:

```{r}
# Load packages & import data
library(tidyverse)
bikes <- read_csv("https://mac-stat.github.io/data/bikeshare.csv") %>% 
  rename(rides = riders_registered)

# Check out the data
head(bikes)
```

Our goal is to explore how registered ridership varies from day to day.
To this end, we'll build various **multiple linear regression** models of `rides` by different combinations of the possible predictors.


\
\
\
\


**EXAMPLE 1: Visualization review**

Let's explore the relationship of `rides` with `windspeed` (quantitative) and `weekend` status (categorical):

```{r}
bikes %>% 
  select(rides, windspeed, weekend) %>% 
  head()
```


a. Write a model statement for this regression model.

E[rides | windspeed, weekend] = $\beta_0$ + $\beta_1$ windspeed + $\beta_2$ ___

b. Construct and discuss a plot of the relationship between these 3 variables. NOTE: Make sure that your description addresses all variables here, as well as the strength of the overall relationship.

```{r}
# Plot of rides vs windspeed & weekend
# Include a representation of the linear regression model
# PRO TIP: Start with a plot of rides by windspeed, then modify it
bikes %>% 
  ggplot(aes(y = rides, x = windspeed)) + 
  geom___() + 
  geom_smooth(method = "lm", se = FALSE)
```



\
\
\
\



**EXAMPLE 2: Model review**

Let's build the linear regression model of this relationship:

```{r}
bike_model_1 <- lm(rides ~ windspeed + weekend, data = bikes)
summary(bike_model_1)
```

Thus estimated model formula is:

E[rides | windspeed, weekend] = 4738.38 - 63.97 windspeed - 925.16 weekendTRUE

This model formula is represented by 2 lines, one corresponding to weekends and the other to weekdays. Simplify the model formula above for weekdays and weekends:       
    
weekdays: rides = ___ - ___ windspeed   

weekends: rides = ___ - ___ windspeed

```{r}
# Record any calculations here

```


\
\
\
\



**EXAMPLE 3: Coefficient interpretation review**


a. The intercept coefficient, 4738.38, represents the *intercept* of the sub-model for weekdays, the reference category. What's its *contextual* interpretation?
    

b. The `windspeed` coefficient, -63.97, represents the *shared slope* of the weekend and weekday sub-models. What's its *contextual* interpretation?   

    

c. The `weekendTRUE` coefficient, -925.16, represents the *change in intercept* for the weekend vs weekday sub-model. What's its *contextual* interpretation?



\
\
\
\


::: {.callout-note title = "Interpreting MLR coefficients (2 predictors)"}

Consider a multiple linear regression model with 2 predictors:

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

Then

- $\beta_0$ ("beta 0") = expected value of $Y$ when $X_1$ and $X_2$ are both 0, ie. when all quantitative predictors are set to 0 and the categorical predictors are set to their reference levels.    

- $\beta_1$ ("beta 1") = change in the expected value of $Y$ for two cases whose values of $X_1$ differ by 1 but who are identical in terms of $X_2$.

:::

\



::: {.callout-note title="OPTIONAL Math Box: Fitting multiple linear regression models"}

Consider a multiple linear regression model with $p$ predictors:

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

We use the *least squares criterion* to estimate $(\beta_0, \beta_1, ..., \beta_p)$, i.e. we pick the coefficients *minimize the sum of squared residuals*.
Just as with simple linear regression, we can write out a formula for these least square estimates.
But for multiple linear regression, this requires linear algebra.

Suppose we have $n$ data points.
Using matrix notation, we can collect our observed response values into a vector $Y$, our predictor values into a matrix $X$, and our regression coefficients into a vector $\beta$.
The column of 1s in $X$ reflects our model's inclusion of an intercept:

$$Y = \left(\begin{array}{c} Y_1 \\ Y_2 \\ \vdots \\ Y_n \end{array} \right) \;\;\;\; \text{ and } \;\;\;\;
X = \left(\begin{array}{ccccc} 
1 & X_{11} & X_{12} & \cdots & X_{1p} \\
1 & X_{21} & X_{22} & \cdots & X_{2p} \\
\vdots & \vdots & \vdots & \cdots & \vdots \\
1 & X_{n1} & X_{n2} & \cdots & X_{np} \\
\end{array} \right) \;\;\;\; \text{ and } \;\;\;\; \beta = \left(\begin{array}{c} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_p \end{array} \right)$$


Then we can express the model $E[Y_i|X_{i1}, X_{i2},...,X_{ip}] = \beta_0 + \beta_1 X_{i1} + \cdots + \beta_p X_{ip}$ for $i \in \{1,...,n\}$ as:

$$E[Y|X] = X\beta$$
Further, let $\hat{\beta}$ denote the vector of sample estimated $\beta$ and $\hat{Y}$ denote the vector of predictions / model values:

$$\hat{Y} = X \hat{\beta}$$

Thus the sum of squared residuals is

$$\sum_{i=1}^n (Y_i - \hat{Y}_i)^2 = (Y - \hat{Y})^T(Y - \hat{Y}) = (Y - X \hat{\beta})^T(Y - X \hat{\beta})$$

and the following formula for sample coefficients $\hat{\beta}$ are the **least squares estimates** of $\beta$, i.e. they *minimize* the sum of squared residuals:

$$\hat{\beta} = (X^TX)^{-1}X^TY$$



:::



\
\
\
\



# Exercises {-}

Thus far, we've explored a couple examples of multiple regression models that have 2 predictors, 1 quantitative and 1 categorical.
Let's explore what happens when *both* predictors are categorical or *both* are quantitative!

\


## Exercise 1: 2 categorical predictors -- visualization {-}

Let's explore the relationship of `rides` with 2 *categorical* predictors, `weekend` status and `season`.
The below code plots `rides` vs just `season`.
*Modify* this code to *also* include information about `weekend`.
HINT: Remember the visualization *principle* that additional categorical predictors require some sort of grouping mechanism that distinguishes between the groups.

```{r}
# Modify this code to plot rides vs season AND weekend
bikes %>% 
  ggplot(aes(y = rides, x = season)) + 
  geom_boxplot()
```



\
\
\
\




## Exercise 2: follow-up {-}

### Part a {-}

Describe, in words, the relationship of ridership with season & weekend status.
In doing so, address the following:

- No matter the weekend status, how is ridership associated with season?
- No matter the season, how is ridership associated with weekend status?
- Does the overall relationship of ridership with season & weekend status appear to be *strong*? 




### Part b {-}

A model of `rides` by `season` alone would be represented by only 4 expected outcomes, 1 for each season.
Considering this and the plot above, how do you *anticipate* a model of `rides` by `season` and `weekend` status will be represented?        

- 2 lines, 1 for each weekend status
- 8 lines, 1 for each possible combination of season & weekend
- 2 expected outcomes, 1 for each weekend status
- 8 expected outcomes, 1 for each possible combination of season & weekend
    


\
\
\
\



## Exercise 3: 2 categorical predictors -- modeling {-}

Build the multiple regression model of `rides` vs `season` and `weekend`:
    
```{r}
bike_model_2 <- lm(rides ~ weekend + season, bikes)
summary(bike_model_2)
```

Thus the estimated model formula is:

E[rides | weekend, season] = 4260.45 - 912.33 weekendTRUE - 116.38 seasonspring + 438.44 seasonsummer - 1719.06 seasonwinter

### Part a {-}

Use this model to make the following predictions:
    
```{r}
# predicted ridership on a fall weekday
4260.45 - 912.33*___ - 116.38*___  + 438.44*___ - 1719.06*___

# predicted ridership on a winter weekday    
4260.45 - 912.33*___ - 116.38*___  + 438.44*___ - 1719.06*___

# predicted ridership on a winter weekend day
4260.45 - 912.33*___ - 116.38*___  + 438.44*___ - 1719.06*___
```

Using the calculations in the chunk above, calculate the *differences* between various predictions:

```{r}
# diff in predicted ridership on a winter weekday vs a fall weekday
# (that is, the diff btwn 2 days that are at the same time of week but DIFFERENT seasons)
# WHERE HAVE YOU OBSERVED THIS NUMBER BEFORE?!
___ - ___
```

```{r}
# diff in predicted ridership on a winter weekend vs a winter weekday
# (that is, the diff btwn 2 days that are in the same season but DIFFERENT times of the week)
# WHERE HAVE YOU OBSERVED THIS NUMBER BEFORE?!
___ - ___
```


### Part b {-}

We only made 3 predictions here.
In general there are *8* possible predictions: 2 weekend categories x 4 season categories.
Is this consistent with your intuition in the previous exercise?



\
\
\
\



## Exercise 4: 2 categorical predictors -- interpret the model {-}

Use your above predictions and visualization to fill in the below interpretations of the model coefficients.
Hint: What is the consequence of plugging in 0 or 1 for the different `weekend` and `season` categories?    


- Interpreting 4260: On average, there are 4260 riders on (weekdays/weekends) during the (fall/spring/summer/winter), i.e. when all predictors are "0" 

- Interpreting -912: On average, *in any season*, there are 912 (more/fewer) riders on weekends than on ___. 

- Interpreting -1719: On average, *on both weekdays and weekends*, we expect there to be 1719 (more/fewer) riders in winter than in ___.





\
\
\
\




## Exercise 5: 2 quantitative predictors -- visualization {-}

Next, consider the relationship of `rides` with 2 *quantitative* predictors: `windspeed` and `temp_feel`.
Check out the plot of this relationship below which reflects the visualization *principle* that quantitative variables require some sort of numerical scaling mechanism -- `rides` and `windspeed` get numerical axes, and `temp_feel` gets a color scale:

![](https://mac-stat.github.io/images/155/bikes_multivar_3.png)

Modify the code below to recreate this plot:

```{r}
bikes %>%
  ggplot(aes(___)) +
  geom_point()
```

OPTIONAL: Check out this *interactive* plot which allows us to explore this point cloud in 3D.

```{r}
library(plotly)
bikes %>% 
  plot_ly(x = ~windspeed, y = ~temp_feel, z = ~rides, 
          type = "scatter3d", 
          mode = "markers",
          marker = list(size = 5, color = ~rides, colorscale = "Viridis"))
```




\
\
\
\






## Exercise 6: follow-up {-}

Describe (in words) the relationship of ridership with windspeed & temperature.
Make sure to comment on all 3 variables, strength of the relationship, and any outliers (if relevant).


    



\
\
\
\



## Exercise 7: 2 quantitative predictors -- modeling {-}

Let's build the multiple regression model of `rides` vs `windspeed` and `temp_feel`:
    
```{r}
bike_model_3 <- lm(rides ~ windspeed + temp_feel, data = bikes)
summary(bike_model_3)
```

Thus the estimated model formula is:

E[rides | windspeed, temp_feel] = -24.06 - 36.54 windspeed + 55.52 temp_feel

### Part a {-}

Interpret the intercept coefficient, -24.06, in context.


### Part b {-}

Interpret the `windspeed` coefficient, -36.54, in context. 


### Part c {-}

Interpret the `temp_feel` coefficient, 55.52, in context. 



    

\
\
\
\



## Exercise 8: Which model is "best"? {-}

We've now observed 3 different models of ridership, each having 2 predictors.
The R-squared values of these models, along with those of the simple linear regression models with each predictor alone, are summarized below.
    
model                     predictors                R-squared
------------------------- ------------------------- ------------
`bike_model_1`            `windspeed` & `weekend`   0.119
`bike_model_2`            `weekend` & `season`      0.349
`bike_model_3`            `windspeed` & `temp_feel` 0.310
`bike_model_4`            `windspeed`               0.047
`bike_model_5`            `temp_feel`               0.296
`bike_model_6`            `weekend`                 0.074
`bike_model_7`            `season`                  0.279


### Part a {-}

Which model does the best job of explaining the variability in ridership from day to day?


### Part b {-}

If you could only pick one predictor, which would it be?


### Part c {-}

What happens to R-squared when we add a second predictor to our model, and why does this make sense?
For example, how does the R-squared for model 1 (with both `windspeed` and `weekend`) compare to those of model 4 (only `windspeed`) and model 6 (only `weekend`)?


### Part d {-}

Are 2 predictors always better than 1?
Provide evidence for your answer and explain why this makes sense.




\
\
\
\




## Exercise 9: Principles of interpretation {-}

These exercises have confirmed the **principles** behind interpreting model coefficients, summarized below.
Review and confirm that these make sense.

\

::: {.callout-note title = "Interpretating MLR coefficients (2+ predictors)"}

Consider a multiple linear regression model:

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

- $\beta_0$ ("beta 0") = expected value of $Y$ when the $X_i$ are *all* 0, ie. when all quantitative predictors are set to 0 and the categorical predictors are set to their reference levels.    

- $\beta_1$ ("beta 1") = change in the expected value of $Y$ for two cases whose values of $X_1$ differ by 1 but who are identical in terms of all other $X$:

    - If $X_1$ is quantitative, $\beta_1$ describes the change in the expected value of $Y$ associated with a 1-unit increase in $X_1$, while at a fixed set of the other $X$.  
    
    - If $X_1$ represents a category of a categorical variable, $\beta_1$ describes the change in the expected value of $Y$ between this category and the reference category, while at a fixed set of the other $X$.    

:::




\
\
\
\



## Exercise 10: Extra Practice 1 {-}

The following exercises provide extra practice.
Consider the relationship of `rides` vs `weekend` and `weather_cat`.    

- Construct a visualization of this relationship.   
- Construct a model of this relationship.    
- Interpret the first 3 model coefficients. 



\
\
\
\

## Exercise 11: Extra Practice 2 {-}

Consider the relationship of `rides` vs `temp_feel` and `humidity`.    

- Construct a visualization of this relationship.   
- Construct a model of this relationship.    
- Interpret the first 3 model coefficients. 


\
\
\
\


## Exercise 12: Extra Practice 3 {-}

Consider the relationship of `rides` vs `temp_feel` and `weather_cat`.    

- Construct a visualization of this relationship.   
- Construct a model of this relationship.    
- Interpret the first 3 model coefficients. 



\
\
\
\

## Exercise 13: Optional CHALLENGE {-}

We've explored models with 2 predictors.
What about 3 predictors?!
Consider the relationship of `rides` vs `temp_feel`, `humidity`, AND `weekend`.

- Construct a visualization of this relationship.   
- Construct a model of this relationship.    
- Interpret *each* model coefficient. 


\
\
\
\


## Exercise 14: Explore more examples {-}

Outside of class, check out the below videos that include extra examples:

- [Interpreting multivariate models (Part I)](https://youtu.be/DnXwu1OWgMM) and [Part II](https://youtu.be/vQ_QZcZCEMU) ([slides](https://drive.google.com/file/d/1fiJL1IbReg6RJmjxRideYzYTyQtWR3Or/view?usp=sharing))






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


