Palmer penguins - back again

Same peng_bal table and rec_base recipe as Tuesday - today we swap models and add recipe steps.

What you already have

• rec_base - step_zv (drop zero-variance columns) → step_dummy (one-hot categoricals) → step_normalize (scale numerics) on peng_bal
• workflow() + decision_tree / rpart tuned on Tuesday
• vfold_cv(..., strata = y) - 5-fold cross-validation (CV) with both species in each fold
• metric_set(accuracy) for Parts B–D (see metrics refresh; richer metrics in Part E)

Minimal skeleton (copy/paste)

rec <- recipe(y ~ ., data = data)
spec <- rand_forest(mtry = tune(), trees = 500) |> set_engine("ranger") |> set_mode("classification")
metrics <- metric_set(accuracy)
wf <- workflow() |> add_recipe(rec) |> add_model(spec)
split <- initial_split(data, prop = 0.8, strata = y)
folds <- vfold_cv(training(split), v = 5, strata = y)
res <- tune_grid(wf, resamples = folds, grid = 10, metrics = metrics)
wf_final <- finalize_workflow(wf, select_best(res, metric = "accuracy"))
collect_metrics(last_fit(wf_final, split, metrics = metrics))

The honest evaluation workflow

You split your data into training and test sets, then use cross-validation on the training set to tune hyperparameters for each model and compute their average performance. For each model type, you select the best-performing hyperparameter combination based on CV results, and then compare these tuned models to choose the best overall model. Finally, you refit that selected model on the full training data and evaluate it once on the untouched test set to get an unbiased estimate of real-world performance.

• Recipe steps (dummy, scale, impute, PCA, upsample, …) must be fit inside each CV fold on training rows only - otherwise information leaks from held-out data.
• fit_resamples() (next) averages CV metrics for a workflow - quick comparisons when specs are fixed or already tuned.
• Chapter 4 — train/test split walks through last_fit() with a real train/test split.

Confusion matrix counts

Fix a positive class (here: Gentoo, event_level = "second" in yardstick). From a holdout confusion matrix:

Most rates below are built from these four counts at one classification threshold.

Key metrics — what each is good at

Recall = sensitivity (two names, same formula). ROC plots sensitivity vs 1 − specificity across thresholds; AUC summarizes that curve (0.5 ≈ random, 1 = perfect ranker).

ROC curve (Gentoo as positive class)

• x-axis: false positive rate = 1 − specificity (Adelie called Gentoo).
• y-axis: true positive rate = sensitivity / recall (Gentoo caught).
• Each point: a different probability threshold on .pred_Gentoo.
• Diagonal: no better than random ranking; higher curve = better separation of Gentoo scores.

Parts B–D use accuracy for fair model shoot-outs; Part E returns to macro recall, F1, and confusion matrices on imbalanced species.

Part B goal — swap model spec only

Keep preprocessing fixed and learn model-family differences fairly: same recipe, same folds, same metric — change add_model(spec) only.

If we change everything at once, we cannot tell what caused a performance difference.

Model catalog and engines

Same rec_base from Tuesday - change spec (and set_engine()); keep workflow() + folds.

Also available in the same pattern: linear_reg(), svm_rbf(), nearest_neighbor().

• mtry: number of predictors considered at each split (RF only).
• tune() marks hyperparameters for tune_grid().
• Lookup: tidymodels - find models · Engines at a glance.

Parsnip: type → engine → mode

Every spec needs: model type → set_engine() → set_mode("classification") (or regression).

rf_tune_spec <- rand_forest(trees = tune(), mtry = tune()) |>
  set_engine("ranger") |>
  set_mode("classification")
extract_parameter_set_dials(rf_tune_spec)

Same recipe, swap the spec

Part A gave wf_tree (rec_base + tree_spec). Below: random forest (bagging) and XGBoost (boosting) on the same folds.

Afternoon lab — Task 4.1

Lab exercises · Solution · Download Rmd

Bagging and boosting (two ensemble ideas)

Bagging (random forests): many trees on bootstrap samples; each split sees a random subset of predictors → lower variance than one deep tree.
rand_forest(mtry = ..., trees = ..., min_n = ...) |> set_engine("ranger")

Boosting (gradient boosted trees): trees added sequentially; each new tree corrects errors left by the ensemble.
boost_tree(trees = ..., tree_depth = ..., learn_rate = ...) |> set_engine("xgboost") — smaller learn_rate usually needs more trees and gives smoother surfaces.

Forest and boosting specs (code)

rf_spec <- rand_forest(mtry = 3, trees = 400, min_n = 2) |>
  set_engine("ranger", importance = "impurity") |>
  set_mode("classification")
xgb_spec <- boost_tree(trees = 150, tree_depth = 3, learn_rate = 0.05) |>
  set_engine("xgboost") |>
  set_mode("classification")

list(tree = tree_spec, rf = rf_spec, xgb = xgb_spec)

$tree
Decision Tree Model Specification (classification)

Main Arguments:
  tree_depth = 4
  min_n = 10

Computational engine: rpart 


$rf
Random Forest Model Specification (classification)

Main Arguments:
  mtry = 3
  trees = 400
  min_n = 2

Engine-Specific Arguments:
  importance = impurity

Computational engine: ranger 


$xgb
Boosted Tree Model Specification (classification)

Main Arguments:
  trees = 150
  tree_depth = 3
  learn_rate = 0.05

Computational engine: xgboost 

wf_rf <- workflow() |> add_recipe(rec_base) |> add_model(rf_spec)
wf_xgb <- workflow() |> add_recipe(rec_base) |> add_model(xgb_spec)

Cross-validated accuracy (fit three specs)

set.seed(7)
rs_rf <- fit_resamples(wf_rf, folds, metrics = metrics_acc)
rs_xgb <- fit_resamples(wf_xgb, folds, metrics = metrics_acc)

cmp_ens <- bind_rows(
  collect_metrics(rs_tree) |> mutate(model = "Single tree"),
  collect_metrics(rs_rf) |> mutate(model = "Random forest"),
  collect_metrics(rs_xgb) |> mutate(model = "XGBoost")
) |>
  filter(.metric == "accuracy")

Cross-validated accuracy (results)

# A tibble: 3 × 3
  model          mean std_err
  <chr>         <dbl>   <dbl>
1 Random forest 1     0      
2 Single tree   0.996 0.00377
3 XGBoost       0.992 0.00462

Afternoon lab — Task 4.2

Lab exercises · Solution · Download Rmd

Fixed settings here for speed - in projects, wrap each spec in tune_grid() like Day 2 (Tuesday) pipeline.

Boosting: watch the boundary sharpen

We fix tree_depth = 2 and learn_rate = 0.3, then increase trees. Each round adds small trees that correct leftover mistakes - the Gentoo probability surface becomes more nonlinear. We show rounds 1, 3, and 6 (the pattern is similar in between).

boost_tree(trees = 6, tree_depth = 2, learn_rate = 0.3) |>
  set_engine("xgboost") |>
  set_mode("classification")

Boosted Tree Model Specification (classification)

Main Arguments:
  trees = 6
  tree_depth = 2
  learn_rate = 0.3

Computational engine: xgboost 

Boosting round 1

One shallow tree: rough split - mostly underfits the Gentoo cloud.

Boosting round 3

Boosting round 6

Bagging vs boosting (bill plane)

• Bagging (RF): many trees on bootstrap samples; predictions average → often smoother boundaries.
• Boosting (XGBoost): trees added sequentially to fix errors → can look sharper after several rounds.

Neural network (small MLP) - schematic

Neural network (small MLP)

• MLP = multilayer perceptron: one hidden layer of neurons; hidden_units sets its width.
• epochs: how many full passes through the training data.
• mlp(hidden_units = ..., penalty = ..., epochs = ...) |> set_engine("nnet") - weight decay via penalty.
• Requires scaling (already in rec_base via step_normalize).
• On small tabular data, trees often match or beat MLPs - still worth knowing the same workflow interface.

MLP spec and workflow (code)

mlp_spec <- mlp(hidden_units = 8, penalty = 0.1, epochs = 150) |>
  set_engine("nnet", trace = FALSE) |>
  set_mode("classification")

wf_mlp <- workflow() |>
  add_recipe(rec_base) |>
  add_model(mlp_spec)

MLP vs forests (cross-validated accuracy)

set.seed(9)
rs_mlp <- fit_resamples(wf_mlp, folds, metrics = metrics_acc)
mlp_cmp <- bind_rows(
  collect_metrics(rs_rf) |> mutate(model = "Random forest"),
  collect_metrics(rs_xgb) |> mutate(model = "XGBoost"),
  collect_metrics(rs_mlp) |> mutate(model = "MLP (nnet)")
) |>
  filter(.metric == "accuracy")

# A tibble: 3 × 4
  model         .metric   mean std_err
  <chr>         <chr>    <dbl>   <dbl>
1 Random forest accuracy 1     0      
2 XGBoost       accuracy 0.992 0.00462
3 MLP (nnet)    accuracy 1     0      

MLP decision boundary (bill plane)

Logistic regression boundary (bill plane)

Part C goal — PCA inside the recipe

PCA combines correlated numeric predictors into orthogonal principal components — fit step_pca() on training folds only, same as step_normalize().
Compare accuracy with/without PCA on the same folds; tradeoff is interpretability (components are not “bill mm” anymore).

Recipe with step_pca() (code)

rec_pca <- recipe(y ~ ., data = peng_pca) |>
  step_zv(all_predictors()) |>
  step_dummy(all_nominal_predictors()) |>
  step_normalize(all_numeric_predictors()) |>
  step_pca(all_numeric_predictors(), num_comp = 2)

Steps: zero-variance removal → dummies → scale numerics → retain 2 principal components.

Variance explained by each component

PCA biplot - scores vs loadings

• Points = scores: where each penguin sits in PC1-PC2 space (compressed view of all numerics).
• Arrows = loadings (preview): which original measurements push in which direction — detail on the next slide.

PCA loadings — what they mean and how to read them

After step_normalize(), each component is a weighted sum of the scaled columns. Loadings = those weights (from prcomp() / step_pca()).

• Sign: positive loading → birds above average on that measurement tend to have higher scores on that PC.
• Magnitude: larger |loading| → that measurement matters more for that component.
• PC1 loads heavily on bill length, flipper length, body mass (bill depth negatively) → an “overall size” axis.
• PC2 is dominated by year → mostly study year, not morphology here.

Each bar below = one original column’s weight on PC1/PC2 (predictors were scaled before PCA).

Takeaway: PCA merged correlated size measures into PC1; the model later sees PC1/PC2, not “mm of bill” directly — harder to explain biologically than raw columns.

XGBoost with PCA vs without (same folds)

set.seed(8)
rs_pca <- fit_resamples(wf_pca, folds_pca, metrics = metrics_acc)
set.seed(8)
rs_no_pca <- fit_resamples(wf_no_pca, folds_pca, metrics = metrics_acc)

XGBoost with PCA vs without (results)

# A tibble: 2 × 3
  model                    mean std_err
  <chr>                   <dbl>   <dbl>
1 XGBoost + PCA (2 comp.) 0.996 0.00377
2 XGBoost, no PCA         0.989 0.00459

PCA tradeoffs and leakage

• Tradeoff: components are not “bill mm” anymore - VIP on PCA columns is harder to explain biologically.
• Leakage check: always fit PCA inside workflow() + resampling, never on the full table before splitting.

Part D goal — impute inside the recipe

Handle missing values without leakage: add step_impute_*() before step_zv() / dummies / normalize — never drop_na() on the full table before splitting.
Dropping rows can bias who remains; imputing before splits leaks information.

Wrong fix: drop rows with drop_na()

tibble(
  rows_before = nrow(peng_na),
  rows_after_drop_na = nrow(tidyr::drop_na(peng_na))
)

# A tibble: 1 × 2
  rows_before rows_after_drop_na
        <int>              <int>
1         265                238

Dropping loses penguins and can bias who remains - never drop on the full table before splitting.

Right fix: impute inside the recipe

• step_impute_median() for numeric bill_length_mm; step_impute_mode() for categorical sex.
• Then the usual step_zv → step_dummy → step_normalize chain from Part A.

rec_impute <- recipe(y ~ ., data = peng_na) |>
  step_impute_median(bill_length_mm) |>
  step_impute_mode(sex) |>
  step_zv(all_predictors()) |>
  step_dummy(all_nominal_predictors()) |>
  step_normalize(all_numeric_predictors())

Before imputation

Rows with missing bill length sit in the left margin (pink band) - we do not know their x-position yet.

After imputation

• Open circles: true bill length (simulation-only reference).
• Gold points: median-imputed bill length for the same penguins.
• Orange arrows: true value → imputed value shift.
• Gray points: other penguins with observed bill length.

sum(is.na(bake(prep_impute, peng_na)$bill_length_mm))

[1] 0

Imputation inside cross-validation

xgb_na_spec <- boost_tree(trees = 80, tree_depth = 3, learn_rate = 0.08) |>
  set_engine("xgboost") |>
  set_mode("classification")
wf_na <- workflow() |>
  add_recipe(rec_impute) |>
  add_model(xgb_na_spec)
folds_na <- vfold_cv(peng_na, v = 5, strata = y)
set.seed(14)
rs_na <- fit_resamples(wf_na, folds_na, metrics = metrics_acc)
collect_metrics(rs_na) |>
  filter(.metric == "accuracy") |>
  select(mean, std_err)

# A tibble: 1 × 2
   mean std_err
  <dbl>   <dbl>
1 0.996 0.00370

The workflow accepts NA in raw rows; each fold learns imputation on analysis data only.

Metrics and deployment (reminders)

• Pick metrics before peeking at holdout data; report subgroup tables when ethically appropriate.
• Data card: document NA handling and class prevalence before modeling.
• Model card: intended use, limits, monitoring - finish the draft you started Monday.

Continue: Day 2 (Tuesday) pipeline.

Part E goal — metrics beyond accuracy

When one class is rare, accuracy can hide failure on the minority class. Practice macro recall, macro F1, and per-class confusion matrices — and know when step_upsample() helps recall.

We now use three-species classification (Adelie, Gentoo, Chinstrap). For a clean imbalance lesson, we keep only 15 Chinstrap rows (deterministic selection), while Adelie and Gentoo stay abundant.

Why accuracy alone can still mislead

majority_label <- peng_imb3 |>
  dplyr::count(y3, sort = TRUE) |>
  dplyr::slice(1) |>
  dplyr::pull(y3) |>
  as.character()

tibble(
  strategy = paste("Always predict", majority_label),
  accuracy_on_peng_imb3 = mean(as.character(peng_imb3$y3) == majority_label)
)

# A tibble: 1 × 2
  strategy              accuracy_on_peng_imb3
  <chr>                                 <dbl>
1 Always predict Adelie                 0.521

• Overall accuracy can look fine while one class (here Chinstrap) is mostly missed.
• So we track macro recall, macro F1, and per-class recall.

Recipes and model (peng_imb3)

mc_tree_spec <- decision_tree(tree_depth = 4, min_n = 20) |>
  set_engine("rpart") |>
  set_mode("classification")

rec_imb3 <- recipe(y3 ~ ., data = peng_imb3) |>
  step_zv(all_predictors()) |>
  step_dummy(all_nominal_predictors()) |>
  step_normalize(all_numeric_predictors())

rec_imb3_up <- recipe(y3 ~ ., data = peng_imb3) |>
  step_upsample(y3) |>
  step_zv(all_predictors()) |>
  step_dummy(all_nominal_predictors()) |>
  step_normalize(all_numeric_predictors())

wf_imb3 <- workflow() |> add_recipe(rec_imb3) |> add_model(mc_tree_spec)
wf_imb3_up <- workflow() |> add_recipe(rec_imb3_up) |> add_model(mc_tree_spec)

Cross-validated metrics (multiclass)

set.seed(11)
rs_imb3 <- fit_resamples(wf_imb3, folds_imb3, metrics = metrics_cls_imb3)
set.seed(12)
rs_imb3_up <- fit_resamples(wf_imb3_up, folds_imb3, metrics = metrics_cls_imb3)

compare_metrics_tbl_multiclass(
  "No upsample" = rs_imb3,
  "Upsample" = rs_imb3_up
)

Holdout confusion matrices (no upsample vs upsample)

set.seed(13)
split_imb3 <- initial_split(peng_imb3, prop = 0.8, strata = y3)
fit_imb3 <- fit(wf_imb3, training(split_imb3))
pred_imb3 <- augment(fit_imb3, testing(split_imb3))

pred_imb3 |>
  conf_mat(truth = y3, estimate = .pred_class)

           Truth
Prediction  Adelie Gentoo Chinstrap
  Adelie        30      0         3
  Gentoo         0     23         0
  Chinstrap      0      1         0

fit_imb3_up <- fit(wf_imb3_up, training(split_imb3))
pred_imb3_up <- augment(fit_imb3_up, testing(split_imb3))

pred_imb3_up |>
  conf_mat(truth = y3, estimate = .pred_class)

           Truth
Prediction  Adelie Gentoo Chinstrap
  Adelie        30      0         2
  Gentoo         0     23         0
  Chinstrap      0      1         1

Per-class recall (holdout)

• In the confusion matrix: rows = true species, columns = predicted species.
• Upsampling is not a free lunch, but it often improves minority-class recall on this engineered imbalance.
• Use macro/per-class metrics when one class is underrepresented.

Part F goal — VIP and local explanations

Interpret fitted models responsibly: VIP ranks influential predictors globally; SHAP explains individual predictions. Neither is a causal claim.

Data the forest sees (bill plane)

Variable importance (VIP)

• VIP = how often splits use each baked predictor (impurity in ranger).
• Correlated bill measures can share importance — not proof that changing one bill mm causes species shift.
• Permutation importance (shuffle one column, watch metric drop) is model-agnostic — see the SHAP notebook.

Afternoon lab — Task 4.3

Lab exercises · Solution · Download Rmd

SHAP - what it means (game theory intuition)

• One penguin, one prediction: the model outputs P(male) (or log-odds). SHAP asks: how much did each measurement move this bird away from a baseline?
• Players and coalitions: treat each feature as a player. A coalition is any subset of features with the others averaged over the data. The model’s output on that coalition is the coalition’s payout.
• Shapley value: for each feature, average its marginal contribution over all orders in which features can join the coalition. That fair share is the SHAP value for feature j (written phi_j).
• Additivity: start from a baseline score, then add one credit per feature. For penguin i:

baseline + phi_i1 + phi_i2 + ...  ≈  P(male) for this bird

(all credits are for the same fitted model on the link/probability scale)

• VIP vs SHAP: VIP ranks predictors globally (how often splits use them). SHAP explains one row (why this penguin scored more male or female).

SHAP task — sex without species

• Outcome: sex (female / male). species is omitted on purpose so the task is not trivial (see sex notebook).
• Model: same recipe + random forest (ranger) pattern as Tuesday; fit on all complete rows here so SHAP curves are stable in class (in production, fit SHAP on training data only).
• Packages: shapviz for plots; kernelshap to compute values on the fitted ranger forest.

pg_sex <- prep_penguins_sex()

rec_sex <- recipe(sex ~ bill_length_mm + bill_depth_mm + flipper_length_mm +
                    body_mass_g + island + year, data = pg_sex) |>
  step_zv(all_predictors()) |>
  step_dummy(all_nominal_predictors()) |>
  step_normalize(all_numeric_predictors())

rf_sex_spec <- rand_forest(trees = 300, mtry = 3, min_n = 2) |>
  set_engine("ranger", probability = TRUE) |>
  set_mode("classification")

wf_sex <- workflow() |>
  add_recipe(rec_sex) |>
  add_model(rf_sex_spec)

fit_sex <- fit(wf_sex, pg_sex)

X_model — bake the forest’s spreadsheet

X_model = bake(rec_prepped, new_data = pg_sex, all_predictors())

• Rows: one per penguin (same birds as pg_sex, after recipe row handling).
• Columns: predictors only — numeric columns ranger actually split on. Outcome sex is omitted.
• kernelshap() edits rows of this object; X_shap = rows we explain; bg_X = background reference flock.

Wrong: pass raw bill_length_mm = 45 while the forest expects a z-scored value. Right: pass the baked row from X_model.

rec_prepped <- prep(extract_recipe(fit_sex, estimated = TRUE), training = pg_sex)
X_model <- bake(rec_prepped, new_data = pg_sex, all_predictors())

head(X_model, 3)

KernelSHAP — loop and code

For one penguin in X_shap: start from a baseline built from bg_X; perturb each baked feature; call pred_fun → P(male); average moves into one SHAP value per feature (they sum ≈ this bird’s final score).

• rf_engine — unwrap ranger from tidymodels via extract_fit_parsnip()
• pred_fun — baked matrix → P(male) (SHAP does not know about recipes)
• X_shap — ~80 random rows to explain; bg_X — ~35 background rows for coalitions
• kernelshap(...) then shapviz(...) for beeswarm / waterfall plots

rf_engine <- extract_fit_parsnip(fit_sex)$fit

pred_fun <- function(object, X_new) {
  as.numeric(predict(object, data = X_new)$predictions[, "male"])
}

set.seed(11)
X_shap <- dplyr::slice_sample(X_model, n = 80)
bg_X <- dplyr::slice_sample(X_model, n = 35)

ks <- kernelshap(rf_engine, X = X_shap, pred_fun = pred_fun, bg_X = bg_X)
shp_sex <- shapviz(ks, X_pred = X_shap)

Afternoon lab — Task 4.4

Lab exercises · Solution · Download Rmd

SHAP beeswarm - global pattern

• y-axis: features; x-axis: SHAP contribution toward P(male).
• Color: feature value (high vs low). Wide spread = this predictor matters for many birds in this model.

SHAP waterfall — two local explanations

Male and female examples from the same fitted model — bars sum from baseline to P(male) for each bird.

Reading SHAP without over-claiming

Reasonable claims

• For this penguin, higher body mass (in the model, holding the coalition story fixed) pushed P(male) up or down.
• Globally, bill depth and mass show the widest SHAP spread in this fitted RF.
• SHAP describes the trained model, not a randomized experiment in the wild.

Avoid

• “Bill length causes males” - observational data; morphometrics are correlated (data card).
• “We should intervene on flipper length” - policy from association, not a causal effect shown by the model.

Humility prompt (assessment): What would you not claim from this SHAP beeswarm?

Full write-up: SHAP notebook

Afternoon lab — Task 4.5

Lab exercises · Solution · Download Rmd

Nested resampling and end-of-day checklist

• Nested CV: an outer loop estimates performance; an inner loop tunes hyperparameters inside each outer training fold — avoid tuning and scoring on the same held-out rows.

Check yourself

• Name three parsnip model types and a valid engine for each
• Why must step_pca() and step_impute_*() live inside the recipe?
• Which metric would you report for a rare-event screening study?
• When might step_upsample() not improve CV scores?
• One reason VIP on correlated bills is not a causal claim
• One sentence you would not say about a SHAP plot on penguin sex