Simple Linear Regression - Model Evaluation

Notes and in-class exercises

You can download the .qmd file for this activity here and open in R-studio. The rendered version is posted in the course website (Activities tab). I often experiment with the class activities (and see it in live!) and make updates, but I always post the final version before class starts. To be sure you have the most up-to-date copy, please download it once you’ve settled in before class begins.

Notes

Learning goals

By the end of this lesson, you should be able to:

• Use residual plots to evaluate the correctness of a model
• Explain the rationale for the R-squared metric of model strength
• Interpret the R-squared metric
• Think about ethical implications of modeling by examining the impacts of biased data, power dynamics, the role of categorization, and the role of emotion and lived experience

Readings and videos

Choose either the reading or the videos to go through before class.

• Reading: Sections 1.7, 3.7, and 3.8 in the STAT 155 Notes
  • Note: You do not need to focus on the “Ladder of Power” in Section 3.8. Transformations in general will be the focus of the next activity we do.
• Videos:
  • Model evaluation: is the model wrong? (slides)
  • Model evaluation: is the model strong? (slides)
  • Model evaluation: is the model fair? (slides)
  • R Code for Evaluating and Using a Linear Model

Model Assumptions

One way to think about model evaluation is to consider whether or not underlying assumptions of our regression models are being met (or not). Asking ourselves if our models are “wrong”, “strong”, and “fair” approaches this from one perspective. To the first question (whether our model is wrong), recall the following four assumptions of linear regression:

1. Linearity
2. Independence
3. Normality
4. Equal Variance

Note that they spell “LINE” (how convenient!).

By assumptions, we mean that the above four “things” are needed mathematically in order for linear regression to “work”.

Whereas we can check some of these assumptions using a residual plot, we need to examine the context of our data collection when checking the Independence assumption. What we mean by independence, is that the residuals in our model do not depend on one another. This may seem like an unsatisfying definition, so here are some examples:

• Suppose I want to understand the association between a person’s high school GPA and their college GPA. I collect data from every graduating senior, at three different high schools. If I have college GPA as my outcome, and high school GPA as my predictor, are my residuals independent? Probably not! It is reasonable to believe that students from the same high school may have similar GPAs, due to resources their high school may have had available, or specific teachers grading differently at one school or another. This is an example of clustering, where we have clusters of students within schools. The independence assumption of our linear regression model would be violated. One way to address this would be to include which high school they went to as an additional covariate in our regression model (we’ll get to this with multiple linear regression), and more advanced methods are covered in a course on Correlated Data.

• Suppose I want to understand the association between a mouse’s weight and their water consumption across time. I collect data for 365 days for ten different mice, recording their weight and water consumption each day of the year. If I have weight as my predictor and water consumption as my outcome, are my residuals independent? Nope! This is an example of correlated data that is longitudinal in nature: I have multiple observations per individual (mouse) across time. A mouse’s weight one day is certainly not independent of it’s weight the following day. The independence assumption of our linear regression model would again be violated. One way to address this would be to include “Mouse ID” as a predictor in our regression model (again, we’ll get to this with multiple linear regression).

All types of data that will violate the independence assumption of linear regression will have some sort of correlation structure (within individual, across time, across space, etc.). Think about clusters. If your observations fall neatly into specific clusters, your data may violate the independence assumption of linear regression.

File organization: Save this file in the “Activities” subfolder of your “STAT155” folder.

Exercises

Exercise 1: Is the model correct?

Let’s revisit the Capital Bikeshare data:

# Load packages and import data
library(readr)
library(ggplot2)
library(dplyr)

bikes <- read_csv("https://mac-stat.github.io/data/bikeshare.csv")

We previously explored a model of daily ridership among registered users as a function of temperature:

# Fit a linear model
bike_model <- lm(riders_registered ~ temp_feel, data = bikes)

# Check it out
summary(bike_model)
## 
## Call:
## lm(formula = riders_registered ~ temp_feel, data = bikes)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3607.1  -959.2  -153.8   998.2  3304.8 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -667.916    251.608  -2.655  0.00811 ** 
## temp_feel     57.892      3.306  17.514  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1310 on 729 degrees of freedom
## Multiple R-squared:  0.2961, Adjusted R-squared:  0.2952 
## F-statistic: 306.7 on 1 and 729 DF,  p-value: < 2.2e-16

Plot this relationship with both a curved and linear trend line. Based on this plot, do you think the model is correct? If not, which of the LINE assumptions does it violate?

# Plot temp_feel vs riders_registered with a model trend
ggplot(bikes, aes(x = temp_feel, y = riders_registered)) + 
    geom_point() + 
    geom_smooth(method = "lm", se = FALSE) +
    geom_smooth(se = FALSE, color = "red")

Exercise 2: Residual plots

Plotting the residuals vs the predictions (also called “fitted values”) for each case can help us assess how wrong our model is. This will be a particularly important tool when evaluating models with multiple predictors. Construct the residual plot for bike_model. As with the scatterplot, this plot indicates that bike_model violates one of the LINE assumptions. Explain which assumption that is and how you can tell that from just the residual plot.

Notes:

• Information about the residuals (.resid) and predictions (.fitted) are stored within our model, thus we start our ggplot() with the model name as opposed to the raw dataset. We will rarely start ggplot() with a model instead of the data.
• We can fix this model by adding a quadratic “transformation term”.

# Check out the residual plot for bike_model
ggplot(bike_model, aes(x = .fitted, y = .resid)) + 
    geom_point() + 
    geom_hline(yintercept = 0) + # Check this first!
    geom_smooth(se = FALSE)

Exercise 3: What’s incorrect about this model?

Consider another example. The mammals data includes data on the average brain weight (g) and body weight (kg) for a variety of mammals:

# Import the data
mammals <- read_csv("https://mac-stat.github.io/data/mammals.csv")

# Check it out
head(mammals)
## # A tibble: 6 × 4
##    ...1 animal            body brain
##   <dbl> <chr>            <dbl> <dbl>
## 1     1 Arctic fox        3.38  44.5
## 2     2 Owl monkey        0.48  15.5
## 3     3 Mountain beaver   1.35   8.1
## 4     4 Cow             465    423  
## 5     5 Grey wolf        36.3  120. 
## 6     6 Goat             27.7  115

Fit a model of brain vs body weight:

# Construct the model
mammal_model <- lm(brain ~ body, mammals)

# Check it out
summary(mammal_model)
## 
## Call:
## lm(formula = brain ~ body, data = mammals)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -810.07  -88.52  -79.64  -13.02 2050.33 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 91.00440   43.55258    2.09   0.0409 *  
## body         0.96650    0.04766   20.28   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 334.7 on 60 degrees of freedom
## Multiple R-squared:  0.8727, Adjusted R-squared:  0.8705 
## F-statistic: 411.2 on 1 and 60 DF,  p-value: < 2.2e-16

1. Construct two plots that will help us evaluate mammal_model:

# Scatterplot of brain weight (y) vs body weight (x)
# Include a model trend line (i.e. a representation of mammal_model)
ggplot(mammals, aes(x = body, y = brain)) + 
    geom_point() + 
    geom_smooth(method = "lm", se = FALSE) +
    geom_smooth(se = FALSE, color = "red")

# Residual plot for mammal_model

2. These two plots confirm that our model is wrong. What is wrong? That is, which of the LINE assumptions are violated? (NOTE: We again can fix this model by “transforming” one or both of the brain and body variables.

Exercise 4: Exploring mammals

Just for fun, let’s dig into the mammals data. Discuss what you observe:

# Label the points by the animal name!
# Discuss: What 2 things are new in this code?
ggplot(mammals, aes(x = body, y = brain, label = animal)) + 
    geom_text() + 
    geom_smooth(method = "lm", se = FALSE) 

# Zoom in
ggplot(mammals, aes(x = body, y = brain, label = animal)) + 
    geom_text() + 
    lims(y = c(0, 1500), x = c(0, 600))

# Zoom in more
ggplot(mammals, aes(x = body, y = brain, label = animal)) + 
    geom_text() + 
    lims(y = c(0, 500), x = c(0, 200))

Exercise 5: Is the model strong? Developing R-squared intuition

The R-squared metric is a way to quantify the strength of a model. It measures how much variation in the outcome/response variable can be explained by the variation in the predictors.

Where does R-squared come from? Well, it turns out that we can partition the variance of the observed response values into the variability that’s explained by the model (the variance of the predictions) and the variability that’s left unexplained by the model (the variance of the residuals):

\[\text{Var(observed) = Var(predicted) + Var(residuals)}\]

Strong models have residuals that don’t deviate far from 0. So the smaller the variance in the residuals (thus larger the variance in the predictions), the stronger the model. Take a look at the picture below and write a few sentences addressing the following:

• The two rows of plots show a stronger and a weaker model. Just by looking at the blue trend line and the dispersion of the points about the line, which row corresponds to the stronger model? How can you tell? Which row would you expect to have a higher correlation?

• What is different about the variance of the residuals from the first to the second row?

• The first row corresponds to the weaker model. We can tell because the points are much more dispersed from the trend line than in the second row. Recall that the correlation metric measures how closely clustered points are about a straight line of best fit, so we would expect the correlation to be lower for the first row than the second row.
• The variance of the residuals is much lower for the second row—the residuals are all quite small. This indicates a stronger model.

Putting this together, the R-squared compares Var(predicted) to Var(response):

\[R^2 = \frac{\text{variance of predicted values}}{\text{variance of observed response values}} = 1 - \frac{\text{variance of residuals}}{\text{variance of observed response values}}\]

R-squared

\[ R^2 = 1 - \frac{SSE}{SSTO} = 1 - \frac{\sum (y_i - \hat{y}_i)^2}{\sum (y_i - \bar{y})^2} \] where \(y_i\) are our observed outcomes, \(i = 1, \dots, n\), \(\hat{y}_i\) are our fitted values/predictions, and \(\bar{y}\) is our observed average outcome.

Exercise 6: R-squared Interpretations

Recall bikemod1 from Exercise 1, where we predicted registered riders by what the temperature felt like on a given day. Use the summary function to look out the model output for bikemod1, and interpret the \(R^2\) value for this model, in the context of the problem. (NOTE: \(R^2\) is reported in output here as “Multiple R-squared”).

# Get R-squared
summary(bike_model)
## 
## Call:
## lm(formula = riders_registered ~ temp_feel, data = bikes)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3607.1  -959.2  -153.8   998.2  3304.8 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -667.916    251.608  -2.655  0.00811 ** 
## temp_feel     57.892      3.306  17.514  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1310 on 729 degrees of freedom
## Multiple R-squared:  0.2961, Adjusted R-squared:  0.2952 
## F-statistic: 306.7 on 1 and 729 DF,  p-value: < 2.2e-16

Multiple R-squared: 0.2961

Interpretation: 29.61% of the variation in number of registered riders on any given day can be explained by the variation in temperature (specifically, what temperature it “feels” like it is).

Exercise 7: Further exploring R-squared

In this exercise, we’ll look at data from a synthetic dataset called Anscombe’s quartet. Load the data in as follows, and look at the first few rows:

data(anscombe)

# Look at the first few rows
head(anscombe)
##   x1 x2 x3 x4   y1   y2    y3   y4
## 1 10 10 10  8 8.04 9.14  7.46 6.58
## 2  8  8  8  8 6.95 8.14  6.77 5.76
## 3 13 13 13  8 7.58 8.74 12.74 7.71
## 4  9  9  9  8 8.81 8.77  7.11 8.84
## 5 11 11 11  8 8.33 9.26  7.81 8.47
## 6 14 14 14  8 9.96 8.10  8.84 7.04

The anscombe data is actually 4 datasets in one: x1 and y1 go together, and so forth. Examine the coefficient estimates (in the “Estimate” column of the “Coefficients:” part) and the “Multiple R-squared” value on the second to last line. What do you notice? How do these models compare?

anscombe_mod1 <- lm(y1 ~ x1, data = anscombe)
anscombe_mod2 <- lm(y2 ~ x2, data = anscombe)
anscombe_mod3 <- lm(y3 ~ x3, data = anscombe)
anscombe_mod4 <- lm(y4 ~ x4, data = anscombe)

summary(anscombe_mod1)
## 
## Call:
## lm(formula = y1 ~ x1, data = anscombe)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.92127 -0.45577 -0.04136  0.70941  1.83882 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.0001     1.1247   2.667  0.02573 * 
## x1            0.5001     0.1179   4.241  0.00217 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.237 on 9 degrees of freedom
## Multiple R-squared:  0.6665, Adjusted R-squared:  0.6295 
## F-statistic: 17.99 on 1 and 9 DF,  p-value: 0.00217
summary(anscombe_mod2)
## 
## Call:
## lm(formula = y2 ~ x2, data = anscombe)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9009 -0.7609  0.1291  0.9491  1.2691 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    3.001      1.125   2.667  0.02576 * 
## x2             0.500      0.118   4.239  0.00218 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.237 on 9 degrees of freedom
## Multiple R-squared:  0.6662, Adjusted R-squared:  0.6292 
## F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002179
summary(anscombe_mod3)
## 
## Call:
## lm(formula = y3 ~ x3, data = anscombe)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1586 -0.6146 -0.2303  0.1540  3.2411 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.0025     1.1245   2.670  0.02562 * 
## x3            0.4997     0.1179   4.239  0.00218 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.236 on 9 degrees of freedom
## Multiple R-squared:  0.6663, Adjusted R-squared:  0.6292 
## F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002176
summary(anscombe_mod4)
## 
## Call:
## lm(formula = y4 ~ x4, data = anscombe)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.751 -0.831  0.000  0.809  1.839 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.0017     1.1239   2.671  0.02559 * 
## x4            0.4999     0.1178   4.243  0.00216 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.236 on 9 degrees of freedom
## Multiple R-squared:  0.6667, Adjusted R-squared:  0.6297 
## F-statistic:    18 on 1 and 9 DF,  p-value: 0.002165

• All of these models have close to the same intercept, slope, and R-squared!

Now take a look at the following scatterplots of the 4 pairs of variables. What do you notice? What takeaway can we draw from this exercise?

ggplot(anscombe, aes(x = x1, y = y1)) +
    geom_point() + 
    geom_smooth(method = "lm", color = "red", se = FALSE)

ggplot(anscombe, aes(x = x2, y = y2)) +
    geom_point() + 
    geom_smooth(method = "lm", color = "red", se = FALSE)

ggplot(anscombe, aes(x = x3, y = y3)) +
    geom_point() + 
    geom_smooth(method = "lm", color = "red", se = FALSE)

ggplot(anscombe, aes(x = x4, y = y4)) +
    geom_point() + 
    geom_smooth(method = "lm", color = "red", se = FALSE)

• x2 and y2: The scatterplot is clearly curved—a “linear” regression model with squared terms, for example, would be more appropriate for this data. (We’ll talk more about ways to handle nonlinear relationships soon!)
• x3 and y3: There is a very clear outlier at about x3 = 13 that we would want to dig into to better understand the context. After that investigation, we might consider removing this outlier and refitting the model.
• x4 and y4: There is clearly something strange going on with most of the cases having an x4 value of exactly 8. We would not want to jump straight into modeling. Instead, we should dig deeper to find out more about this data.

Complete Exercise 8-9 after the class.

Exercise 8: Biased data, biased results: example 1

In the above exercises, we focused on exploring our first 2 model evaluation questions: Is it correct? Is it strong? In the next exercises, let’s explore the third question: Is it fair?

Data are not neutral. Data can reflect personal biases, institutional biases, power dynamics, societal biases, the limits of our knowledge, etc. And biased data can lead to biased analyses. Consider the example of a large company that developed a model / algorithm to review the résumés of applicants for software developer & other tech positions. The model then gave each applicant a score indicating their hire-ability or potential for success at the company. You can think of this model as something like:

\[E[\text{potential | résumé features}] = \beta_0 + \beta_1 (\text{résumé features})\]

Skim this Reuter’s article about the company’s résumé model.

• Explain why the data used by this model are not neutral.

• What are the potential implications, personal or societal, of the results produced from this biased data?

Exercise 9: Biased data, biased results: example 2

When working with categorical variables, we’ve seen that our units of observation fall into neat groups. Reality isn’t so discrete. For example, check out questions 6 and 9 on page 2 of the 2020 US Census. With your group, discuss the following:

• What are a couple of issues you see with these questions?

• What impact might this type of data collection have on a subsequent analysis of the census responses and the policies it might inform?

• Can you think of a better way to write these questions while still preserving the privacy of respondents?

FOR A DEEPER DISCUSSION: Read Chapter 4 of Data Feminism on “What gets counted counts”.

Additional practice

Exercise 10: Data & visualization drill

1. Let’s practice some data wrangling with our peaks data. In addition to some basic functions (e.g. head()), use the tidyverse functions we’ve been accumulating: select(), summarize(), filter().

# How many hikes are in the data set?


# Show just the peak's name and its elevation for the first 6 hikes
___ %>% 
  ___ %>% 
  head()

# Calculate the average hiking time


# Show data for the hikes that are *at least* 17 miles long


# Calculate the average hiking time for hikes that are at least 17 miles long
# HINT: You'll need to use 2 tidyverse functions!
___ %>% 
  ___ %>% 
  ___

# Calculate the average hiking time for hikes with an "easy" rating
# HINT: You'll need to use 2 tidyverse functions again!

## Error in parse(text = input): <text>:5:2: unexpected input
## 4: # Show just the peak's name and its elevation for the first 6 hikes
## 5: __
##     ^

2. Let’s practice some visualization!

# Construct a visualization of hike rating


# Construct a visualization of hike length


# Construct a visualization of the relationship between hiking time and length (distance)

3. Summarize, in words, what you learned from each plot above.

Exercise 10: Types of models

1. Is it possible to build a model that’s correct but weak? If not, explain. If yes, sketch an example of what this might look like.

2. Is it possible to build a model that’s wrong but strong? If not, explain. If yes, sketch an example of what this might look like.

Exercise 11: Model practice

Let’s explore the temperature data for Hobart, including the temperature in celsius at 3pm and 9am:

temps <- read_csv("https://mac-stat.github.io/data/weather_3_locations.csv") %>% 
  filter(location == "Hobart") %>% 
  select(date, temp3pm, temp9am)

head(temps)
## # A tibble: 6 × 3
##   date       temp3pm temp9am
##   <date>       <dbl>   <dbl>
## 1 2020-01-01    22.1    17.4
## 2 2020-01-02    19      19  
## 3 2020-01-03    19.7    16.2
## 4 2020-01-04    20.6    19.9
## 5 2020-01-05    21.9    14.2
## 6 2020-01-06    20.8    15.3

Part a

Between temp3pm and temp9am, which is the response variable and which is the predictor?

Part b

Visualize the relationship between these 2 variables, including representations of the raw data and the simple linear regression model. Describe your observations.

Part c

Build the linear regression model. Write out the estimated model formula and interpret the 2 coefficients.

Solutions

Exercise 1: Is the model correct?

The red curved trend line shows a clear downward trend around 85 degrees, which contextually makes plenty of sense—extremely hot days would naturally see less riders. Overall the combination of the upward trend and downward trend makes for a curved relationship that is not captured well by a straight line of best fit. Specifically, a simple linear regression model would violate the Linearity assumption.

# Load packages and import data
library(readr)
library(ggplot2)
library(dplyr)

bikes <- read_csv("https://mac-stat.github.io/data/bikeshare.csv")

ggplot(bikes, aes(x = temp_feel, y = riders_registered)) + 
    geom_point() + 
    geom_smooth(se = FALSE, color = "red") +
    geom_smooth(method = "lm", se = FALSE)

Exercise 2: Residual plots

The residual plot shows a lingering trend in the residuals—the blue curve traces the trend in the residuals, and it does not lie flat on the y = 0 line. This again suggests that the Linearity assumption is violated.

bike_model <- lm(riders_registered ~ temp_feel, data = bikes)

# Check out the residual plot for bike_model
ggplot(bike_model, aes(x = .fitted, y = .resid)) + 
    geom_point() + 
    geom_hline(yintercept = 0) +
    geom_smooth(se = FALSE)

Exercise 3: What’s incorrect about this model?

# Import the data
mammals <- read_csv("https://mac-stat.github.io/data/mammals.csv")

# Check it out
head(mammals)
## # A tibble: 6 × 4
##    ...1 animal            body brain
##   <dbl> <chr>            <dbl> <dbl>
## 1     1 Arctic fox        3.38  44.5
## 2     2 Owl monkey        0.48  15.5
## 3     3 Mountain beaver   1.35   8.1
## 4     4 Cow             465    423  
## 5     5 Grey wolf        36.3  120. 
## 6     6 Goat             27.7  115

# Construct the model
mammal_model <- lm(brain ~ body, mammals)

# Check it out
summary(mammal_model)
## 
## Call:
## lm(formula = brain ~ body, data = mammals)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -810.07  -88.52  -79.64  -13.02 2050.33 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 91.00440   43.55258    2.09   0.0409 *  
## body         0.96650    0.04766   20.28   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 334.7 on 60 degrees of freedom
## Multiple R-squared:  0.8727, Adjusted R-squared:  0.8705 
## F-statistic: 411.2 on 1 and 60 DF,  p-value: < 2.2e-16

# Scatterplot of brain weight (y) vs body weight (x)
# Include a model trend line (i.e. a representation of mammal_model)
ggplot(mammals, aes(y = brain, x = body)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE)

# Residual plot for mammal_model
ggplot(mammal_model, aes(x = .fitted, y = .resid)) + 
    geom_point() + 
    geom_hline(yintercept = 0) +
    geom_smooth(se = FALSE)

2. The biggest issue here is that the assumption of equal variance is violated. There’s much greater variability in the residuals as the predictions increase. This is because there’s much greater variability in the brain weights (y) as body weights (x) increase.

Exercise 4: Exploring mammals

Answers will vary.

Exercise 5: Is the model strong? Developing R-squared intuition

The R-squared metric is a way to quantify the strength of a model. It measures how much variation in the outcome/response variable can be explained by the model.

Where does R-squared come from? Well, it turns out that we can partition the variance of the observed response values into the variability that’s explained by the model (the variance of the predictions) and the variability that’s left unexplained by the model (the variance of the residuals):

\[\text{Var(observed) = Var(predicted) + Var(residuals)}\]

“Good” models have residuals that don’t deviate far from 0. So the smaller the variance in the residuals (thus larger the variance in the predictions), the stronger the model. Take a look at the picture below and write a few sentences addressing the following:

• The first row corresponds to the weaker model. We can tell because the points are much more dispersed from the trend line than in the second row. Recall that the correlation metric measures how closely clustered points are about a straight line of best fit, so we would expect the correlation to be lower for the first row than the second row.
• The variance of the residuals is much lower for the second row—the residuals are all quite small. This indicates a stronger model.

Exercise 6: R-squared Interpretations

summary(bike_model)
## 
## Call:
## lm(formula = riders_registered ~ temp_feel, data = bikes)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3607.1  -959.2  -153.8   998.2  3304.8 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -667.916    251.608  -2.655  0.00811 ** 
## temp_feel     57.892      3.306  17.514  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1310 on 729 degrees of freedom
## Multiple R-squared:  0.2961, Adjusted R-squared:  0.2952 
## F-statistic: 306.7 on 1 and 729 DF,  p-value: < 2.2e-16

Multiple R-squared: 0.2961

Interpretation: 29.61% of the variation in number of registered riders on any given day can be explained by the variation in temperature (specifically, what temperature it “feels” like it is).

Exercise 7: Further exploring R-squared

In this exercise, we’ll look at data from a synthetic dataset called Anscombe’s quartet. Load the data in as follows, and look at the first few rows:

data(anscombe)

# Look at the first few rows
head(anscombe)
##   x1 x2 x3 x4   y1   y2    y3   y4
## 1 10 10 10  8 8.04 9.14  7.46 6.58
## 2  8  8  8  8 6.95 8.14  6.77 5.76
## 3 13 13 13  8 7.58 8.74 12.74 7.71
## 4  9  9  9  8 8.81 8.77  7.11 8.84
## 5 11 11 11  8 8.33 9.26  7.81 8.47
## 6 14 14 14  8 9.96 8.10  8.84 7.04

All of these models have close to the same intercept, slope, and R-squared!

anscombe_mod1 <- lm(y1 ~ x1, data = anscombe)
anscombe_mod2 <- lm(y2 ~ x2, data = anscombe)
anscombe_mod3 <- lm(y3 ~ x3, data = anscombe)
anscombe_mod4 <- lm(y4 ~ x4, data = anscombe)

summary(anscombe_mod1)
## 
## Call:
## lm(formula = y1 ~ x1, data = anscombe)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.92127 -0.45577 -0.04136  0.70941  1.83882 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.0001     1.1247   2.667  0.02573 * 
## x1            0.5001     0.1179   4.241  0.00217 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.237 on 9 degrees of freedom
## Multiple R-squared:  0.6665, Adjusted R-squared:  0.6295 
## F-statistic: 17.99 on 1 and 9 DF,  p-value: 0.00217
summary(anscombe_mod2)
## 
## Call:
## lm(formula = y2 ~ x2, data = anscombe)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9009 -0.7609  0.1291  0.9491  1.2691 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    3.001      1.125   2.667  0.02576 * 
## x2             0.500      0.118   4.239  0.00218 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.237 on 9 degrees of freedom
## Multiple R-squared:  0.6662, Adjusted R-squared:  0.6292 
## F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002179
summary(anscombe_mod3)
## 
## Call:
## lm(formula = y3 ~ x3, data = anscombe)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1586 -0.6146 -0.2303  0.1540  3.2411 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.0025     1.1245   2.670  0.02562 * 
## x3            0.4997     0.1179   4.239  0.00218 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.236 on 9 degrees of freedom
## Multiple R-squared:  0.6663, Adjusted R-squared:  0.6292 
## F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002176
summary(anscombe_mod4)
## 
## Call:
## lm(formula = y4 ~ x4, data = anscombe)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.751 -0.831  0.000  0.809  1.839 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.0017     1.1239   2.671  0.02559 * 
## x4            0.4999     0.1178   4.243  0.00216 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.236 on 9 degrees of freedom
## Multiple R-squared:  0.6667, Adjusted R-squared:  0.6297 
## F-statistic:    18 on 1 and 9 DF,  p-value: 0.002165

But when we look at the scatterplots, they all look substantially different, and we would want to approach our modeling differently for each one:

• x1 and y1: A linear model seems appropriate for this data.
• x2 and y2: The scatterplot is clearly curved—a “linear” regression model with squared terms, for example, would be more appropriate for this data. (We’ll talk more about ways to handle nonlinear relationships soon!)
• x3 and y3: There is a very clear outlier at about x3 = 13 that we would want to dig into to better understand the context. After that investigation, we might consider removing this outlier and refitting the model.
• x4 and y4: There is clearly something strange going on with most of the cases having an x4 value of exactly 8. We would not want to jump straight into modeling. Instead, we should dig deeper to find out more about this data.

ggplot(anscombe, aes(x = x1, y = y1)) +
    geom_point() + 
    geom_smooth(method = "lm", color = "red", se = FALSE)

ggplot(anscombe, aes(x = x2, y = y2)) +
    geom_point() + 
    geom_smooth(method = "lm", color = "red", se = FALSE)

ggplot(anscombe, aes(x = x3, y = y3)) +
    geom_point() + 
    geom_smooth(method = "lm", color = "red", se = FALSE)

ggplot(anscombe, aes(x = x4, y = y4)) +
    geom_point() + 
    geom_smooth(method = "lm", color = "red", se = FALSE)

Exercises 8 - 9

No solutions for these exercises. These require longer discussions, not discrete answers.

Exercise 10: Data & visualization drill

# STEP 1: Load packages and data, and examine data structure
library(tidyverse)
peaks <- read_csv("https://mac-stat.github.io/data/high_peaks.csv")
head(peaks)
## # A tibble: 6 × 7
##   peak           elevation difficulty ascent length  time rating   
##   <chr>              <dbl>      <dbl>  <dbl>  <dbl> <dbl> <chr>    
## 1 Mt. Marcy           5344          5   3166   14.8  10   moderate 
## 2 Algonquin Peak      5114          5   2936    9.6   9   moderate 
## 3 Mt. Haystack        4960          7   3570   17.8  12   difficult
## 4 Mt. Skylight        4926          7   4265   17.9  15   difficult
## 5 Whiteface Mtn.      4867          4   2535   10.4   8.5 easy     
## 6 Dix Mtn.            4857          5   2800   13.2  10   moderate
dim(peaks)
## [1] 46  7

# How many hikes are in the data set?
nrow(peaks)
## [1] 46

# Show just the peak's name and its elevation for the first 6 hikes
peaks %>% 
  select(peak, elevation) %>% 
  head()
## # A tibble: 6 × 2
##   peak           elevation
##   <chr>              <dbl>
## 1 Mt. Marcy           5344
## 2 Algonquin Peak      5114
## 3 Mt. Haystack        4960
## 4 Mt. Skylight        4926
## 5 Whiteface Mtn.      4867
## 6 Dix Mtn.            4857

# Calculate the average hiking time
peaks %>% 
  summarize(mean(time))
## # A tibble: 1 × 1
##   `mean(time)`
##          <dbl>
## 1         10.7

# Show data for the hikes that are at least 17 miles long
peaks %>% 
  filter(length >= 17)
## # A tibble: 7 × 7
##   peak          elevation difficulty ascent length  time rating   
##   <chr>             <dbl>      <dbl>  <dbl>  <dbl> <dbl> <chr>    
## 1 Mt. Haystack       4960          7   3570   17.8  12   difficult
## 2 Mt. Skylight       4926          7   4265   17.9  15   difficult
## 3 Mt. Redfield       4606          7   3225   17.5  14   difficult
## 4 Panther Peak       4442          6   3762   17.6  13.5 moderate 
## 5 Mt. Donaldson      4140          7   3490   17    17   difficult
## 6 Mt. Emmons         4040          7   3490   18    18   difficult
## 7 Cliff Mtn.         3960          6   2160   17.2  12   moderate

# Calculate the average hiking time for hikes that are more than 14 miles long
peaks %>% 
  filter(length >= 17) %>% 
  summarize(mean(time))
## # A tibble: 1 × 1
##   `mean(time)`
##          <dbl>
## 1         14.5

# Calculate the average hiking time for hikes with an "easy" rating
peaks %>% 
  filter(rating == "easy") %>% 
  summarize(mean(time))
## # A tibble: 1 × 1
##   `mean(time)`
##          <dbl>
## 1            8

# Construct a visualization of hike rating
peaks %>% 
  ggplot(aes(x = rating)) + 
  geom_bar()

# Construct a visualization of hike length
# a histogram or boxplot would also work
# BUT a bar plot would NOT be appropriate
peaks %>% 
  ggplot(aes(x = length)) + 
  geom_density()

# Construct a visualization of the relationship of hiking time with length (distance)
peaks %>% 
  ggplot(aes(y = time, x = length)) + 
  geom_point()

3. observations
   • The majority of hikes are rated as “moderate”, with “difficult” hikes being the least numerous.
   • Hike lengths are roughly normally distributed, with a slight left skew, around an average hike length of 12 miles and ranging from as short as ~6 miles to as long as ~18 miles.
   • There’s a strong, positive association between hiking time and hike length – the longer the hike, the longer it tends to take to complete the hike.

Exercise 11: Types of models

1. Yes. Below is an example:

2. No. It’s possible to have a strong relationship between Y and X, but we model that relationship incorrectly (e.g. if they have a quadratic relationship but we model it with a line). But it’s impossible for a wrong model to be strong.

Exercise 11: Model practice

temps <- read_csv("https://mac-stat.github.io/data/weather_3_locations.csv") %>% 
  filter(location == "Hobart") %>% 
  select(date, temp3pm, temp9am)

head(temps)
## # A tibble: 6 × 3
##   date       temp3pm temp9am
##   <date>       <dbl>   <dbl>
## 1 2020-01-01    22.1    17.4
## 2 2020-01-02    19      19  
## 3 2020-01-03    19.7    16.2
## 4 2020-01-04    20.6    19.9
## 5 2020-01-05    21.9    14.2
## 6 2020-01-06    20.8    15.3

Part a

temp3pm = response, temp9am = predictor. it doesn’t make contextual sense to predict 9am temperature from 3pm temperature (without time travel).

Part b

There is a relatively strong, positive, linear assocation between these variables – the warmer it is at 9am, the warmer we expect it to be at 3pm.

temps %>% 
  ggplot(aes(y = temp3pm, x = temp9am)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE)

Part c

E[temp3pm | temp9am] = 5.21 + 0.88 temp9am

• When the 9am temperature is 0 degrees, the expected 3pm temperature is 5.21 degrees. Or, among days that are 0 degrees at 9am, the average 3pm temperature is 5.21 degrees.

• For every 1 degree increase in 3pm temperature, the average / expected 3pm temperature increases by 0.88 degrees.

temp_model <- lm(temp3pm ~ temp9am, temps)
summary(temp_model)
## 
## Call:
## lm(formula = temp3pm ~ temp9am, data = temps)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.2361  -1.8142  -0.1131   1.6880  17.0775 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.20618    0.32983   15.78   <2e-16 ***
## temp9am      0.88022    0.02448   35.95   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.963 on 785 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.6222, Adjusted R-squared:  0.6217 
## F-statistic:  1293 on 1 and 785 DF,  p-value: < 2.2e-16