1 Overview

Following our previous practical session, we noticed discrepancies between train and test performance (overfitting). In this session, we will use k-fold Cross-Validation (CV) to get an averaged, more realistic performance metric over the folds on the training data. We will also expand our toolkit to include Decision Trees, Random Forests, and XGBoost.

2 Install packages (if you haven’t done so)

install.packages(c("tidymodels", "glmnet", "rpart", "ranger", "xgboost", "mlbench", "knitr"))

3 Load Packages

library(tidymodels)
library(glmnet)
library(rpart)
library(ranger)
library(xgboost)
library(mlbench)

4 Load and Prepare the Data

We load the dataset and prepare it by releveling the outcome variable and removing missing values.

data("PimaIndiansDiabetes2")
pima <- PimaIndiansDiabetes2

Relevel so the positive class (“pos”) is the first level for yardstick metrics.

pima$diabetes <- relevel(pima$diabetes, ref = "pos")

Remove observations with NA (missing data). (Note: this reduces N from 768 to 392. In practice, imputation might be preferred).

pima <- na.omit(pima)

5 Data Splitting & Resampling

We split the data into training (70%) and testing (30%) sets, stratified by the diabetes outcome to maintain class proportions.

set.seed(123)
pima_split <- initial_split(pima, prop = 0.70, strata = diabetes)
pima_train <- training(pima_split)

We create 5-fold cross-validation folds from the training data.

set.seed(234)
cv_folds <- vfold_cv(pima_train, v = 10, strata = diabetes)

Create a recipe

pima_rec <- recipe(diabetes ~ ., data = pima_train) |> 
  step_normalize(all_numeric_predictors()) 

Normalizing numeric predictors is crucial for engines like glmnet.

6 Model Specifications

In tidymodels, we define the “blueprint” of the model before passing any data. Common elements used across these blueprints include:

  • tune(): A placeholder indicating that the best hyperparameter value will be found later during cross-validation.
  • set_engine(): Specifies the underlying R package to do the math (e.g., glmnet, xgboost).
  • set_mode("classification"): Specifies we are predicting a categorical outcome (“pos” or “neg”).

Standard Logistic Regression (Baseline)

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

Penalized Logistic Regression (Elastic Net)

  • penalty: Controls the total amount of shrinkage (regularization) applied to prevent overfitting.
  • mixture: Balances between Ridge regression (shrinks all coefficients) and Lasso regression (forces some coefficients to exactly 0).
penal_spec <- logistic_reg(penalty = tune(), mixture = tune()) |> 
  set_engine("glmnet") |> 
  set_mode("classification")

Decision Tree

  • cost_complexity: The pruning parameter. Penalizes the tree for getting too complicated, preventing it from fitting to noise.
  • tree_depth: Limits how many levels deep the single tree is allowed to grow.
dt_spec    <- decision_tree(cost_complexity = tune(), tree_depth = tune()) |> 
  set_engine("rpart") |> 
  set_mode("classification")

Random Forest

  • mtry: The number of random predictors (columns) the model is allowed to consider at each split.
  • min_n: The minimum number of data points required in a node before it is allowed to split further.
  • trees = 500: Hardcodes the model to build exactly 500 trees in the ensemble (a solid default).
rf_spec    <- rand_forest(mtry = tune(), min_n = tune(), trees = 500) |> 
  set_engine("ranger") |> 
  set_mode("classification")

XGBoost (Gradient Boosted Trees)

  • tree_depth: Limits how deep each individual boosted tree can grow.
  • learn_rate: Shrinkage. Controls how much each sequential tree contributes to the final answer.
  • trees = 500: Builds 500 trees sequentially, where each new tree tries to fix the mistakes of the previous one.
xgb_spec   <- boost_tree(tree_depth = tune(), learn_rate = tune(), trees = 500) |> 
  set_engine("xgboost") |> 
  set_mode("classification")

7 Workflows

Next, we bundle the preprocessing recipe and model specifications into workflows.

log_wf   <- workflow() |> add_recipe(pima_rec) |> add_model(log_spec)
penal_wf <- workflow() |> add_recipe(pima_rec) |> add_model(penal_spec)
dt_wf    <- workflow() |> add_recipe(pima_rec) |> add_model(dt_spec)
rf_wf    <- workflow() |> add_recipe(pima_rec) |> add_model(rf_spec)
xgb_wf   <- workflow() |> add_recipe(pima_rec) |> add_model(xgb_spec)

8 Define Evaluation Metrics

We specify our desired yardstick metrics.

eval_metrics <- metric_set(accuracy, f_meas, roc_auc, brier_class)

9 Master Helper Function

To streamline tuning, CV evaluation, and test set evaluation across multiple models without repeating code, we define a helper function evaluate_model().

evaluate_model <- function(workflow, model_name, requires_tuning = TRUE, seed = 345) {
  set.seed(seed)
  
  # Step A: Tune and Finalize Workflow (if required)
  if (requires_tuning) {
    tune_res <- tune_grid(
      workflow,
      resamples = cv_folds,
      grid = 20,                      # Fixed: Matched to manual code
      metrics = metric_set(roc_auc)   # Fixed: Use only AUC for tuning to prevent f_meas warnings
    )
    best_params <- select_best(tune_res, metric = "roc_auc")
    final_wf <- finalize_workflow(workflow, best_params)
  } else {
    final_wf <- workflow
  }
  
  # Step B: Calculate 5-Fold CV Performance (on Training Data)
  set.seed(seed)
  cv_resamples <- fit_resamples(
    final_wf,
    resamples = cv_folds,
    metrics = eval_metrics            # Evaluate all metrics on the final best model
  )
  
  cv_metrics <- collect_metrics(cv_resamples) |> 
    dplyr::mutate(model = paste0(model_name, " (CV Train)")) |> 
    dplyr::select(model, .metric, mean) |> 
    dplyr::rename(.estimate = mean)
  
  # Step C: Calculate Final Test Set Performance
  set.seed(seed)
  test_fit <- last_fit(
    final_wf,
    split = pima_split,
    metrics = eval_metrics            # Evaluate all metrics on the test set
  )
  
  test_metrics <- collect_metrics(test_fit) |> 
    dplyr::mutate(model = paste0(model_name, " (Test)")) |> 
    dplyr::select(model, .metric, .estimate)
  
  # Step D: Combine and return
  dplyr::bind_rows(cv_metrics, test_metrics)
}

10 Run Evaluations

We now run our evaluation pipeline on all models. (Note: Standard logistic regression does not have hyperparameters to tune, so we bypass tuning.)

log_results   <- evaluate_model(log_wf, "Standard Logistic", requires_tuning = FALSE)
penal_results <- evaluate_model(penal_wf, "Penalized Logistic")
dt_results    <- evaluate_model(dt_wf, "Decision Tree")
rf_results    <- evaluate_model(rf_wf, "Random Forest")
#> i Creating pre-processing data to finalize 1 unknown parameter: "mtry"
xgb_results   <- evaluate_model(xgb_wf, "XGBoost")

11 Compare Results

Finally, we bring all the results together to compare the CV metrics against the completely unseen Test Set metrics to assess generalization.

final_comparison <- dplyr::bind_rows(
  log_results, penal_results, dt_results, rf_results, xgb_results
) |> 
  # Pivot to match ml-classification1.R format
  tidyr::pivot_wider(names_from = .metric, values_from = .estimate) |>
  # Reorder columns explicitly to match the first practical
  dplyr::select(model, accuracy, f_meas, roc_auc, brier_class)

# Print the final output
final_comparison |> 
  knitr::kable(
    col.names = c("Model", "Accuracy", "F1 Score", "AUC", "Brier Score"),
    digits = 3, 
    align = "lcccc"
  )
Model Accuracy F1 Score AUC Brier Score
Standard Logistic (CV Train) 0.815 0.699 0.860 0.142
Standard Logistic (Test) 0.729 0.543 0.821 0.169
Penalized Logistic (CV Train) 0.815 0.693 0.862 0.141
Penalized Logistic (Test) 0.729 0.543 0.821 0.168
Decision Tree (CV Train) 0.788 0.703 0.841 0.146
Decision Tree (Test) 0.797 0.714 0.812 0.147
Random Forest (CV Train) 0.789 0.661 0.861 0.144
Random Forest (Test) 0.754 0.592 0.821 0.157
XGBoost (CV Train) 0.789 0.647 0.871 0.140
XGBoost (Test) 0.780 0.629 0.814 0.156