Gradient Boosting Tree Home Equity Loan
Pham Dinh Khanh
July 7, 2018
1 Giới thiệu tập dữ liệu
Giải thích dữ liệu hmeq. HMEQ là bộ dữ liệu về các khoản vay có tài sản đảm bảo là nhà. Dữ liệu bao gồm 5960 khoản vay được phân loại thành các nhóm vỡ nợ hoặc không vỡ nợ. Trong đó bao gồm 2 kiểu dữ liệu dạng numeric và dạng category.
Các biến dữ liệu numeric:
- BAD: 1 = Khoản vay bị vỡ nợ; 0 = ứng viên đã trả nợ.
- LOAN: Giá trị khoản vay.
- MORTDUE: Giá trị tài sản đảm bảo được định giá.
- VALUE: Giá trị tài sản hiện tại.
- YOJ: Số năm trong nghề.
- DEROG:Số lượng báo cáo vỡ nợ.
- DELINQ: Giá trị tín dụng chưa thanh toán.
- CLAGE: Độ tuổi tín dụng lâu nhất trong tháng.
- NINQ: Số lượng yêu cầu tín dụng gần đây.
- CLNO:Số lượng hạn mức tín dụng.
- DEBTINC:Tỷ lệ nợ/thu nhập.
Các biến dữ liệu category:
- REASON: DebtCon = debt consolidation (ghép nợ); HomeImp = home improvement (sửa sang nhà cửa).
- JOB: loại công việc, bao gồm các nhóm “Office”, “Sales”, “Mananger”, “Professional Executive”, “Self business” và các công việc khác.
2 Data preparation
Kiểm tra tỷ lệ các hợp đồng nợ xấu:
hmeq %>% group_by(BAD) %>%
summarise(n = n()/nrow(.)*100)Ta có thể thấy dữ liệu không có mất cân bằng giữa các nhóm nợ xấu và thông thường. Do đó chúng ta có thể sử dụng accuracy và đồ thị ROI làm tiêu chuẩn đánh giá mức độ dự báo chính xác của mô hình. Thống kê số lượng các quan sát NA.
library(ggplot2)
dfNA <- hmeq %>% apply(2, is.na) %>% apply(2, sum) %>% as.matrix() %>% as.data.frame()
dfNA$Var <- rownames(dfNA)
dfNA$Percent <- dfNA$V1/nrow(hmeq)*100
dfNA %>% ggplot(aes(x = reorder(Var, -desc(Percent)), y = Percent)) +
geom_col(fill = 'blue') +
labs(
x = 'Percent',
y = 'Variable',
title = 'Missing value percent'
) +
coord_flip()
Có nhiều cách để xử lý missing value như điền giá trị missing bằng mean, median, giá trị có tần xuất xuất hiện cao nhất, 0, hoặc loại bỏ missing value. Các được sử dụng phổ biến là fill missing value bằng giá trị mean vì đây là giá trị đại diện cho biến nhất. Ta sẽ thực hiện
f_fillNA <- function(x){
replace(x, is.na(x), mean(x, na.rm = TRUE))
}
hmeq <- replace(hmeq, TRUE, lapply(hmeq, f_fillNA))## Warning in mean.default(x, na.rm = TRUE): argument is not numeric or
## logical: returning NA
## Warning in mean.default(x, na.rm = TRUE): argument is not numeric or
## logical: returning NAThống kê lại số lượng các giá trị NA
is.na(hmeq) %>% sum()## [1] 0Đối với các biến category sẽ không thể đưa trực tiếp vào hồi qui mà phải chuyển sang biến numeric. Chúng ta sẽ có 2 cách để tạo biến mới. Đó là:
- Label coding: Gán giá trị cho các category một số nguyên.
- One-hot coding: Tạo thành các biến mới sao cho mỗi một category sẽ tương ứng với một biến và khi giá trị rơi vào category nào thì biến ứng với category đó sẽ nhận giá trị bằng 1 và các biến còn lại nhận giá trị 0.
Cách 2 có ưu điểm hơn cách 1 bởi ta không biết trước nhóm nào sẽ có ảnh hưởng tới biến target mạnh hơn. Do đó nếu gán các giá trị nguyên một cách ngẫu nhiên cho các category thì sẽ không chính xác. Theo cách 2 sẽ đảm bảo được sự công bằng giữa các category vì các category này đều được gán là 1 khi xảy ra và không xảy ra thì gán là 0. Độ lớn của sự tác động của từng category sẽ được thể hiện qua hệ số hồi qui khác nhau trong phương trình hồi qui. Do đó ta sẽ thực hiện tạo biến one-hot coding (rơi vào 1 thì hot) theo cách 2.
Trước tiên ta thống kê số lượng của từng nhóm
hmeq %>% group_by(REASON) %>%
summarise(n = 100*n()/nrow(hmeq))Do số lượng các REASON bị blank khá nhỏ nên ta có thể gán giá trị có tần xuất lớn nhất là DebtCon vào nhóm này. Sau đó tiến hành tạo biến one-hot coding
library(tidyr)
library(onehot)
hmeq$REASON[hmeq$REASON == ""] <- "DebtCon"
hmeq %>% group_by(REASON) %>%
summarise(n = 100*n()/nrow(hmeq))hmeq$REASON1 <- 0
hmeq$REASON1[hmeq$REASON == 'DebtCon'] <- 1
hmeq$REASON2 <- 0
hmeq$REASON2[hmeq$REASON == 'HomeImp'] <- 1Biểu diễn phân phối GOOD, BAD theo từng category của JOB
library(reshape2)##
## Attaching package: 'reshape2'## The following objects are masked from 'package:data.table':
##
## dcast, melt## The following object is masked from 'package:tidyr':
##
## smithsdfJOB <- hmeq %>% group_by(JOB, BAD) %>%
summarise(n = n()) %>% acast(JOB ~ BAD) %>%
as.data.frame(row.names(rownames(.))) %>%
mutate(PercentBAD = `1`/(`0`+`1`),
JOB = rownames(.)) %>% arrange(desc(PercentBAD)) ## Using n as value column: use value.var to override.dfJOBdfJOB %>%
ggplot(aes(x = reorder(JOB,PercentBAD), y = PercentBAD)) +
geom_col(fill = 'blue') +
labs(
x = 'JOB',
y = 'percentageBAD',
title = 'PercentageBAD in each group JOB'
) +
coord_flip()
Ta thấy tỷ lệ %BAD tại một số nhóm ngành nghề là gần bằng nhau. Trong đó có một số nhóm khá nhỏ như Sales, Self, Office, blank ta có thể gộp chung vào các nhóm lớn còn lại có tỷ lệ này là gần nhất để các nhóm không bị mất cân bằng số lượng. Cụ thể có thể gộp: Sales và Self; Mgr và Other; ProfExe, Office và blank.
hmeq$JOB1 <- 0
hmeq$JOB1[hmeq$JOB %in% c('Sales','Self')] <- 1
hmeq$JOB2 <- 0
hmeq$JOB2[hmeq$JOB %in% c('Mgr','Other')] <- 1
hmeq$JOB3 <- 0
hmeq$JOB3[hmeq$JOB %in% c('ProfExe','Office','')] <- 1Sau khi hoàn thiện xử lý mẫu ta sẽ tiến hành feature rescaling các biến. Thông thường có 3 phương pháp feature rescaling:
- Standardization. Biến được rescaling theo kì vọng và độ lệch chuẩn như sau: \[\mathbf{x} = \frac{\mathbf{x}-\overline{\mathbf{x}}}{\sigma(\mathbf{x})}\]
- Rescaling theo range. Biến được đưa về các range [0,1] hoặc [-1,1] mà vẫn giữa nguyên tỷ lệ khoảng cách giữa các biến. Chẳng hạn khi feature theo khoảng [0,1]: \[\mathbf{x'} = \frac{\mathbf{x}-\min({\mathbf{x})}}{\max({\mathbf{x}}) - \min({\mathbf{x}})}\] Khi biến được feature theo khoảng [-1, 1] công thức sẽ như sau:
\[\mathbf{x'} = \frac{2(\mathbf{x}-{\overline{\mathbf{x}}})}{\max({\mathbf{x}}) - \min({\mathbf{x}})}\]
- Rescaling về unit length. Khi đó độ dài khoảng cách norm2 bằng 1 (tức trong không gian Euclid vecor sẽ có độ lớn bằng 1).
\[\mathbf{x'} = \frac{\mathbf{x}}{||\mathbf{x}||_2}\]
Cách 1 sẽ có nhược điểm là đối với các data point là oulier thì giá trị của nó sẽ cách range scaling rất xa. Khi đó rescaling sẽ không đảm bảo các biến nằm trọn trong 1 range đủ nhỏ để thuật toán gradient hội tụ nhanh hơn. Cách 2 theo sẽ khác phục được hạn chế này, tuy nhiên nếu lấy range nằm trong khoảng [0, 1] thì không giữ được sự biến đổi âm dương của biến. Do đó thông thường khi chọn cách 2 ta sẽ lấy range thuộc khoảng [-1, 1]. Cách 3 sẽ tạo ra biến mới có tổng bình phương các biến bằng 1 trong khi cách 2 chỉ giới hạn các biến nằm trong khoảng [-1, 1] (trong trường hợp range là [-1, 1]) nên sự biến động của các biến của cách 3 nhỏ hơn rất nhiều so với cách 2 dẫn tới kết quả phân loại của mô hình cũng thấp hơn.
Tổng kết lại chúng ta sẽ sử dụng cách 2 theo những tiêu chí biến động theo chiều âm/dương trong range [-1,1] và sự biến động của biến đủ lớn để phân loại biến target. Trong caret ta sử dụng hàm preProcess để scale các biến:
library(caret)## Loading required package: lattice##
## Attaching package: 'caret'## The following object is masked from 'package:purrr':
##
## lift#Lọc riêng biến target
y <- hmeq$BAD
#Lấy những biến numeric và lọc biến target
hmeq <- hmeq %>% select_if(is.numeric) %>% select(-BAD)
#Rescale biến
preProc <- preProcess(hmeq, method = 'range', rangeBounds = c(-1, 1))
hmeq <- predict(preProc, hmeq)
head(hmeq)Để việc training data được nhanh ta chuyển bảng sang dạng data.matrix
hmeq <- hmeq %>% data.matrix()3 Training model
Chúng ta sẽ sử dụng phương pháp Gradient Boosting để training model. Phương pháp này hiện nay đang dẫn đầu về mức độ chính xác trong các thuật toán supervised learning và thường được sử dụng trong rất nhiều các cuộc thi phân tích dữ liệu như kaggle, KDnuggets, DrivenData và cho kết qủa thuộc top đầu. Về lý thuyết thì thuật toán là xây dựng một decision tree và tìm cách để cải thiện độ chính xác tại mỗi nhánh của decision tree bằng các thuật toán gradient boosting. Quá trình này được lặp lại nhiều lần vào sau mỗi một vòng lặp thì metric của model (có thể là loss function hoặc accuracy) sẽ được cải thiện.
Đầu tiên ta tiến hành chia mẫu train, test. Sử dụng hàm createDataPartition để chia mẫu dự báo theo tỷ lệ 9:1 với điều kiện cân bằng tỷ lệ này tại các group của biến target.
library(xgboost)##
## Attaching package: 'xgboost'## The following object is masked from 'package:dplyr':
##
## sliceset.seed(1)
id_train <- caret::createDataPartition(y, p = 0.9, list = F) %>% c()
#Khởi tạo data.matrix train và test
dtrain <- xgb.DMatrix(data = hmeq[id_train, ], label = y[id_train])
dtest <- xgb.DMatrix(data = hmeq[-id_train, ], label = y[-id_train])Khởi tạo list các parameter như bên dưới:
p <- list(objective = "binary:logistic",
booster = "gbtree",
eval_metric = "auc",
nthread = 7,
eta = 0.05,
max_depth = 6,
min_child_weight = 30,
gamma = 0,
subsample = 0.85,
colsample_bytree = 0.7,
colsample_bylevel = 0.632,
alpha = 0,
lambda = 0,
nrounds = 2000)Ở đây ta quan tâm đến các thông số chính:
objective: Phương thức được sử dụng để dự báo model. Do model thuộc lớp phân loại binary nên chúng ta chọn binary:logistic.
booster: Thuật toán sử dụng để tăng cường độ chính xác cho model (boost model). Có 2 thuật toán chính là gbtree (gradient boosting tree), gblinear (gradient boosting linear). Mỗi thuật toán sẽ sử dụng một phương pháp khác nhau để tối ưu hàm loss function. Tìm hiểu sâu hơn về thuật toán chúng ta có thể xem ở link gradient boosting tree.
eval_metric: evaluation metric sử dụng để tối ưu model. Chúng ta chọn trong trường hợp này là auc (area under curve) là diện tích của đường cong auc sử dụng để đánh giá sức mạnh phân loại model.
nthread: xgboost cho phép thuật toán xử lý song song và số core được qui định trong tham số nthread.
max_depth: Chiều sâu của decision tree.
min_child_weight: Tổng trọng số nhỏ nhất của một child node. Trong quá trình processing model nếu tổng trọng số của child node nhỏ hơn min_child_weight thì node sẽ không bị phân chia nhỏ hơn nữa.
eta, gamma, alpha, lambda: Các tham số để tối ưu hóa hàm loss function.
nrounds: Số vòng lặp tối đa của thuật toán
subsample: Tỷ lệ mẫu được lấy ra từ tập training để phát triển decision tree.
colsample_bytree, colsample_bylevel: Tỷ lệ lấy mẫu của mỗi biến khi khởi tạo một tree mới.
Cuối cùng ta sẽ training thuật toán dựa trên các parmeter vừa thiết lập. Ở đây để tránh hiện tượng overfiting ta sẽ train model trên tập dtrain và đánh giá metric auc trên tập dtest. Chúng ta không đánh giá mức độ dự báo trên chính tập dtrain vì mục tiêu phân loại là nâng cao khả năng dự báo trên tập dữ liệu chưa biết là dtest.
set.seed(0)
m_xgb <- xgb.train(p, dtrain, p$nrounds, list(val = dtest),
#print metric sau mỗi 50 vòng lặp
print_every_n = 50,
#dừng lại sau 1000 rounds nếu model vẫn không được cải thiện để tránh hiện tượng vượt quá global optimal value.
early_stopping_rounds = 1000)## [1] val-auc:0.832698
## Will train until val_auc hasn't improved in 1000 rounds.
##
## [51] val-auc:0.883292
## [101] val-auc:0.891688
## [151] val-auc:0.898708
## [201] val-auc:0.903479
## [251] val-auc:0.908375
## [301] val-auc:0.911208
## [351] val-auc:0.914062
## [401] val-auc:0.916333
## [451] val-auc:0.917771
## [501] val-auc:0.919354
## [551] val-auc:0.920833
## [601] val-auc:0.921562
## [651] val-auc:0.921125
## [701] val-auc:0.923146
## [751] val-auc:0.923063
## [801] val-auc:0.923375
## [851] val-auc:0.924104
## [901] val-auc:0.924667
## [951] val-auc:0.926188
## [1001] val-auc:0.926854
## [1051] val-auc:0.927542
## [1101] val-auc:0.928000
## [1151] val-auc:0.928521
## [1201] val-auc:0.928771
## [1251] val-auc:0.929521
## [1301] val-auc:0.929708
## [1351] val-auc:0.929875
## [1401] val-auc:0.930000
## [1451] val-auc:0.930979
## [1501] val-auc:0.931021
## [1551] val-auc:0.930979
## [1601] val-auc:0.931479
## [1651] val-auc:0.931938
## [1701] val-auc:0.932021
## [1751] val-auc:0.932271
## [1801] val-auc:0.932458
## [1851] val-auc:0.932292
## [1901] val-auc:0.932375
## [1951] val-auc:0.932458
## [2000] val-auc:0.932417Ta thấy tỷ lệ auc của model được cải thiện dần qua thời gian. Ở một số giai đoạn ngắn auc giảm dần nhưng sau đó lại tiếp tục giữ được đà tăng do thuật toán gradient boosting đi qua những điểm local optimal value, đà tăng sẽ giúp thuật toán đi đến được đích cuối cùng là global optimal value. Do training model khá tốn thời gian nên chúng ta có thể save model dưới dạng binary để sử dụng cho các lần sau:
#save model into binary file m_xgb
xgb.save(m_xgb,'m_xgb')## [1] TRUE#load model m_xgb file
m_xgb <- xgb.load('m_xgb')Đánh giá mức độ quan trọng của các biến
cols <- colnames(hmeq)
xgb.importance(cols, model=m_xgb) %>%
xgb.plot.importance(top_n = 30)
Từ biểu đồ trên ta có thể thấy tỷ lệ nợ/thu nhập (DEBTINCOM) vẫn là biến có ảnh hưởng nhất đến khả năng vỡ nợ của khách hàng. Điều này rất phù hợp với thực tế. Sau đó lần lượt là các biến CLAGE (tuổi giao dịch dài nhất trong tháng) và DEBLINQ (giá trị tín dụng chưa thanh toán).
Dự báo xác xuất vỡ nợ trên tập dtest
Probability <- predict(m_xgb, dtest)
Pre_label <- ifelse(Probability < 0.5, 0, 1)
Act_label <- y[-id_train]
confusionMatrix(Pre_label %>% as.factor(), Act_label %>% as.factor())## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 483 30
## 1 17 66
##
## Accuracy : 0.9211
## 95% CI : (0.8965, 0.9415)
## No Information Rate : 0.8389
## P-Value [Acc > NIR] : 2.085e-09
##
## Kappa : 0.6913
## Mcnemar's Test P-Value : 0.08005
##
## Sensitivity : 0.9660
## Specificity : 0.6875
## Pos Pred Value : 0.9415
## Neg Pred Value : 0.7952
## Prevalence : 0.8389
## Detection Rate : 0.8104
## Detection Prevalence : 0.8607
## Balanced Accuracy : 0.8267
##
## 'Positive' Class : 0
## Kiểm tra tỷ lệ dự báo chính xác này trên tập train
Probability_train <- predict(m_xgb, dtrain)
Pre_label_train <- ifelse(Probability_train < 0.5, 0, 1)
Act_label_train <- y[id_train]
confusionMatrix(Pre_label_train %>% as.factor(), Act_label_train %>% as.factor())## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 4229 130
## 1 42 963
##
## Accuracy : 0.9679
## 95% CI : (0.9629, 0.9725)
## No Information Rate : 0.7962
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.8981
## Mcnemar's Test P-Value : 3.274e-11
##
## Sensitivity : 0.9902
## Specificity : 0.8811
## Pos Pred Value : 0.9702
## Neg Pred Value : 0.9582
## Prevalence : 0.7962
## Detection Rate : 0.7884
## Detection Prevalence : 0.8126
## Balanced Accuracy : 0.9356
##
## 'Positive' Class : 0
## Ta thấy trên tập train thì kết quả dự báo khả quan hơn tập test khi accuracy = 96.79%. Tuy nhiên tỷ lệ dự báo chính xác ở tập test cũng khá cao (92.11%) nên có thể khẳng định rằng không có hiện tượng overfiting. Ta hoàn toàn có thể áp dụng model vào thực tế.
4 Tunning Probability Threshold
Quan trọng nhất trong dự báo binary classification là phải dự báo được xác xuất để các sự kiện xảy ra. Xác xuất sẽ làm nền tảng để ta triển khai các phân tích khác tùy theo mục đích nghiên cứu. Đối với ngân hàng nếu xem lợi nhuận là mục tiêu thì có thể tunning lợi nhuận dựa trên sự thay đổi Probability Threshold. Trong khuôn khổ bài này chúng ta chỉ quan tâm đến accuracy. Do đó mục đích các nghiên cứu là tìm ra Probability Threshold mang lại accuracy là cao nhất.
eval_thres <- function(thres, prob, actual){
prob <- case_when(prob >= thres ~ 1,
prob < thres ~ 0)
#Summarise results
cm <- confusionMatrix(prob %>% as.factor(), actual %>% as.factor())
cross_tbl <- cm$table %>%
as.vector() %>%
matrix(ncol = 4) %>%
data.frame()
names(cross_tbl) <- c('Good-Good','Bad-Good','Good-Bad','Bad-Bad')
ten <- c(Threshold = thres, cm$byClass, cm$overall) %>% names()
kq <- c(Threshold = thres,
cm$byClass,
cm$overall) %>%
matrix(ncol = 19) %>%
data.frame()
names(kq) <- ten
kq <- bind_cols(kq, cross_tbl)
kq
}
df_tunningThreshold <- sapply(seq(0.2, 0.8, by = 0.05), eval_thres,
prob = Probability, actual = Act_label) %>% t() %>%
as.data.frame()
library(ggplot2)
theme_set(theme_minimal())
library(hrbrthemes)## NOTE: Either Arial Narrow or Roboto Condensed fonts are *required* to use these themes.## Please use hrbrthemes::import_roboto_condensed() to install Roboto Condensed and## if Arial Narrow is not on your system, please see http://bit.ly/arialnarrowdf_tunningThreshold %>% transmute(Threshold = Threshold %>% unlist(),
Accuracy = Accuracy %>% unlist()) %>%
ggplot(aes(x = Threshold, y = Accuracy)) +
geom_line(color = 'red') +
geom_point(size = 3, color = 'green') +
scale_y_percent() +
scale_x_continuous(breaks = seq(0.05, 0.95, by = 0.05)) +
labs(y = "Accuracy Rate",
title = "Variation of Accuracy Rate Metrics by Threshold",
subtitle = "Data Used: Home Equity from\nhttps://www.kaggle.com/ajay1735/hmeq-data#hmeq.csv")
Ngưỡng xác xuất mang lại accuracy cao nhất là 0.65. Khi đó Chúng ta phân biệt được 60 hợp đồng nợ xấu và 490 hợp đồng tốt chính xác. Kết quả được thể hiện ở bên dưới
df_tunningThreshold %>% filter(Threshold == 0.65) %>% as.matrix()## Threshold Sensitivity Specificity Pos Pred Value Neg Pred Value
## [1,] 0.65 0.98 0.625 0.9315589 0.8571429
## Precision Recall F1 Prevalence Detection Rate
## [1,] 0.9315589 0.98 0.9551657 0.8389262 0.8221477
## Detection Prevalence Balanced Accuracy Accuracy Kappa
## [1,] 0.8825503 0.8025 0.9228188 0.67933
## AccuracyLower AccuracyUpper AccuracyNull AccuracyPValue McnemarPValue
## [1,] 0.8983912 0.9429422 0.8389262 9.119745e-10 0.0002277626
## Good-Good Bad-Good Good-Bad Bad-Bad
## [1,] 490 10 36 60Để thể hiện rằng ngưỡng mà ta lựa chọn là đáng tin cậy trong mọi trường hợp thì chúng ta nên tiến hành cross validation model nhiều lần. Tuy nhiên việc này khá tốn thời gian và xin dành cho bạn đọc. Về ý tưởng của cross validation các bạn có thể tham khảo ở link crosss validation.