Notes and in-class exercises
By the end of this lesson, you should be able to:
Today is a day to discover ideas, so no readings or videos to go through before class, but if you want to see today’s ideas presented in a different way, you can take a look at the following after class:
Context: Last time we talked about different types of research questions (descriptive, predictive, and causal) and worked through how to translate these questions to data explorations and modeling. Also recall that when we talked about model evaluation previously, we asked if the model was correct, strong, and fair.
Today we’ll revisit predictive research questions today and some considerations when trying to build strong models. In particular, we’ll consider nuances and limitations to indiscriminately adding more predictors to our model. To explore these ideas, we’ll explore data on physical characteristics of different penguin species:
# Load packages
library(ggplot2)
library(dplyr)
Attaching package: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, union# Load data
library(palmerpenguins)
data(penguins)You can find a codebook for these data by typing ?penguins in your console (not Rmd). Our goal throughout will be to build a model of bill lengths (in mm):
(Art by @allison_horst)
To get started, the flipper_length_mm variable currently measures flipper length in mm. Let’s create and save a new variable named flipper_length_cm which measures flipper length in cm. NOTE: There are 10mm in a cm.
penguins <- penguins %>%
mutate(flipper_length_cm = flipper_length_mm / 10)Run the code chunk below to build a bunch of models that you’ll be exploring in the exercises:
penguin_model_1 <- lm(bill_length_mm ~ flipper_length_mm, penguins)
penguin_model_2 <- lm(bill_length_mm ~ flipper_length_cm, penguins)
penguin_model_3 <- lm(bill_length_mm ~ flipper_length_mm + flipper_length_cm, penguins)
penguin_model_4 <- lm(bill_length_mm ~ body_mass_g, penguins)
penguin_model_5 <- lm(bill_length_mm ~ flipper_length_mm + body_mass_g, penguins)What can a penguin’s flipper (arm) length tell us about their bill length? To answer this question, we’ll consider 3 of our models:
| model | predictors |
|---|---|
penguin_model_1 |
flipper_length_mm |
penguin_model_2 |
flipper_length_cm |
penguin_model_3 |
flipper_length_mm + flipper_length_cm |
Plots of the first two models are below:
ggplot(penguins, aes(y = bill_length_mm, x = flipper_length_mm)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE)
ggplot(penguins, aes(y = bill_length_mm, x = flipper_length_cm)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE)Before examining the model summaries, check your intuition. Do you think the penguin_model_2 R-squared will be less than, equal to, or more than that of penguin_model_1? Similarly, how do you think the penguin_model_3 R-squared will compare to that of penguin_model_1?
Check your intuition: Examine the R-squared values for the three penguin models and summarize how these compare.
summary(penguin_model_1)$r.squared
summary(penguin_model_2)$r.squared
summary(penguin_model_3)$r.squaredflipper_length_mm vs flipper_length_cm.summary(penguin_model_3), the flipper_length_cm coefficient is NA. Explain why this makes sense. HINT: Thinking about what you learned about controlling for covariates, why wouldn’t it make sense to interpret this coefficient? BONUS: For those of you that have taken MATH 236, this has to do with matrices that are not of full rank!body_mass_gIn this exercise you’ll consider 3 models of bill_length_mm:
| model | predictors |
|---|---|
penguin_model_1 |
flipper_length_mm |
penguin_model_4 |
body_mass_g |
penguin_model_5 |
flipper_length_mm + body_mass_g |
Which is the better predictor of bill_length_mm: flipper_length_mm or body_mass_g? Provide some numerical evidence.
penguin_model_5 incorporates both flipper_length_mm and body_mass_g as predictors. Before examining a model summary, ask your gut: Will the penguin_model_5 R-squared be close to 0.35, close to 0.43, or greater than 0.6?
Check your intuition. Report the penguin_model_5 R-squared and summarize how this compares to that of penguin_model_1 and penguin_model_4.
Explain why your observation in part c makes sense. Support your reasoning with a plot of the 2 predictors: flipper_length_mm vs body_mass_g.
The exercises above have illustrated special phenomena in multivariate modeling:
Recall that we examined 5 models:
| model | predictors |
|---|---|
penguin_model_1 |
flipper_length_mm |
penguin_model_2 |
flipper_length_cm |
penguin_model_3 |
flipper_length_mm + flipper_length_cm |
penguin_model_4 |
body_mass_g |
penguin_model_5 |
flipper_length_mm + body_mass_g |
Let’s dive deeper into important considerations when building a strong model. We’ll use a subset of the penguins data for exploring these ideas.
# For illustration purposes only, take a sample of 10 penguins.
# We'll discuss this code later in the course!
set.seed(155)
penguins_small <- sample_n(penguins, size = 10) %>%
mutate(flipper_length_mm = jitter(flipper_length_mm))Consider 3 models of bill length:
# A model with one predictor (flipper_length_mm)
poly_mod_1 <- lm(bill_length_mm ~ flipper_length_mm, penguins_small)
# A model with two predictors (flipper_length_mm and flipper_length_mm^2)
poly_mod_2 <- lm(bill_length_mm ~ poly(flipper_length_mm, 2), penguins_small)
# A model with nine predictors (flipper_length_mm, flipper_length_mm^2, ... on up to flipper_length_mm^9)
poly_mod_9 <- lm(bill_length_mm ~ poly(flipper_length_mm, 9), penguins_small)summary(poly_mod_1)$r.squared
summary(poly_mod_2)$r.squared
summary(poly_mod_9)$r.squared# A plot of model 1
ggplot(penguins_small, aes(y = bill_length_mm, x = flipper_length_mm)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE)# A plot of model 2
ggplot(penguins_small, aes(y = bill_length_mm, x = flipper_length_mm)) +
geom_point() +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = FALSE)# A plot of model 9
ggplot(penguins_small, aes(y = bill_length_mm, x = flipper_length_mm)) +
geom_point() +
geom_smooth(method = "lm", formula = y ~ poly(x, 9), se = FALSE)bill_length_mm and flipper_length_mm?penguins_small sample? For example, would you expect these new penguins to fall on the wiggly model 9 curve?Model 9 provides an example of a model that is overfit to our sample data. That is, it picks up the tiny details of our data at the cost of losing the more general trends of the relationship of interest. Check out the following xkcd comic. Which plot pokes fun at overfitting?
Some other goodies:
Zooming out, explain some limitations of relying on R-squared to measure the strength / usefulness of a model.
We’ve seen that, unless a predictor is redundant with another, R-squared will increase. Even if that predictor is strongly multicollinear with another. Even if that predictor isn’t a good predictor! Thus if we only look at R-squared we might get overly greedy. We can check our greedy impulses a few ways. We take a more in depth approach in STAT 253, but one quick alternative is reported right in our model summary() tables. Adjusted R-squared includes a penalty for incorporating more and more predictors. Mathematically (where \(n\) is the sample size and \(p\) is the number of non-intercept coefficients):
\[ \text{Adjusted } R^2 = 1 - (1 - R^2) \left( \frac{n-1}{n-p-1} \right) \]
Thus unlike R-squared, Adjusted R-squared can decrease when the information that a predictor contributes to a model isn’t enough to offset the complexity it adds to that model. Consider two models:
example_1 <- lm(bill_length_mm ~ species, penguins)
example_2 <- lm(bill_length_mm ~ species + island, penguins)island that resulted in a decreased Adjusted R-squared. Note: it’s not necessarily the case that island is a bad predictor on its own!Today we looked at some cautions surrounding indiscriminately adding variables to a model. Summarize key takeaways.
Response: Put your response here.
| model | predictors |
|---|---|
penguin_model_1 |
flipper_length_mm |
penguin_model_2 |
flipper_length_cm |
penguin_model_3 |
flipper_length_mm + flipper_length_cm |
Plots of the first two models are below:
ggplot(penguins, aes(y = bill_length_mm, x = flipper_length_mm)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE)`geom_smooth()` using formula = 'y ~ x'Warning: Removed 2 rows containing non-finite outside the scale range
(`stat_smooth()`).Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).
ggplot(penguins, aes(y = bill_length_mm, x = flipper_length_cm)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE)`geom_smooth()` using formula = 'y ~ x'Warning: Removed 2 rows containing non-finite outside the scale range (`stat_smooth()`).
Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).
Your intuition–answers will vary
The R-squared values are all the same!
summary(penguin_model_1)$r.squared[1] 0.430574summary(penguin_model_2)$r.squared[1] 0.430574summary(penguin_model_3)$r.squared[1] 0.430574ggplot(penguins, aes(x = flipper_length_cm, y = flipper_length_mm)) +
geom_point()Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).
NA means that the coefficient couldn’t be estimated. In penguin_model_3, the interpretation of the flipper_length_cm coefficient is the average change in bill length per centimeter change in flipper length, while holding flipper length in millimeters constant…this is impossible! We can’t hold flipper length in millimeters fixed while varying flipper length in centimeters—if one changes the other must. (In linear algebra terms, the matrix underlying our data is not of full rank.)body_mass_gIn this exercise you’ll consider 3 models of bill_length_mm:
| model | predictors |
|---|---|
penguin_model_1 |
flipper_length_mm |
penguin_model_4 |
body_mass_g |
penguin_model_5 |
flipper_length_mm + body_mass_g |
flipper_length_mm is a better predictor than body_mass_g because penguin_model_1 has an R-squared value of 0.4306 vs 0.3542 for penguin_model_4.
Intuition check–answers will vary
R-squared is for penguin_model_5 which is slightly higher than that of penguin_model_1 and penguin_model_4.
d.flipper_length_mm and body_mass_g are positively correlated and thus contain related information, but not completely redundant information. There’s some information in flipper length in explaining bill length that isn’t captured by body mass, and vice-versa.
ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g)) +
geom_point()Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).
| model | predictors |
|---|---|
penguin_model_1 |
flipper_length_mm |
penguin_model_2 |
flipper_length_cm |
penguin_model_3 |
flipper_length_mm + flipper_length_cm |
penguin_model_4 |
body_mass_g |
penguin_model_5 |
flipper_length_mm + body_mass_g |
penguin_model_3 had redundant predictors: `flipper_length_mm and flipper_length_cmpenguin_model_5 had multicollinear predictors: flipper_length_mm and body_mass_g were related but not redundantpoly_mod_9 might seem to be the best. BUT we’ll see where this reasoning is flawed soon!The bottom left plot pokes fun at overfitting.
It measures how well our model explains / predicts our sample data, not how well it explains / predicts the broader population. It also has the feature that any non-redundant predictor added to a model will increase the R-squared.
example_1 <- lm(bill_length_mm ~ species, penguins)
example_2 <- lm(bill_length_mm ~ species + island, penguins)
summary(example_1)
Call:
lm(formula = bill_length_mm ~ species, data = penguins)
Residuals:
Min 1Q Median 3Q Max
-7.9338 -2.2049 0.0086 2.0662 12.0951
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 38.7914 0.2409 161.05 <2e-16 ***
speciesChinstrap 10.0424 0.4323 23.23 <2e-16 ***
speciesGentoo 8.7135 0.3595 24.24 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.96 on 339 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.7078, Adjusted R-squared: 0.7061
F-statistic: 410.6 on 2 and 339 DF, p-value: < 2.2e-16summary(example_2)
Call:
lm(formula = bill_length_mm ~ species + island, data = penguins)
Residuals:
Min 1Q Median 3Q Max
-7.9338 -2.2049 -0.0049 2.0951 12.0951
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 38.97500 0.44697 87.198 <2e-16 ***
speciesChinstrap 10.33204 0.53502 19.312 <2e-16 ***
speciesGentoo 8.52988 0.52082 16.378 <2e-16 ***
islandDream -0.47321 0.59729 -0.792 0.429
islandTorgersen -0.02402 0.61004 -0.039 0.969
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.965 on 337 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.7085, Adjusted R-squared: 0.7051
F-statistic: 204.8 on 4 and 337 DF, p-value: < 2.2e-16island didn’t provide useful information about bill length beyond what was already provided by species.