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", "DT"))
install.packages("remotes")
# For Windows, install RTools from https://cran.r-project.org/bin/windows/Rtools/
install.packages("parameters") # for building vignettes
remotes::install_github("wnarifin/MJMSTable", build_vignettes = TRUE) # for easy descriptive statisticslibrary(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 <- PimaIndiansDiabetes2Relevel 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 10-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(). Updated: This function now returns
the aggregated metrics, raw test predictions (for CIs), and the best
tuned hyperparameters.
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,
metrics = metric_set(roc_auc)
)
best_params <- select_best(tune_res, metric = "roc_auc")
final_wf <- finalize_workflow(workflow, best_params)
} else {
final_wf <- workflow
best_params <- NULL
}
# Step B: Calculate CV Performance (on Training Data)
set.seed(seed)
cv_resamples <- fit_resamples(
final_wf,
resamples = cv_folds,
metrics = eval_metrics
)
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 & Save Predictions
set.seed(seed)
test_fit <- last_fit(
final_wf,
split = pima_split,
metrics = eval_metrics
)
test_metrics <- collect_metrics(test_fit) |>
dplyr::mutate(model = paste0(model_name, " (Test)")) |>
dplyr::select(model, .metric, .estimate)
test_preds <- collect_predictions(test_fit) |>
dplyr::mutate(model = model_name)
# Step D: Combine metrics and return as a list
combined_metrics <- dplyr::bind_rows(cv_metrics, test_metrics)
return(list(
metrics = combined_metrics,
test_preds = test_preds,
best_params = best_params
))
}We now run our evaluation pipeline on all models.
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$metrics,
penal_results$metrics,
dt_results$metrics,
rf_results$metrics,
xgb_results$metrics
) |>
# 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 as an interactive table
DT::datatable(
final_comparison,
colnames = c("Model", "Accuracy", "F1 Score", "AUC", "Brier Score"),
rownames = FALSE,
options = list(
pageLength = 10, # Shows 10 rows (perfect for Train & Test for 5 models)
autoWidth = TRUE,
dom = 't', # Hides the search bar and pagination for a cleaner look
columnDefs = list(
list(className = 'dt-center', targets = 1:4) # Centers the metric columns
)
)
) |>
DT::formatRound(columns = 2:5, digits = 3) # Rounds the 4 metric columns to 3 decimal placesdiabetes, a binary
classification target indicating the presence (“pos”) or absence (“neg”)
of diabetes.MJMSTable::descriptive_tbl(
pima,
group_var = "diabetes",
included_var = names(pima[,-9]) # diabetes is the 9th variable in the dataset
)| Variables | pos ( n = 130) | neg ( n = 262) | Total n (%) |
|---|---|---|---|
| pregnant, Mean (SD)a | 4.5 (3.9) | 2.7 (2.6) | 3.3 (3.2) |
| glucose, Mean (SD)a | 145.2 (29.8) | 111.4 (24.6) | 122.6 (30.9) |
| pressure, Mean (SD)a | 74.1 (13.0) | 69.0 (11.9) | 70.7 (12.5) |
| triceps, Mean (SD)a | 33.0 (9.6) | 27.3 (10.4) | 29.1 (10.5) |
| insulin, Mean (SD)a | 206.8 (132.7) | 130.9 (102.6) | 156.1 (118.8) |
| mass, Mean (SD)a | 35.8 (6.7) | 31.8 (6.8) | 33.1 (7.0) |
| pedigree, Mean (SD)a | 0.6 (0.4) | 0.5 (0.3) | 0.5 (0.3) |
| age, Mean (SD)a | 35.9 (10.6) | 28.3 (9.0) | 30.9 (10.2) |
| a Mean (SD) | |||
We evaluated five distinct machine learning algorithms. Where applicable, hyperparameter tuning was conducted using a grid search (testing 20 pseudo-random combinations) optimized for the Area Under the Receiver Operating Characteristic Curve (ROC AUC).
Table: Algorithms and Hyperparameter Settings
hyper_df <- data.frame(
Algorithm = c(
"Standard Logistic Regression",
"Penalized Logistic (Elastic Net)",
"Decision Tree",
"Random Forest",
"XGBoost"
),
Engine = c("glm", "glmnet", "rpart", "ranger", "xgboost"),
`Tuned Hyperparameters (Used Values)` = c(
"None (Baseline)",
sprintf("penalty = %.4f, mixture = %.4f", penal_results$best_params$penalty, penal_results$best_params$mixture),
sprintf("cost_complexity = %.4f, tree_depth = %d", dt_results$best_params$cost_complexity, dt_results$best_params$tree_depth),
sprintf("mtry = %d, min_n = %d", rf_results$best_params$mtry, rf_results$best_params$min_n),
sprintf("tree_depth = %d, learn_rate = %.4f", xgb_results$best_params$tree_depth, xgb_results$best_params$learn_rate)
),
`Fixed Hyperparameters` = c("None", "None", "None", "trees = 500", "trees = 500"),
check.names = FALSE
)
knitr::kable(hyper_df, align = "llll")| Algorithm | Engine | Tuned Hyperparameters (Used Values) | Fixed Hyperparameters |
|---|---|---|---|
| Standard Logistic Regression | glm | None (Baseline) | None |
| Penalized Logistic (Elastic Net) | glmnet | penalty = 0.0078, mixture = 0.1000 | None |
| Decision Tree | rpart | cost_complexity = 0.0000, tree_depth = 3 | None |
| Random Forest | ranger | mtry = 2, min_n = 18 | trees = 500 |
| XGBoost | xgboost | tree_depth = 1, learn_rate = 0.0280 | trees = 500 |
To fulfill the TRIPOD+AI reporting guidelines for model evaluation, we present the final unseen test set performance metrics alongside their 95% confidence intervals (CIs). We use bootstrapping (1,000 resamples) on the test set predictions to calculate these CIs robustly.
# 1. Define the specific TRIPOD+AI metric set (added precision and recall)
tripod_metrics <- metric_set(accuracy, precision, recall, f_meas, roc_auc, brier_class)
# 2. Combine the raw test set predictions extracted from our evaluation function
all_test_preds <- dplyr::bind_rows(
log_results$test_preds,
penal_results$test_preds,
dt_results$test_preds,
rf_results$test_preds,
xgb_results$test_preds
)
# 3. Define a bootstrap helper function for CIs
calc_tripod_ci <- function(preds_data, seed = 123) {
set.seed(seed)
# Create 1,000 bootstrap resamples of the test set predictions
boots <- rsample::bootstraps(preds_data, times = 1000)
# Calculate metrics for each resample
boot_metrics <- boots |>
dplyr::mutate(
metrics = purrr::map(splits, function(split) {
split_data <- rsample::analysis(split)
# Suppress warnings for edge-case resamples that lack certain classes
suppressWarnings(
tripod_metrics(split_data, truth = diabetes, estimate = .pred_class, .pred_pos)
)
})
) |>
tidyr::unnest(metrics)
# Summarize the distribution to find the 2.5th and 97.5th percentiles (95% CI)
boot_metrics |>
dplyr::group_by(.metric) |>
dplyr::summarize(
estimate = mean(.estimate, na.rm = TRUE),
lower_ci = quantile(.estimate, 0.025, na.rm = TRUE),
upper_ci = quantile(.estimate, 0.975, na.rm = TRUE),
.groups = "drop"
)
}
# 4. Apply the function to each model's test predictions
tripod_results <- all_test_preds |>
dplyr::group_by(model) |>
tidyr::nest() |>
dplyr::mutate(ci_data = purrr::map(data, calc_tripod_ci)) |>
dplyr::select(model, ci_data) |>
tidyr::unnest(ci_data)
# 5. Format and print the final table
tripod_table <- tripod_results |>
dplyr::mutate(
# Format exactly as Estimate (Lower CI, Upper CI)
formatted = sprintf("%.3f (%.3f, %.3f)", estimate, lower_ci, upper_ci)
) |>
dplyr::select(model, .metric, formatted) |>
tidyr::pivot_wider(names_from = .metric, values_from = formatted) |>
# Reorder columns logically
dplyr::select(model, accuracy, precision, recall, f_meas, roc_auc, brier_class)
tripod_table |>
knitr::kable(
col.names = c("Model", "Accuracy", "Precision", "Recall", "F1 Score", "AUC", "Brier Score"),
align = "lcccccc"
)| Model | Accuracy | Precision | Recall | F1 Score | AUC | Brier Score |
|---|---|---|---|---|---|---|
| Standard Logistic | 0.728 (0.644, 0.814) | 0.609 (0.423, 0.783) | 0.487 (0.324, 0.639) | 0.538 (0.381, 0.676) | 0.820 (0.739, 0.894) | 0.169 (0.125, 0.216) |
| Penalized Logistic | 0.728 (0.644, 0.814) | 0.609 (0.423, 0.783) | 0.487 (0.324, 0.639) | 0.538 (0.381, 0.676) | 0.820 (0.741, 0.893) | 0.169 (0.125, 0.214) |
| Decision Tree | 0.796 (0.720, 0.864) | 0.663 (0.513, 0.796) | 0.773 (0.632, 0.900) | 0.711 (0.590, 0.819) | 0.814 (0.723, 0.894) | 0.147 (0.104, 0.191) |
| Random Forest | 0.754 (0.669, 0.831) | 0.652 (0.481, 0.815) | 0.538 (0.375, 0.700) | 0.586 (0.440, 0.727) | 0.821 (0.734, 0.899) | 0.157 (0.120, 0.193) |
| XGBoost | 0.780 (0.703, 0.856) | 0.708 (0.531, 0.867) | 0.565 (0.395, 0.725) | 0.624 (0.476, 0.761) | 0.815 (0.720, 0.897) | 0.156 (0.116, 0.198) |
Once we have thoroughly evaluated our models and selected the best candidate for deployment (for this demonstration, we will use the Decision Tree), we can use it to predict the diabetes status of new, unseen patients.
First, we finalize our Decision Tree workflow using the exact tuned hyperparameters we extracted during cross-validation. Then, we fit (train) this finalized model on our full training dataset.
# 1. Finalize the workflow with the best tuned parameters
final_dt_wf <- finalize_workflow(dt_wf, dt_results$best_params)
# 2. Fit the finalized model to the full training set
final_dt_model <- fit(final_dt_wf, data = pima_train)Now, let’s imagine a new patient comes into the clinic. We measure their vitals and record their clinical metrics into a new data frame. (Note: The column names and data types must perfectly match the original training data).
# 3. Create a new observation (mock patient data)
new_patient <- tibble::tibble(
pregnant = 2,
glucose = 155,
pressure = 72,
triceps = 30,
insulin = 125,
mass = 33.5,
pedigree = 0.55,
age = 35
)
# Display the new patient's profile
knitr::kable(new_patient, caption = "New Patient Clinical Profile")| pregnant | glucose | pressure | triceps | insulin | mass | pedigree | age |
|---|---|---|---|---|---|---|---|
| 2 | 155 | 72 | 30 | 125 | 33.5 | 0.55 | 35 |
Finally, we pass this new patient data into our fitted model using
the predict() function. We can ask the model for both the
hard classification (Positive/Negative) and the underlying probability
(how confident the model is in its decision).
# 4. Generate predictions
pred_class <- predict(final_dt_model, new_data = new_patient)
pred_probs <- predict(final_dt_model, new_data = new_patient, type = "prob")
# Combine the patient data with the predictions into a single, clean table
final_prediction <- dplyr::bind_cols(
new_patient,
pred_class,
pred_probs
)
# Display the final prediction
knitr::kable(
final_prediction,
caption = "Model Prediction Results",
digits = 3
)| pregnant | glucose | pressure | triceps | insulin | mass | pedigree | age | .pred_class | .pred_pos | .pred_neg |
|---|---|---|---|---|---|---|---|---|---|---|
| 2 | 155 | 72 | 30 | 125 | 33.5 | 0.55 | 35 | pos | 0.603 | 0.397 |