Week 08: Multiple Regression I

Date: October 15, 2026

Today…

One predictor was never the plan.

  • Why a second predictor changes the first one’s meaning
  • “Holding constant” — what it actually does
  • A slope that flips sign when we add the right variable
  • Did the new predictor earn its place? (\(\Delta R^2\), nested \(F\))
  • b or β: which one do you report?
  • When predictors step on each other (multicollinearity)
library(here)
library(tidyverse)
library(rio)
library(broom)
library(psych)
library(patchwork)

options(scipen = 999)

Today…

https://bsky.app/profile/andrew.heiss.phd/post/3ly6hqr3mtk2z

Where we have been

Week Model What was new
5 \(Y = b_0 + b_1X\) one continuous predictor
6 \(Y = b_0 + b_1D\) the predictor can be a group
7 \(SS_{total} = SS_{model} + SS_{res}\) how well the line fits

Everything so far has been one predictor at a time. That is almost never the model you want.

Regression

Overall, we are providing a model to give us a “best guess” on predicting our outcome

\[ Y_i = b_0 + b_1X_i + e_i \]

Explaining Variance

Explaining Variance

Explaining Variance

Model Interpretation

Once we have a model, we will be able to interpret the coefficients. For a bivariate regression (simple regression), this was fairly straightforward

We would look to the beta (effect size) and interpret it as: for every 1 unit increase in our X variable, there will be beta units increase in our Y variable.

The only equation we need

\[Y_i = b_0 + b_1X_{1i} + b_2X_{2i} + \dots + b_kX_{ki} + e_i\]

  • \(b_0\) — predicted \(Y\) when every predictor is 0
  • \(b_1\) — change in \(Y\) per 1-unit change in \(X_1\), holding the other predictors constant
  • \(e_i\) — still the part we miss

Important

That phrase — holding the other predictors constant — is the entire week. It is not decoration. It changes what the number means.

Explaining Variance

Part 1: Why Bother?

You have been hired by the county 🚒

The fire department wants to know what drives the cost of a fire.

They hand you 300 fires: how many trucks were dispatched, how many acres burned, and the damage in thousands of dollars.

Code 
set.seed(640)
n <- 300

acres  <- round(runif(n, 0.5, 40), 1)
trucks <- round(1 + 0.22 * acres + rnorm(n, 0, 1.2))
trucks[trucks < 1] <- 1
damage <- 10 + 4.2 * acres - 2.8 * trucks + rnorm(n, 0, 12)

fire_data <- tibble(acres, trucks, damage)
glimpse(fire_data)
Rows: 300
Columns: 3
$ acres  <dbl> 10.3, 8.0, 29.4, 31.7, 8.0, 25.9, 12.5, 38.6, 29.3, 14.0, 32.1,…
$ trucks <dbl> 1, 3, 8, 7, 2, 7, 2, 10, 8, 3, 6, 6, 8, 8, 1, 4, 1, 6, 7, 11, 1…
$ damage <dbl> 68.4656354, 29.8639841, 98.0352116, 123.1168377, 42.6368658, 10…

Start where we always start

Code 
fire_data %>%
  ggplot(aes(x = trucks, y = damage)) +
  geom_point(alpha = .5) +
  geom_smooth(method = "lm", se = FALSE, color = "#B3261E") +
  labs(x = "Trucks dispatched", y = "Damage ($1,000s)") +
  theme_bw(base_size = 16)

Bivariate relationship

This relationship can also be represented by a correlation, or a regression

cor.test(fire_data$trucks, fire_data$damage)

    Pearson's product-moment correlation

data:  fire_data$trucks and fire_data$damage
t = 24.728, df = 298, p-value < 0.00000000000000022
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.7790614 0.8539156
sample estimates:
      cor 
0.8199638 
## OR
fire.lm1 <- lm(damage ~ trucks, data = fire_data)
summary(fire.lm1)

Call:
lm(formula = damage ~ trucks, data = fire_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-73.661 -16.168   1.857  15.619  76.910 

Coefficients:
            Estimate Std. Error t value             Pr(>|t|)    
(Intercept)   8.8397     3.2784   2.696              0.00741 ** 
trucks       12.7299     0.5148  24.728 < 0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 25.26 on 298 degrees of freedom
Multiple R-squared:  0.6723,    Adjusted R-squared:  0.6712 
F-statistic: 611.5 on 1 and 298 DF,  p-value: < 0.00000000000000022

Bivariate relationship

Based on the visualization, the correlation and the regression, what would we conclude from this analysis?

Bivariate relationship

A linear regression was used to examine the impact of number of trucks on the scene and the amount of damage that was done by a fire. The number of trucks positively impacted the damage done by the fire (b = 11.67, p < .001), suggesting that as more trucks are present, the fire does more damage.

“Trucks strongly predict damage. Send fewer trucks.”

Warning

Every number in that sentence is correct. The recommendation would get people killed.

🧠 Side quest: draw the thing

Before you fit anything, sketch what you think is going on. Arrows for “causes.”

  • Does trucks → damage?
  • Does something cause both?

A variable that causes both your predictor and your outcome is a confound. Leave it out and your slope absorbs its effect.

These sketches have a name. Remember when we talked about Directed Acyclic Graphs (DAG’s)?

🧠 Side quest: draw the thing

The variable we ignored

Code 
fire_data %>%
  ggplot(aes(x = trucks, y = damage, color = acres)) +
  geom_point(alpha = .8) +
  scale_color_viridis_c(option = "inferno", direction = -1) +
  labs(x = "Trucks dispatched", y = "Damage ($1,000s)",
       color = "Acres\nburned") +
  theme_bw(base_size = 16)

Big fires get more trucks and do more damage. The slope we estimated was mostly fire size wearing a truck costume.

Put acres in the model

fire.lm2 <- lm(damage ~ trucks + acres,
               data = fire_data)
summary(fire.lm2)

Call:
lm(formula = damage ~ trucks + acres, data = fire_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-34.190  -9.226   0.103   8.002  40.010 

Coefficients:
            Estimate Std. Error t value             Pr(>|t|)    
(Intercept)  10.9421     1.6520   6.624       0.000000000164 ***
trucks       -2.6306     0.5794  -4.541       0.000008163552 ***
acres         4.0861     0.1378  29.645 < 0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.72 on 297 degrees of freedom
Multiple R-squared:  0.9172,    Adjusted R-squared:  0.9167 
F-statistic:  1646 on 2 and 297 DF,  p-value: < 0.00000000000000022
  • trucks: +12.73 → −2.63. The sign flipped.
  • acres: +4.09 per acre
  • \(R^2\): .672 → .917

Each additional truck is associated with $2,630 less damage — among fires of the same size.

What just happened

Before

We compared a 2-truck fire to a 9-truck fire. Those are not the same fire. One was a kitchen; one was a neighborhood.

After

We compare 2-truck and 3-truck fires of the same acreage. Now the extra truck can only help.

“Holding constant” is a comparison, not a magic word

\(b_1\) answers: among cases that are identical on every other predictor in the model, what is the difference in \(Y\) per unit of \(X_1\)?

Identical on the predictors you included. Nothing else.

Fire Data: Formal Writeup

“The current set of analyses sought to examine the influence of the number of firetrucks on the scene and the amount of damage that a fire caused, controlling for the size of the fire. The damage amount (in \$1,000s) was regressed on the number of firetrucks present and the size of the fire (in acres burned). The overall multiple regression was statistically significant ( \(R^2\) = .92, p < .001), meaning the two variables (Number of Trucks and Size of the Fire) accounted for 92% of the variance in damage done. Each of the two independent variables also had a statistically significant effect on damage (p’s < .001). The number of trucks was negatively related to the amount of damage done (b = -2.63), meaning that for each additional truck that was on the scene, there was a $2,630 reduction in damage, controlling for the size of the fire. Additionally, the size of the fire increased the amount of damage by approximately $4,000 for each additional acre burned.”

Omitted variables do not politely disappear

When you leave out a variable that is related to both \(X\) and \(Y\), its effect does not vanish — it gets folded into the coefficients you did estimate.

  • Sometimes the slope is inflated (the trucks case)
  • Sometimes it is shrunk (next week’s tutoring example)
  • Sometimes it appears out of nowhere (Lab 8’s suppression)
  • Sometimes nothing moves — and that is a finding too

You cannot tell which from the output. You can only tell from what you chose to measure.

Goldilocks Problem

Models that are underfit, overfit or just-right

Another Example: Computer time & Frustration 💻 😠

Note

This example will focus mostly on the model comparisons

Code 
# Step 1: Simulate the Data 
n <- 150 # Number of students

# Simulate total daily computer use in hours
computer_use <- rnorm(n, mean = 6, sd = 1.5)

# Simulate hours spent on stats homework, now as a whole number
stats_homework <- round(0.7 * computer_use + rnorm(n, 0, 1))
# Ensure no negative values after rounding
stats_homework[stats_homework < 0] <- 0
stats_homework[stats_homework > computer_use] <- 1

# Simulate frustration level (driven mostly by stats homework)
# Scale of 1-100
frustration <- 10 + (2 * computer_use) + (8 * stats_homework) + rnorm(n, 0, 10)

# Create a dataframe and clean up any impossible values
student_data <- tibble(
  frustration = frustration,
  computer_use = computer_use,
  stats_homework = stats_homework
) %>%
  filter(frustration > 0, computer_use > 0)

# Bin the 'stats_homework' variable
student_data <- student_data %>%
  mutate(
    stats_homework_bin = cut(stats_homework,
                           breaks = quantile(stats_homework, probs = seq(0, 1, by = 1/3)),
                           include.lowest = TRUE,
                           labels = c("Low", "Medium", "High"))
  )

We have collected data on graduate students and their weekly frustration levels. Data that were also collected included the amount of time they were spending on their computer as well as the amount of statistics homework that they had.

Frustration Model

Taking a bivariate look at things, we might have an initial model like this:

\[ Frustration=b_0 + b_1Comp.Hours + e_i \]

Frustration Model

Code 
p1 <- student_data %>%
  ggplot(aes(x = computer_use, y = frustration)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm",
              se = FALSE,
              color = "darkorange") +
  labs(
    title = "Frustration & Computer Use",
    subtitle = "Linear Fit",
    x = "Daily Computer Use (Hours)",
    y = "Frustration Level (1-100)"
  ) +
  theme_minimal()

p2 <- student_data %>%
  ggplot(aes(x = computer_use, y = frustration)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess",
              se = FALSE,
              color = "darkorange") +
  labs(
    title = "Frustration & Computer Use",
    subtitle = "Loess Fit",
    x = "Daily Computer Use (Hours)",
    y = "Frustration Level (1-100)"
  ) +
  theme_minimal()

p1 + p2

What are the correlations?

student_data %>% select(-stats_homework_bin) %>% cor() %>% round(3)
               frustration computer_use stats_homework
frustration          1.000        0.717          0.845
computer_use         0.717        1.000          0.751
stats_homework       0.845        0.751          1.000

Well that doesn’t look nice…

What are the correlations? 💅

https://strengejacke.github.io/sjPlot/reference/tab_corr.html

student_data %>% 
  select(-stats_homework_bin) %>% 
  sjPlot::tab_corr(var.labels = c("Frustration", "Computer Use", 
                                  "Statistics Homework"),
                   triangle = "lower")
  Frustration Computer Use Statistics Homework
Frustration      
Computer Use 0.717***    
Statistics Homework 0.845*** 0.751***  
Computed correlation used pearson-method with listwise-deletion.

Side Quest: Model Assumptions

  1. Linearity - A linear relationship exists between predictors and outcome
  2. Multicollinearity - The variables are minimally related to one another (VIF)
  3. Independence of Observations - Each observation in the model is independent of one another
  4. Homoscedasticity of Residuals - Residuals have constant variance across all points of the model
  5. Multivariate Normality - Residuals of the model are normally distributed (QQ Plots)

All of these are checked when using check_model() within the easystats library

Frustration Model: lm()

frus.lm1 <- lm(frustration ~ computer_use, data = student_data)
summary(frus.lm1)

Call:
lm(formula = frustration ~ computer_use, data = student_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-30.562  -8.140   1.604   7.439  33.613 

Coefficients:
             Estimate Std. Error t value            Pr(>|t|)    
(Intercept)   -1.1840     4.5220  -0.262               0.794    
computer_use   9.2552     0.7427  12.462 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.14 on 147 degrees of freedom
Multiple R-squared:  0.5137,    Adjusted R-squared:  0.5104 
F-statistic: 155.3 on 1 and 147 DF,  p-value: < 0.00000000000000022

Frustration Model: Multivariate

Code 
ggplot(student_data, aes(x = computer_use, y = frustration, color = stats_homework_bin)) +
  geom_point(alpha = 0.8, size = 2) +
  labs(
    title = "Frustration & Computer Use by Stats",
    x = "Daily Computer Use (Hours)",
    y = "Frustration Level (1-100)",
    color = "Stats HW (Hours)"
  ) +
  theme_minimal()

Frustration Model: Multivariate

Code 
student_data %>% 
  ggplot(aes(x = computer_use, y = frustration, 
             color = stats_homework_bin)) +
    geom_point(alpha = 0.5) +
    # Fit a line FOR EACH GROUP to show the true relationship
    geom_smooth(method = "lm", se = FALSE) +
    labs(
      title = "Frustration & Computer Use by Stats",
      x = "Daily Computer Use (Hours)",
      y = "Frustration Level (1-100)",
      color = "Stats HW (Hours)"
    ) +
    theme_bw()

Frustration Model

frus.lm2 <- lm(frustration ~ computer_use + stats_homework,
               data = student_data)
summary(frus.lm2)

Call:
lm(formula = frustration ~ computer_use + stats_homework, data = student_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-19.5242  -7.7787   0.5553   6.0824  23.6036 

Coefficients:
               Estimate Std. Error t value             Pr(>|t|)    
(Intercept)      6.0628     3.4473   1.759              0.08073 .  
computer_use     2.4187     0.8416   2.874              0.00466 ** 
stats_homework   8.3376     0.7712  10.811 < 0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.828 on 146 degrees of freedom
Multiple R-squared:  0.7299,    Adjusted R-squared:  0.7262 
F-statistic: 197.3 on 2 and 146 DF,  p-value: < 0.00000000000000022

Comparing Models

Maybe we would like to test to see if including the Stats Homework variable improved the model

Since the first regression is nested within the second, we can run an ANOVA

anova(frus.lm1, frus.lm2)
Analysis of Variance Table

Model 1: frustration ~ computer_use
Model 2: frustration ~ computer_use + stats_homework
  Res.Df   RSS Df Sum of Sq      F                Pr(>F)    
1    147 25391                                              
2    146 14102  1     11289 116.88 < 0.00000000000000022 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Statistic Meaning Interpretation
Res.Df (Residual Degrees of Freedom) Remaining degrees of freedom after estimating parameters 148 for Model 1, 147 for Model 2 (one more parameter estimated)
RSS (Residual Sum of Squares) Total unexplained variance by the model Decreases substantially from 24,468 → 13,620 with the new predictor, showing improved fit
Df (Difference in df) Number of parameters added 1 (adding stats_homework)
Sum of Sq (Difference in RSS) Variance explained by the added predictor 10,848 units of variance explained by stats_homework
F F-statistic testing whether this variance reduction is significant 117.07 — very large, meaning substantial improvement
Pr(>F) p-value associated with F

Model Assumptions: Frustration

performance::check_model(frus.lm2)

Model Assumptions

  1. Linearity - A linear relationship exists between predictors and outcome
    1. Apply Transformation (square root or log); Drop the variable
  2. Multicollinearity - The variables are minimally related to one another (VIF)
    1. Remove variable; different type of regression (lasso, partial least squares)
  3. Independence of Observations - Each observation in the model is independent of one another
    1. Dummy coding, lags; Use MLM
  4. Homoscedasticity of Residuals - Residuals have constant variance across all points of the model
    1. Apply Transformation (square root or log);
  5. Multivariate Normality - Residuals of the model are normally distributed (QQ Plots)
    1. Check for outliers; Apply Transformation (square root or log)

Model Selection - Stepwise

  • Backward elimination

    • Start with the full model (all potential predictors)

    • Predictors eliminated one-at-a-time from the model until we cannot improve the model any further

  • Forward selection

    • Start with nothing in the model

    • Continue to add predictors one-at-a-time until we cannot find any predictors that improve the model any further

Model Selection

ALWAYS start with theory and previous work

You should NEVER just go into data with just a lot of predictors and backward selection in mind. 🧨

One More Example - Bringing it together

520 students, one survey 😵‍💫

A student-stress survey: everything rated 1–5. The column names are, charitably, a disaster — so we rename by position.

stress <- import(here("files", "data", "student_stress.xlsx")) %>%
  set_names(c("sleep_quality", "headaches", "academic_perf",
              "study_load", "extracurricular", "stress"))

stress %>% describe(fast = TRUE)
                vars   n mean   sd median min max range  skew kurtosis   se
sleep_quality      1 520 3.12 1.10    3.0   1   5     4  0.01    -0.66 0.05
headaches          2 520 2.18 1.25    2.0   1   5     4  0.75    -0.48 0.05
academic_perf      3 520 3.33 1.06    3.0   1   5     4 -0.34    -0.28 0.05
study_load         4 520 2.75 1.37    2.5   1   5     4  0.23    -1.25 0.06
extracurricular    5 520 2.68 1.47    3.0   1   5     4  0.32    -1.26 0.06
stress             6 520 2.88 1.36    3.0   1   5     4  0.09    -1.18 0.06

Who correlates with stress?

stress %>% 
  sjPlot::tab_corr(var.labels = c("Sleep Quality", "Headaches", 
                                  "Academic Performance", "Study Load", 
                                  "Extracurricular", "Stress"),
                   triangle = "lower")
  Sleep Quality Headaches Academic Performance Study Load Extracurricular Stress
Sleep Quality            
Headaches -0.059          
Academic Performance 0.254*** -0.205***        
Study Load 0.065 -0.007 0.096*      
Extracurricular 0.001 -0.073 0.153*** 0.175***    
Stress 0.165*** -0.072 0.055 0.392*** 0.052  
Computed correlation used pearson-method with listwise-deletion.
  • study_load is the obvious candidate: r = .39
  • sleep_quality is next: r = .17
  • And the two predictors barely know each other: r = .065

Model 1: study load alone

s_one <- lm(stress ~ study_load, data = stress)
summary(s_one)

Call:
lm(formula = stress ~ study_load, data = stress)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.7468 -0.9719 -0.1969  1.0844  2.8031 

Coefficients:
            Estimate Std. Error t value            Pr(>|t|)    
(Intercept)   1.8095     0.1229  14.722 <0.0000000000000002 ***
study_load    0.3875     0.0400   9.687 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.251 on 518 degrees of freedom
Multiple R-squared:  0.1534,    Adjusted R-squared:  0.1517 
F-statistic: 93.84 on 1 and 518 DF,  p-value: < 0.00000000000000022

Each 1-point increase in study load → 0.39 points more stress. \(R^2\) = .153.

Model 2: add sleep quality

s_two <- lm(stress ~ study_load + sleep_quality, data = stress)
tidy(s_two)
# A tibble: 3 × 5
  term          estimate std.error statistic  p.value
  <chr>            <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)      1.29     0.191       6.76 3.78e-11
2 study_load       0.378    0.0397      9.54 5.55e-20
3 sleep_quality    0.173    0.0495      3.50 5.00e- 4
Predictor Alone Together
study_load 0.388 0.378
sleep_quality — 0.174

The study-load slope moved by 0.01. Essentially nothing.

Why nothing moved

A coefficient only changes when the new predictor is related to both the outcome and the old predictor.

cor.test(stress$study_load, stress$sleep_quality)

    Pearson's product-moment correlation

data:  stress$study_load and stress$sleep_quality
t = 1.4933, df = 518, p-value = 0.136
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.0206320  0.1506088
sample estimates:
      cor 
0.0654704
  • r = .065 → the two predictors carry separate information
  • So each one keeps the effect it had on its own
  • Their contributions simply add

Tip

“I added a predictor and the first one didn’t budge” is a real result, and it is the common one. It means your two predictors are not competing for the same variance. Say that, instead of apologizing for it.

Interpreting Model 2 out loud

\[\widehat{Stress} = 1.29 + 0.378(StudyLoad) + 0.174(SleepQuality)\]

  • 1.29 — predicted stress for a student who rated study load and sleep quality as 0. Neither is possible (the scale starts at 1), so this intercept is a mathematical anchor, not a person.
  • 0.378 — comparing two students with the same sleep quality, one point more study load goes with 0.38 points more stress
  • 0.174 — comparing two students with the same study load, one point more sleep quality goes with 0.17 points more stress

Wait. Read that last one again. 🤨

Better sleep quality predicts more stress (b = +0.174, p < .001), holding study load constant.

Three possibilities, and you cannot choose between them from this output:

  1. The item is scored backwards — maybe 5 means “I rate my sleep quality as a problem”
  2. The relationship is real and reflects something else (students who sleep well take on more)
  3. The survey question was ambiguous and students interpreted it inconsistently

Important

This is the midterm’s lesson arriving in a lecture. Check the direction of every coefficient against what you expected. When it disagrees, the first suspect is always scoring — not the universe.

Part 3: Did It Earn Its Place?

Three ways to ask the same question

glance(s_one) %>% select(r.squared, adj.r.squared, sigma)
# A tibble: 1 × 3
  r.squared adj.r.squared sigma
      <dbl>         <dbl> <dbl>
1     0.153         0.152  1.25
glance(s_two) %>% select(r.squared, adj.r.squared, sigma)
# A tibble: 1 × 3
  r.squared adj.r.squared sigma
      <dbl>         <dbl> <dbl>
1     0.173         0.170  1.24
Model \(R^2\) adj. \(R^2\) RSE
study load .1534 .1517 1.2506
+ sleep quality .1730 .1698 1.2372
  • \(\Delta R^2\) = .0196 — about 2% more of stress accounted for
  • Adjusted \(R^2\) went up, so the predictor paid its own df
  • Typical miss shrank from 1.251 to 1.237 (on a variable with SD = 1.358)

The formal test: nested models

Model 1 is nested inside Model 2 — same data, same outcome, Model 1 is Model 2 with one coefficient forced to zero. That lets us test the difference directly.

anova(s_one, s_two)
Analysis of Variance Table

Model 1: stress ~ study_load
Model 2: stress ~ study_load + sleep_quality
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    518 810.12                                  
2    517 791.34  1    18.779 12.269 0.0005005 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Reading that table

Column What it is Here
Res.Df df left after estimating 518 → 517
RSS unexplained SS 810.12 → 791.34
Df predictors added 1
Sum of Sq SS the new predictor picked up 18.78
F is that drop bigger than noise? 12.27
Pr(>F) the p-value .0005

With one added predictor, this \(F\) is just the square of its \(t\): \(3.503^2 = 12.27\). Same test, different table. ✅

Now let’s add three useless things

s_junk <- lm(stress ~ study_load + sleep_quality +
               extracurricular + headaches + academic_perf,
             data = stress)
tidy(s_junk)
# A tibble: 6 × 5
  term            estimate std.error statistic  p.value
  <chr>              <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)       1.59      0.275      5.80  1.13e- 8
2 study_load        0.383     0.0404     9.50  8.00e-20
3 sleep_quality     0.178     0.0512     3.47  5.68e- 4
4 extracurricular  -0.0148    0.0380    -0.389 6.97e- 1
5 headaches        -0.0740    0.0445    -1.66  9.73e- 2
6 academic_perf    -0.0384    0.0547    -0.703 4.83e- 1

Extracurriculars (p = .70), headaches (p = .10), academic performance (p = .48).

Watch the two R²s diverge

Model \(R^2\) adj. \(R^2\) RSE
2 predictors .1730 .1698 1.2372
5 predictors .1779 ⬆️ .1699 ➡️ 1.2371
  • \(R^2\) rose by .005, because it always rises
  • Adjusted \(R^2\) moved by .0001 — three predictors bought you nothing
  • RSE: unchanged to four decimals

Important

Week 7’s warning, now with consequences. \(R^2\) cannot tell you whether a predictor belongs in the model. Adjusted \(R^2\), the nested \(F\), and your theory can.

Part 4: Reporting It

b or β?

stress %>%
  mutate(across(everything(), ~ as.numeric(scale(.x)))) %>%
  lm(stress ~ study_load + sleep_quality, data = .) %>%
  tidy() %>%
  select(term, estimate)
# A tibble: 3 × 2
  term          estimate
  <chr>            <dbl>
1 (Intercept)   1.26e-16
2 study_load    3.82e- 1
3 sleep_quality 1.40e- 1

Report b when

  • the units mean something (dollars, hours, mg)
  • you want to compare predictors within this sample

Report β when

  • the units are arbitrary (survey points)
  • you want to compare across studies or scales

Here: β = .38 vs .14. Study load matters roughly 2.7× as much as sleep quality for stress.

When predictors step on each other

If two predictors are nearly the same variable, the model cannot tell them apart. Coefficients get unstable and standard errors balloon.

VIF (variance inflation factor) measures it: \(VIF_j = \frac{1}{1 - R^2_j}\), where \(R^2_j\) is from regressing predictor \(j\) on all the others.

With two predictors you can compute it by hand — regress one predictor on the other:

# stress model
1 / (1 - summary(lm(study_load ~ sleep_quality, data = stress))$r.squared)
[1] 1.004305
# fire model
1 / (1 - summary(lm(trucks ~ acres, data = fire_data))$r.squared)
[1] 4.997326
  • Stress model: VIF = 1.004 — these predictors are strangers
  • Fire model: VIF = 5.00 — acres and trucks are nearly redundant (r = .89), and the SE on trucks rose from 0.515 to 0.579 as a result

Rules of thumb: > 5 is worth a sentence, > 10 is a problem. Note that the fire model still worked — high VIF widens your SEs, it does not invalidate the model.

Checking the Model - Easy

There are a lot of ways to take a look at the model

One of the easiest ways is to use performance and their check_model() function

performance::check_model(s_junk)

The write-up

Stress was regressed on self-reported study load and sleep quality in a sample of 520 students. The overall model was significant, F(2, 517) = 54.07, p < .001, accounting for 17.3% of the variance in stress (adj. R² = .170). Holding sleep quality constant, each one-point increase in study load was associated with a 0.38-point increase in stress (b = 0.38, SE = 0.04, t = 9.54, p < .001). Holding study load constant, each one-point increase in rated sleep quality was associated with a 0.17-point increase in stress (b = 0.17, SE = 0.05, t = 3.50, p < .001). Adding sleep quality improved fit over study load alone, ΔR² = .020, F(1, 517) = 12.27, p < .001. The positive coefficient for sleep quality was contrary to expectation and may reflect the direction in which that item was scored.

Tip

Everything a reader needs: the overall test, each coefficient with its holding-constant clause, the model comparison, and an honest flag on the odd result.

🛠️ Activity: Does Spending More Make More?

Tip

Week 8 Exercises · ~45 min · Pairs

1,763 films, 1970–2013: budget, international gross, year, and whether the film passes the Bechdel test.

  • Budget alone explains 35% of gross
  • Add one more variable and the budget slope gets bigger
  • You will find out why, and write it up like a consultant

Open the activity

✍️ Wrap-Up

Finish these at the bottom of your activity file:

  1. In your own words: what does “holding X constant” let you compare?
  2. A colleague says “my \(R^2\) went up, so the new predictor helped.” What do you ask them for?
  3. One thing that is still muddy → journal.

For next time

📂 Midterm Project — due Tuesday 10/27 at 11:59pm

Part 6 is now unlocked. You have everything you need: the two-predictor model, the “holding constant” interpretation, \(\Delta R^2\), and the model-comparison table.

Target for this weekend: Parts 1–5 done, Part 6 drafted. That leaves next week for writing instead of panicking.

No lab is due this week or next — the midterm is the assignment.

📖 Chapter 14.3 — ST · Chapter 8.1 — MSR

Next week: what happens when the groups from Week 6 meet the predictors from this week — and what to do when a slope is different for different groups.