reshape2 package today: install this if you haven’t already.
Highlight: Collaborative Learning
Take 5 minutes to reflect upon your work throughout Unit 1, particularly with respect to collaboration.
Reflect upon your strengths and what you might change in the next unit:
If you did not finish Exercises #9 and #10 from last class, please take time after class today to do so.
world = supervised learning
We want to model some output variable \(y\) using a set of potential predictors (\(x_1, x_2, ..., x_p\)).
task = regression
\(y\) is quantitative
model = linear regression
We’ll assume that the relationship between \(y\) and (\(x_1, x_2, ..., x_p\)) can be represented by
\[y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_p x_p + \varepsilon\]
In model building, the decision of which predictors to use depends upon our goal.
Inferential models
Predictive models
Model selection algorithms can help us build a predictive model of \(y\) using a set of potential predictors (\(x_1, x_2, ..., x_p\)).
There are 3 general approaches to this task:
Variable selection (today)
Identify a subset of predictors to use in our model of \(y\).
Shrinkage / regularization (next class)
Shrink / regularize the coefficients of all predictors toward or to 0.
Dimension reduction (later in the semester)
Combine the predictors into a smaller set of new predictors.
Open the Part 1 QMD. Scroll down to the # Exercises section.
As a group, you’ll design a variable selection algorithm to pick which predictors to use in a predictive model of height. Specifically, you will:
NOTE: This will NOT be perfect! Our goals are to:
Let’s build a predictive model of height in inches using one or more of 12 possible predictors. Other than age and weight, these are circumferences measured in cm.
# Load packages
library(tidyverse)
library(tidymodels)
# Load data
humans <- read.csv("https://kegrinde.github.io/stat253_coursenotes/data/bodyfat1.csv")
names(humans) [1] "age" "weight" "neck" "chest" "abdomen" "hip" "thigh"
[8] "knee" "ankle" "biceps" "forearm" "wrist" "height" A heat map displays correlations for each pair of variables in our dataset. Not only is height correlated with multiple predictors, the predictors are correlated with one another (mulicollinear)! We don’t need all of them in our model.
# Get the correlation matrix
library(reshape2)
cor_matrix <- cor(humans)
cor_matrix[lower.tri(cor_matrix)] <- NA
cor_matrix <- cor_matrix %>%
melt() %>%
na.omit() %>%
rename(correlation = value)
# Visualize the correlation for each pair of variables
ggplot(cor_matrix, aes(x = Var1, y = Var2, fill = correlation)) +
geom_tile(color = "white") +
scale_fill_gradient2(
low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1)) +
labs(x = "", y = "") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
coord_fixed()
# STEP 1: model specification
lm_spec <- linear_reg() %>%
set_mode("regression") %>%
set_engine("lm")
# STEP 2: model estimation
height_model_1 <- lm_spec %>%
fit(height ~ age + weight + neck + chest + abdomen + hip + thigh + knee + ankle + biceps + forearm + wrist, data = humans)
# Check it out
height_model_1 %>%
tidy()
# CV MAE
set.seed(253)
lm_spec %>%
fit_resamples(
height ~ age + weight + neck + chest + abdomen + hip + thigh + knee + ankle + biceps + forearm + wrist,
resamples = vfold_cv(humans, v = 10),
metrics = metric_set(mae)
) %>%
collect_metrics()Open the Part 2 QMD to take notes.
Let’s consider three existing variable selection algorithms.
Heads up: these algorithms are important to building intuition for the questions and challenges in model selection, BUT have major drawbacks.
Suppose we used this algorithm for our height model with 12 possible predictors. What’s the main drawback?
It’s computationally expensive. For our humans example, we’d need to build 4096 models:
2^12[1] 4096. . .
Let’s try out the first few steps!
. . .
# Load packages and data
library(tidyverse)
library(tidymodels)
humans <- read.csv("https://kegrinde.github.io/stat253_coursenotes/data/bodyfat1.csv")# STEP 1: model specifications
lm_spec <- linear_reg() %>%
set_mode("regression") %>%
set_engine("lm")# STEP 2: model estimate (using all 12 predictors to start)
# Pick apart this code and make it easier to identify the least "significant" predictor!!!
lm_spec %>%
fit(height ~ age + weight + neck + chest + abdomen + hip + thigh + knee + ankle + biceps + forearm + wrist,
data = humans) %>%
tidy() %>%
filter(term != "(Intercept)") %>%
mutate(p.value = round(p.value, 4))# 11 predictors (tweak the code)
lm_spec %>%
fit(height ~ age + weight + neck + chest + abdomen + hip + thigh + knee + ankle + biceps + forearm + wrist,
data = humans) %>%
tidy() %>%
filter(term != "(Intercept)") %>%
mutate(p.value = round(p.value, 4))# 10 predictors (tweak the code)
lm_spec %>%
fit(height ~ age + weight + neck + chest + abdomen + hip + thigh + knee + ankle + biceps + forearm + wrist,
data = humans) %>%
tidy() %>%
filter(term != "(Intercept)") %>%
mutate(p.value = round(p.value, 4))# All 12 predictors
lm_spec %>%
fit(height ~ age + weight + neck + chest + abdomen + hip + thigh + knee + ankle + biceps + forearm + wrist,
data = humans) %>%
tidy() %>% # use tidy to get p-values for each coefficient
filter(term != "(Intercept)") %>% # exclude the intercept
mutate(p.value = round(p.value, 4)) %>% # round the p-values for easier viewing
arrange(desc(p.value)) # added this line to arrange from largest to smallest p-value# A tibble: 12 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 biceps -0.0808 0.746 -0.108 0.915
2 neck 0.139 1.17 0.119 0.906
3 knee 0.151 0.941 0.160 0.874
4 wrist 0.836 2.32 0.361 0.721
5 ankle -0.888 1.28 -0.693 0.494
6 abdomen 0.283 0.354 0.798 0.432
7 age -0.112 0.132 -0.847 0.405
8 chest -0.459 0.473 -0.971 0.340
9 forearm 2.25 1.80 1.25 0.223
10 weight 0.379 0.213 1.78 0.0864
11 hip -0.921 0.510 -1.81 0.0822
12 thigh -1.24 0.646 -1.92 0.066 # 11 predictors (got rid of biceps)
lm_spec %>%
fit(height ~ age + weight + neck + chest + abdomen + hip + thigh + knee + ankle + forearm + wrist,
data = humans) %>%
tidy() %>%
filter(term != "(Intercept)") %>%
mutate(p.value = round(p.value, 4)) %>%
arrange(desc(p.value))# A tibble: 11 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 neck 0.161 1.14 0.142 0.888
2 knee 0.180 0.886 0.203 0.841
3 wrist 0.907 2.18 0.416 0.681
4 ankle -0.878 1.26 -0.699 0.490
5 abdomen 0.281 0.348 0.809 0.425
6 age -0.111 0.130 -0.858 0.398
7 chest -0.453 0.461 -0.982 0.334
8 forearm 2.17 1.62 1.34 0.192
9 hip -0.902 0.470 -1.92 0.0652
10 weight 0.369 0.190 1.94 0.0623
11 thigh -1.26 0.602 -2.09 0.0454# 10 predictors (got rid of neck)
lm_spec %>%
fit(height ~ age + weight + chest + abdomen + hip + thigh + knee + ankle + forearm + wrist,
data = humans) %>%
tidy() %>%
filter(term != "(Intercept)") %>%
mutate(p.value = round(p.value, 4)) %>%
arrange(desc(p.value))# A tibble: 10 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 knee 0.166 0.866 0.191 0.850
2 wrist 0.985 2.07 0.475 0.639
3 ankle -0.884 1.23 -0.716 0.480
4 age -0.111 0.127 -0.869 0.392
5 abdomen 0.298 0.322 0.924 0.363
6 chest -0.460 0.451 -1.02 0.316
7 forearm 2.29 1.37 1.66 0.107
8 weight 0.377 0.179 2.11 0.0435
9 thigh -1.26 0.591 -2.14 0.0409
10 hip -0.931 0.416 -2.24 0.0331Below is the complete model sequence along with 10-fold CV MAE for each model (using set.seed(253)).
| pred | CV MAE | predictor list |
|---|---|---|
| 12 | 5.728 | weight, hip, forearm, thigh, chest, abdomen, age, ankle, wrist, knee, neck, biceps |
| 11 | 5.523 | weight, hip, forearm, thigh, chest, abdomen, age, ankle, wrist, knee, neck |
| 10 | 5.413 | weight, hip, forearm, thigh, chest, abdomen, age, ankle, wrist, knee |
| 9 | 5.368 | weight, hip, forearm, thigh, chest, abdomen, age, ankle, wrist |
| 8 | 5.047 | weight, hip, forearm, thigh, chest, abdomen, age, ankle |
| 7 | 5.013 | weight, hip, forearm, thigh, chest, abdomen, age |
| 6 | 4.684 | weight, hip, forearm, thigh, chest, abdomen |
| 5 | 4.460 | weight, hip, forearm, thigh, chest |
| 4 | 4.386 | weight, hip, forearm, thigh |
| 3 | 4.091 | weight, hip, forearm |
| 2 | 3.733 | weight, hip |
| 1 | 3.658 | weight |
DISCUSS:
(Review) Interpret the CV MAE for the model of height by weight alone.
Is this algorithm more or less computationally expensive than the best subset algorithm?
The predictors neck and wrist, in that order, are the most strongly correlated with height. Where do these appear in the backward sequence and what does this mean?
cor(humans)[,'height'] %>%
sort() thigh hip age abdomen knee chest
-0.11301249 -0.10648937 -0.05853538 -0.02173587 0.02345904 0.05838830
biceps ankle weight forearm wrist neck
0.07441696 0.07920867 0.11228791 0.16968040 0.28967468 0.29147610
height
1.00000000 We have to pick just 1 of the 12 models as our final model.
That is, we have to pick a value for our tuning parameter, the number of predictors.
It helps to plot the CV MAE for each model in the sequence.
Here’s what we saw above:
data.frame(
predictors = c(12:1),
mae = c(5.728, 5.523, 5.413, 5.368, 5.047, 5.013, 4.684, 4.460, 4.386, 4.091, 3.733, 3.658)) %>%
ggplot(aes(x = predictors, y = mae)) +
geom_point() +
geom_line() +
scale_x_continuous(breaks = c(1:12))
Here’s another example from a different subset of these data:
When tuning or finalizing a model building algorithm, we (humans!) have our own choices to make. For one, we need to decide what we prefer:
In deciding, here are some human considerations:
For each scenario below, which model would you pick: (1) the model with the lowest prediction errors; or (2) a parsimonious model with slightly worse predictions?
Google asks us to re-build their search algorithm.
A small non-profit hires us to help them build a predictive model of the donation dollars they’ll receive throughout the year.
Start with 0 predictors. Add the predictor with the smallest p-value. To this model, add a second predictor with the smallest p-value. Continue until all predictors are in the model.
more. For 12 predictors, we’d have to build 12 models in step 1, 11 models in step 2, etc. Thus 12 + 11 + … + 1 = 78 models total.
Variable selection algorithms are a nice, intuitive place to start our discussion of model selection techniques.
BUT we will not use them.
They are frowned upon in the broader ML community, so much so that tidymodels doesn’t even implement them! Why?
# Exercises section in the Part 2 QMD.The video for today introduced the concepts of recipes and workflows in the tidymodels framework. These concepts will become important to our new modeling algorithms. Though they aren’t necessary to linear regression models, let’s explore them in this familiar setting.
Run through the following discussion and code one step at a time. Take note of the general process, concepts, and questions you have.
STEP 1: model specification
This specifies the structure or general modeling algorithm we plan to use.
It does not specify anything about the variables of interest or our data.
lm_spec <- linear_reg() %>%
set_mode("regression") %>%
set_engine("lm")
# Check it out
lm_specSTEP 2: recipe specification
Just as a cooking recipe specifies the ingredients and how to prepare them, a tidymodels recipe specifies:
It does not specify anything about the model structure we’ll use to explore this relationship.
# A simple recipe with NO pre-processing
data_recipe <- recipe(height ~ wrist + ankle, data = humans)
# Check it out
data_recipeSTEP 3: workflow creation (model + recipe)
This specifies the general workflow of our modeling process, including our model structure and our variable recipe.
model_workflow <- workflow() %>%
add_recipe(data_recipe) %>%
add_model(lm_spec)
# Check it out
model_workflowSTEP 4: Model estimation
This step estimates or fits our model of interest using our entire sample data.
The model (lm_spec) and variable details (here just height ~ wrist + ankle) are specified in the workflow, so we do not need to give that information again!
my_model <- model_workflow %>%
fit(data = humans)STEPS 5: Model evaluation
To get in-sample metrics, use my_model like normal.
# example: calculate of in-sample metrics
my_model %>%
glance()To get CV metrics, pass the workflow to fit_resamples along with information about how to randomly create folds.
set.seed(253)
my_model_cv <- model_workflow %>%
fit_resamples(resamples = vfold_cv(humans, v = 10),
metrics = metric_set(rsq))Then, proceed as usual… (my_model_cv %>% collect_metrics(), etc. )