Linear Regression II

Grayson White

Stat 243
Week 3 | Fall 2026

Goals for this week (and next time)

Today:

• Regression diagnostics
• Diagnostics activity

This Week:

• More linear regression:
  • Inference for regression (W)
  • categorical predictors (W)
  • multiple linear regression (W)
  • linear regression extensions (F)
• Quiz on Lab 1 (F)

(Parametric) Statistical Learning

\[ y = f(x) + \epsilon \]

Goal:

• Determine a reasonable form for \(f()\). (Ex: Line, curve, …)

• Estimate \(f()\) with \(\hat{f}()\) using the data.

• Generate predicted values: \(\hat y = \hat{f}(x)\).

Simple Linear Regression Model

\[ y = \beta_o + \beta_1 x + \epsilon \]

Consider this model when:

• Response variable \((y)\): quantitative

• Explanatory variable \((x)\): quantitative

  • Have only ONE explanatory variable.

• AND, \(f()\) can be approximated by a line.

• Need to determine the best estimates of \(\beta_o\) and \(\beta_1\).

Distinguishing between the population and the sample

\[ y = \beta_o + \beta_1 x + \epsilon \]

• Parameters:
  • Based on the population
  • Unknown then if don’t have data on the whole population
  • EX: \(\beta_o\) and \(\beta_1\)

\[ \hat{y} = \hat{ \beta}_o + \hat{\beta}_1 x \]

• Statistics:
  • Based on the sample data
  • Known
  • Usually estimate a population parameter
  • EX: \(\hat{\beta}_o\) and \(\hat{\beta}_1\)

New dataset: The Ultimate Halloween Candy Power Ranking

“The social contract of Halloween is simple: Provide adequate treats to costumed masses, or be prepared for late-night tricks from those dissatisfied with your offer. To help you avoid that type of vengeance, and to help you make good decisions at the supermarket this weekend, we wanted to figure out what Halloween candy people most prefer. So we devised an experiment: Pit dozens of fun-sized candy varietals against one another, and let the wisdom of the crowd decide which one was best.” – Walt Hickey

“While we don’t know who exactly voted, we do know this: 8,371 different IP addresses voted on about 269,000 randomly generated matchups. So, not a scientific survey or anything, but a good sample of what candy people like.”

Example: The Ultimate Halloween Candy Power Ranking

Example: The Ultimate Halloween Candy Power Ranking

candy <- read_csv("https://raw.githubusercontent.com/fivethirtyeight/data/master/candy-power-ranking/candy-data.csv") %>%
  mutate(pricepercent = pricepercent*100)

glimpse(candy)

Rows: 85
Columns: 13
$ competitorname   <chr> "100 Grand", "3 Musketeers", "One dime", "One quarter…
$ chocolate        <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,…
$ fruity           <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1,…
$ caramel          <dbl> 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,…
$ peanutyalmondy   <dbl> 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ nougat           <dbl> 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,…
$ crispedricewafer <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ hard             <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1,…
$ bar              <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,…
$ pluribus         <dbl> 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1,…
$ sugarpercent     <dbl> 0.732, 0.604, 0.011, 0.011, 0.906, 0.465, 0.604, 0.31…
$ pricepercent     <dbl> 86.0, 51.1, 11.6, 51.1, 51.1, 76.7, 76.7, 51.1, 32.5,…
$ winpercent       <dbl> 66.97173, 67.60294, 32.26109, 46.11650, 52.34146, 50.…

New Modeling Goal

New modeling goal: predict win percentage by using the price percentage variable.

ggplot(data = candy, 
       mapping = aes(x = pricepercent,
                     y = winpercent)) +
  geom_point(alpha = 0.6, size = 4, 
             color = "chocolate4") +
  geom_smooth(method = "lm", se = FALSE,
              color = "deeppink2")

Constructing the Simple Linear Regression Model in R

We can use the lm() function to construct the simple linear regression model in R.

mod <- lm(winpercent ~ pricepercent, data = candy)

library(moderndive)
get_regression_table(mod)

# A tibble: 2 × 7
  term         estimate std_error statistic p_value lower_ci upper_ci
  <chr>           <dbl>     <dbl>     <dbl>   <dbl>    <dbl>    <dbl>
1 intercept      41.979     2.908    14.435   0       36.195   47.763
2 pricepercent    0.178     0.053     3.352   0.001    0.072    0.283

What is the fitted model form?

\[\begin{align*} \hat{y} &= \hat{\beta_o} + \hat{\beta_1} \times x_{pricepercent} \\ &= 41.979 + 0.178 \times x_{pricepercent} \end{align*}\]

How do we interpret the coefficients?

Prediction

new_cases <- data.frame(pricepercent = c(25, 85, 150))
predict(mod, newdata = new_cases)

       1        2        3 
46.42443 57.09409 68.65289 

• We didn’t have any treats in our sample with a price percentage of 85%. Can we still make this prediction?
  • Called interpolation

• We didn’t have any treats in our sample with a price percentage of 150%. Can we still make this prediction?
  • Called extrapolation

Cautions

• Careful to only predict values within the range of \(x\) values in the sample.

• Make sure to investigate outliers: observations that fall far from the cloud of points.

  • High leverage points

A closer look at our model

get_regression_table(mod)

# A tibble: 2 × 7
  term         estimate std_error statistic p_value lower_ci upper_ci
  <chr>           <dbl>     <dbl>     <dbl>   <dbl>    <dbl>    <dbl>
1 intercept      41.979     2.908    14.435   0       36.195   47.763
2 pricepercent    0.178     0.053     3.352   0.001    0.072    0.283

\[\begin{align*} \hat{y} &= \hat{\beta_o} + \hat{\beta_1} \times x_{pricepercent} \\ &= 41.979 + 0.178 \times x_{pricepercent} \end{align*}\]

What assumptions have we made?

Linear Regression Assumptions

We can always find the line of best fit to explore data, but…

To make accurate predictions or inferences, certain conditions should be met.

\[\begin{align*} y &= \beta_o + \beta_1 + x + \epsilon \\ \epsilon &\overset{\mbox{ind}}{\sim} \mathcal{N}(0,~ \sigma^2) \end{align*}\]

• Mathematically, all linear regression assumptions are encoded in the above equation.

• However, we can break them down into four easy to remember mnemonics.

Linear Regression Assumptions

\[\begin{align*} y &= \beta_o + \beta_1 + x + \epsilon \\ \epsilon &\overset{\mbox{ind}}{\sim} \mathcal{N}(0,~ \sigma^2) \end{align*}\]

1. Linearity: The relationship between explanatory and response variables must be approximately linear

   • Mathematically: \(y = \beta_0 + \beta_1 x + \epsilon\) where \(\mathbb{E}(\epsilon \ | \ x) = 0\)

   • Check using scatterplot of data, or residual plot

Linear Regression Assumptions

\[\begin{align*} y &= \beta_o + \beta_1 + x + \epsilon \\ \epsilon &\overset{\mbox{ind}}{\sim} \mathcal{N}(0,~ \sigma^2) \end{align*}\]

2. Independence: The observations should be independent of one another.

   • Mathematically: \(\epsilon_i\) and \(\epsilon_j\) are independent (given \(x\))
   • Can check by considering data context, and
   • by looking at residual scatterplots too

Linear Regression Assumptions

\[\begin{align*} y &= \beta_o + \beta_1 + x + \epsilon \\ \epsilon &\overset{\mbox{ind}}{\sim} \mathcal{N}(0,~ \sigma^2) \end{align*}\]

3. Normality: The distribution of residuals should be approximately bell-shaped, unimodal, symmetric, and centered at 0 at every “slice” of the explanatory variable

   • Mathematically: \(\epsilon \ | \ x \sim N(0, \sigma^2)\)
   • Simple check: look at histogram of residuals
   • Better to use a “Q-Q plot”

Linear Regression Assumptions

\[\begin{align*} y &= \beta_o + \beta_1 + x + \epsilon \\ \epsilon &\overset{\mbox{ind}}{\sim} \mathcal{N}(0,~ \sigma^2) \end{align*}\]

4. Equal Variability: Variance of residuals should be roughly constant across data set. Also called “homoscedasticity”. Models that violate this assumption are sometimes called “heteroscedastic”

   • Mathematically: \(\mbox{Var}(\epsilon \ | \ x) = \sigma^2\) for constant \(\sigma^2\)
   • Check using residual plot.

A cute way to remember this: “LINE”

1. Linearity: The relationship between explanatory and response variables must be approximately linear

2. Independence: The observations should be independent of one another.

3. Normality: The distribution of residuals should be approximately bell-shaped, unimodal, symmetric, and centered at 0 at every “slice” of the explanatory variable

4. Equal Variability: Variance of residuals should be roughly constant across data set. Also called “homoscedasticity”. Models that violate this assumption are sometimes called “heteroscedastic”

Assessing conditions: diagnostics plots

In order to assess if we’ve met the described conditions, we will utilize four common diagnostics plots:

• Residual vs. fitted plot (i.e. ‘residual plot’)
  • for assessing linearity, independence, and equal variability
• Normal Q-Q plot
  • for assessing normality
• Residual histogram
  • for assessing normality
• Residual vs. leverage plot
  • for assessing the influence of outliers

In order to make these plots, we’ll use the gglm (grammer of graphics for linear model diagnostics) package.

# Run the below code (without the '#') once
# install.packages("gglm") 
# Then load with:
library(gglm)

Recall out simple linear regression model

mod <- lm(winpercent ~ pricepercent, data = candy)

Let’s check if this meets the LINE assumptions

Residual vs. fitted plot

ggplot(mod) + 
  stat_fitted_resid()

• Linearity: ✅

• Independence: ✅

• Equal Variability: ✅

Checking normality: residual histogram

ggplot(mod) + 
  stat_resid_hist(bins = 15)

• Normality: 🤔
  • Looks pretty good!
  • We’ll look at a Q-Q plot to further assess!

Checking normality: Q-Q plot

ggplot(mod) + 
  stat_normal_qq()

• Normality: ✅

Checking outliers: Residuals vs. leverage plot

ggplot(mod) + 
  stat_resid_leverage()

• Looks pretty good! 🤩

Diagnostic shortcut: gglm::gglm()

gglm(mod)

Now let’s simulate some data…

A valid model

set.seed(794)
n <- 80
x <- runif(n)
beta_1 <- 3.5
sigma <- .6
y <- beta_1 * x + rnorm(n, sd = sigma)
df <- data.frame(x, y)
m1 <- lm(y ~ x, data = df)

A valid model

ggplot(df, aes(x = x, y = y)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE)

A valid model

gglm(m1)

What can break?

1. Unbiased estimation of \(\hat\beta\)

2. Accurate assessment of \(\mbox{SE}(\hat\beta)\)

   • p-values and hypothesis tests
   • CIs

Checking for breaks

it <- 5000
beta_hats <- rep(NA, it)
capture <- rep(FALSE, it)
for(i in 1:it) {
  y <- beta_1 * x + rnorm(n, sd = sigma)
  m <- lm(y ~ x)
  beta_hats[i] <- m$coef[2]
  ci <- confint(m)[2, ]
  capture[i] <- (ci[1] < beta_1 & beta_1 < ci[2])
}

Checking for breaks

ggplot(data.frame(beta_hat = beta_hats), aes(x = beta_hat)) + 
  geom_histogram() + 
  geom_vline(xintercept = beta_1, color = "tomato")

• Bias in \(\hat\beta\)’s?

mean(beta_hats) - beta_1

[1] 0.000283525

• Accurate assessment of \(\mbox{SE}(\hat\beta)\)

mean(capture)

[1] 0.9482

Let’s break something

Non-constant variance

set.seed(794)
n <- 80
x <- runif(n)
y <- beta_1 * x + rnorm(n, sd = .93 * x)
df <- data.frame(x, y)
m1 <- lm(y ~ x, data = df)

Non-constant variance

ggplot(df, aes(x = x, y = y)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE)

Non-constant variance

gglm(m1)

Checking for breaks

it <- 5000
beta_hats <- rep(NA, it)
capture <- rep(FALSE, it)
for(i in 1:it) {
  y <- beta_1 * x + rnorm(n, sd = .93 * x)
  m <- lm(y ~ x)
  beta_hats[i] <- m$coef[2]
  ci <- confint(m)[2, ]
  capture[i] <- (ci[1] < beta_1 & beta_1 < ci[2])
}

Checking for breaks

ggplot(data.frame(beta_hat = beta_hats), 
       aes(x = beta_hat)) + 
  geom_histogram() + 
  geom_vline(xintercept = beta_1, color = "tomato")

• Bias in \(\hat\beta\)’s?

mean(beta_hats) - beta_1

[1] -0.004152782

• Accurate assessment of \(\mbox{SE}(\hat\beta)\)

mean(capture)

[1] 0.9106

Activity

With a partner (there needs to be one computer between you), choose one of the following components of the linear model to break, then assess what goes wrong.

1. Normality of errors
2. Independence of errors
3. Linear trend in mean function

You can build off the code used in these slides.

Next time

• Inference for regression
• categorical predictors
• multiple linear regression