This practical session uses the Pima Indians Diabetes dataset. We will use a 70:30 data split, 5-fold cross-validation for hyperparameter tuning, and evaluate the models using Accuracy, F1-score, AUC, and Brier’s score.
install.packages(c("tidymodels", "mlbench", "knitr"))
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(mlbench)We load the dataset
?PimaIndiansDiabetes2#> _P_i_m_a _I_n_d_i_a_n_s _D_i_a_b_e_t_e_s _D_a_t_a_b_a_s_e
#>
#> _D_e_s_c_r_i_p_t_i_o_n:
#>
#> A data frame with 768 observations on 9 variables.
#>
#> _U_s_a_g_e:
#>
#> data("PimaIndiansDiabetes", package = "mlbench")
#> data("PimaIndiansDiabetes2", package = "mlbench")
#>
#> _F_o_r_m_a_t:
#>
#> pregnant Number of times pregnant
#> glucose Plasma glucose concentration (glucose tolerance test)
#> pressure Diastolic blood pressure (mm Hg)
#> triceps Triceps skin fold thickness (mm)
#> insulin 2-Hour serum insulin (mu U/ml)
#> mass Body mass index
#> pedigree Diabetes pedigree function
#> age Age (years)
#> diabetes Class variable (test for diabetes)
#>
#> _D_e_t_a_i_l_s:
#>
#> The data set 'PimaIndiansDiabetes2' contains a corrected version
#> of the original data set. While the UCI repository index claims
#> that there are no missing values, closer inspection of the data
#> shows several physical impossibilities, e.g., blood pressure or
#> body mass index of 0. In 'PimaIndiansDiabetes2', all zero values
#> of 'glucose', 'pressure', 'triceps', 'insulin' and 'mass' have
#> been set to 'NA', see also Wahba, Gu, Wang, and Chappell (1995)
#> and Ripley (1996).
#>
#> _S_o_u_r_c_e:
#>
#> • Original owners: National Institute of Diabetes and Digestive
#> and Kidney Diseases
#>
#> • Donor of database: Vincent Sigillito (vgs@aplcen.apl.jhu.edu)
#>
#> These data have been taken from the UCI Repository Of Machine
#> Learning Databases (Blake and Merz 1998) and were converted to R
#> format by Friedrich Leisch in the late 1990s.
#>
#> The data no longer seems to be available from the UC Irvine
#> Machine Learning Repository (now at
#> <https://archive.ics.uci.edu/>).
#>
#> _R_e_f_e_r_e_n_c_e_s:
#>
#> Blake CL, Merz CJ (1998). "UCI Repository of Machine Learning
#> Databases." University of California, Irvine, Department of
#> Information and Computer Science. Formerly available from
#> 'http://www.ics.uci.edu/~mlearn/MLRepository.html'.
#>
#> Ripley BD (1996). _Pattern Recognition and Neural Networks_.
#> Cambridge University Press, Cambridge.
#> doi:10.1017/CBO9780511812651
#> <https://doi.org/10.1017/CBO9780511812651>.
#>
#> Wahba G, Gu C, Wang Y, Chappell R (1995). "Chapter Soft
#> Classification, a.k.a. Risk Estimation, via Penalized Log
#> Likelihood and Smoothing Spline Analysis of Variance." In Wolpert
#> DH (ed.), _The Mathematics of Generalization_, 331-359.
#> Addison-Wesley, Reading, MA.
#>
#> _E_x_a_m_p_l_e_s:
#>
#> data("PimaIndiansDiabetes", package = "mlbench")
#> summary(PimaIndiansDiabetes)
#>
#> data("PimaIndiansDiabetes2", package = "mlbench")
#> summary(PimaIndiansDiabetes2)
#> data("PimaIndiansDiabetes2")
pima <- PimaIndiansDiabetes2
str(pima)#> 'data.frame': 768 obs. of 9 variables:
#> $ pregnant: num 6 1 8 1 0 5 3 10 2 8 ...
#> $ glucose : num 148 85 183 89 137 116 78 115 197 125 ...
#> $ pressure: num 72 66 64 66 40 74 50 NA 70 96 ...
#> $ triceps : num 35 29 NA 23 35 NA 32 NA 45 NA ...
#> $ insulin : num NA NA NA 94 168 NA 88 NA 543 NA ...
#> $ mass : num 33.6 26.6 23.3 28.1 43.1 25.6 31 35.3 30.5 NA ...
#> $ pedigree: num 0.627 0.351 0.672 0.167 2.288 ...
#> $ age : num 50 31 32 21 33 30 26 29 53 54 ...
#> $ diabetes: Factor w/ 2 levels "neg","pos": 2 1 2 1 2 1 2 1 2 2 ...Relevel so the positive class (“pos”) is the first level for yardstick metrics.
pima$diabetes <- relevel(pima$diabetes, ref = "pos")
levels(pima$diabetes)#> [1] "pos" "neg"Remove observations with NA (missing data). In practice we can save them with imputation.
pima <- na.omit(pima)
str(pima)#> 'data.frame': 392 obs. of 9 variables:
#> $ pregnant: num 1 0 3 2 1 5 0 1 1 3 ...
#> $ glucose : num 89 137 78 197 189 166 118 103 115 126 ...
#> $ pressure: num 66 40 50 70 60 72 84 30 70 88 ...
#> $ triceps : num 23 35 32 45 23 19 47 38 30 41 ...
#> $ insulin : num 94 168 88 543 846 175 230 83 96 235 ...
#> $ mass : num 28.1 43.1 31 30.5 30.1 25.8 45.8 43.3 34.6 39.3 ...
#> $ pedigree: num 0.167 2.288 0.248 0.158 0.398 ...
#> $ age : num 21 33 26 53 59 51 31 33 32 27 ...
#> $ diabetes: Factor w/ 2 levels "pos","neg": 2 1 1 1 1 1 1 2 1 2 ...
#> - attr(*, "na.action")= 'omit' Named int [1:376] 1 2 3 6 8 10 11 12 13 16 ...
#> ..- attr(*, "names")= chr [1:376] "1" "2" "3" "6" ...summary(pima)#> pregnant glucose pressure triceps
#> Min. : 0.000 Min. : 56.0 Min. : 24.00 Min. : 7.00
#> 1st Qu.: 1.000 1st Qu.: 99.0 1st Qu.: 62.00 1st Qu.:21.00
#> Median : 2.000 Median :119.0 Median : 70.00 Median :29.00
#> Mean : 3.301 Mean :122.6 Mean : 70.66 Mean :29.15
#> 3rd Qu.: 5.000 3rd Qu.:143.0 3rd Qu.: 78.00 3rd Qu.:37.00
#> Max. :17.000 Max. :198.0 Max. :110.00 Max. :63.00
#> insulin mass pedigree age diabetes
#> Min. : 14.00 Min. :18.20 Min. :0.0850 Min. :21.00 pos:130
#> 1st Qu.: 76.75 1st Qu.:28.40 1st Qu.:0.2697 1st Qu.:23.00 neg:262
#> Median :125.50 Median :33.20 Median :0.4495 Median :27.00
#> Mean :156.06 Mean :33.09 Mean :0.5230 Mean :30.86
#> 3rd Qu.:190.00 3rd Qu.:37.10 3rd Qu.:0.6870 3rd Qu.:36.00
#> Max. :846.00 Max. :67.10 Max. :2.4200 Max. :81.00Descriptive statistics
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 split the data into training (70%) and testing (30%) sets, stratified by the diabetes outcome to maintain class proportions. We also create 5-fold cross-validation folds from the training data.
set.seed(123)
# 70:30 split, stratified by outcome
pima_split <- initial_split(pima, prop = 0.70, strata = diabetes)
pima_train <- training(pima_split)
pima_test <- testing(pima_split)
# 5-fold Cross-Validation
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()) # data pre-processing
# Normalizing numeric predictors is crucial for glmnetWe specify two models: standard logistic regression using the
glm engine, and penalized logistic regression using the
glmnet engine. For the penalized model, we tag
penalty and mixture for tuning.
Standard Logistic Regression (No tuning needed)
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")Workflows bundle the model specification and the formula together
log_wf <- workflow() |> add_recipe(pima_rec) |> add_model(log_spec)
penal_wf <- workflow() |> add_recipe(pima_rec) |> add_model(penal_spec)We fit the standard model directly to the training data.
# Fit standard model
log_fit <- fit(log_wf, data = pima_train)For the penalized model, we use the 5-fold CV to test 20 random hyperparameter combinations, optimizing for AUC.
# Tune penalized model
set.seed(345)
tune_res <- tune_grid(
penal_wf,
resamples = cv_folds,
grid = 20,
metrics = metric_set(roc_auc)
)
# Extract best hyperparameters and finalize the model
best_params <- select_best(tune_res, metric = "roc_auc")
final_penal_wf <- finalize_workflow(penal_wf, best_params)
penal_fit <- fit(final_penal_wf, data = pima_train)Standard Logistic
predictions <- pima_train |>
bind_cols(
predict(log_fit, new_data = pima_train),
predict(log_fit, new_data = pima_train, type = "prob")
)
predictions |> conf_mat(truth = diabetes, estimate = .pred_class)#> Truth
#> Prediction pos neg
#> pos 58 18
#> neg 33 165predictions |> accuracy(truth = diabetes, estimate = .pred_class)#> # A tibble: 1 × 3
#> .metric .estimator .estimate
#> <chr> <chr> <dbl>
#> 1 accuracy binary 0.814predictions |> f_meas(truth = diabetes, estimate = .pred_class)#> # A tibble: 1 × 3
#> .metric .estimator .estimate
#> <chr> <chr> <dbl>
#> 1 f_meas binary 0.695predictions |> roc_auc(truth = diabetes, .pred_pos)#> # A tibble: 1 × 3
#> .metric .estimator .estimate
#> <chr> <chr> <dbl>
#> 1 roc_auc binary 0.874predictions |> brier_class(truth = diabetes, .pred_pos)#> # A tibble: 1 × 3
#> .metric .estimator .estimate
#> <chr> <chr> <dbl>
#> 1 brier_class binary 0.133To be more concise, we define all metrics: Accuracy, F1-score (f_meas), AUC (roc_auc), Brier score (brier_class).
eval_metrics <- metric_set(accuracy, f_meas, roc_auc, brier_class)
eval_metrics(
predictions,
truth = diabetes,
estimate = .pred_class,
.pred_pos
) -> log_metric0; log_metric0#> # A tibble: 4 × 3
#> .metric .estimator .estimate
#> <chr> <chr> <dbl>
#> 1 accuracy binary 0.814
#> 2 f_meas binary 0.695
#> 3 roc_auc binary 0.874
#> 4 brier_class binary 0.133Penalised Logistic
predictions <- pima_train |>
bind_cols(
predict(penal_fit, new_data = pima_train),
predict(penal_fit, new_data = pima_train, type = "prob")
)
eval_metrics(
predictions,
truth = diabetes,
estimate = .pred_class,
.pred_pos
) -> penal_metric0; penal_metric0#> # A tibble: 4 × 3
#> .metric .estimator .estimate
#> <chr> <chr> <dbl>
#> 1 accuracy binary 0.818
#> 2 f_meas binary 0.699
#> 3 roc_auc binary 0.873
#> 4 brier_class binary 0.133Standard Logistic
predictions <- pima_test |>
bind_cols(
predict(log_fit, new_data = pima_test),
predict(log_fit, new_data = pima_test, type = "prob")
)
eval_metrics(
predictions,
truth = diabetes,
estimate = .pred_class,
.pred_pos
) -> log_metric; log_metric#> # A tibble: 4 × 3
#> .metric .estimator .estimate
#> <chr> <chr> <dbl>
#> 1 accuracy binary 0.729
#> 2 f_meas binary 0.543
#> 3 roc_auc binary 0.821
#> 4 brier_class binary 0.169Penalised Logistic
predictions <- pima_test |>
bind_cols(
predict(penal_fit, new_data = pima_test),
predict(penal_fit, new_data = pima_test, type = "prob")
)
eval_metrics(
predictions,
truth = diabetes,
estimate = .pred_class,
.pred_pos
) -> penal_metric; penal_metric#> # A tibble: 4 × 3
#> .metric .estimator .estimate
#> <chr> <chr> <dbl>
#> 1 accuracy binary 0.729
#> 2 f_meas binary 0.543
#> 3 roc_auc binary 0.821
#> 4 brier_class binary 0.168Finally, we compare the metrics for all models on train & test set and format them into a comparison table.
log_metric0 |>
dplyr::mutate(model = "Standard Logistic (Train)") |>
dplyr::select(model, .metric, .estimate) -> log_perf0
penal_metric0 |>
dplyr::mutate(model = "Penalized Logistic (Train)") |>
dplyr::select(model, .metric, .estimate) -> penal_perf0
log_metric |>
dplyr::mutate(model = "Standard Logistic (Test)") |>
dplyr::select(model, .metric, .estimate) -> log_perf
penal_metric |>
dplyr::mutate(model = "Penalized Logistic (Test)") |>
dplyr::select(model, .metric, .estimate) -> penal_perf
# Print the final output
comparison <- dplyr::bind_rows(
log_perf0, log_perf,
penal_perf0, penal_perf) |>
tidyr::pivot_wider(names_from = .metric, values_from = .estimate)
comparison |>
knitr::kable(
col.names = c("Model", "Accuracy", "F1 Score", "AUC", "Brier Score"),
digits = 3, # Optional: rounds the decimals to make the table cleaner
align = "lcccc" # Optional: Left-align the first column, center the rest
)| Model | Accuracy | F1 Score | AUC | Brier Score |
|---|---|---|---|---|
| Standard Logistic (Train) | 0.814 | 0.695 | 0.874 | 0.133 |
| Standard Logistic (Test) | 0.729 | 0.543 | 0.821 | 0.169 |
| Penalized Logistic (Train) | 0.818 | 0.699 | 0.873 | 0.133 |
| Penalized Logistic (Test) | 0.729 | 0.543 | 0.821 | 0.168 |
Notice that there are discrepancies between train and test performance. Train performance is higher than when it is tested on the unseen test set. This is overfitting, the model unable to generalize to unseen data. How to get better idea of the real performance of the model in the training data? We will use k-fold CV in the next training session to get averaged performance over the folds.