Tò mò tìm hiểu coi những thuật toán nào thường được sử dụng bởi các đội thuộc Top 5 trên sân chơi Kaggle về khoa học dữ liệu thì vấp phải cụm từ Primary ML software used by top-5 teams on Kaggle với cái hình như trên. Đại ý là Keras Deep Learning xếp đầu bảng, LightGBM và Xgboost lần lượt xếp vị trị thứ hai và ba. Tuy vậy quan điểm của mình là Using the Right Tool for the Right Job: không có một mô hình nào chiếm ưu thế, với mỗi bộ dữ liệu, một vấn đề và hoàn cảnh khác nhau thì một mô hình được cho là “yếu” lại có thể hiệu quả hơn.
Cụ thể hơn chúng ta xét một bài toán rất cụ thể: bài toán Binary Classification với bộ dữ liệu IBM Watson Telco Customer Churn Data Set. Mục tiêu của bài toán là xây dựng mô hình dự báo nhân viên sẽ rời bỏ công ti và tiêu chuẩn lựa chọn mô hình tốt, giả sử, là ROC-AUC. Chúng ta sẽ so sánh khả năng dự báo của Artificial Neural Network và Xgboost căn cứ vào tiêu chí này.
Trước hết chúng ta thực hiện việc tiền xử lí số liệu - chuẩn bị số liệu như đã trình bày trong post trước:
# Load data:
library(tidyverse)
churn_data_raw <- read_csv("WA_Fn-UseC_-Telco-Customer-Churn.csv")
# Remove unnecessary data:
churn_data_tbl <- churn_data_raw %>%
select(-customerID) %>%
drop_na()
# Conduct one-hot encoding process (https://rdrr.io/rforge/caret/man/dummyVars.html):
library(caret)
dummies <- dummyVars("~ .", data = churn_data_tbl %>% select(-Churn))
predict(dummies, churn_data_tbl %>% select(-Churn)) %>% as.data.frame() -> df_features
# Normalize 0-1 for data:
df_features %>%
mutate_all(function(x) {(x - min(x)) / (max(x) - min(x))}) %>%
mutate(Churn = case_when(churn_data_tbl$Churn == "Yes" ~ 1, TRUE ~ 0)) -> df_final
# Split data (train data for testing and test data for evaluation/comparing):
set.seed(1)
id <- createDataPartition(df_final$Churn, p = 0.7, list = FALSE)
train <- df_final[id, ]
test <- df_final[-id, ]
x_train_tbl <- train %>%
select(-Churn) %>%
as.matrix()
n <- ncol(x_train_tbl)
y_train_vec <- train %>% pull(Churn)
x_test_tbl <- test %>%
select(-Churn) %>%
as.matrix()
y_test_vec <- test %>% pull(Churn)Trước hết chúng ta xây dựng một Artificial Neural Network (ANN) hai layers trong đó tỉ lệ drop-out lần lượt là 20%, 20% với learning rate = 0.1:
#--------------------------------------
# Train Artificial Neural Network
# with some parameters selected
#--------------------------------------
# Building our Artificial Neural Network:
library(keras)
model_keras <- keras_model_sequential()
model_keras %>%
# First hidden layer:
layer_dense(units = 16,
kernel_initializer = "uniform",
activation = "relu",
input_shape = n) %>%
# Dropout rate to prevent overfitting:
layer_dropout(rate = 0.2, seed = 29) %>%
# Second hidden layer:
layer_dense(units = 16,
kernel_initializer = "uniform",
activation = "relu") %>%
# Dropout rate to prevent overfitting again:
layer_dropout(rate = 0.2, seed = 29) %>%
# Output layer:
layer_dense(units = 1,
kernel_initializer = "uniform",
activation = "sigmoid") %>%
# Compile ANN:
compile(optimizer = optimizer_adamax(lr = 0.1),
loss = "binary_crossentropy",
metrics = "accuracy")Rồi thực hiện huấn luyện ANN trên train data:
# Train the ANN keras model to the training data:
fit_keras <- fit(object = model_keras,
x = x_train_tbl,
y = y_train_vec,
batch_size = 2^6,
epochs = 30,
validation_split = 0.30)Chúng ta có thể xem quá trình huấn luyện ANN bằng công cụ hình ảnh như sau:
# Training process by plot:
plot(fit_keras) +
labs(title = "Figure 1: Deep Learning Training Results",
subtitle = "Data Source: IBM Watson Telco Customer Churn Data Set")OK. Với ANN đã có chúng ta có thể đánh giá sơ bộ chất lượng dự báo của nó trên test data. Trước hết là ma trận nhầm lẫn với ngưỡng mặc định 0.5:
# Predicted Probability:
prob_keras <- predict_proba(model_keras, x_test_tbl) %>% as.vector()
# Function calculates confusion matrix:
my_cm <- function(prediction, actual, cutoff) {
pred <- case_when(prediction >= cutoff ~ "Yes", TRUE ~ "No") %>% as.factor()
thuc_te <- case_when(actual == 1 ~ "Yes", TRUE ~ "No") %>% as.factor()
confusionMatrix(pred, thuc_te, positive = "Yes")
}
# Confusion matrix with default cutoff = 0.5:
my_cm(prob_keras, y_test_vec, cutoff = 0.5)## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 1333 228
## Yes 200 348
##
## Accuracy : 0.7971
## 95% CI : (0.7793, 0.814)
## No Information Rate : 0.7269
## P-Value [Acc > NIR] : 5.496e-14
##
## Kappa : 0.481
##
## Mcnemar's Test P-Value : 0.1919
##
## Sensitivity : 0.6042
## Specificity : 0.8695
## Pos Pred Value : 0.6350
## Neg Pred Value : 0.8539
## Prevalence : 0.2731
## Detection Rate : 0.1650
## Detection Prevalence : 0.2598
## Balanced Accuracy : 0.7369
##
## 'Positive' Class : Yes
##
Và ROC-AUC - tiêu chuẩn được chọn để đánh giá chất lượng dự báo của mô hình là:
# Function calculates ROC-AUC:
library(pROC)
getAUC <- function(prediction, actual) {
auc <- roc(actual, prediction)$auc %>% as.numeric()
return(auc)
}
# AUC on Test data by Artificial Neural Network:
auc_ann <- getAUC(prob_keras, y_test_vec)
auc_ann ## [1] 0.8454012
AUC = 0.8454012 khá cao theo một số tiêu chuẩn thường thấy để đánh giá chất lượng của mô hình phân loại.
Chúng ta huấn luyện Xgboost mặc định không tinh chỉnh bất cứ tham số nào như sau:
#------------------------------------------
# Train XGBoost with default parameters
#------------------------------------------
library(xgboost)
dtrain <- xgb.DMatrix(data = x_train_tbl, label = y_train_vec)
# Train a default XGBoost:
xgb_default <- xgboost(data = dtrain,
eval_metric = "auc",
objective = "binary:logistic",
verbose = 0,
nround = 30)Rồi đánh giá khả năng phân loại của Xgboost trên test data qua ma trận nhầm lẫn và AUC như đã làm với Keras ANN:
# Use Xgboost for predicting:
pd_xgb <- predict(xgb_default, x_test_tbl)
# COnfusion matrix by Xgboost:
my_cm(pd_xgb, y_test_vec, cutoff = 0.5)## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 1388 267
## Yes 145 309
##
## Accuracy : 0.8046
## 95% CI : (0.7871, 0.8214)
## No Information Rate : 0.7269
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.4732
##
## Mcnemar's Test P-Value : 2.503e-09
##
## Sensitivity : 0.5365
## Specificity : 0.9054
## Pos Pred Value : 0.6806
## Neg Pred Value : 0.8387
## Prevalence : 0.2731
## Detection Rate : 0.1465
## Detection Prevalence : 0.2153
## Balanced Accuracy : 0.7209
##
## 'Positive' Class : Yes
##
## [1] 0.8534515
Như vậy Xgboost mặc định là tốt hơn Keras KNN vì AUC = 0.8534515 cao hơn 0.8454012. Không những thế các metrics đánh giá model performance trên ma trận nhầm lẫn của Xgboost cũng ngon hơn Keras ANN.
Phe “yêu thích” Keras chưa chịu thua. Trong một nỗ lực đánh bại Xgboost họ sử dụng Bayesian Optimization để tinh chỉnh và tìm kiếm tham số tối ưu cho Keras KNN. Giả định răng có ba tham số dropout1, dropout2 và learning_rate (trong số vô vàn tham số của Keras KNN) sẽ được tinh chỉnh. Vậy trước hết là hàm mục tiêu (objective function) trả về thước đo cần tối đa hóa:
# Define objective function:
keras_fit <- function(dropout1, dropout2, learning_rate){
model_keras <- keras_model_sequential()
model_keras %>%
layer_dense(units = 16,
kernel_initializer = "uniform",
activation = "relu",
input_shape = n) %>%
layer_dropout(rate = dropout1) %>%
layer_dense(units = 16,
kernel_initializer = "uniform",
activation = "relu") %>%
layer_dropout(rate = dropout2) %>%
layer_dense(units = 1,
kernel_initializer = "uniform",
activation = "sigmoid") %>%
compile(optimizer = optimizer_adamax(lr = learning_rate),
loss = "binary_crossentropy",
metrics = "accuracy")
fit_keras <- fit(object = model_keras,
x = x_train_tbl,
y = y_train_vec,
verbose = 0,
batch_size = 2^6,
epochs = 30,
validation_split = 0.30)
prob_keras <- predict_proba(model_keras, x_test_tbl) %>% as.vector()
auc <- getAUC(prob_keras, y_test_vec)
result <- list(Score = auc, Pred = NULL)
return(result)
}Rồi thiết lập không gian tìm kiếm và không gian search ban đầu:
# Define domain space:
search_bound_keras <- list(dropout1 = c(0.1, 0.3),
dropout2 = c(0.1, 0.3),
learning_rate = c(0.1, 0.5))
# Define initial search:
set.seed(29)
search_grid_keras <- data.frame(dropout1 = runif(10, 0.1, 0.3),
dropout2 = runif(10, 0.1, 0.3),
learning_rate = runif(10, 0.1, 0.5))Rồi thực hiện tinh chỉnh sử dụng Bayesian Optimization:
# Search optimal hyperparameter by Bayesian Optimization:
library(rBayesianOptimization)
set.seed(29)
system.time(bayes_keras <- BayesianOptimization(keras_fit,
bounds = search_bound_keras,
init_points = 0,
init_grid_dt = search_grid_keras,
n_iter = 10,
kappa = 2.576,
verbose = FALSE,
eps = 0.0,
acq = "ucb"))##
## Best Parameters Found:
## Round = 19 dropout1 = 0.3000 dropout2 = 0.1360 learning_rate = 0.1000 Value = 0.8435
## user system elapsed
## 868.931 640.457 220.980
Quá trình tinh chỉnh có thể mất nhiều thời gian nên chúng ta sử dụng hàm system.time() để tính luôn thời gian tinh chỉnh. Và kết quả cao nhất của AUC cho Keras ANN là:
## [1] 0.843528
Tương ứng với các tham số tối ưu là:
## dropout1 dropout2 learning_rate
## 0.3000000 0.1359839 0.1000000
Thật không may là AUC tương ứng với tham số tối ưu tìm được vẫn còn thua AUC của Xgboost mặc định. Và giả định rằng AUC của Keras KNN tinh chỉnh có cao hơn Xgboost mặc định thì “phe Xgboost” họ vẫn có thể tinh chỉnh các tham số của Xgboost.
Cũng cần phải lưu ý rằng huấn luyện Keras ANN mất rất nhiều thời gian, ít nhất là trên CPu. Trên GPU thì có thể đỡ hơn. Một nhược điểm của Keras ANN chính là AUC đặc biệt nhạy cảm với tham số. Thật vậy, AUC có thể nằm đâu đó trong một khoảng rất rộng từ dưới 0.5 đến 0.839 như ta có thể thấy:
# Some statistics about ROC-AUC:
bayes_keras$History %>%
as.data.frame() %>%
pull(Value) %>%
summary()## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.5000 0.7655 0.8250 0.7547 0.8361 0.8435
Không có mô hình nào chiếm ưu thế tuyệt đối. Vay vậy việc xây dựng, đánh giá và lựa chọn mô hình / tool sử dụng hay rộng hơn và cách tiếp cận nó sẽ nên là Right Tool for the Right Job.
## R version 3.6.2 (2019-12-12)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=vi_VN LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=vi_VN LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=vi_VN LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=vi_VN LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rBayesianOptimization_1.1.0 xgboost_0.90.0.2
## [3] pROC_1.14.0 keras_2.2.5.0
## [5] caret_6.0-84 lattice_0.20-40
## [7] forcats_0.4.0 stringr_1.4.0
## [9] dplyr_0.8.0.1 purrr_0.3.3
## [11] readr_1.3.1 tidyr_0.8.3
## [13] tibble_2.1.3 ggplot2_3.2.1
## [15] tidyverse_1.2.1
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.0 jsonlite_1.6 splines_3.6.2
## [4] foreach_1.4.4 prodlim_2018.04.18 modelr_0.1.4
## [7] assertthat_0.2.1 stats4_3.6.2 GPfit_1.0-8
## [10] cellranger_1.1.0 yaml_2.2.0 ipred_0.9-9
## [13] pillar_1.4.3 backports_1.1.5 glue_1.3.1
## [16] reticulate_1.14 digest_0.6.23 rvest_0.3.3
## [19] colorspace_1.4-1 recipes_0.1.5 htmltools_0.3.6
## [22] Matrix_1.2-18 plyr_1.8.5 timeDate_3043.102
## [25] pkgconfig_2.0.3 lhs_1.0.1 broom_0.5.2
## [28] haven_2.1.0 scales_1.1.0 whisker_0.3-2
## [31] gower_0.2.0 lava_1.6.5 generics_0.0.2
## [34] farver_2.0.1 withr_2.1.2 nnet_7.3-12
## [37] lazyeval_0.2.2 cli_2.0.0 survival_3.1-8
## [40] magrittr_1.5 crayon_1.3.4 readxl_1.3.1
## [43] evaluate_0.13 fansi_0.4.0 nlme_3.1-144
## [46] MASS_7.3-51.5 xml2_1.2.0 class_7.3-15
## [49] tools_3.6.2 data.table_1.12.2 hms_0.4.2
## [52] lifecycle_0.1.0 munsell_0.5.0 e1071_1.7-1
## [55] compiler_3.6.2 rlang_0.4.2 grid_3.6.2
## [58] iterators_1.0.10 rstudioapi_0.10 rappdirs_0.3.1
## [61] base64enc_0.1-3 labeling_0.3 rmarkdown_1.12
## [64] gtable_0.3.0 ModelMetrics_1.2.2 codetools_0.2-16
## [67] reshape2_1.4.3 R6_2.4.1 tfruns_1.4
## [70] lubridate_1.7.4 knitr_1.22 tensorflow_2.0.0
## [73] zeallot_0.1.0 stringi_1.4.3 Rcpp_1.0.3
## [76] rpart_4.1-15 tidyselect_0.2.5 xfun_0.6