LogisticModelClassifyHomeEquityBadLoan

Logistic Model

Logistic model là một trong những phương pháp phân loại model phổ biến nhất trong các model phân loại thống kê. Giả định của model về hàm phân phối xác xuất dựa trên hàm sigmoid.

Giả sử Y là một sự kiện ngẫu nhiên nhận các giá trị là [0,1]. Trong đó X là tập không gian các vector ảnh hưởng đến Y. Khi đó xác xuất để P(Y=1|X) được thể hiện qua hàm sigmoid như sau:

\[ P(Y=1|X;\theta) = \frac{1}{1+e^{{\space\theta}^{T}X}} \] Từ phương trình của hàm sigmoid ta suy ra miền giá trị của nó nằm trong (0,1). Trong đó \(\theta ^{T}\) là vector của hệ số ước lượng đã được transpose và X là matrix biến độc lập. Mục đích cuối cùng của ta khi ước lượng vector \(\theta\) là giá trị của hàm sigmoid càng gần 1 với những điểm được ước lượng là 1 và giá trị của hàm sigmoid càng gần 0 với những điểm được ước lượng là 0. Tức là hàm số sau đây phải đạt giá trị là lớn nhất: \[ P(Y_{i}|X_{i};\theta) = p_{i}^{Y_{i}}(1-p_{i})^{1-Y_{i}}\] Với mọi \(i \in {\over{1,n}}\).\(p_{i}\) là xác xuất ước lượng của sự kiện Y=1 xảy ra.

Trên toàn bộ tập training ta có thể tìm hàm giá trị maximum của hàm số:

\[ P(Y|X;\theta) = \prod_{i=1}^{n}\space {p_{i}^{Y_{i}}(1-p_{i})^{1-Y_{i}}}\]

Để thuận lợi cho tính toán đạo hàm. Ta logarit hai vế và đổi dấu để tìm minimum hàm Loss có đạng:

\[ J(\theta) = -\frac{1}{n} * \sum_{i=1}^{n}{(log(1-Y_{i})*(1-p_{i})+log(Y_{i})*p_{i}}) \]

Trong R sẽ xây dựng các thuật toán di chuyển theo Gradient descent ứng với mỗi \(\theta\) như bên dưới:

\[\theta_{i_{1}} := \theta_{i_{0}} - \alpha \frac{\delta J(\theta)}{\delta \theta_{i_{0}}} \] \[...\] \[\theta_{i_{n}} := \theta_{i_{n-1}} - \alpha \frac{\delta J(\theta)}{\delta \theta_{i}} \] Tương tự với các hệ số \(\theta_{j}\) khác.n là số vòng lặp.

Trong đó \(\alpha\) là Learning Rate được sử dụng để xác định giá trị thay đổi của mối bước nhảy trong vòng lặp của chúng ta nhằm mục đích di chuyển đến điểm cực trị địa phương. Learning Rate càng lớn thì càng nhanh hội tụ đến cực trị địa phương, Learning Rate càng nhỏ tốc độ hội tụ càng chậm. Tuy nhiên nếu Learning Rate quá lớn sẽ khiến chúng ta step out khỏi cực trị địa phương sau một vài vòng lặp và hàm Loss sẽ không hội tụ về minimum. Quá trình di chuyển theo Gradient descent của các \(\theta\) được thực hiện đồng thời trong thuật toán. Nghĩa là trong khi \(\theta_{i}\) thay đổi thì \(\theta_{j}\) cũng sẽ thay đổi đồng thời.

Trong R có những package chứa các thuật toán tính toán \(\theta\) chẳng hạn như caret, stats,…. Tuy nhiên các thuật toán này đều không dựa trên cách tiếp cận gradient descent mà \(\theta\) đã được công thức hóa sau khi chứng minh được nghiệm tối ưu hàm maximum loglikelihood. Những tính toán này rất phức tạp nên không trình bày ở đây. Còn đối với thuật toán gradient descent thì đơn giản hơn. Các bạn cũng có thể tự xây dựng cho mình những thuật toán tương tự như thế để giải một bài toán tối ưu bất kì miễn là xác định được hàm Loss.

Model phân loại nợ xấu

Tiếp tục bài trước chúng ta sẽ thực hiện model phân loại nợ xấu và đo lường tỷ lệ dự báo chính xác các hợp đồng nợ xấu đối với model Logistic trên tập dữ liệu hmeq.

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ợ.

• 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.
• REASON: DebtCon = debt consolidation (ghép nợ); HomeImp = home improvement (sửa sang nhà cửa).
• JOB: Công việc.
• 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.

Sử dụng tập training và testing như ở bài phân tích trước.

library(rio)
library(dplyr)
testing_chuanhoa <- import("D:/Rcode/RPub/BankingData/testing_chuanhoa.rds")
training_chuanhoa <- import("D:/Rcode/RPub/BankingData/training_chuanhoa.rds")
testing_logistic <- select(testing_chuanhoa,-STT)
training_logistic <- select(training_chuanhoa,-STT)

Hồi qui model logistic trên tập training. Có nhiều package hỗ trợ hồi qui logistic như stats, carret, …. Tuy nhiên package nguyên thủy và lâu đời nhất là stats, các package sau này sẽ kế thừa phương pháp hồi qui logistic từ package này. Vì vậy chúng ta sẽ tìm hiểu package stats trước:

library(stats)
glm.fit = glm(BAD ~ LOAN+MORTDUE+VALUE+YOJ+DEROG+DELINQ+CLAGE+NINQ     +CLNO+DEBTINC+REASON.LV+JOB.LV, data = training_logistic, family = binomial)

summary(glm.fit)

## 
## Call:
## glm(formula = BAD ~ LOAN + MORTDUE + VALUE + YOJ + DEROG + DELINQ + 
##     CLAGE + NINQ + CLNO + DEBTINC + REASON.LV + JOB.LV, family = binomial, 
##     data = training_logistic)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3156  -0.4987  -0.3486  -0.2147   4.9974  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.87649    0.05213 -35.995  < 2e-16 ***
## LOAN        -0.05386    0.05108  -1.055   0.2916    
## MORTDUE      0.17537    0.07647   2.293   0.0218 *  
## VALUE       -0.03848    0.07750  -0.496   0.6196    
## YOJ         -0.05727    0.04883  -1.173   0.2409    
## DEROG        0.44444    0.04645   9.569  < 2e-16 ***
## DELINQ       0.70222    0.04894  14.349  < 2e-16 ***
## CLAGE       -0.56077    0.05919  -9.474  < 2e-16 ***
## NINQ         0.26907    0.04163   6.464 1.02e-10 ***
## CLNO         0.04130    0.05492   0.752   0.4521    
## DEBTINC     -0.91487    0.04249 -21.532  < 2e-16 ***
## REASON.LV   -0.10324    0.04605  -2.242   0.0250 *  
## JOB.LV       0.22140    0.04851   4.564 5.03e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 4764.2  on 4766  degrees of freedom
## Residual deviance: 3307.4  on 4754  degrees of freedom
## AIC: 3333.4
## 
## Number of Fisher Scoring iterations: 5

Từ hệ số ước lượng của phương trình hồi qui logistic ta có thể thấy các biến có ý nghĩa thống kê trong việc dự báo model logistic bao gồm: MORTDUE,DEROG, DELINQ, CLAGE, NINQ, DEBTINC, REASON.LV, JOB.LV ở mức ý nghĩa 0.05%.

Do đó hồi qui lại model logistic theo các biến này:

glm.fit = glm(BAD ~ MORTDUE+DEROG+DELINQ+CLAGE+NINQ+DEBTINC+REASON.LV+ JOB.LV, data = training_logistic, family = binomial)

summary(glm.fit)

## 
## Call:
## glm(formula = BAD ~ MORTDUE + DEROG + DELINQ + CLAGE + NINQ + 
##     DEBTINC + REASON.LV + JOB.LV, family = binomial, data = training_logistic)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3038  -0.5018  -0.3488  -0.2183   5.0051  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.87249    0.05194 -36.052  < 2e-16 ***
## MORTDUE      0.15144    0.04956   3.056  0.00224 ** 
## DEROG        0.45011    0.04614   9.756  < 2e-16 ***
## DELINQ       0.71128    0.04814  14.776  < 2e-16 ***
## CLAGE       -0.56392    0.05151 -10.947  < 2e-16 ***
## NINQ         0.27179    0.04107   6.617 3.66e-11 ***
## DEBTINC     -0.91662    0.04232 -21.660  < 2e-16 ***
## REASON.LV   -0.10120    0.04463  -2.268  0.02335 *  
## JOB.LV       0.21835    0.04818   4.532 5.83e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 4764.2  on 4766  degrees of freedom
## Residual deviance: 3311.5  on 4758  degrees of freedom
## AIC: 3329.5
## 
## Number of Fisher Scoring iterations: 5

Biến DEROG có mức ý nghĩa thống kê nhỏ nhất trong tập dữ liệu của chúng ta. Biến này có giá trị dương có ý nghĩa rằng khi DEROG (số lượng báo cáo vỡ nợ) tăng thì khả năng khách hàng tiếp tục vỡ nợ sẽ tăng cao.Tương tự là biến DELINQ (giá trị tín dụng chưa thanh toán) cũng cho thấy khi lượng tiền chưa thanh toán càng lớn thì xác xuất tiếp tục vỡ nợ càng tăng. Trái lại đối với NINQ là số lượng yêu cầu tín dụng gần đây càng nhỏ thì khả năng vỡ nợ càng thấp.

Sử dụng hàm predict với type = “response” để trả về giá trị xác xuất được tính toán từ hàm sigmoid.

glm.probs <- predict(glm.fit, type = "response") 
glm.probs %>% head()

##          1          2          3          4          5          6 
## 0.48040833 0.74128726 0.38848297 0.48912279 0.09438458 0.95726376

Nếu coi 0.5 là ngưỡng phân biệt hợp đồng vỡ nợ hay không vỡ nợ (>=0.5 BAD, < 0.5 GOOD) ta sẽ gán nhãn cho glm.pred dựa trên giá trị của gml.props.

glm.pred <- rep(1,length(glm.probs))
glm.pred[glm.probs < 0.5] <- 0

Tỷ lệ dự báo chính xác:

library(gmodels)
CrossTable(glm.pred,training_logistic$BAD)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## | Chi-square contribution |
## |           N / Row Total |
## |           N / Col Total |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  4767 
## 
##  
##              | training_logistic$BAD 
##     glm.pred |         0 |         1 | Row Total | 
## -------------|-----------|-----------|-----------|
##            0 |      3686 |       481 |      4167 | 
##              |    36.787 |   147.614 |           | 
##              |     0.885 |     0.115 |     0.874 | 
##              |     0.966 |     0.506 |           | 
##              |     0.773 |     0.101 |           | 
## -------------|-----------|-----------|-----------|
##            1 |       130 |       470 |       600 | 
##              |   255.488 |  1025.177 |           | 
##              |     0.217 |     0.783 |     0.126 | 
##              |     0.034 |     0.494 |           | 
##              |     0.027 |     0.099 |           | 
## -------------|-----------|-----------|-----------|
## Column Total |      3816 |       951 |      4767 | 
##              |     0.801 |     0.199 |           | 
## -------------|-----------|-----------|-----------|
## 
## 

Với ngưỡng phân chia hợp đồng nợ xấu là 0.5 thì tỷ lệ nợ xấu phân loại đúng là 49.4%. Luôn có một sự đánh đổi giữa tỷ lệ phân loại chính xác hợp đồng nợ xấu và hợp đồng tốt. Nếu muốn phân loại chính xác 100% hợp đồng nợ xấu ta chỉ cần hạ ngưỡng phân chia xuống 0 và tất cả các hợp đồng đều được phân loại là xấu. Khi đó chỉ xét riêng các hợp đồng nợ xấu thì 100% sẽ phân loại chính xác. Tuy nhiên đối với các hợp đồng tốt tỷ lệ chính xác sẽ trở về 0% do tất cả bị phân loại thành xấu. Như vậy luôn có sự đánh đổi giữa 2 tỷ lệ này. Đây được gọi là tỷ lệ đánh đổi giữa TPR (tỷ lệ phân loại đúng trường hợp tốt) và TNR (tỷ lệ phân loại đúng trường hợp xấu). Do vậy để xác định đâu là ngưỡng phân chia hợp lý nhất thì ta sẽ căn cứ trên tổng tỷ lệ dự báo chính xác cả 2 loại này. Dưới đây ta sẽ xây dựng vòng lặp xác định ngưỡng tối ưu:

df_Accuracy <- data.frame()
for(i in 1:19){
  Threshold <- 0.05*i
  glm.pred <- rep(1,length(glm.probs))
  glm.pred[glm.probs < Threshold] <- 0  
  tbl <- table(glm.pred,training_logistic$BAD)/length(glm.pred)
  Accuracy <- tbl[1,1]+tbl[2,2]
  df_Accuracy <- rbind(df_Accuracy,c(Threshold,Accuracy))
  colnames(df_Accuracy) <- c("Threshold","Accuracy")
}
df_Accuracy

##    Threshold  Accuracy
## 1       0.05 0.4172435
## 2       0.10 0.6612125
## 3       0.15 0.7740717
## 4       0.20 0.8132998
## 5       0.25 0.8313405
## 6       0.30 0.8487518
## 7       0.35 0.8588211
## 8       0.40 0.8659534
## 9       0.45 0.8688903
## 10      0.50 0.8718271
## 11      0.55 0.8682610
## 12      0.60 0.8623872
## 13      0.65 0.8539962
## 14      0.70 0.8466541
## 15      0.75 0.8426683
## 16      0.80 0.8395217
## 17      0.85 0.8319698
## 18      0.90 0.8248374
## 19      0.95 0.8164464

Ta thấy ngưỡng phân loại có mức độ chính xác cao nhất là̀ 0.5. Giá trị của Threshold và Accuracy có mối quan hệ dạng hàm lõm.

Visualization:

library(ggplot2)
library(dplyr)
df_Accuracy %>% ggplot()+
  geom_line(aes(Threshold,Accuracy)) +
  labs(
    title="Accuracy vs Threshold",
    subtitle="Created by Pham Dinh Khanh"
  ) 

Sử dụng ngưỡng phân chia là 0.5 cho phân loại trên tập test.

glm.probs <- predict(glm.fit,testing_logistic,type = "response")
glm.probs %>% head()

##          1          2          3          4          5          6 
## 0.44022028 0.50722019 0.41353407 0.01830806 0.46107861 0.26859001

df_Accuracy_Testing <- data.frame()
for(i in 1:19){
  Threshold <- 0.05*i
  glm.pred <- rep(1,length(glm.probs))
  glm.pred[glm.probs < Threshold] <- 0  
  tbl <- table(glm.pred,testing_logistic$BAD)/length(glm.pred)
  Accuracy <- tbl[1,1]+tbl[2,2]
  df_Accuracy_Testing <- rbind(df_Accuracy_Testing ,c(Threshold,Accuracy))
  colnames(df_Accuracy_Testing) <- c("Threshold","Accuracy")
}
df_Accuracy_Testing

##    Threshold  Accuracy
## 1       0.05 0.4266555
## 2       0.10 0.6831517
## 3       0.15 0.7862531
## 4       0.20 0.8197821
## 5       0.25 0.8407376
## 6       0.30 0.8642079
## 7       0.35 0.8725901
## 8       0.40 0.8717519
## 9       0.45 0.8717519
## 10      0.50 0.8633697
## 11      0.55 0.8650461
## 12      0.60 0.8600168
## 13      0.65 0.8549874
## 14      0.70 0.8466052
## 15      0.75 0.8424141
## 16      0.80 0.8382230
## 17      0.85 0.8281643
## 18      0.90 0.8248114
## 19      0.95 0.8147527

Sử dụng kiểm định t-test kiểm tra sự khác biệt về Accuracy giữa tập train và test ứng với mỗi threshold:

t.test(df_Accuracy$Accuracy,df_Accuracy_Testing$Accuracy)

## 
##  Welch Two Sample t-test
## 
## data:  df_Accuracy$Accuracy and df_Accuracy_Testing$Accuracy
## t = -0.11906, df = 35.969, p-value = 0.9059
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.07320037  0.06508263
## sample estimates:
## mean of x mean of y 
## 0.8104292 0.8144880

P-value > 0.05 cho thấy không có sự khác biệt về mức độ chính xác của 2 chuỗi này. Điều đó cho thấy tỷ lệ chính xác trên tập train và test khá khớp nhau. Vẽ đồ thị đường Threshold-Accuracy của tập train và test

library(data.table)
df <- merge(df_Accuracy,df_Accuracy_Testing, by="Threshold")
colnames(df) <- c("Threshold","training","testing")
df <- df %>% melt(id=c("Threshold"))
df %>% ggplot()+
  geom_line(aes(Threshold,value,color = variable)) +
  labs(
    title="Accuracy vs Threshold",
    subtitle="Created by Pham Dinh Khanh"
  ) 

Kết quả của model dự báo trên tập test với ngưỡng phân chia là 0.5.

glm.pred <- rep(1,length(glm.probs))
  glm.pred[glm.probs < 0.5] <- 0  
  tbl <- table(glm.pred,testing_logistic$BAD)/length(glm.pred)
  tbl

##         
## glm.pred          0          1
##        0 0.77451802 0.11064543
##        1 0.02598491 0.08885163

So với phương pháp Knn với mức độ chính xác dự báo của cả 2 loại gồm hợp đồng GOOD và BAD là 87% thì Logistic có độ chính xác ngang bằng ở ngưỡng phân loại là 0.5. Tuy nhiên mức độ dự báo chính xác hợp đồng nợ xấu của Knn lại cao hơn còn mức độ dự báo chính xác hợp đồng tốt của Logistic lại cao hơn. Điều này xuất phát từ 2 nguyên nhân:

• cách tiếp cận và lựa chọn biến của 2 phương pháp là khác nhau. Knn dựa trên mức độ tương quan tuyến tính còn Logistic dựa trên ý nghĩa thống kê của hệ số ước lượng.

• Knn là phương pháp non-parametric dựa trên tính toán khoảng cách mà không đánh giá tác động của các biến lên sự phân loại của model trong khi đó Logistic là phương pháp parametric thể hiện sự tường mình hơn qua phương trình sigmoid ước lượng xác xuất. Các biến dự báo của Logistic có thể hiện tác động rõ ràng lên kết quả dự báo.

Thông thường thì các phương pháp non-parametric dễ gặp phải hiện tượng overfiting hơn (tức là trên tập train thì dự báo tốt nhưng trên tập test thì dự báo không chuẩn xác) so với các phương pháp parametric nên trong dự báo chúng ta cần phải backtest để kiểm tra tính ổn định của model.