Monotonic Optimal Binning for Feature Selection: Empirical Evidence from Corporate Bankruptcy Prediction, Kaggle 2021 Competition

R Data Science Series

Motivation

Post 1 cho thấy việc sử dụng AUC để lựa chọn biến số cho mô hình Logistic để đạt được AUC = 0.903 trên Test Data chỉ bằng sử dụng hồi quy Logistic là một kết quả tương đối tốt. Nếu sử dụng danh sách các biến được lựa chọn để huấn luyện Automated Machine Learning thì chúng ta có thể đạt kết quả còn tốt hơn với AUC = 0.9181.

Post 2 chỉ ra rằng việc lựa chọn biến số dựa trên Information Value có thể đạt được kết quả còn cao hơn nữa (AUC = 0.937021 trên Test Data). Kết quả này thua kém không nhiều so với AUC đạt được của Team xếp thứ nhất (Team KotaShimomura có AUC = 0.93798).

Với cố gắng vượt qua được kết quả của Team KotaShimomura, post này sẽ sử dụng Monotonic Optimal Binning - MOB như là một kĩ thuật cho Feature Selection. Cụ thể từ 94 biến số ban đầu sẽ chỉ lựa chọn ra 41 biến số tiềm năng (đã được “rời rạc hóa”) dựa trên IV từ MOB. Sau đó từ danh sách các biến tiềm năng này lại chọn ra sự kết hợp tối ưu nhất các biến bằng Stepwise Selection (bạn đọc quan tâm có thể tìm hiểu thêm về thuật toán này tại trang 229 cuốn An Introduction to Statistical Learning). Kết quả thực nghiệm cho thấy chúng ta có thể đạt AUC = 0.9383908 trên Test Data. Danh sách Team dẫn đầu ở Private Score là Team KotaShimomura với AUC = 0.93798.

Vẫn như ở hai post trước, post này vẫn sử dụng bộ dữ liệu đã được sử dụng cho cuộc thi Corporate Bankruptcy Prediction 2021 trên Kaggle (có thể download tại đây).

Short Introduction to Monotonic Optimal Binning

MOB là một kĩ thuật “rời rạc hóa” được sử dụng phổ biến khi mô hình hóa rủi ro và phân loại-xếp hạng tín dụng (như PD, LGD, EAD). Ưu điểm của nó là làm cho mô hình dễ giải thích và nhất quán - một yêu cầu thường giữ vị trí quan trọng tại các tổ chức tài chính như ngân hàng.

MOB có thể được thực hiện bằng nhiều R package khác nhau. Trong post này chúng ta sử dụng thư viện mob của Liu. Dưới đây là R codes thực hiện MOB cho biến OperatingGrossMargin:

# Clear our R environment: 
rm(list = ls())

# Load tidyverse package: 
library(tidyverse)

# Load data: 
read_csv("F:/data.csv/data.csv") -> data

# Rename for all columns: 

old_names <- names(data)

old_names %>% str_replace_all("[^a-z|^A-Z]", "") -> new_names

names(data) <- new_names

# MOB for OperatingProfitRate: 

library(mob)

bin_result <- iso_bin(data$OperatingGrossMargin, data$Bankrupt)

# Decreasing event rates (or Bankrupt rate): 

bin_result$tbl %>%
  select(-rule) %>%
  knitr::kable() 

bin	freq	miss	bads	rate	woe	iv	ks
1	15	0	4	0.2667	2.3894	0.0395	1.65
2	9	0	2	0.2222	2.1483	0.0173	2.45
3	99	0	19	0.1919	1.9635	0.1458	9.88
4	90	0	15	0.1667	1.7916	0.1018	15.56
5	130	0	19	0.1462	1.6360	0.1138	22.51
6	54	0	6	0.1111	1.3216	0.0264	24.51
7	71	0	7	0.0986	1.1881	0.0263	26.73
8	13	0	1	0.0769	0.9161	0.0025	27.00
9	934	0	50	0.0535	0.5286	0.0493	36.33
10	277	0	12	0.0433	0.3062	0.0044	37.77
11	576	0	21	0.0365	0.1266	0.0014	38.90
12	55	0	2	0.0364	0.1239	0.0001	39.01
13	804	0	22	0.0274	-0.1698	0.0031	37.16
14	1863	0	21	0.0113	-1.0730	0.1971	18.79
15	1379	0	15	0.0109	-1.1091	0.1536	4.94
16	309	0	3	0.0097	-1.2239	0.0401	1.67
17	141	0	1	0.0071	-1.5406	0.0257	0.00

Biến OperatingGrossMargin sẽ được “rời rạc hóa” thành 17 khoảng khác nhau sao cho tỉ lệ Bankrupt là đơn điệu giảm. Ví dụ, ở khoảng OperatingGrossMargin <= 0.5408913 thì Bankrupt Rate = 0.2667 là cao nhất. Tại các khoảng tăng dần của OperatingGrossMargin thì Bankrupt Rate cũng giảm dần một cách đơn điệu.

IV from Monotonic Binning for Feature Selection

Monotonic Binning có thể được thực hiện bằng nhiều thuật toán. Từ đơn giản như K-mean Clustering cho đến việc sử dụng các thuật toán Machine Learning như XGBoost. Post này sử dụng cách tiếp cận isotonic regression để thực hiện MOB với hàm iso_bin(). Chúng ta khảo sát AUC trên Test Data nếu mô hình Logistic chỉ sử dụng chỉ một biến số sau khi đã được rời rạc hóa. Dưới đây là R codes:

# Remove NetIncomeFlag column: 
data %>% select(-NetIncomeFlag) -> df

# Set response and predictors: 

response <- "Bankrupt"

predictors <- names(df %>% select(-Bankrupt))

# Split our data: 

library(caret)

set.seed(1)

id <- createDataPartition(y = df %>% pull(response), p = 0.7, list = FALSE)

# 70% data for training và validation: 
train_valid <- df[id, ] 

# 30% data will be used for evaluating model performance: 

df_test <- df[-id, ] # Test data.  

set.seed(1)

id_new <- createDataPartition(y = train_valid %>% pull(response), p = nrow(df_test) / nrow(train_valid), list = FALSE)

df_train <- train_valid[-id_new, ] # Train data. 

df_valid <- train_valid[id_new, ] # Validation data. 

# Discretize selected variables: 

n_var <- length(predictors)

auc_space <- NULL

iv_space <- NULL

library(pROC)

for (j in 1:n_var) {
  
  var_j <- predictors[j]
  
  var_j_new <- str_c(var_j, "_woe")
  
  x <- df_train %>% pull(var_j)
  
  y <- df_train$Bankrupt
  
  bin_result_j <- iso_bin(x, y)
  
  # Calculate IV: 
  
  bin_result_j$tbl %>% pull(iv) %>% sum() -> iv_j
  
  iv_space <- c(iv_space, iv_j)
  
  # Discretize variable: 
  
  cal_woe(x, bin_result_j) -> var_woe_j
  
  #------------------------------------
  # Train Logistic on Discretized Data
  #------------------------------------
  
  df_train %>% 
    mutate(Bankrupt = case_when(Bankrupt == 1 ~ "Bankrupt", TRUE ~ "NonBankrupt") %>% as.factor()) %>% 
    pull(Bankrupt) -> label_train
  
  df_j <- data.frame(Bankrupt = label_train, var_woe = var_woe_j)
  
  logit <- glm(Bankrupt ~., family = "binomial", data = df_j)
  
  #-----------------------
  # PD on Validation data
  #-----------------------
  
  df_valid %>% 
    mutate(Bankrupt = case_when(Bankrupt == 1 ~ "Bankrupt", TRUE ~ "NonBankrupt") %>% as.factor()) %>% 
    pull(Bankrupt) -> label_valid
  
  df_valid_j <- data.frame(var_woe =  cal_woe(df_valid %>% pull(var_j), bin_result_j))
  
  prob_pred <- predict(logit, df_valid_j, type = "response")
  
  # AUC on Validation Data: 
  
  my_auc <- roc(label_valid, prob_pred)$auc %>% as.numeric()
  
  auc_space <- c(auc_space, my_auc)
  
}


df_iv_single_var <- data.frame(variable = predictors, iv = iv_space, auc = auc_space)

# Some Results: 

df_iv_single_var %>% 
  arrange(-iv) %>% 
  head() %>% 
  knitr::kable()

variable	iv	auc
PersistentEPSintheLastFourSeasons	3.2601	0.8932600
NetprofitbeforetaxPaidincapital	3.1445	0.8844458
PerShareNetprofitbeforetaxYuan	3.0799	0.8849536
NetIncometoTotalAssets	3.0300	0.8577858
RetainedEarningstoTotalAssets	2.9773	0.8638168
ROACbeforeinterestanddepreciationbeforeinterest	2.8790	0.8646053

Ở đây IV dựa trên MOB được xếp theo thứ tự giảm dần. Biến PersistentEPSintheLastFourSeasons có IV cao nhất là 3.2601 và nếu ta chỉ sử dụng biến số này cho mô hình Logistic thì AUC đạt được trên Test Data sẽ là 0.8932600.

Figure 1 dưới đây chỉ ra ngưỡng IV được chọn để có AUC tối ưu trên Test Data (điểm màu đỏ):

# Function discretizes variable (for Train, Validation and Test Data): 

discretize_fun <- function(variable_selected) {
  
  x <- df_train %>% pull(variable_selected)
  
  y <- df_train$Bankrupt
  
  bin_result <- iso_bin(x, y)
  
  cal_woe(x, bin_result) -> var_woe_train
  
  cal_woe(df_valid %>% pull(variable_selected), bin_result) -> var_woe_valid
  
  cal_woe(df_test %>% pull(variable_selected), bin_result) -> var_woe_test
  
  var_woe <- c(var_woe_train, var_woe_valid, var_woe_test)
  
  data_set <- c(rep("train", nrow(df_train)), rep("validation", nrow(df_valid)), rep("test", nrow(df_test)))
  
  data.frame(var_woe = var_woe, data_set = data_set) -> df_result
  
  return(df_result %>% mutate(original_var = variable_selected))
  
}


do.call("bind_rows", lapply(predictors, discretize_fun)) -> df_woe

df_woe %>% 
  filter(data_set == "train") %>% 
  pull(var_woe) %>% 
  matrix(ncol = length(predictors), byrow = FALSE) %>% 
  as.data.frame() %>% 
  set_names(predictors) %>% 
  mutate(Bankrupt = case_when(df_train$Bankrupt == 1 ~ "Bankrupt", TRUE ~ "NonBankrupt") %>% as.factor()) -> df_train_woe


df_woe %>% 
  filter(data_set == "validation") %>% 
  pull(var_woe) %>% 
  matrix(ncol = length(predictors), byrow = FALSE) %>% 
  as.data.frame() %>% 
  set_names(predictors) %>% 
  mutate(Bankrupt = case_when(df_valid$Bankrupt == 1 ~ "Bankrupt", TRUE ~ "NonBankrupt") %>% as.factor()) -> df_validation_woe


df_woe %>% 
  filter(data_set == "test") %>% 
  pull(var_woe) %>% 
  matrix(ncol = length(predictors), byrow = FALSE) %>% 
  as.data.frame() %>% 
  set_names(predictors) %>% 
  mutate(Bankrupt = case_when(df_test$Bankrupt == 1 ~ "Bankrupt", TRUE ~ "NonBankrupt") %>% as.factor()) -> df_test_woe


returnROC_AUC_ValidData <- function(predictor_selected) {
  
  f <- as.formula(paste0(response, " ~ ", paste(predictor_selected, collapse = " + ")))
  
  logit <- glm(f, family = "binomial", data = df_train_woe)
  
  prob_pred <- predict(logit, df_validation_woe, type = "response")
  
  my_auc <- roc(df_validation_woe$Bankrupt, prob_pred)$auc %>% as.numeric()
  
  return(my_auc)
  
}

# Set a sequence of thresholds: 

iv_thresholds <- seq(min(df_iv_single_var$iv), max(df_iv_single_var$iv), 0.05)

auc_space <- NULL

# AUC by threshold: 

for (j in iv_thresholds) {
  
  df_iv_single_var %>% 
    filter(iv >= j) %>% 
    pull(variable) -> predictors_for_modelling
  
  returnROC_AUC_ValidData(predictors_for_modelling) -> my_auc
  
  auc_space <- c(auc_space, my_auc)
  
}

tibble(iv_thresholds = iv_thresholds, auc = auc_space) -> df_auc_threshold

# Features create max AUC on validation data: 

auc_max <- df_auc_threshold %>% slice(which.max(auc))

df_auc_threshold %>% 
  ggplot(aes(iv_thresholds, auc)) + 
  geom_line(size = 1, color = "blue") + 
  geom_point(data = auc_max, aes(iv_thresholds, auc), color = "red", size = 2) + 
  labs(x = "IV Threshold", 
       y = "AUC on Valid Data", 
       title = "Figure 1: AUC on Validation Data by Information Value Threshold")

Như vậy nếu chọn các biến số mà có IV >= 1.14 thì AUC trên Validation Data sẽ là 0.938:

auc_max

## # A tibble: 1 x 2
##   iv_thresholds   auc
##           <dbl> <dbl>
## 1          1.14 0.938

Dưới đây là danh sách các biến thỏa mãn IV >= 1.14:

df_iv_single_var %>% 
  filter(iv >= auc_max$iv_thresholds) %>% 
  pull(variable) -> var_auc_938

df_iv_single_var %>% 
  filter(iv >= auc_max$iv_thresholds) %>% 
  select(-auc) %>% 
  knitr::kable()

variable	iv
ROACbeforeinterestanddepreciationbeforeinterest	2.8790
ROAAbeforeinterestandaftertax	2.4156
ROABbeforeinterestanddepreciationaftertax	2.7952
OperatingProfitRate	1.4413
PretaxnetInterestRate	2.1303
AftertaxnetInterestRate	1.9496
Nonindustryincomeandexpenditurerevenue	2.0579
Continuousinterestrateaftertax	2.4234
TaxrateA	1.3241
NetValuePerShareB	2.2337
NetValuePerShareA	2.2157
NetValuePerShareC	2.1619
PersistentEPSintheLastFourSeasons	3.2601
OperatingProfitPerShareYuan	1.8248
PerShareNetprofitbeforetaxYuan	3.0799
NetValueGrowthRate	2.3407
CurrentRatio	1.4081
QuickRatio	1.8383
InterestExpenseRatio	1.3197
TotaldebtTotalnetworth	2.0893
Debtratio	2.0640
NetworthAssets	2.0640
Borrowingdependency	2.3064
OperatingprofitPaidincapital	1.8207
NetprofitbeforetaxPaidincapital	3.1445
Operatingprofitperperson	1.5079
WorkingCapitaltoTotalAssets	1.2208
CashTotalAssets	1.1772
QuickAssetsCurrentLiability	1.4090
OperatingFundstoLiability	1.2125
CurrentLiabilitiesEquity	1.5337
RetainedEarningstoTotalAssets	2.9773
TotalincomeTotalexpense	2.4484
CurrentLiabilitytoEquity	1.5337
CurrentLiabilitytoCurrentAssets	1.4016
NetIncometoTotalAssets	3.0300
NetIncometoStockholdersEquity	2.8223
LiabilitytoEquity	1.8260
DegreeofFinancialLeverageDFL	1.4161
InterestCoverageRatioInterestexpensetoEBIT	1.3253
EquitytoLiability	2.0658

Nếu sử dụng 41 các biến số này cho mô hình Logistic thì chúng ta sẽ đạt AUC = 0.9365559 trên Test Data:

f938 <- as.formula(paste0(response, " ~ ", paste(var_auc_938, collapse = " + ")))

logit938 <- glm(f938, family = "binomial", data = df_train_woe)

prob_pred938 <- predict(logit938, df_test_woe, type = "response")

my_auc <- roc(df_test_woe$Bankrupt, prob_pred938)$auc %>% as.numeric()

my_auc

## [1] 0.9365559

Chúng ta có thể cải thiện kết quả này với Stepwise Selection.

Stepwise Selection

R codes dưới đây thực hiện Stepwise Selection cho mô hình Logistic và tính AUC trên Test Data:

# Conduct Stepwise Selection: 
logit_back <- step(logit938, steps = 1000, direction = "backward", trace = 0) 

# AUC on Test Data
roc(df_test_woe$Bankrupt, predict(logit_back, df_test_woe, type = "response"))$auc %>% as.numeric() 

## [1] 0.9383908

Kết quả 0.9383908 chỉ ra rằng bằng cách sử dụng MOB kết hợp với Stepwise Selection chúng ta có thể đạt được kết quả tốt hơn KotaShimomura là Team đang xếp vị trí thứ nhất với AUC = 0.93798.

Summary

1. Sử dụng MOB kết hợp với Stepwise Selection cho mô hình Logistic chúng ta có thể đạt được thứ hạng cao hơn vị trí của Team đang xếp vị trí thứ nhất là KotaShimomura.

2. Đây là dữ liệu bất cân bằng rất cao (chỉ có 3.2% các quan sát là Bankrupt) nên có thể cần xem xét đến khả năng sử dụng các giải pháp resampling dữ liệu như SMOTE, upsampling - downsampling để đạt kết quả tốt hơn nữa.

References

1. Monotone optimal binning algorithm for credit risk modeling.
2. Purposeful selection of variables in logistic regression.
3. Monotonic binning using XGBOOST.
4. An Introduction to Statistical Learning.
5. Optimal binning: mathematical programming formulation.
6. Step away from stepwise.

R Environment and OS

sessionInfo()

## R version 4.1.2 (2021-11-01)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19043)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] pROC_1.18.0     caret_6.0-90    lattice_0.20-45 mob_0.4.2      
##  [5] forcats_0.5.1   stringr_1.4.0   dplyr_1.0.7     purrr_0.3.4    
##  [9] readr_2.1.0     tidyr_1.1.4     tibble_3.1.6    ggplot2_3.3.5  
## [13] tidyverse_1.3.1
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-153         fs_1.5.0             lubridate_1.8.0     
##  [4] httr_1.4.2           tools_4.1.2          backports_1.3.0     
##  [7] bslib_0.3.1          utf8_1.2.2           R6_2.5.1            
## [10] rpart_4.1-15         DBI_1.1.1            colorspace_2.0-2    
## [13] nnet_7.3-16          withr_2.4.3          gbm_2.1.8           
## [16] tidyselect_1.1.1     compiler_4.1.2       cli_3.1.0           
## [19] rvest_1.0.2          xml2_1.3.2           sass_0.4.0          
## [22] scales_1.1.1         digest_0.6.28        rmarkdown_2.11      
## [25] pkgconfig_2.0.3      htmltools_0.5.2      parallelly_1.30.0   
## [28] highr_0.9            dbplyr_2.1.1         fastmap_1.1.0       
## [31] rlang_0.4.12         readxl_1.3.1         rstudioapi_0.13     
## [34] jquerylib_0.1.4      generics_0.1.1       jsonlite_1.7.3      
## [37] ModelMetrics_1.2.2.2 magrittr_2.0.1       Matrix_1.3-4        
## [40] Rcpp_1.0.7           munsell_0.5.0        fansi_0.5.0         
## [43] lifecycle_1.0.1      stringi_1.7.6        yaml_2.2.1          
## [46] MASS_7.3-54          plyr_1.8.6           recipes_0.1.17      
## [49] grid_4.1.2           parallel_4.1.2       listenv_0.8.0       
## [52] crayon_1.4.2         haven_2.4.3          splines_4.1.2       
## [55] hms_1.1.1            knitr_1.36           pillar_1.6.4        
## [58] stats4_4.1.2         reshape2_1.4.4       future.apply_1.8.1  
## [61] codetools_0.2-18     reprex_2.0.1         glue_1.5.0          
## [64] evaluate_0.14        data.table_1.14.2    modelr_0.1.8        
## [67] vctrs_0.3.8          tzdb_0.2.0           foreach_1.5.1       
## [70] cellranger_1.1.0     Rborist_0.2-3        gtable_0.3.0        
## [73] future_1.23.0        assertthat_0.2.1     xfun_0.28           
## [76] gower_0.2.2          prodlim_2019.11.13   broom_0.7.10        
## [79] class_7.3-19         survival_3.2-13      timeDate_3043.102   
## [82] iterators_1.0.13     lava_1.6.10          globals_0.14.0      
## [85] ellipsis_0.3.2       ipred_0.9-12