Kaggle|Give Me Some Credit信用卡評分建模分析(R語言)
1.目的
本文是基於Kaggle|Give Me Some Credit專案(資料地址:https://www.kaggle.com/c/GiveMeSomeCredit),通過對消費者的人口特徵、信用歷史記錄、交易記錄等大量資料進行系統的分析、挖掘資料蘊含的行為模式、信用特徵，發展出預測行的模式，結合信用卡評分的構建原理，採用R語言完成資料的清洗，主要包括缺失資料的填充、異常的刪除和資料的分箱；呼叫Logistic迴歸模型建立信用卡評分的基礎模型，藉助自變數的證據權重轉換（WOE）建立信用卡評分卡，並開發一個簡單的信用評分系統。
2.背景
信用評分卡模型在國外是一種成熟的預測方法，尤其在信用風險評估以及金融風險控制領域更是得到了比較廣泛的使用，其原理是將模型變數WOE編碼方式離散化之後運用logistic迴歸模型進行的一種二分類變數的廣義線性模型。
客戶申請評分卡由一系列特徵項組成，每個特徵項相當於申請表上的一個問題（例如，年齡、銀行流水、收入等）。每一個特徵項都有一系列可能的屬性，相當於每一個問題的一系列可能答案（例如，對於年齡這個問題，答案可能就有30歲以下、30到45等）。在開發評分卡系統模型中，先確定屬性與申請人未來信用表現之間的相互關係，然後給屬性分配適當的分數權重，分配的分數權重要反映這種相互關係。分數權重越大，說明該屬性表示的信用表現越好。一個申請的得分是其屬性分值的簡單求和。如果申請人的信用評分大於等於金融放款機構所設定的界限分數，此申請處於可接受的風險水平並將被批准；低於界限分數的申請人將被拒絕或給予標示以便進一步審查。
3.資料處理
3.1資料說明
資料來自Kaggle的give me some credit專案，其中cs-training.csv是有15萬條的樣本資料,包含了12個變數,大致情況如下:
每個變數的解釋如下:
SeriousDlqin2yrs：超過90天或更糟的逾期拖欠
RevolvingUtilization Of UnsecuredLines：無擔保放款的迴圈利用：除了不動產和像車 貸那樣除以信用額度總和的無分期付款債務的信用卡和個人信用額度總額
Age：借款人年齡
NumberOfTime30-59DaysPastDueNotWorse：30-59天逾期次數
DebtRatio：負債比例
MonthlyIncome：月收入
Number Of OpenCreditLinesAndLoans：貸款數量NumberOfTimes90DaysLate：90天逾期次數：借款者有90天或更高逾期的次數
NumberReal Estate Loans Or Lines：不動產貸款或額度數量：抵押貸款和不動產放款包括房屋淨值信貸額度
Number Of Time 60-89Days PastDue Not Worse：60-89天逾期次數NumberOfDependents：家屬數量,不包括本人在內的家屬數量
3.2匯入資料
setwd('E:/data for give me some credit')
credit_data <- read.csv('cs-training.csv',stringsAsFactors = FALSE)
str(credit_data) #顯示資料結構
'data.frame': 150000 obs. of 12 variables:
$ X : int 1 2 3 4 5 6 7 8 9 10 ...
$ SeriousDlqin2yrs : int 1 0 0 0 0 0 0 0 0 0 ...
$ RevolvingUtilizationOfUnsecuredLines: num 0.766 0.957 0.658 0.234 0.907 ...
$ age : int 45 40 38 30 49 74 57 39 27 57 ...
$ NumberOfTime30.59DaysPastDueNotWorse: int 2 0 1 0 1 0 0 0 0 0 ...
$ DebtRatio : num 0.803 0.1219 0.0851 0.036 0.0249 ...
$ MonthlyIncome : int 9120 2600 3042 3300 63588 3500 NA 3500 NA 23684 ...
$ NumberOfOpenCreditLinesAndLoans : int 13 4 2 5 7 3 8 8 2 9 ...
$ NumberOfTimes90DaysLate : int 0 0 1 0 0 0 0 0 0 0 ...
$ NumberRealEstateLoansOrLines : int 6 0 0 0 1 1 3 0 0 4 ...
$ NumberOfTime60.89DaysPastDueNotWorse: int 0 0 0 0 0 0 0 0 0 0 ...
$ NumberOfDependents : int 2 1 0 0 0 1 0 0 NA 2 ...3.3更換變數名
第一列為資料的索引,沒什麼意義,因此可以去除,需要預測的變數是SeriousDlqin2yrs,所以將這個變數重新命名為y,其餘變數重新命名為x1~x10
> cr_data <- credit_data[,2:12]
> names(cr_data) <- c('y','x1','x2','x3','x4','x5','x6','x7','x8','x9','x10')
> head(cr_data)
y x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
1 1 0.7661266 45 2 0.80298213 9120 13 0 6 0 2
2 0 0.9571510 40 0 0.12187620 2600 4 0 0 0 1
3 0 0.6581801 38 1 0.08511338 3042 2 1 0 0 0
4 0 0.2338098 30 0 0.03604968 3300 5 0 0 0 0
5 0 0.9072394 49 1 0.02492570 63588 7 0 1 0 0
6 0 0.2131787 74 0 0.37560697 3500 3 0 1 0 1cr_data$y <- 1-cr_data$y #0表示違約,1表示不違約3.4處理缺失值
這裡需要用到VIM包和mice包,首先識別出缺失資料,VIM包中的matrixplot函式可以將資料集中的缺失資料以視覺化的方式顯示出來,在本例中紅色表示缺失值,顏色越深表示缺失值越多,然後再用mice包中的mice.pattern()函式生成一個數據框來展示缺失值的個數.我自己的電腦VIM打死也沒安裝成功,暫時不演示如何用圖表去看缺失值了,程式碼就是matrixplot(cr_data).
從圖中可以看到x5變數(月收入)缺失值為29731個,x10(家屬數量)缺失值為3924個
> library(lattice)
> library(mice)
> md.pattern(cr_data)
y x1 x2 x3 x4 x6 x7 x8 x9 x10 x5
120269 1 1 1 1 1 1 1 1 1 1 1 0
25807 1 1 1 1 1 1 1 1 1 1 0 1
3924 1 1 1 1 1 1 1 1 1 0 0 2
0 0 0 0 0 0 0 0 0 3924 29731 33655可以看到x5的缺失值較多,直接刪除可能會影響結果,可以採用中位數/均值/眾數等方式填充缺失資料,由於x5呈明顯的正態分佈,因此用中位數去填充x5的缺失資料
> library(ggplot2)
> #x5(月收入)的概率密度圖
> ggplot(data = cr_data,aes(x5))+geom_density(fill="lightskyblue")+xlim(0,25000)+geom_vline(aes(xintercept=median(cr_data$x5[!is.na(cr_data$x5)])),colour="red",linetype="dashed",lwd=1)
> #中位數填充x5缺失值
>cr_data$x5[is.na(cr_data$x5)]=median(cr_data$x5[!is.na(cr_data$x5)])#x10的缺失值數量所佔比列較小,因此,直接將缺失值刪除
> cr_data <- na.omit(cr_data)3.5處理異常值
關於異常值的檢測,常見的檢測方法有:
（1）單變數異常值檢測：呼叫R語言中的boxplot.stats()函式可以實現單變數異常值的檢測，該函式基於資料生成的箱體圖進行異常值的檢測。在函式的返回結果中，引數out是由異常值組成的列表。
（2）使用LOF（local outlier factor，區域性異常因子）進行異常值檢測：LOF（區域性異常因子）是用於識別基於密度的區域性異常值的演算法。使用LOF，一個點的區域性密度會與它的鄰居進行比較。如果前者明顯低於後者（有一個大於1 的LOF值），該點位於一個稀疏區域，對於它的鄰居而言，這就表明，該點是一個異常值。LOF的缺點就是它只對數值資料有效。R語言中DMwR和dprep包的lofactor()函式使用LOF演算法計算區域性異常因子。
（3）聚類進行異常值的檢測：通過把資料聚成類，將那些不屬於任務一類的資料作為異常值。比如，使用k-means演算法來檢測異常。使用k-means演算法，資料被分成k組，通過把它們分配到最近的聚類中心。然後，我們能夠計算每個物件到聚類中心的距離（或相似性），並且選擇最大的距離作為異常值。
3.5.1 處理age異常值
> boxplot(cr_data$x2)
結合現實中的信用卡申請條件,申請人的年齡必須在18~65歲之間,因此將小於18和大於65的資料進行刪除,並用boxplot.stats()檢查刪除後的異常值為0
> cr_data <- cr_data[cr_data$x2 <= 65 & cr_data$x2 >= 18,]
> boxplot.stats(cr_data$x2)
$`stats`
[1] 21 39 48 56 65
$n
[1] 119063
$conf
[1] 47.92216 48.07784
$out
integer(0)3.5.2 處理x3,x7,x9的異常值
> boxplot(cr_data$x3,cr_data$x7,cr_data$x9)
可以看到異常值均為接近100的值,均刪除.
> boxplot(cr_data$x3,cr_data$x7,cr_data$x9)
> unique(cr_data$x3)
[1] 2 0 1 3 4 5 7 10 6 98 12 8 9 96 13 11
> unique(cr_data$x7)
[1] 0 1 3 2 5 4 98 10 9 6 7 8 15 96 11 13 14 17 12
> unique(cr_data$x9)
[1] 0 1 2 5 3 98 4 6 7 8 96 11 9
> cr_data <- cr_data[-which(cr_data$x3==96),]
> cr_data <- cr_data[-which(cr_data$x3==98),]
> unique(cr_data$x7)
[1] 0 1 3 2 5 4 10 9 6 7 8 15 11 13 14 17 12
> unique(cr_data$x9)
[1] 0 1 2 5 3 4 6 7 8 11 9
> #可以看出已無異常值4.建模分析
4.1 變數的相關性分析
建模之前需要先檢驗變數之間的相關性,,如果變數之間具有強相關性,則會影響模型的準確性.呼叫R中的cor()函式來計算不同變數之間的相關係數,同時,呼叫corrplot包中的corrplot()函式來將相關係數視覺化
> library(corrplot)
corrplot 0.84 loaded
> cor1 <- cor(cr_data)
> corrplot(cor1,method = 'number')
由上圖可知:各個變數之間的相關係數較小,相關性較弱,不存在明顯的多重共線問題,採用logistic迴歸需要考慮多重共線問題,在此處由於個變數之間的相關性較小,可以初步判斷不存在多重共線問題.在建模之後也可以通過VIF(方差膨脹因子)來檢驗多重共線問題.當存在多維共線問題時,需要進行降維或剔除處理.
4.2 切分資料集
> table(cr_data$y)
0 1
9044 109785 由上表可知：響應變數y存在著明顯的類失衡問題，y等於0的觀測值為9044，僅佔所有觀測值的7.6%，因此，我們需要對非平衡資料進行處理，基於smote演算法，呼叫R語言中caret包中的createDataPartition對稀有資料進行超級取樣。
> library(caret)
> set.seed(500)
> splitindex <- createDataPartition(cr_data$y,times = 1,p=0.5,list = FALSE)
> traindata <- cr_data[splitindex,]
> testdata <- cr_data[-splitindex,]
> prop.table(table(traindata$y))
0 1
0.0771354 0.9228646
> prop.table(table(testdata$y))
0 1
0.07508331 0.92491669 由上表可知：兩者分類後的結果是平衡的，y等於1的概率均為7.6%左右，處於良好的水平，因此，可以採用切割後的資料進行建模和預測分析。
4.3建模分析
Logistic迴歸模型在信用卡評分開發的基礎模型，由於其自身的特點以及對自變數進行證據權重轉換（WOE），logistic迴歸的結果可以直接轉換為一個彙總表，即所謂的標準評分卡格式。呼叫R語言中glm（）函式對所有變數進行logistic迴歸建模.
> fit=glm(y~.,data=traindata,family = "binomial")
Warning message:
glm.fit:擬合機率算出來是數值零或一
> summary(fit)
Call:
glm(formula = y ~ ., family = "binomial", data = traindata)
Deviance Residuals:
Min 1Q Median 3Q Max
-5.2021 0.2552 0.3014 0.3581 4.4247
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.769e+00 7.395e-02 23.919 < 2e-16 ***
x1 -6.035e-06 7.113e-05 -0.085 0.93239
x2 2.466e-02 1.668e-03 14.785 < 2e-16 ***
x3 -5.599e-01 1.628e-02 -34.399 < 2e-16 ***
x4 5.078e-05 1.778e-05 2.856 0.00429 **
x5 4.136e-05 5.055e-06 8.181 2.81e-16 ***
x6 8.239e-03 4.079e-03 2.020 0.04339 *
x7 -8.260e-01 2.462e-02 -33.550 < 2e-16 ***
x8 -1.148e-01 1.663e-02 -6.899 5.23e-12 ***
x9 -7.321e-01 3.339e-02 -21.927 < 2e-16 ***
x10 -5.823e-02 1.394e-02 -4.176 2.96e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 32288 on 59414 degrees of freedom
Residual deviance: 26558 on 59404 degrees of freedom
AIC: 26580
Number of Fisher Scoring iterations: 6由上述結果可知,變數x1,x4,x6對響應變數y的貢獻不顯著,因此直接刪除這三個變數,利用剩下的變數來進行logistic迴歸
> fit2 <- glm(y~x2+x3+x5+x7+x8+x9+x10,data=traindata,family = "binomial")
Warning message:
glm.fit:擬合機率算出來是數值零或一
> summary(fit2)
Call:
glm(formula = y ~ x2 + x3 + x5 + x7 + x8 + x9 + x10, family = "binomial",
data = traindata)
Deviance Residuals:
Min 1Q Median 3Q Max
-5.1559 0.2554 0.3019 0.3583 4.4472
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.786e+00 7.358e-02 24.278 < 2e-16 ***
x2 2.578e-02 1.633e-03 15.788 < 2e-16 ***
x3 -5.551e-01 1.609e-02 -34.497 < 2e-16 ***
x5 4.037e-05 4.976e-06 8.114 4.91e-16 ***
x7 -8.352e-01 2.435e-02 -34.294 < 2e-16 ***
x8 -9.317e-02 1.532e-02 -6.081 1.19e-09 ***
x9 -7.341e-01 3.342e-02 -21.967 < 2e-16 ***
x10 -6.084e-02 1.388e-02 -4.383 1.17e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 32288 on 59414 degrees of freedom
Residual deviance: 26571 on 59407 degrees of freedom
AIC: 26587
Number of Fisher Scoring iterations: 6
4.4 模型評估
通常一個二值分類器可以通過ROC曲線和AUC值進行評價，由於很多二元分類器會產生一個概率值，而非僅僅的0~1預測值。我們可以使用臨界點，當概率預測值大於該臨界點時，劃分為1，概率預測值小於該臨界點時，劃分為0。得到二元預測值後，可以構建一個混淆矩陣來評價二元分類器的預測效果。所有的訓練資料都會落入這個矩陣中，而對角線上的數字代表了預測正確的數目，即true positive + true nagetive。同時，計算出TPR（真正率或稱為靈敏度）和TNR（真負率或稱為特異度），其中，偽陽性率（FPR）=1-TNR。除了分類器的訓練引數，臨界點的選擇，也會大大的影響TPR和TNR。有時可以根據具體問題和需要，來選擇具體的臨界點。
如果我們選擇一系列的臨界點，就會得到一系列的TPR和TNR，將這些值對應的點連線起來，就構成了ROC曲線。ROC曲線可以幫助我們清楚的瞭解到這個分類器的效能表現，還能方便比較不同分類器的效能。在繪製ROC曲線的時候，習慣上是使用1-TNR作為橫座標即FPR（false positive rate），TPR作為縱座標。這是就形成了ROC曲線。
而AUC（Area Under Curve）被定義為ROC曲線下的面積，顯然這個面積的數值不會大於1。又由於ROC曲線一般都處於y=x這條直線的上方，所以AUC的取值範圍在0.5和1之間。使用AUC值作為評價標準是因為很多時候ROC曲線並不能清晰的說明哪個分類器的效果更好，而作為一個數值，對應AUC更大的分類器效果更好。
呼叫R語言中pROC包中的roc函式計算分類器的AUC值，可以方便的比較兩個分類器，並且自動標註出最優的臨界點。如下圖所示：最優點FPR=1-TNR=0.840，TPR=0.636，AUC值為0.796，說明該模型的預測效果不錯，正確率較高。
> pre=predict(fit2,testdata)
> library(pROC)
Type 'citation("pROC")' for a citation.
載入程輯包：‘pROC’
The following objects are masked from ‘package:stats’:
cov, smooth, var
> modelroc=roc(testdata$y,pre)
> plot(modelroc, print.auc=TRUE, auc.polygon=TRUE, grid=c(0.1, 0.2),
+ grid.col=c("green", "red"), max.auc.polygon=TRUE,auc.polygon.col="skyblue", print.thres=TRUE)5.WOE轉換
證據權重（Weight of Evidence,WOE）轉換可以將Logistic迴歸模型轉化為標準評分卡格式，引入WOE轉換的目的並不是為了提高模型質量，而是由於一些變數不應該被納入模型，或者是因為它們不能增加模型值，或者是因為與其模型相關係數有關的誤差較大，其實建立標準信用評分卡也可以不採用WOE轉換。這種情況下，Logistic迴歸模型需要處理更大數量的自變數。儘管這樣會增加建模程式的複雜性，但最終得到的評分卡都是一樣的。
用WOE(x)替換變數x，WOE()=ln[(違約/總違約)/(正常/總正常)]。由於模型中剔除x1,x4,x6三個變數，因此對剩下的變數進行WOE轉換。
5.1 資料分箱
WOE分箱原則：
1.分箱數量適中，不宜過多和過少。
2.各個分箱內的記錄數應該合理，不應過多或者或過少。
3.結合目標變數，分箱應該表現出明顯的趨勢。
4.相鄰分箱的目標變數分佈差異儘可能大。
x2變數（age分箱）
> cutx2=c(-Inf,25,30,35,40,45,50,55,60,65)
> plot(cut(traindata$x2,cutx2))
x3變數分箱
> cutx3=c(-Inf,0,1,3,5,Inf)
> plot(cut(traindata$x3,cutx3))
x5變數分箱
> cutx5=c(-Inf,2000,3000,4000,5000,6000,7500,9500,12000,Inf)
> plot(cut(traindata$x5,cutx5))
x7變數分箱
> cutx7=c(-Inf,0,1,3,5,Inf)
> plot(cut(traindata$x7,cutx7))
x8分箱
> cutx8= c(-Inf,0,1,2,3,5,Inf)
> plot(cut(traindata$x8,cutx8))
x9分箱
> cutx9 = c(-Inf,0,1,3,5,Inf)
> plot(cut(traindata$x9,cutx9))
x10分箱
> cutx10 = c(-Inf,0,1,2,3,5,Inf)
> plot(cut(traindata$x10,cutx10))
5.2WOE數值計算
> #x2的WOE值
> agelessthan30=getwoe(traindata$x2,-Inf,30)
> age30to40=getwoe(traindata$x2,30,40)
> age40to50=getwoe(traindata$x2,40,50)
> age50to60=getwoe(traindata$x2,50,60)
> age60to65=getwoe(traindata$x2,60,65)
> age.woe=c(agelessthan30,age30to40,age40to50,age50to60,age60to65)
> age.woe
[1] 0.43555387 0.27427494 0.06637011 -0.20620943 -0.74500267同理，計算出其他變數的WOE值，分別如下：
> x3_woe
[1] -0.5207976 0.8374273 1.6003370 2.2152915 2.7804131
> x5_woe
[1] 0.22048191 0.38990261 0.37583957 0.18627954 -0.02563007 -0.15179000 -0.44092717
[8] -0.54503261 -0.51895581
> x7_woe
[1] -0.3789825 1.8726850 2.5803602 3.0893005
> x8_woe
[1] 0.3073539 -0.2468977 -0.2640414 -0.1193822 0.3895990 1.0255582
> x9_woe
[1] -0.2714009 1.6945721 2.6061341 2.9212868 3.3982109
> x10_woe
[1] -0.09893059 0.06819497 0.07729239 0.15846434 0.28712719 0.52653127
待續。。。