DEV Community

Cover image for Get started with tidymodels and the Palmer penguins 🐧
Julia Silge
Julia Silge

Posted on • Originally published at juliasilge.com on

Get started with tidymodels and the Palmer penguins 🐧

#machinelearning #datascience #rstats #tutorial

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>
Enter fullscreen mode Exit fullscreen mode

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)
Enter fullscreen mode Exit fullscreen mode

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)
Enter fullscreen mode Exit fullscreen mode

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)
Enter fullscreen mode Exit fullscreen mode

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
Enter fullscreen mode Exit fullscreen mode

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
Enter fullscreen mode Exit fullscreen mode

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 ~ .
Enter fullscreen mode Exit fullscreen mode

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
Enter fullscreen mode Exit fullscreen mode

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
Enter fullscreen mode Exit fullscreen mode

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
Enter fullscreen mode Exit fullscreen mode

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
Enter fullscreen mode Exit fullscreen mode

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
Enter fullscreen mode Exit fullscreen mode

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()
Enter fullscreen mode Exit fullscreen mode

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…
Enter fullscreen mode Exit fullscreen mode

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
Enter fullscreen mode Exit fullscreen mode

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
Enter fullscreen mode Exit fullscreen mode
  • 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)
Enter fullscreen mode Exit fullscreen mode

Yes, the male and female penguins are much more separated now.

Top comments (0)

Read next

mikeyoung44 profile image

Reinforcement Learning Tunes Code LLMs With Execution Feedback

Mike Young -

techsourabh profile image

5 Best ETL Tools: A Comprehensive Comparison Guide

Sourabh Gupta -

sebconejo profile image

How to create a newsletter signup form in just 20 minutes with shadcn/UI and Manifest

Sébastien Conejo -

mikeyoung44 profile image

Universal Text Segmentation Model Outperforms Specialized Systems Across Diverse Domains

Mike Young -