AS3-3 Predicting Loan Repayment

rm(list=ls(all=T))
options(digits=4, scipen=12)
library(magrittr)

---

Introduction

議題：議題：使用貸款人的資料，預測他會不會還款

---

1 資料整理 Preparing the Dataset

【1.1 基礎機率】What proportion of the loans in the dataset were not paid in full?

loans = read.csv("data/loans.csv")
str(loans)

'data.frame':   9578 obs. of  14 variables:
 $ credit.policy    : int  1 1 1 1 1 1 1 1 1 1 ...
 $ purpose          : Factor w/ 7 levels "all_other","credit_card",..: 3 2 3 3 2 2 3 1 5 3 ...
 $ int.rate         : num  0.119 0.107 0.136 0.101 0.143 ...
 $ installment      : num  829 228 367 162 103 ...
 $ log.annual.inc   : num  11.4 11.1 10.4 11.4 11.3 ...
 $ dti              : num  19.5 14.3 11.6 8.1 15 ...
 $ fico             : int  737 707 682 712 667 727 667 722 682 707 ...
 $ days.with.cr.line: num  5640 2760 4710 2700 4066 ...
 $ revol.bal        : int  28854 33623 3511 33667 4740 50807 3839 24220 69909 5630 ...
 $ revol.util       : num  52.1 76.7 25.6 73.2 39.5 51 76.8 68.6 51.1 23 ...
 $ inq.last.6mths   : int  0 0 1 1 0 0 0 0 1 1 ...
 $ delinq.2yrs      : int  0 0 0 0 1 0 0 0 0 0 ...
 $ pub.rec          : int  0 0 0 0 0 0 1 0 0 0 ...
 $ not.fully.paid   : int  0 0 0 0 0 0 1 1 0 0 ...

summary(loans)

 credit.policy                 purpose        int.rate      installment   
 Min.   :0.000   all_other         :2331   Min.   :0.060   Min.   : 15.7  
 1st Qu.:1.000   credit_card       :1262   1st Qu.:0.104   1st Qu.:163.8  
 Median :1.000   debt_consolidation:3957   Median :0.122   Median :268.9  
 Mean   :0.805   educational       : 343   Mean   :0.123   Mean   :319.1  
 3rd Qu.:1.000   home_improvement  : 629   3rd Qu.:0.141   3rd Qu.:432.8  
 Max.   :1.000   major_purchase    : 437   Max.   :0.216   Max.   :940.1  
                 small_business    : 619                                  
 log.annual.inc       dti             fico     days.with.cr.line   revol.bal      
 Min.   : 7.55   Min.   : 0.00   Min.   :612   Min.   :  179     Min.   :      0  
 1st Qu.:10.56   1st Qu.: 7.21   1st Qu.:682   1st Qu.: 2820     1st Qu.:   3187  
 Median :10.93   Median :12.66   Median :707   Median : 4140     Median :   8596  
 Mean   :10.93   Mean   :12.61   Mean   :711   Mean   : 4562     Mean   :  16914  
 3rd Qu.:11.29   3rd Qu.:17.95   3rd Qu.:737   3rd Qu.: 5730     3rd Qu.:  18250  
 Max.   :14.53   Max.   :29.96   Max.   :827   Max.   :17640     Max.   :1207359  
 NA's   :4                                     NA's   :29                         
   revol.util    inq.last.6mths   delinq.2yrs        pub.rec      not.fully.paid
 Min.   :  0.0   Min.   : 0.00   Min.   : 0.000   Min.   :0.000   Min.   :0.00  
 1st Qu.: 22.7   1st Qu.: 0.00   1st Qu.: 0.000   1st Qu.:0.000   1st Qu.:0.00  
 Median : 46.4   Median : 1.00   Median : 0.000   Median :0.000   Median :0.00  
 Mean   : 46.9   Mean   : 1.57   Mean   : 0.164   Mean   :0.062   Mean   :0.16  
 3rd Qu.: 71.0   3rd Qu.: 2.00   3rd Qu.: 0.000   3rd Qu.:0.000   3rd Qu.:0.00  
 Max.   :119.0   Max.   :33.00   Max.   :13.000   Max.   :5.000   Max.   :1.00  
 NA's   :62      NA's   :29      NA's   :29       NA's   :29                    

Calculate the proportion of the loans in the dataset were not paid infull.

table(loans$not.fully.paid)

   0    1 
8045 1533 

print("Propotion = "); 1533/(8045+1533)

[1] "Propotion = "
[1] 0.1601

【1.2 檢查缺項】Which of the following variables has at least one missing observation?

print("Check summary(loans) will get : log.annual.inc, days.with.cr.line, revol.util, inq.last.6mths, delinq.2yrs, pub.rec")

[1] "Check summary(loans) will get : log.annual.inc, days.with.cr.line, revol.util, inq.last.6mths, delinq.2yrs, pub.rec"

【1.3 決定是否要補缺項】Which of the following is the best reason to fill in the missing values for these variables instead of removing observations with missing data?

print("第一個選項:移除過多資料會造成Overfitting，然而有任一缺陷資料的row其實只有62筆，不太會發生。")

[1] "第一個選項:移除過多資料會造成Overfitting，然而有任一缺陷資料的row其實只有62筆，不太會發生。"

print("第三個選項:移除NA資料會造成confusing matrix比率失衡，然而實際計算後發現比率其實相差不多。")

[1] "第三個選項:移除NA資料會造成confusing matrix比率失衡，然而實際計算後發現比率其實相差不多。"

print("所以較正確的選項應該是第二個，此模型也須包含預測NA的使用者(然而使用Mice套件補齊資料後，難道就能夠預測NA的資料嗎?")

[1] "所以較正確的選項應該是第二個，此模型也須包含預測NA的使用者(然而使用Mice套件補齊資料後，難道就能夠預測NA的資料嗎?"

【1.4 補缺項工具】What best describes the process we just used to handle missing values?

library(mice)

Loading required package: lattice

Attaching package: 'mice'

The following objects are masked from 'package:base':

    cbind, rbind

set.seed(144)
vars.for.imputation = setdiff(names(loans), "not.fully.paid")
imputed = complete(mice(loans[vars.for.imputation]))

 iter imp variable
  1   1  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  1   2  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  1   3  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  1   4  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  1   5  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  2   1  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  2   2  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  2   3  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  2   4  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  2   5  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  3   1  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  3   2  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  3   3  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  3   4  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  3   5  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  4   1  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  4   2  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  4   3  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  4   4  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  4   5  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  5   1  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  5   2  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  5   3  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  5   4  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec
  5   5  log.annual.inc  days.with.cr.line  revol.util  inq.last.6mths  delinq.2yrs  pub.rec

loans[vars.for.imputation] = imputed
# We predicted missing variable values using the available independent variables for each observation. 

---

2 建立模型 Prediction Models

【2.1 顯著性】 Now, use logistic regression trained on the training set to predict the dependent variable not.fully.paid using all the independent variables.

Which independent variables are significant in our model?

library(caTools)
#Split data
set.seed(144)
spl = sample.split(loans$not.fully.paid, 0.7)
train = subset(loans, spl == TRUE)
test = subset(loans, spl == FALSE)
mod = glm(not.fully.paid~., data=train, family="binomial")
summary(mod)

Call:
glm(formula = not.fully.paid ~ ., family = "binomial", data = train)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-2.187  -0.621  -0.495  -0.361   2.639  

Coefficients:
                             Estimate  Std. Error z value     Pr(>|z|)    
(Intercept)                9.15800765  1.55514067    5.89 0.0000000039 ***
credit.policy             -0.34924658  0.10082598   -3.46      0.00053 ***
purposecredit_card        -0.61439777  0.13440351   -4.57 0.0000048472 ***
purposedebt_consolidation -0.32165576  0.09182845   -3.50      0.00046 ***
purposeeducational         0.13578079  0.17527403    0.77      0.43853    
purposehome_improvement    0.17435349  0.14791842    1.18      0.23851    
purposemajor_purchase     -0.48141854  0.20079306   -2.40      0.01650 *  
purposesmall_business      0.41343357  0.14183237    2.91      0.00356 ** 
int.rate                   0.62211379  2.08490225    0.30      0.76541    
installment                0.00127297  0.00020918    6.09 0.0000000012 ***
log.annual.inc            -0.43129368  0.07145262   -6.04 0.0000000016 ***
dti                        0.00462735  0.00549967    0.84      0.40013    
fico                      -0.00929381  0.00170831   -5.44 0.0000000532 ***
days.with.cr.line          0.00000219  0.00001588    0.14      0.89044    
revol.bal                  0.00000303  0.00000117    2.60      0.00928 ** 
revol.util                 0.00191601  0.00153293    1.25      0.21134    
inq.last.6mths             0.08073987  0.01586849    5.09 0.0000003617 ***
delinq.2yrs               -0.08338381  0.06554209   -1.27      0.20330    
pub.rec                    0.33104300  0.11375811    2.91      0.00361 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5896.6  on 6704  degrees of freedom
Residual deviance: 5487.4  on 6686  degrees of freedom
AIC: 5525

Number of Fisher Scoring iterations: 5

【2.2 從回歸係數估計邊際效用】Consider two loan applications, which are identical other than the fact that the borrower in Application A has FICO credit score 700 while the borrower in Application B has FICO credit score 710. What is the value of Logit(A) - Logit(B)? What is the value of O(A)/O(B)?

# FICO of A is 10 less than B, So the logit will more 10 * 0.009293 = 0.09293 => Logit(A) - Logit(B)
# Logit = log(odd)  -> odd = exp(Logit) 
# odd(A) / odd(b) = exp(Logit A) / exp(Logit B) = exp(Logit A - Logit B) =>
print("the difference of logits is 0.09293")

[1] "the difference of logits is 0.09293"

print("the ratio of odds is exp(0.09293) = ")

[1] "the ratio of odds is exp(0.09293) = "

exp(0.09293)

[1] 1.097

【2.3 混淆矩陣、正確率 vs 底線機率】What is the accuracy of the logistic regression model? What is the accuracy of the baseline model?

predicted.risk = predict(mod,newdata = test, type = "response")
test$predicted.risk = predicted.risk
table(test$not.fully.paid,test$predicted.risk >= 0.5)

   
    FALSE TRUE
  0  2400   13
  1   457    3

# test accuracy
(2400+3) / (2400+13+457+3)

[1] 0.8364

# baseline accuracy
table(test$not.fully.paid)

   0    1 
2413  460 

2413 / (2413+460)

[1] 0.8399

【2.4 ROC & AUC】Use the ROCR package to compute the test set AUC.

library(ROCR)
ROCRpred = prediction(predicted.risk , test$not.fully.paid)
as.numeric(performance(ROCRpred, "auc")@y.values)

[1] 0.6719

---

3 提高底線 Smart Baseline

【3.1 高底線模型】 Using the training set, build a bivariate logistic regression model (aka a logistic regression model with a single independent variable) that predicts the dependent variable not.fully.paid using only the variable int.rate.

The variable int.rate is highly significant in the bivariate model, but it is not significant at the 0.05 level in the model trained with all the independent variables. What is the most likely explanation for this difference?

md2 = glm(not.fully.paid ~ int.rate, data = train, family = "binomial")
summary(md2)

Call:
glm(formula = not.fully.paid ~ int.rate, family = "binomial", 
    data = train)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.055  -0.627  -0.544  -0.436   2.291  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -3.673      0.169   -21.8   <2e-16 ***
int.rate      15.921      1.270    12.5   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5896.6  on 6704  degrees of freedom
Residual deviance: 5734.8  on 6703  degrees of freedom
AIC: 5739

Number of Fisher Scoring iterations: 4

#int.rate在前面的model中，coefficient相當不顯著，然而當移除所有變數後，反而出現了顯著的係數，可能發生了共線性問題 =>此變數與其他變數有很高的相關性。

【3.2 高底線模型的預測值】What is the highest predicted probability of a loan not being paid in full on the testing set? With a logistic regression cutoff of 0.5, how many loans would be predicted as not being paid in full on the testing set?

pretest2 = predict(md2,newdata = test, type = "response")
table(pretest2 >= 0.5)

FALSE 
 2873 

#以模型2預測的結果來說，50%以上機率未還款的資料筆數為0。

【3.3 高底線模型的辨識率】What is the test set AUC of the bivariate model?

library(ROCR)
ROCRpred2 = prediction(pretest2 , test$not.fully.paid)
as.numeric(performance(ROCRpred2, "auc")@y.values)

[1] 0.6239

---

4 預估投資獲利 Computing the Profitability of an Investment

【4.1 投資價值的算法】How much does a $10 investment with an annual interest rate of 6% pay back after 3 years, using continuous compounding of interest?

10 * exp(0.06*3)

[1] 11.97

【4.2 投資獲利的算法，合約完成】While the investment has value c * exp(rt) dollars after collecting interest, the investor had to pay $c for the investment. What is the profit to the investor if the investment is paid back in full?

# C * exp(r*t) - C = C (exp(r*t) - 1)

【4.3 投資獲利的算法，違約】Now, consider the case where the investor made a $c investment, but it was not paid back in full. Assume, conservatively, that no money was received from the borrower (often a lender will receive some but not all of the value of the loan, making this a pessimistic assumption of how much is received). What is the profit to the investor in this scenario?

# -c 

---

5 簡單投資策略 A Simple Investment Strategy

【5.1 計算測試資料的實際投報率】What is the maximum profit of a $10 investment in any loan in the testing set?

test$profit = 10 * exp(test$int.rate*3) - 10
test$profit[test$not.fully.paid == 1] = -10
max(test$profit)

[1] 8.895

---

6 面對不確定性的投資策略 An Investment Strategy Based on Risk

A simple investment strategy of equally investing in all the loans would yield profit $20.94 for a $100 investment. But this simple investment strategy does not leverage the prediction model we built earlier in this problem.

【6.1 高利率、高風險】What is the average profit of a $1 investment in one of these high-interest loans (do not include the $ sign in your answer)? What proportion of the high-interest loans were not paid back in full?

higherInterset = subset(test,int.rate>=0.15)
higherInterset$profit = exp(higherInterset$int.rate*3) -1
higherInterset$profit[higherInterset$not.fully.paid == 1] = -1
mean(higherInterset$profit) * 10

[1] 2.251

table(higherInterset$not.fully.paid)

  0   1 
327 110 

110 / (327+110)

[1] 0.2517

【6.2 高利率之中的低風險】What is the profit of the investor, who invested $1 in each of these 100 loans? How many of 100 selected loans were not paid back in full?

cutoff = sort(higherInterset$predicted.risk, decreasing=FALSE)[100]
selectedLoan = subset(higherInterset,predicted.risk<=cutoff)
selectedLoan$profit = exp(selectedLoan$int.rate*3) -1
selectedLoan$profit[selectedLoan$not.fully.paid == 1] = -1
mean(selectedLoan$profit) * 10 #Answer1

[1] 3.128

table(selectedLoan$not.fully.paid)

 0  1 
81 19 

---

【Q】利用我們建好的模型，你可以設計出比上述的方法獲利更高的投資方法嗎？請詳述你的作法？

trySet = subset(test,int.rate>=0.145)
trySet$profit = exp(trySet$int.rate*3) -1
trySet$profit[trySet$not.fully.paid == 1] = -1
cutoff2 = sort(trySet$predicted.risk , decreasing = F)[100]
selectedLoan2 = subset(trySet, predicted.risk <= cutoff2)
selectedLoan2$profit = exp(selectedLoan2$int.rate*3) -1
selectedLoan2$profit[selectedLoan2$not.fully.paid == 1] = -1
mean(selectedLoan2$profit)

[1] 0.3737

# 不考量風險，直接進行隨機抽樣。

---