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.
install.packages(c("tidymodels", "glmnet", "rpart", "ranger", "xgboost", "mlbench", "knitr"))
library(tidymodels)
library(glmnet)
library(rpart)
library(ranger)
library(xgboost)
library(mlbench)
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)
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.
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")
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)
We specify our desired yardstick metrics.
eval_metrics <- metric_set(accuracy, f_meas, roc_auc, brier_class)
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)
}
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")
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 |