Elements of Data Science
SDS 322E

H. Sherry Zhang
Department of Statistics and Data Sciences
University of Texas at Austin

Get the repository for today:
library(usethis)
create_from_github("SDS322E-26FALL/1302-compare-models", fork = FALSE)

Learning objectives

• Compare several machine learning models for the same question fairly, all in tidymodels
• Choose the pre-processing that each model needs
• Choose a model with cross-validation, and evaluate it once on the testing data
• Use this as a template for Project 2

The question

Which model predicts penguin sex best, from bill_len and bill_dep?

We have seen three models for a binary response:

• Logistic regression: logistic_reg(), no hyperparameters
• KNN: nearest_neighbor(), tune neighbors
• Random forest: rand_forest(), tune mtry and min_n

Most of the code is what you have already seen in each lecture. What’s new is putting them together fairly:

• the same training-testing split and the same cross-validation folds for all models
• a different recipe for each model
• compare the models on the cross-validation results, and use the testing data once, for the chosen model

The plan

• Step 0: Training-testing split, shared by all models
• Steps 1-3, for each model: the model, the recipe (depends on the model), the workflow
• Step 4: Cross-validation, for all models:
  • 4a: the cross-validation folds, shared by all models
  • 4b: the hyperparameter grid for each model
  • 4c: fit each model on the folds
  • 4d: average the metrics over the folds
  • 4e: the best version of each model
• Step 4.5 (new): Compare the models on their cross-validation results, and choose one
• Step 5/6: Refit the chosen model on the training data, and evaluate it once on the testing data

Step 0: split once, for all models

penguins_clean <- as_tibble(datasets::penguins) |> na.omit() |> select(-year)

set.seed(5)
penguins_split <- initial_split(penguins_clean, prop = 0.75, strata = sex)
penguins_train <- training(penguins_split)
penguins_test  <- testing(penguins_split)

• strata = sex: the training and the testing data keep the same proportion of females and males
  • The penguins are almost an equal split (49.5% female, 50.5% male), so a random split is already close to balanced: we didn’t need strata before
  • With an uneven response (e.g. 10% vs 90%), use strata: otherwise some folds may have very few of the rare class
• All three models use the same split: they are trained on the same penguins, and the chosen one is tested on the same penguins

Step 1: specify the models

log_spec <- logistic_reg() |> set_engine("glm")

knn_spec <- nearest_neighbor(mode = "classification", neighbors = tune()) |> set_engine("kknn")

rf_spec <- rand_forest(mode = "classification", trees = 500, mtry = tune(), min_n = tune()) |> set_engine("ranger")

• Same as in each lecture, all in classification mode: the response sex is binary
• Logistic regression has nothing to tune; KNN and random forest have hyperparameters left as tune()

Step 2: a recipe for each model

Different models need different pre-processing:

log_rec <- recipe(sex ~ bill_len + bill_dep, data = penguins_train)

knn_rec <- recipe(sex ~ bill_len + bill_dep, data = penguins_train) |>
  step_normalize(all_numeric_predictors())

rf_rec <- recipe(sex ~ bill_len + bill_dep, data = penguins_train)

• Logistic regression and random forest: nothing to do, both predictors are numeric. A tree split doesn’t care about the scale of a predictor
• KNN: step_normalize(), because distances need all predictors on the same scale
• With categorical predictors (e.g. species), logistic regression and KNN would also need step_dummy(); the random forest can use them directly

Step 2: a recipe for each model (new)

What does the KNN recipe do? Use prep() and bake() to look at the processed training data:

knn_rec |> prep() |> bake(new_data = NULL) |> head(5)

# A tibble: 5 × 3
  bill_len bill_dep sex   
     <dbl>    <dbl> <fct> 
1   -0.748    0.434 female
2   -1.41     1.08  female
3   -1.01     0.335 female
4   -0.601    0.237 female
5   -1.43     0.335 female

• bill_len and bill_dep now have mean 0 and standard deviation 1
• You don’t need prep() and bake() to fit a model, the workflow does it for you. But they are useful to check what your recipe does

Step 3: a workflow for each model

log_wf <- workflow() |> add_model(log_spec) |> add_recipe(log_rec)
knn_wf <- workflow() |> add_model(knn_spec) |> add_recipe(knn_rec)
rf_wf  <- workflow() |> add_model(rf_spec)  |> add_recipe(rf_rec)

Each workflow keeps a model with its own recipe, so the right pre-processing always goes with the right model.

Step 4a: folds, shared by all models

set.seed(5)
penguins_folds <- vfold_cv(penguins_train, v = 10, strata = sex)
penguins_folds

#  10-fold cross-validation using stratification 
# A tibble: 10 × 2
   splits           id    
   <list>           <chr> 
 1 <split [223/26]> Fold01
 2 <split [223/26]> Fold02
 3 <split [223/26]> Fold03
 4 <split [224/25]> Fold04
 5 <split [224/25]> Fold05
 6 <split [224/25]> Fold06
 7 <split [225/24]> Fold07
 8 <split [225/24]> Fold08
 9 <split [225/24]> Fold09
10 <split [225/24]> Fold10

• All three models use the same penguins_folds: each model is validated on exactly the same penguins, so their cross-validation results can be compared directly
• strata = sex again: every fold keeps the same proportion of females and males

Step 4b: hyperparameter grids

knn_grid <- grid_regular(neighbors(range = c(1, 30)), levels = 30)
rf_grid <- grid_regular(mtry(range = c(1, 2)), min_n(range = c(2, 20)), levels = c(2, 4))
rf_grid

# A tibble: 8 × 2
   mtry min_n
  <int> <int>
1     1     2
2     2     2
3     1     8
4     2     8
5     1    14
6     2    14
7     1    20
8     2    20

• Logistic regression: no hyperparameters, no grid
• KNN: 30 values of \(k\), from 1 to 30
• Random forest: with only two predictors, mtry can only be 1 or 2: 2 \(\times\) 4 values of min_n = 8 combinations

Step 4c: fit each model on the folds

Choose the metrics with metric_set() (new), and use the same ones for all models:

cls_metrics <- metric_set(accuracy, roc_auc)

log_res <- fit_resamples(log_wf, resamples = penguins_folds, metrics = cls_metrics,
                         control = control_resamples(save_pred = TRUE))
knn_res <- tune_grid(knn_wf, resamples = penguins_folds, grid = knn_grid, metrics = cls_metrics,
                     control = control_grid(save_pred = TRUE))
set.seed(5)
rf_res  <- tune_grid(rf_wf, resamples = penguins_folds, grid = rf_grid, metrics = cls_metrics,
                     control = control_grid(save_pred = TRUE))

• fit_resamples() (new): like tune_grid(), but for a model with nothing to tune: it fits the workflow on each fold and computes the validation metrics
• save_pred = TRUE: keep the predictions on each validation fold, so we can draw ROC curves later
• How many models? Logistic regression: 10 fits. KNN: 30 \(\times\) 10 = 300 fits. Random forest: 8 \(\times\) 10 = 80 forests of 500 trees, the slowest to tune

Step 4d: average over the folds

collect_metrics(log_res)

# A tibble: 2 × 6
  .metric  .estimator  mean     n std_err .config        
  <chr>    <chr>      <dbl> <int>   <dbl> <chr>          
1 accuracy binary     0.773    10  0.0292 pre0_mod0_post0
2 roc_auc  binary     0.848    10  0.0294 pre0_mod0_post0

collect_metrics(knn_res) |> filter(.metric == "accuracy") |> head(3)

# A tibble: 3 × 7
  neighbors .metric  .estimator  mean     n std_err .config         
      <int> <chr>    <chr>      <dbl> <int>   <dbl> <chr>           
1         1 accuracy binary     0.786    10  0.0196 pre0_mod01_post0
2         2 accuracy binary     0.786    10  0.0196 pre0_mod02_post0
3         3 accuracy binary     0.786    10  0.0196 pre0_mod03_post0

• mean: the average over the n = 10 folds; std_err: how much it varies from fold to fold. We will need both to compare models
• One row per metric for logistic regression; one row per metric and per value of \(k\) for KNN

Step 4d: average over the folds

autoplot(knn_res, metric = "accuracy")

Step 4e: best version of each model

show_best(knn_res, metric = "accuracy", n = 3)

# A tibble: 3 × 7
  neighbors .metric  .estimator  mean     n std_err .config         
      <int> <chr>    <chr>      <dbl> <int>   <dbl> <chr>           
1        12 accuracy binary     0.852    10  0.0221 pre0_mod12_post0
2        11 accuracy binary     0.848    10  0.0255 pre0_mod11_post0
3        13 accuracy binary     0.848    10  0.0227 pre0_mod13_post0

show_best(rf_res, metric = "accuracy", n = 3)

# A tibble: 3 × 8
   mtry min_n .metric  .estimator  mean     n std_err .config        
  <int> <int> <chr>    <chr>      <dbl> <int>   <dbl> <chr>          
1     1     2 accuracy binary     0.828    10  0.0292 pre0_mod1_post0
2     2     2 accuracy binary     0.828    10  0.0225 pre0_mod5_post0
3     1    20 accuracy binary     0.824    10  0.0241 pre0_mod4_post0

• For KNN and random forest, keep the best hyperparameters: this is the version of each model we compare
• Logistic regression has only one version

Step 4.5: compare the models

Take the best version of each model, and compare their cross-validation accuracy:

cv_compare <- bind_rows(
  collect_metrics(log_res) |> mutate(model = "Logistic regression"),
  show_best(knn_res, metric = "accuracy", n = 1) |> mutate(model = "KNN"),
  show_best(rf_res, metric = "accuracy", n = 1) |> mutate(model = "Random forest")
) |>
  filter(.metric == "accuracy") |>
  select(model, mean, std_err, n)
cv_compare

# A tibble: 3 × 4
  model                mean std_err     n
  <chr>               <dbl>   <dbl> <int>
1 Logistic regression 0.773  0.0292    10
2 KNN                 0.852  0.0221    10
3 Random forest       0.828  0.0292    10

cv_compare |>
  ggplot(aes(x = mean, y = fct_reorder(model, mean))) +
  geom_pointrange(aes(xmin = mean - std_err, xmax = mean + std_err)) +
  labs(x = "Cross-validation accuracy, mean ± 1 standard error", y = NULL)

Step 4.5: compare the ROC curves

Accuracy depends on the 0.5 cutoff. The ROC curve compares the models over all cutoffs, using the predictions on the validation folds:

• The closer a curve gets to the top-left corner, the better the model separates females and males (male = positive class)

• The same conclusion as accuracy: KNN \(\approx\) random forest, both clearly better than logistic regression

Why does logistic regression lose?

• Within each species, males have bigger bills. But the species differ so much that a male Adelie and a female Chinstrap can have similar bills
• So the boundary between female and male is not one straight line: it is in a different place for each species. Logistic regression can only draw one straight line; KNN and random forest can draw a flexible boundary

Step 5/6: evaluate the chosen model once

Plug the best \(k\) into the KNN workflow, refit on all the training data, and evaluate once on the testing data. last_fit() does the refit and the evaluation in one step:

best_k <- select_best(knn_res, metric = "accuracy")
knn_final <- last_fit(finalize_workflow(knn_wf, best_k), split = penguins_split, metrics = cls_metrics)
collect_metrics(knn_final)

# A tibble: 2 × 4
  .metric  .estimator .estimate .config        
  <chr>    <chr>          <dbl> <chr>          
1 accuracy binary         0.857 pre0_mod0_post0
2 roc_auc  binary         0.921 pre0_mod0_post0

• This is the number to report: how well the chosen model, KNN with \(k\) = 12, predicts new penguins

What about a regression problem?

If your response is numeric (e.g. predict body_mass), the same steps apply, with a few changes:

• Models: linear_reg() instead of logistic_reg(), and mode = "regression" for KNN and random forest

• Split: strata is optional: for a numeric response, strata = body_mass splits it into bins so every fold covers light and heavy penguins

• Metrics: RMSE (in the units of the response) and \(R^2\), with metric_set(rmse, rsq): for RMSE, lower is better. Use select_best(metric = "rmse")

• Recipes: the same rules. Dummy variables for categorical predictors in linear regression and KNN, normalization for KNN, nothing for the random forest

• Don’t assume the flexible model wins: if the response is close to a straight-line function of the predictors (as body mass is of flipper length), linear regression can do as well as KNN and random forest, and is much easier to explain. Compare, and let the cross-validation results decide

For Project 2

Use this lecture as a template:

1. The question: what is one row, what is the response, and which predictors (no leakage!)?
2. Step 0: one training-testing split and one set of folds, shared by all models
3. Steps 1-4, for each model: a model, a recipe that suits it, a workflow, and cross-validation
4. Step 4.5: compare the models on the cross-validation metric, with its standard error: are the differences big enough to matter? Then choose a model, and explain why (performance, simplicity, interpretability)
5. Step 5/6: evaluate the chosen model once on the testing data, and report that number