Outline

Train a classification model to predict flower species based on the four available features

The model formula will have the form species ~ . where . represents all explanatory variables in the data.

Splitting the data into training and test sets

By default, initial split() provides a 75:25 split for our train and test sets respectively. Since our dataset is small to begin with, we are going to make bootstrap resamples from the training data. The function bootstraps() will split the data into training and test sets, then repeat this process with replacement a specified number of times (25 is the default).

set.seed(123)
tidy_split <- initial_split(iris_df)
tidy_split

## <Analysis/Assess/Total>
## <112/38/150>

iris_train <- training(tidy_split)
iris_test <- testing(tidy_split)

iris_boots <- bootstraps(iris_train, times = 30) 
iris_boots

## # Bootstrap sampling 
## # A tibble: 30 × 2
##    splits           id         
##    <list>           <chr>      
##  1 <split [112/45]> Bootstrap01
##  2 <split [112/43]> Bootstrap02
##  3 <split [112/39]> Bootstrap03
##  4 <split [112/40]> Bootstrap04
##  5 <split [112/39]> Bootstrap05
##  6 <split [112/41]> Bootstrap06
##  7 <split [112/35]> Bootstrap07
##  8 <split [112/37]> Bootstrap08
##  9 <split [112/42]> Bootstrap09
## 10 <split [112/37]> Bootstrap10
## # … with 20 more rows

Recipes

Recipes is a powerful tool with functions for a wide range of feature engineering tasks designed to prepare data for modeling. Typing recipes:: into the Rstudio console is a great way to browse the available functions in the package.

Let’s create a simple recipe to demonstrate optional feature engineering steps for our numeric data.

iris_rec <- recipe(species ~., data = iris_train) %>%
  step_pca(all_predictors()) %>%
  step_normalize(all_predictors())

prep <-  prep(iris_rec)

kable(head(iris_juice <- juice(prep)))

Creating models with Parsnip

Let’s set up two different models: first, a generalized linear model or glmnet. In this step we will create the model, workflow and fit the bootstraps. Let’s take a look at the output from each step.

# set seed
set.seed(1234)

# generate the glmnet model with parsnip
glmnet_mod <- multinom_reg(penalty = 0) %>% 
  set_engine("glmnet") %>% 
  set_mode("classification")
glmnet_mod

## Multinomial Regression Model Specification (classification)
## 
## Main Arguments:
##   penalty = 0
## 
## Computational engine: glmnet

# create a workflow
glmnet_wf <- workflow() %>%
  add_formula(species ~ .) 
glmnet_wf

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: None
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## species ~ .

# add the model to the workflow and use iris_boots to fit our model 25 times
glmnet_results <- glmnet_wf %>%
  add_model(glmnet_mod) %>% 
  fit_resamples(
    resamples = iris_boots,
    control = control_resamples(extract = extract_model,
                             save_pred = TRUE)
    )
glmnet_results

## # Resampling results
## # Bootstrap sampling 
## # A tibble: 30 × 6
##    splits           id          .metrics         .notes  .extracts  .predictions
##    <list>           <chr>       <list>           <list>  <list>     <list>      
##  1 <split [112/45]> Bootstrap01 <tibble [2 × 4]> <tibbl… <tibble [… <tibble [45…
##  2 <split [112/43]> Bootstrap02 <tibble [2 × 4]> <tibbl… <tibble [… <tibble [43…
##  3 <split [112/39]> Bootstrap03 <tibble [2 × 4]> <tibbl… <tibble [… <tibble [39…
##  4 <split [112/40]> Bootstrap04 <tibble [2 × 4]> <tibbl… <tibble [… <tibble [40…
##  5 <split [112/39]> Bootstrap05 <tibble [2 × 4]> <tibbl… <tibble [… <tibble [39…
##  6 <split [112/41]> Bootstrap06 <tibble [2 × 4]> <tibbl… <tibble [… <tibble [41…
##  7 <split [112/35]> Bootstrap07 <tibble [2 × 4]> <tibbl… <tibble [… <tibble [35…
##  8 <split [112/37]> Bootstrap08 <tibble [2 × 4]> <tibbl… <tibble [… <tibble [37…
##  9 <split [112/42]> Bootstrap09 <tibble [2 × 4]> <tibbl… <tibble [… <tibble [42…
## 10 <split [112/37]> Bootstrap10 <tibble [2 × 4]> <tibbl… <tibble [… <tibble [37…
## # … with 20 more rows

# look at the model metrics
collect_metrics(glmnet_results)

## # A tibble: 2 × 6
##   .metric  .estimator  mean     n std_err .config             
##   <chr>    <chr>      <dbl> <int>   <dbl> <chr>               
## 1 accuracy multiclass 0.958    30 0.00507 Preprocessor1_Model1
## 2 roc_auc  hand_till  0.994    30 0.00119 Preprocessor1_Model1

Now for a random forest model. We only need to change a few things and walah!

set.seed(1234)
rf_mod <- rand_forest() %>%
  set_engine("ranger") %>%
  set_mode("classification")

# We set up a workflow and add the parts of our model together like legos
rf_wf <- workflow() %>%
  add_formula(species ~ .)

# Here we fit our 25 resampled datasets 
rf_results <- rf_wf %>%
  add_model(rf_mod) %>% 
  fit_resamples(
    resamples = iris_boots,
    control = control_resamples(save_pred = TRUE)
    )
collect_metrics(rf_results)

## # A tibble: 2 × 6
##   .metric  .estimator  mean     n  std_err .config             
##   <chr>    <chr>      <dbl> <int>    <dbl> <chr>               
## 1 accuracy multiclass 0.953    30 0.00449  Preprocessor1_Model1
## 2 roc_auc  hand_till  0.995    30 0.000800 Preprocessor1_Model1

Here’s a look at the confusion matrix summaries for both models. The confusion matrix let’s us see the correct and incorrect predictions of our models in a single table.

glmnet_results %>%
  conf_mat_resampled() 

## # A tibble: 9 × 3
##   Prediction Truth        Freq
##   <fct>      <fct>       <dbl>
## 1 setosa     setosa     14    
## 2 setosa     versicolor  0    
## 3 setosa     virginica   0    
## 4 versicolor setosa      0    
## 5 versicolor versicolor 10.2  
## 6 versicolor virginica   0.833
## 7 virginica  setosa      0    
## 8 virginica  versicolor  0.867
## 9 virginica  virginica  14.2

rf_results %>%
  conf_mat_resampled() 

## # A tibble: 9 × 3
##   Prediction Truth        Freq
##   <fct>      <fct>       <dbl>
## 1 setosa     setosa     14    
## 2 setosa     versicolor  0    
## 3 setosa     virginica   0    
## 4 versicolor setosa      0    
## 5 versicolor versicolor 10.2  
## 6 versicolor virginica   1.03 
## 7 virginica  setosa      0    
## 8 virginica  versicolor  0.867
## 9 virginica  virginica  14

The ROC curve helps us visually interpret our model performance at every threshold.

glmnet_results %>%
  collect_predictions() %>%
  group_by(id) %>%
  roc_curve(species, .pred_setosa:.pred_virginica) %>%
  autoplot()

rf_results %>%
  collect_predictions() %>%
  group_by(id) %>%
  roc_curve(species, .pred_setosa:.pred_virginica) %>%
  autoplot() +
  theme(legend.position = 'none')

Final fit

By using the last_fit(tidy_split), we are able to train our model on the training set and test the model on the testing set in one fell swoop! Note, this is the only time we use the test set.

final_glmnet <- glmnet_wf %>%
    add_model(glmnet_mod) %>%
    last_fit(tidy_split)

final_glmnet

## # Resampling results
## # Manual resampling 
## # A tibble: 1 × 6
##   splits           id               .metrics   .notes    .predictions  .workflow
##   <list>           <chr>            <list>     <list>    <list>        <list>   
## 1 <split [112/38]> train/test split <tibble [… <tibble … <tibble [38 … <workflo…

final_rf <- rf_wf %>%
    add_model(rf_mod) %>%
    last_fit(tidy_split)

final_rf

## # Resampling results
## # Manual resampling 
## # A tibble: 1 × 6
##   splits           id               .metrics   .notes    .predictions  .workflow
##   <list>           <chr>            <list>     <list>    <list>        <list>   
## 1 <split [112/38]> train/test split <tibble [… <tibble … <tibble [38 … <workflo…

Confusion Matrices

Finally, let’s generate a multiclass confusion matrix with the results from our test data. The confusion matrix provides a count of each outcome for all possible outcomes. The columns contain the true values and the predictions are assigned to the rows.

collect_metrics(final_glmnet)

## # A tibble: 2 × 4
##   .metric  .estimator .estimate .config             
##   <chr>    <chr>          <dbl> <chr>               
## 1 accuracy multiclass     0.974 Preprocessor1_Model1
## 2 roc_auc  hand_till      0.991 Preprocessor1_Model1

collect_predictions(final_glmnet) %>%
  conf_mat(species, .pred_class) %>% 
  autoplot(type = 'heatmap') 

collect_metrics(final_rf)

## # A tibble: 2 × 4
##   .metric  .estimator .estimate .config             
##   <chr>    <chr>          <dbl> <chr>               
## 1 accuracy multiclass     0.974 Preprocessor1_Model1
## 2 roc_auc  hand_till      0.998 Preprocessor1_Model1

collect_predictions(final_rf) %>%
  conf_mat(species, .pred_class) %>% 
  autoplot(type = 'heatmap')

Final thoughts

Both models exhibit near perfect predictive power but are they really that good? From our visual analysis, we can confidently say that the combination of explanatory features provide for a clean separation of species. So yes, our toy model is that good!

Special thanks to Julia Silge, David Robinson and Andrew Couch for creating and sharing many amazing learning resources for mastering the tidyverse and tidymodels data science packages.