Lately I’ve been publishing screencasts demonstrating how to use the tidymodels framework, from first steps in modeling to how to evaluate complex models. Today’s screencast is good for folks just getting started with tidymodels, using this week’s #TidyTuesday dataset on penguins. 🐧
Here is the code I used in the video, for those who prefer reading instead of or in addition to video.
Explore the data
This week’s #TidyTuesday dataset is from palmerpenguins, observations of Antarctic penguins who live on the Palmer Archipelago. You can read more about how this dataset came to be in this post on the RStudio Education blog. Our modeling goal here is to predict the sex of the penguins using a classification model, based on other observations in the dataset.
library(tidyverse)
library(palmerpenguins)
penguins
## # A tibble: 344 x 8
## species island bill_length_mm bill_depth_mm flipper_length_…
## <fct> <fct> <dbl> <dbl> <int>
## 1 Adelie Torge… 39.1 18.7 181
## 2 Adelie Torge… 39.5 17.4 186
## 3 Adelie Torge… 40.3 18 195
## 4 Adelie Torge… NA NA NA
## 5 Adelie Torge… 36.7 19.3 193
## 6 Adelie Torge… 39.3 20.6 190
## 7 Adelie Torge… 38.9 17.8 181
## 8 Adelie Torge… 39.2 19.6 195
## 9 Adelie Torge… 34.1 18.1 193
## 10 Adelie Torge… 42 20.2 190
## # … with 334 more rows, and 3 more variables: body_mass_g <int>,
## # sex <fct>, year <int>
If you try building a classification model for species, you will likely find an almost perfect fit, because these kinds of observations are actually what distinguish different species. Sex, on the other hand, is a little messier.
penguins %>%
filter(!is.na(sex)) %>%
ggplot(aes(flipper_length_mm, bill_length_mm, color = sex, size = body_mass_g)) +
geom_point(alpha = 0.5) +
facet_wrap(~species)
It looks like female penguins are smaller with different bills, but let’s get ready for modeling to find out more! We will not use the island or year information in our model.
penguins_df <- penguins %>%
filter(!is.na(sex)) %>%
select(-year, -island)
Build a model
We can start by loading the tidymodels metapackage, and splitting our data into training and testing sets.
library(tidymodels)
set.seed(123)
penguin_split <- initial_split(penguins_df, strata = sex)
penguin_train <- training(penguin_split)
penguin_test <- testing(penguin_split)
Next, let’s create bootstrap resamples of the training data, to evaluate our models.
set.seed(123)
penguin_boot <- bootstraps(penguin_train)
penguin_boot
## # Bootstrap sampling
## # A tibble: 25 x 2
## splits id
## <list> <chr>
## 1 <split [250/93]> Bootstrap01
## 2 <split [250/92]> Bootstrap02
## 3 <split [250/90]> Bootstrap03
## 4 <split [250/92]> Bootstrap04
## 5 <split [250/86]> Bootstrap05
## 6 <split [250/88]> Bootstrap06
## 7 <split [250/96]> Bootstrap07
## 8 <split [250/89]> Bootstrap08
## 9 <split [250/96]> Bootstrap09
## 10 <split [250/90]> Bootstrap10
## # … with 15 more rows
Let’s compare two different models, a logistic regression model and a random forest model. We start by creating the model specifications.
glm_spec <- logistic_reg() %>%
set_engine("glm")
glm_spec
## Logistic Regression Model Specification (classification)
##
## Computational engine: glm
rf_spec <- rand_forest() %>%
set_mode("classification") %>%
set_engine("ranger")
rf_spec
## Random Forest Model Specification (classification)
##
## Computational engine: ranger
Next let’s start putting together a tidymodels workflow(), a helper object to help manage modeling pipelines with pieces that fit together like Lego blocks. Notice that there is no model yet: Model: None.
penguin_wf <- workflow() %>%
add_formula(sex ~ .)
penguin_wf
## ══ Workflow ══════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: None
##
## ── Preprocessor ──────────────────────────────────────────────────────
## sex ~ .
Now we can add a model, and the fit to each of the resamples. First, we can fit the logistic regression model.
glm_rs <- penguin_wf %>%
add_model(glm_spec) %>%
fit_resamples(
resamples = penguin_boot,
control = control_resamples(save_pred = TRUE)
)
glm_rs
## # Resampling results
## # Bootstrap sampling
## # A tibble: 25 x 5
## splits id .metrics .notes .predictions
## <list> <chr> <list> <list> <list>
## 1 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [93 × …
## 2 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [92 × …
## 3 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [90 × …
## 4 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [92 × …
## 5 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [86 × …
## 6 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [88 × …
## 7 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [96 × …
## 8 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [89 × …
## 9 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [96 × …
## 10 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [90 × …
## # … with 15 more rows
Second, we can fit the random forest model.
rf_rs <- penguin_wf %>%
add_model(rf_spec) %>%
fit_resamples(
resamples = penguin_boot,
control = control_resamples(save_pred = TRUE)
)
rf_rs
## # Resampling results
## # Bootstrap sampling
## # A tibble: 25 x 5
## splits id .metrics .notes .predictions
## <list> <chr> <list> <list> <list>
## 1 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [93 × …
## 2 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [92 × …
## 3 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [90 × …
## 4 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [92 × …
## 5 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [86 × …
## 6 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [88 × …
## 7 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [96 × …
## 8 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [89 × …
## 9 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [96 × …
## 10 <split [250/… Bootstra… <tibble [2 × … <tibble [0 … <tibble [90 × …
## # … with 15 more rows
We have fit each of our candidate models to our resampled training set!
Evaluate model
Now let’s check out how we did.
collect_metrics(rf_rs)
## # A tibble: 2 x 5
## .metric .estimator mean n std_err
## <chr> <chr> <dbl> <int> <dbl>
## 1 accuracy binary 0.893 25 0.00691
## 2 roc_auc binary 0.958 25 0.00366
Pretty nice! The function collect_metrics() extracts and formats the .metrics column from resampling results like the ones we have here.
collect_metrics(glm_rs)
## # A tibble: 2 x 5
## .metric .estimator mean n std_err
## <chr> <chr> <dbl> <int> <dbl>
## 1 accuracy binary 0.897 25 0.00631
## 2 roc_auc binary 0.964 25 0.00368
So… also great! If I am in a situation where a more complex model like a random forest performs the same as a simpler model like logistic regression, then I will choose the simpler model. Let’s dig deeper into how it is doing. For example, how is it predicting the two classes?
glm_rs %>%
conf_mat_resampled()
## # A tibble: 4 x 3
## Prediction Truth Freq
## <fct> <fct> <dbl>
## 1 female female 40.6
## 2 female male 4.48
## 3 male female 4.92
## 4 male male 41.4
About the same, which is good. We can also make an ROC curve.
glm_rs %>%
collect_predictions() %>%
group_by(id) %>%
roc_curve(sex, .pred_female) %>%
ggplot(aes(1 - specificity, sensitivity, color = id)) +
geom_abline(lty = 2, color = "gray80", size = 1.5) +
geom_path(show.legend = FALSE, alpha = 0.6, size = 1.2) +
coord_equal()
This ROC curve is more jagged than others you may have seen because the dataset is small.
It is finally time for us to return to the testing set. Notice that we have not used the testing set yet during this whole analysis; the testing set is precious and can only be used to estimate performance on new data. Let’s fit one more time to the training data and evaluate on the testing data using the function last_fit().
penguin_final <- penguin_wf %>%
add_model(glm_spec) %>%
last_fit(penguin_split)
penguin_final
## # Resampling results
## # Monte Carlo cross-validation (0.75/0.25) with 1 resamples
## # A tibble: 1 x 6
## splits id .metrics .notes .predictions .workflow
## <list> <chr> <list> <list> <list> <list>
## 1 <split [2… train/tes… <tibble [2 … <tibble … <tibble [83 … <workflo…
The metrics and predictions here are on the testing data.
collect_metrics(penguin_final)
## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.940
## 2 roc_auc binary 0.991
collect_predictions(penguin_final) %>%
conf_mat(sex, .pred_class)
## Truth
## Prediction female male
## female 39 3
## male 2 39
The coefficients (which we can get out using tidy()) have been estimated using the training data. If we use exponentiate = TRUE, we have odds ratios.
penguin_final$.workflow[[1]] %>%
tidy(exponentiate = TRUE)
## # A tibble: 7 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 3.12e-35 13.5 -5.90 0.00000000369
## 2 speciesChinstrap 1.34e- 3 1.70 -3.89 0.000101
## 3 speciesGentoo 1.08e- 4 2.89 -3.16 0.00159
## 4 bill_length_mm 1.78e+ 0 0.137 4.20 0.0000268
## 5 bill_depth_mm 3.89e+ 0 0.373 3.64 0.000273
## 6 flipper_length_mm 1.07e+ 0 0.0538 1.31 0.189
## 7 body_mass_g 1.01e+ 0 0.00108 4.70 0.00000260
- The largest odds ratio is for bill depth, with the second largest for bill length. An increase of 1 mm in bill depth corresponds to almost 4x higher odds of being male. The characteristics of a penguin’s bill must be associated with their sex.
- We don’t have strong evidence that flipper length is different between male and female penguins, controlling for the other measures; maybe we should explore that by changing that first plot!
penguins %>%
filter(!is.na(sex)) %>%
ggplot(aes(bill_depth_mm, bill_length_mm, color = sex, size = body_mass_g)) +
geom_point(alpha = 0.5) +
facet_wrap(~species)
Yes, the male and female penguins are much more separated now.
Top comments (0)