Kelebihan tidymodels

Berikut beberapa kelebihan tidymodels:

1. Konsisten dan mudah digunakan: Menggunakan prinsip tidyverse, sehingga kode lebih mudah dibaca dan ditulis. Semua alat dalam tidymodels memiliki antarmuka yang seragam, membuat penggunaannya lebih intuitif.

2. Komprehensif: Menyediakan alat untuk semua tahap pemodelan, mulai dari pra-pemrosesan data, resampling, pemilihan model, tuning hiperparameter, hingga evaluasi model.

3. Fleksibel: Mendukung berbagai jenis model seperti regresi, klasifikasi, dan clustering, serta berbagai engine seperti glm untuk regresi, ranger untuk Random Forest, xgboost untuk Gradient Boosting, dan banyak lagi.

4. Dukungan tuning hiperparameter: Mendukung proses tuning hiperparameter langsung, memungkinkan pencarian kombinasi hiperparameter optimal secara otomatis.

Pemodelan dengan tidymodels

Secara umum, pemodelan dengan tidymodels dibagi menjadi lima tahap:

1. Pra-pemrosesan data: Meliputi penanganan missing value dan outliers, normalisasi data, serta encoding variabel kategorik. Paket recipes dalam tidymodels membantu mengotomatisasi langkah-langkah ini dan memastikan semua proses pada data dilakukan secara konsisten baik saat pelatihan maupun pengujian, atau bahkan ketika terdapat data baru.

2. Memilih model: Dengan paket parsnip, kita bisa menentukan model tanpa terikat pada library tertentu. Misalnya, membuat model Random Forest atau Gradient Boosting hanya dengan beberapa baris kode. parsnip menyediakan antarmuka konsisten untuk berbagai jenis model meskipun menggunakan engine berbeda.

3. Membangun workflow: Paket workflows memungkinkan penggabungan seluruh proses pra-pemrosesan data dan model menjadi satu objek workflow yang terintegrasi. Workflow ini mencakup semua langkah dari pengolahan data hingga fitting model, memudahkan pengelolaan dan pengujian berbagai skenario pemodelan, memastikan setiap langkah terkoordinasi dengan baik.

4. Validasi dan tuning hiperparameter: Paket rsample digunakan untuk membagi data ke dalam set latihan dan pengujian serta melakukan validasi silang untuk mengevaluasi kinerja model dengan lebih robust. Sementara itu, paket tune membantu mencari kombinasi hiperparameter optimal untuk model.

5. Evaluasi model: Paket yardstick menyediakan berbagai fungsi untuk mengevaluasi kinerja model, seperti akurasi, presisi, recall, F1-score, dan menampilkan confusion matrix.

Dengan tidymodels, seluruh proses pemodelan data menjadi lebih mudah dan terstruktur, memungkinkan kita untuk fokus pada analisis dan interpretasi hasil.

Pada bagian selanjutnya kita akan membahas contoh penggunaan tidymodels untuk membangun model random forest dan gradient boosting. Data yang akan digunakan adalah dataset Heart Failure yang diunduh pada halaman berikut: Heart Failure Prediction Dataset.

# loading dataset
data <- read.csv('https://raw.githubusercontent.com/sainsdataid/dataset/main/heart.csv')

data$HeartDisease <- as.factor(data$HeartDisease)

head(data)

##   Age Sex ChestPainType RestingBP Cholesterol FastingBS RestingECG MaxHR
## 1  40   M           ATA       140         289         0     Normal   172
## 2  49   F           NAP       160         180         0     Normal   156
## 3  37   M           ATA       130         283         0         ST    98
## 4  48   F           ASY       138         214         0     Normal   108
## 5  54   M           NAP       150         195         0     Normal   122
## 6  39   M           NAP       120         339         0     Normal   170
##   ExerciseAngina Oldpeak ST_Slope HeartDisease
## 1              N     0.0       Up            0
## 2              N     1.0     Flat            1
## 3              N     0.0       Up            0
## 4              Y     1.5     Flat            1
## 5              N     0.0       Up            0
## 6              N     0.0       Up            0

Data Preprocessing

Tahapan ini mencakup pembagian data latih dan data uji, imputasi nilai prediktor yang terdapat missing value serta transformasi nilai menggunakan normalisasi atau teknik lainnya. Selain tiga hal tersebut masih banyak hal lain yang dapat juga dilakukan sesuai keperluan dan kondisi data. Untuk mendalaminya lebih lanjut silahkan cek dokumentasi paket recipes.

Pada contoh ini data dibagi sebanyak 70% sebagai data latih dan 30% untuk data uji. Selanjutnya kita lakukan proses imputasi (jika diperlukan), misalnya mengganti missing value dengan nilai mediannya. Imputasi untuk data nominal atau kategorik juga dapat dilakukan dengan nilai modus (kategori yang paling banyak muncul pada data).

Tahapan berikutnya yaitu melakukan normalisasi data. (Sebagai catatan, model-model berbasis pohon sebenarnya tidak terlalu terpengaruh dengan perubahan skala data), contoh di bawah hanya untuk menunjukkan bahwa proses tersebut dapat otomatisasi.

Tahapan preprocessing hanya dilakukan berdasarkan data latih saja. Perlakuan ini penting agar model dapat diuji dengan data yang sama sekali tidak terkait dengan proses pelatihan, baik secara langsung maupun tidak langsung. konteks secara langsung tentunya mudah dipahami, yaitu data yang digunakan untuk pengujian juga digunakan dalam proses pelatihan. Sementara itu, keterkaitan tidak langsung dapat terjadi pada tahapan preprocessing ini seperti imputasi maupun normalisasi. Misalkan, jika missing value dihitung berdasarkan median atau rata-rata seluruh data, maka tentu saja hasil imputasi pada data latih secara tidak langsung mengandung informasi dari data uji (data leakage). Begitu pula saat melakukan transformasi data, baik normalisasi maupun transformasi lainnya.

# install.packages("tidymodels")
library(tidymodels)

## ── Attaching packages ────────────────────────────────────── tidymodels 1.2.0 ──

## ✔ broom        1.0.6      ✔ recipes      1.0.10
## ✔ dials        1.2.1      ✔ rsample      1.2.1 
## ✔ dplyr        1.1.4      ✔ tibble       3.2.1 
## ✔ ggplot2      3.5.1      ✔ tidyr        1.3.1 
## ✔ infer        1.0.7      ✔ tune         1.2.1 
## ✔ modeldata    1.3.0      ✔ workflows    1.1.4 
## ✔ parsnip      1.2.1      ✔ workflowsets 1.1.0 
## ✔ purrr        1.0.2      ✔ yardstick    1.3.1

## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ purrr::discard() masks scales::discard()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ recipes::step()  masks stats::step()
## • Learn how to get started at https://www.tidymodels.org/start/

set.seed(111)
# Splitting data for training dan testing
data_split <- initial_split(data, prop = 0.7)
data_train <- training(data_split)
data_test <- testing(data_split)

# preprocessing
# note : tree-based model generally can perform well without any scaling
# (this following syntax for ilustration purposed only)

data_recipe <- recipe(HeartDisease ~ ., data = data_train) %>%
  step_impute_median(all_numeric_predictors()) %>%     # missing value imputation with median (numeric)
  step_impute_mode(all_nominal_predictors()) %>%       # missing values imputation with mode (nominal)
  step_normalize(all_numeric_predictors())             # features normalization

Pemodelan

Fungsi rand_forest memiliki beberapa parameter utama yang dapat diatur yaitu:

• mtry: jumlah prediktor yang akan diambil secara acak pada setiap pemisahan ketika membuat model pohon.

• trees: jumlah pohon yang dibentuk dalam ensemble.

• min_n: jumlah minimum titik data dalam sebuah node pada pohon agar node tersebut dapat dibagi lebih lanjut.

Kita juga perlu mengatur engine serta mode yang digunakan. Dalam hal ini menggunakan engine ranger dan mode classification. Misal, kita mengatur jumlah pohon sebanyak 100, sementara untuk parameter lainnya tidak ditentukan dan menggunakan nilai defaultnya.

Setelah inisiasi model, kita dapat membuat workflow dengan menentukan model dan tahapan preprocessing yang akan dilakukan dan dilanjutkan dengan proses pelatihan (fitting) model.

# model initialization
rf_model <- rand_forest(trees = 100) %>%
  set_engine("ranger") %>%                  
  set_mode("classification")               

# creating a workflow
rf_workflow <- workflow() %>%
  add_model(rf_model) %>%
  add_recipe(data_recipe)

# model fitting
rf_fit <- rf_workflow %>%
  fit(data = data_train)

rf_fit

## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: rand_forest()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 3 Recipe Steps
## 
## • step_impute_median()
## • step_impute_mode()
## • step_normalize()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Ranger result
## 
## Call:
##  ranger::ranger(x = maybe_data_frame(x), y = y, num.trees = ~100,      num.threads = 1, verbose = FALSE, seed = sample.int(10^5,          1), probability = TRUE) 
## 
## Type:                             Probability estimation 
## Number of trees:                  100 
## Sample size:                      642 
## Number of independent variables:  11 
## Mtry:                             3 
## Target node size:                 10 
## Variable importance mode:         none 
## Splitrule:                        gini 
## OOB prediction error (Brier s.):  0.1069484

Evaluasi dan Prediksi

Model yang sudah dilatih perlu dievaluasi menggunakan data uji. Pertama-tama buat prediksi pada data uji menggunakan fungsi predict. Selanjutnya kita dapat menghitung akurasi model menggunakan fungsi metrics. Berdasarkan hasil di bawah ini, nilai akurasi model sebesar 0,870. Selain nilai akurasi kita juga dapat menampilkan confusion matrix dengan fungsi con_mat. Jika memerlukan nilai metrik lainnya (i.e precision, recall, f1-score, balanced accuracy, sensitiviy, specificity dan sebagainya) dapat juga dihitung menggunakan berbagai fungsi yang ada.

# Predicting on test data
rf_preds <- predict(rf_fit, data_test, type = "class") %>%
  bind_cols(data_test)

print(rf_preds)

## # A tibble: 276 × 13
##    .pred_class   Age Sex   ChestPainType RestingBP Cholesterol FastingBS
##    <fct>       <int> <chr> <chr>             <int>       <int>     <int>
##  1 0              40 M     ATA                 140         289         0
##  2 0              49 F     NAP                 160         180         0
##  3 0              45 F     ATA                 130         237         0
##  4 1              38 M     ASY                 110         196         0
##  5 0              43 F     TA                  100         223         0
##  6 0              49 F     ATA                 124         201         0
##  7 1              44 M     ATA                 150         288         0
##  8 0              52 M     ATA                 120         284         0
##  9 0              53 F     ATA                 113         468         0
## 10 1              54 M     ASY                 125         224         0
## # ℹ 266 more rows
## # ℹ 6 more variables: RestingECG <chr>, MaxHR <int>, ExerciseAngina <chr>,
## #   Oldpeak <dbl>, ST_Slope <chr>, HeartDisease <fct>

# Calculating accuracy value
rf_metrics <- rf_preds %>%
  metrics(truth = HeartDisease, estimate = .pred_class)

print(rf_metrics)

## # A tibble: 2 × 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.870
## 2 kap      binary         0.734

# Calculating confusion matrix
conf_mat <- rf_preds %>%
  conf_mat(truth = HeartDisease, estimate = .pred_class)

print(conf_mat)

##           Truth
## Prediction   0   1
##          0 100  19
##          1  17 140

# Calculating other metrics 
# (balanced accuracy, precision, recall, F1-score)
rf_metrics_oth <- rf_preds %>%
  summarise(
    Precision = precision_vec(truth = HeartDisease, estimate = .pred_class),
    Recall = recall_vec(truth = HeartDisease, estimate = .pred_class),
    Bal_Accuracy = bal_accuracy_vec(truth = HeartDisease, estimate = .pred_class),
    F1_Score = f_meas_vec(truth = HeartDisease, estimate = .pred_class)
  )

print(rf_metrics_oth)

## # A tibble: 1 × 4
##   Precision Recall Bal_Accuracy F1_Score
##       <dbl>  <dbl>        <dbl>    <dbl>
## 1     0.840  0.855        0.868    0.847

Tuning Hiperparameter

Model sebelumnya dibangun dengan menetapkan nilai hiperparameternya secara langsung. Nilai-nilai ini dapat juga kita tuning dengan fungsi tune() untuk mencari nilai terbaik berdasarkan ruang yang diberikan. Terdapat beberapa teknik yang didukung yaitu grid_regular, grid_random, grid_entropy dan bayes_opt. Pada contoh di bawah ini kita mengunakan fungsi grid_random dan mengatur nilai-nilai hiperparameter sebagai berikut:

• mtry: 2, 3, 4

• min_n: 2, 3, 4, …9, 10

• trees: 100, 101, …500

Dari seluruh kombinasi tersebut, pencarian hanya akan dilakukan sebanyak 100 kali (size=100) berdasarkan kombinasi yang ditentukan secara acak. Nilai hiperparameter terbaik selanjutnya ditentukan berdasarkan metrik evaluasi yang digunakan misalkan accuracy.

Model terbaik yang kita peroleh memiliki pengaturan mtry=2, min_n=2 dan tress=149. Hasil tuning tersebut memiliki nilai akurasi 0,884 dan lebih baik dibandingkan akurasi model pada bagian sebelumnya.

set.seed(111)  # Set the seed for reproducibility

# Random forest model with tunable hyperparameters
rf_model <- rand_forest(
  mtry = tune(),
  min_n = tune(),
  trees = tune()
) %>%
  set_engine("ranger") %>%
  set_mode("classification")

# Creating a workflow
rf_workflow <- workflow() %>%
  add_model(rf_model) %>%
  add_recipe(data_recipe)

# Determining hyperparameter space for random search
rf_grid <- grid_random(
  mtry(range = c(2,4)),
  min_n(range = c(2, 10)),
  trees(range = c(100, 500)),
  size = 100                   # Perform random search on 100 combinations
)

# Creating folds for training/validating data
cv_folds <- vfold_cv(data_train, v = 5)

# Performing hyperparameter tuning
rf_tune <- tune_grid(
  rf_workflow,
  resamples = cv_folds,
  grid = rf_grid
)

rf_best <- select_best(rf_tune, metric="accuracy")

print(rf_best)

## # A tibble: 1 × 4
##    mtry trees min_n .config               
##   <int> <int> <int> <chr>                 
## 1     2   129     9 Preprocessor1_Model032

rf_final <- finalize_workflow(rf_workflow, rf_best)

rf_final_fit <- rf_final %>%
  fit(data = data_train)

rf_final_preds <- predict(rf_final_fit, data_test, type = "class") %>%
  bind_cols(data_test)

# Calculating and printing final metrics
rf_final_metrics <- rf_final_preds %>%
  metrics(truth = HeartDisease, estimate = .pred_class)

print(rf_final_metrics)

## # A tibble: 2 × 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.877
## 2 kap      binary         0.747

Menyimpan dan Memuat Model

Model yang sudah diperoleh selanjutnya dapat disimpan menjadi file rds. menyimpan model memungkinkan kita untuk menggunakan kembali model tersebut tanpa perlu melakukan pelatihan kembali. Untuk menyimpan model, kita dapat menggunakan fungsi saveRDS dan untuk memuat model dengan fungsi readRDS. Model yang dimuat melalui fungsi tersebut dapat digunakan misalnya untuk melakukan prediksi.

# saving model into RDS file
saveRDS(rf_final_fit, file = "model_ok.rds")

# loading  model
loaded_model <- readRDS("model_ok.rds")

data_new <- data.frame(Age=41, 
                       Sex="M",
                       ChestPainType="NAP",
                       RestingBP=130,
                       Cholesterol=320,
                       FastingBS=0,
                       RestingECG="Normal",
                       MaxHR=170,
                       ExerciseAngina="N",
                       Oldpeak=0.0,
                       ST_Slope="Up")

# prediksi data baru
preds <- predict(loaded_model, data_new)

# melihat hasil prediksi
print(preds)

## # A tibble: 1 × 1
##   .pred_class
##   <fct>      
## 1 0

Hiperparameter Model

Hiperparameter pada model gradient boosting adalah sebagai berikut:

• mtry: jumlah (atau proporsi) prediktor yang akan diambil secara acak pada setiap pemisahan ketika membuat model pohon.

• trees: jumlah pohon yang dibuat.

• min_n:jumlah minimum titik data dalam sebuah node yang diperlukan agar node tersebut dapat dibagi lebih lanjut.

• tree_depth: kedalaman maksimum pohon (jumlah pemisahan).

• learn_rate: angka yang mengatur laju di mana algoritma boosting beradaptasi dari iterasi ke iterasi.

• loss_reduction: pengurangan dalam loss function yang diperlukan agar proses pemisahan dapat dilanjutkan.

• sample_size: angka untuk jumlah (atau proporsi) data yang diekspos ke tahapan pelatihan. Untuk xgboost, pengambilan sampel dilakukan pada setiap iterasi sementara pada C5.0 mengambil sampel sekali selama pelatihan.

• stop_iter: jumlah iterasi maksimum tanpa peningkatan yang signifikan dalam performa sebelum pelatihan dihentikan.

Tuning Hiperparameter

Kode di bawah ini menyajikan alur pemodelan dengan gradient boosting menggunakan engine xgboost. Secara umum prosesnya sama dengan pada model random forest. Beberapa perbedaan yang ditambahkan/diubah mencakup tiga hal yaitu:

• Preprocessing dengan fungsi set_dummy untuk membuat peubah dummy dari prediktor kategorik. Ha ini diperlukan karena engine xgboost hanya dapat bekerja dengan data numerik saja.

• Tuning hiperparameter dilakukan menggunakan fungsi grid_max_entropy (bukan suatu keharusan, untuk memperkaya variasi contoh saja).

• Daftar hiperparameter yang tentunya disesuaikan untuk model gradient boosting.

Hasil tuning menunjukkan model terbaik memiliki nilai akurasi 0,884 dengan rincian hiperparameter optimal seperti yang tertera pada output kode di bawah ini.

set.seed(111)  # Set the seed for reproducibility

data_recipe <- recipe(HeartDisease ~ ., data = data_train) %>%
  step_impute_median(all_numeric_predictors()) %>%    # missing value imputation with median (numeric)
  step_impute_mode(all_nominal_predictors()) %>%      # missing values imputation with mode 
  step_dummy(all_nominal_predictors()) %>%            # dummy var. for categorical predictor
  step_normalize(all_numeric_predictors())            # features normalization

# Define the XGBoost model with tunable hyperparameters
xgb_model <- boost_tree(
  trees = tune(),
  tree_depth = tune(),
  learn_rate = tune(),
  loss_reduction = tune(),
  min_n = tune(),
  mtry = tune()
) %>%
  set_engine("xgboost") %>%
  set_mode("classification")

# Create a workflow
xgb_workflow <- workflow() %>%
  add_model(xgb_model) %>%
  add_recipe(data_recipe)

# Define the hyperparameter 
xgb_grid <- grid_max_entropy(
  trees(range = c(100, 500)),
  tree_depth(range = c(2, 10)),
  learn_rate(range = c(0.01, 0.2)),
  loss_reduction(range = c(0, 10)),
  min_n(range = c(1, 20)),
  mtry(range = c(2, 4)),
  size = 100 
)

# Create folds for cross-validation
cv_folds <- vfold_cv(data_train, v = 5)

# Perform hyperparameter tuning
xgb_tune <- tune_grid(
  xgb_workflow,
  resamples = cv_folds,
  grid = xgb_grid
)

# Select the best hyperparameters based on accuracy
xgb_best <- select_best(xgb_tune, metric="accuracy")

print(xgb_best)

## # A tibble: 1 × 7
##    mtry trees min_n tree_depth learn_rate loss_reduction .config               
##   <int> <int> <int>      <int>      <dbl>          <dbl> <chr>                 
## 1     3   320     2          2       1.08           8.55 Preprocessor1_Model003

# Finalize the workflow with the best hyperparameters
xgb_final <- finalize_workflow(xgb_workflow, xgb_best)

# Fit the final model on the training data
xgb_final_fit <- xgb_final %>%
  fit(data = data_train)

# Predict on the test data
xgb_final_preds <- predict(xgb_final_fit, data_test, type = "class") %>%
  bind_cols(data_test)

# Calculate the final metrics
xgb_final_metrics <- xgb_final_preds %>%
  metrics(truth = HeartDisease, estimate = .pred_class)

print(xgb_final_metrics)

## # A tibble: 2 × 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.884
## 2 kap      binary         0.764