信用评分卡

数据描述：数据属于个人消费类贷款，只考虑评分卡最终实施时能够使用到的数据应从如下一些方面获取数据：

– 基本属性：包括了借款人当时的年龄。

– 偿债能力：包括了借款人的月收入、负债比率。

– 信用往来：两年内35-59天逾期次数、两年内60-89天逾期次数、两年内90天或高于90天逾期的次数。

– 财产状况：包括了开放式信贷和贷款数量、不动产贷款或额度数量。

– 贷款属性：暂无。

– 其他因素：包括了借款人的家属数量（不包括本人在内）。

变量名	变量类型	变量描述
SeriousDlqin2yrs	Y/N	超过90天或更糟的逾期拖欠
Revolving Utilization OfUnsecuredLines	percentage	无担保放款的循环利用：除了不动产和像车贷那样除以信用额度总和的无分期付款债务的信用卡和个人信用额度总额
age	integer	借款人当时的年龄
NumberOfTime30-59DaysPastDueNotWorse	integer	35-59天逾期但不糟糕次数
DebtRatio	percentage	负债比率
MonthlyIncome	real	月收入
NumberOfOpenCreditLinesAndLoans	integer	开放式信贷和贷款数量，开放式贷款（分期付款如汽车贷款或抵押贷款）和信贷（如信用卡）的数量
NumberOfTimes90DaysLate	integer	90天逾期次数：借款者有90天或更高逾期的次数
NumberRealEstateLoansOrLines	integer	不动产贷款或额度数量：抵押贷款和不动产放款包括房屋净值信贷额度
NumberOfTime60-89DaysPastDueNotWorse	integer	60-89天逾期但不糟糕次数：借款人在在过去两年内有60-89天逾期还款但不糟糕的次数
NumberOfDependents	integer	家属数量：不包括本人在内的家属数量

导入数据

rm(list = ls())
setwd('C:\\Users\\18073619\\Desktop\\R_workspace\\Machine_Learning')
traindata = read.csv('cs-training.csv')
library(VIM)
library(mice)
library(DMwR)
library(ggplot2)
library(corrplot)
library(caret)
library(pROC)

Part I 缺失值的分析和处理

可视化缺失值(默认缺失值为红色)

matrixplot(traindata)
md.pattern(traindata)

使用KNN方法对缺失值进行填补

traindata<-knnImputation(traindata,k=10,meth = "weighAvg")

traindata <- read.csv('cs-training_nona.csv',header = T)
traindata = traindata[,-1]

Part II 异常值处理和分析

异常值处理的方法

单变量异常值检测：在R语言中使用函数boxplot.stats()可以实现单变量检测，该函数根据返回的统计数据生成箱线图。在上述函数的返回结果中，有一个参数out，它是由异常值组成的列表。更明确的说就是里面列出了箱线图中箱须线外面的数据点。比如我们可以查看月收入分布，第一幅图为没有删除异常值的箱线图。第二幅箱线图删除异常值后，可以发现月收入主要集中分布在3000-8000之间。但是在这份分析报告中，因为我们对业务尚不熟悉，不好将大于8000的数据直接归为异常值，因此对该变量未做处理。
使用LOF（局部异常因子）检测异常值：LOF（局部异常因子）是一种基于密度识别异常值的算法。算法实现是：将一个点的局部密度与分布在它周围的点的密度相比较，如果前者明显的比后者小，那么这个点相对于周围的点来说就处于一个相对比较稀疏的区域，这就表明该点事一个异常值。LOF算法的缺点是它只对数值型数据有效。包‘DMwR’和包‘dprep’中的lofactor()可以计算LOF算法中的局部异常因子。
通过聚类检测异常值：检测异常值的另外一种方式就是聚类。先把数据聚成不同的类，选择不属于任何类的数据作为异常值。例如，基于密度的聚类DBSCAN算法的实现就是将与数据稠密区域紧密相连的数据对象划分为一个类，因此与其他对象分离的数据就会作为异常值。也可以使用K均值算法实现异常值的检测。首先通过把数据划分为k组，划分方式是选择距离各自簇中心最近的点为一组；然后计算每个对象和对应的簇中心的距离（或者相似度），并挑出拥有最大的距离的点作为异常值。

unique(traindata$age) #unique函数去重

可以看到年龄中存在0值，显然是异常值，予以剔除。

traindata <- traindata[-which(traindata$age==0),]
traindata <- traindata[-which(traindata$NumberOfTime30.59DaysPastDueNotWorse %in% c(96,98)),]

Part III 变量分析

单变量分析

我们可以简单地看下部分变量的分布，比如对于age变量，如下图
：

colnames(traindata) = c('y','X1','X2','X3','X4','X5','X6','X7','X8','X9','X10')
ggplot(traindata, aes(x = X2, y = ..density..)) + geom_histogram(fill = "blue", colour = "grey60", size = 0.2, alpha = 0.2) + geom_density()

可以看到年龄变量大致呈正态分布，符合统计分析的假设。再比如月收入变量，也可以做图观察观察，如下：

ggplot(traindata, aes(x = X5, y = ..density..)) + geom_histogram(fill = "blue", colour = "grey60", size = 0.2, alpha = 0.2) + geom_density() + xlim(1, 20000)

月收入也大致呈正态分布，符合统计分析的需要。

变量之间的相关性

建模之前首先得检验变量之间的相关性，如果变量之间相关性显著，会影响模型的预测效果。下面通过corrplot函数，画出各变量之间，包括响应变量与自变量的相关性。

cor1<-cor(traindata[,1:11])
corrplot(cor1)
corrplot(cor1,method = "number")

由上图可以看出，各变量之间的相关性是非常小的。其实Logistic回归同样需要检验多重共线性问题，不过此处由于各变量之间的相关性较小，可以初步判断不存在多重共线性问题，当然我们在建模后还可以用VIF（方差膨胀因子）来检验多重共线性问题。如果存在多重共线性，即有可能存在两个变量高度相关，需要降维或剔除处理。

切分数据集

table(traindata$y)

由上表看出，对于响应变量SeriousDlqin2yrs，存在明显的类失衡问题，SeriousDlqin2yrs等于1的观测为9879，仅为所有观测值的6.6%。因此我们需要对非平衡数据进行处理，在这里可以采用SMOTE算法，用R对稀有事件进行超级采样。

我们利用caret包中的createDataPartition（数据分割功能）函数将数据随机分成相同的两份。

set.seed(1234)
splitIndex<-createDataPartition(traindata$y,time=1,p=0.5,list=FALSE)
train<-traindata[splitIndex,]
test<-traindata[-splitIndex,]

对于分割后的训练集和测试集均有74865个数据，分类结果的平衡性如下：

prop.table(table(train$y))
prop.table(table(test$y))

两者的分类结果是平衡的，仍然有6.6%左右的代表，我们仍然处于良好的水平。因此可以采用这份切割的数据进行建模及预测。

Logistic 回归

Logistic 回归在信用评分卡开发中起到核心作用。由于其特点，以及对自变量进行了证据权重转换（WOE），Logistics回归的结果可以直接转换为一个汇总表，即所谓的标准评分卡格式。

基本公式

Logistic回归模型本事是一个非线性回归模型，经过logit转换（连接函数）将相应变量Y和线性自变量相联系，可以得到一个线性的形式，使用线性回归模型对参数进行估计，所以说logistic回归模型是一个广义线性模型。

下面简单介绍下Logistics回归模型。考虑具有n个独立变量的向量 $x = (x_1,x_2,...,x_n)$ ,设条件概率 $P(y=1|x) = p$ 为根据观测相对于某时间x发生的概率。那么，Logistic回归模型可以表示为：
$P(y=1|x) = \pi(x) = \frac{1}{1+e^{-g(x)}}$
其中， $g(x) = w_0+w_1x_1+w_2x_2+...+w_nx_n$ ,那么在 $x$ 条件下 $y$ 不发生的概率为：
$P(y=0|x) = 1-P(y=1|x) = 1-\frac{1}{1+e^{-g(x)}} = \frac{1}{1+e^{g(x)}}$
这个比值成为事件的发生比（the odds of experiencing an event），简称odds.对odds取对数：
$ln(\frac{P}{1-P}) = g(x) = w_0+w_1x_1+w_2x_2+...+w_nx_n$

fit<-glm(y~.,train,family = "binomial")
summary(fit)
fit2<-glm(y~X2+X3+X5+X7+X8+X9+X10,train,family = "binomial")
summary(fit2)

利用全变量进行回归，模型拟合效果并不是很好，其中 $x_1,x_4,x_6$ 三个变量的p值未能通过检验，在此直接剔除这三个变量，利用剩余的变量对y进行回归。

第二个回归模型所有变量都通过了检验，甚至AIC值（赤池信息准则）更小，所有模型的拟合效果更好些。

模型评估

通常一个二值分类器可以通过ROC(Receiver Operating Characteristic)曲线和AUC值来评价优劣。

很多二元分类器会产生一个概率预测值，而非仅仅是0-1预测值。我们可以使用某个临界点（例如0.5），以划分哪些预测为1，哪些预测为0。得到二元预测值后，可以构建一个混淆矩阵来评价二元分类器的预测效果。所有的训练数据都会落入这个矩阵中，而对角线上的数字代表了预测正确的数目，即true positive + true nagetive。同时可以相应算出TPR（真正率或称为灵敏度）和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更大的分类器效果更好。

| | 真实情况|
|---|---|-----|
|| P| N|
|预测情况|True positive,TP|False Positive,FP|

	False Negative,FN	True Negative,TN

$fp　rate = \frac{FP}{P}$
$tp　rate = \frac{TP}{P}$
$precision = \frac{TP}{TP+FP}$
$recall = \frac{TP}{P}$
$accuracy = \frac{TP+TN}{P+N}$

下面首先利用模型对test数据进行预测，生成概率预测值。

pre <- predict(fit2,test)

在R中，可以利用pROC包，它能方便比较两个分类器，还能自动标注出最优的临界点，图看起来也比较漂亮。在下图中最优点FPR=1-TNR=0.845，TPR=0.638，AUC值为0.8102，说明该模型的预测效果还是不错的，正确较高。

modelroc <- roc(test$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)

WOE转换

证据权重（Weight of Evidence,WOE）转换可以将Logistic回归模型转变为标准评分卡格式。引入WOE转换的目的并不是为了提高模型质量，只是一些变量不应该被纳入模型，这或者是因为它们不能增加模型值，或者是因为与其模型相关系数有关的误差较大，其实建立标准信用评分卡也可以不采用WOE转换。这种情况下，Logistic回归模型需要处理更大数量的自变量。尽管这样会增加建模程序的复杂性，但最终得到的评分卡都是一样的。

用WOE(x)替换变量x。WOE(x)=ln[(违约/总违约)/(正常/总正常)]。

通过上述的Logistic回归，剔除 $x_1$ , $x_4$ , $x_6$ 三个变量，对剩下的变量进行WOE转换。

1、进行分箱
age变量( $x_2$ )

cutx2= c(-Inf,30,35,40,45,50,55,60,65,75,Inf)  
plot(cut(train$X2,cutx2))

NumberOfTime30-59DaysPastDueNotWorse变量(x3)

cutx3 = c(-Inf,0,1,3,5,Inf)  
plot(cut(train$X3,cutx3))

MonthlyIncome变量(x5)：

cutx5 = c(-Inf,1000,2000,3000,4000,5000,6000,7500,9500,12000,Inf)  
plot(cut(train$X5,cutx5))

NumberOfTimes90DaysLate变量(x7)：

cutx7 = c(-Inf,0,1,3,5,10,Inf)  
plot(cut(train$X7,cutx7))

NumberRealEstateLoansOrLines变量(x8)：

cutx8= c(-Inf,0,1,2,3,5,Inf)  
plot(cut(train$X8,cutx8))

NumberOfTime60-89DaysPastDueNotWorse变量(x9)：

cutx9 = c(-Inf,0,1,3,5,Inf)  
plot(cut(train$X9,cutx9))

NumberOfDependents变量(x10)：

cutx10 = c(-Inf,0,1,2,3,5,Inf)  
plot(cut(train$X10,cutx10))

计算WOE值

age变量(x2)：

cutx2= c(-Inf,30,35,40,45,50,55,60,65,75,Inf)  
plot(cut(train$X2,cutx2))

NumberOfTime30-59DaysPastDueNotWorse变量(x3)

cutx3 = c(-Inf,0,1,3,5,Inf)
plot(cut(train$X3,cutx3))

MonthlyIncome变量(x5)：

cutx5 = c(-Inf,1000,2000,3000,4000,5000,6000,7500,9500,12000,Inf)  
plot(cut(train$X5,cutx5))

NumberOfTimes90DaysLate变量(x7)：

cutx7 = c(-Inf,0,1,3,5,10,Inf)  
plot(cut(train$X7,cutx7))

NumberRealEstateLoansOrLines变量(x8)：

cutx8= c(-Inf,0,1,2,3,5,Inf)  
plot(cut(train$X8,cutx8))

NumberOfTime60-89DaysPastDueNotWorse变量(x9)：

cutx9 = c(-Inf,0,1,3,5,Inf)  
plot(cut(train$X9,cutx9))

NumberOfDependents变量(x10)：

cutx10 = c(-Inf,0,1,2,3,5,Inf)  
plot(cut(train$X10,cutx10))

计算WOE

定义WOE转换函数：

getWOE1 <- function(a,p,q,y = train$y){
  totalgood = as.numeric(table(y))[1]  
  totalbad = as.numeric(table(y))[2]
  Good <- as.numeric(table(y[a > p & a <= q]))[1]  
  Bad <- as.numeric(table(y[a > p & a <= q]))[2]  
  WOE <- log((Bad/totalbad)/(Good/totalgood),base = exp(1))  
  return(WOE)  
}
getWOE <- function(aa,cutt,y = train$y){
  woe = rep(NA,length(aa))
  dd = list()
  cutt = cbind(cutt[1:(length(cutt)-1)],cutt[2:(length(cutt))])
  dwoe = NULL
  for( j in 1:nrow(cutt)){
    woe[aa>cutt[j,1] & aa <= cutt[j,2]] = getWOE1(aa,cutt[j,1],cutt[j,2])
    dwoe = c(dwoe,getWOE1(aa,cutt[j,1],cutt[j,2]))
  }
  dd$woe = woe
  names(dwoe) = paste(cutt[,1],'~',cutt[,2])
  dd$dwoe = dwoe
  return(dd)  
}

age.WOE = getWOE(train$X2,cutx2)
NumberOfTime30_59DaysPastDueNotWorse.WOE = getWOE(train$X3,cutx3)
MonthlyIncome.WOE = getWOE(train$X5,cutx5)
NumberOfTimes90DaysLate.WOE = getWOE(train$X7,cutx7)
NumberRealEstateLoansOrLine.WOE = getWOE(train$X8,cutx8)
NumberOfTime60_89DaysPastDueNotWorse.WOE = getWOE(train$X9,cutx9)
NumberOfDependents.WOE = getWOE(train$X10,cutx10)

trainWOE = cbind.data.frame(age.WOE[[1]],NumberOfTime30_59DaysPastDueNotWorse.WOE[[1]],
                            MonthlyIncome.WOE[[1]],NumberOfTime60_89DaysPastDueNotWorse.WOE[[1]],
                            NumberOfTimes90DaysLate.WOE[[1]],NumberRealEstateLoansOrLine.WOE[[1]],
                            NumberOfDependents.WOE[[1]])

评分卡的创建和实施

标准评分卡采用的格式是评分卡中的每一个变量都遵循一系列IF-THEN法则，变量的值决定了该变量所分配的分值，总分就是各变量分值的和。

评分卡设定的分值刻度可以通过将分值表示为违约和正常概率比对数(log(odds))的线性表达式来定义：
$Score = A - B \times log(odds),odds = \frac{p}{1-p}$

在前面的流程中，我们得到了每个变量经WOE变换后各类别的WOE值和各变量的logistic回归模型的系数 $(\beta_0,\beta_1,...)$ ，那么每条记录的违约和正常概率比的对数(log(odds))就可以得到。
$Score = A-B\times(\beta_0+\beta_1 x_1+...+\beta_p x_p)$

$Score = A - B \times \pmatrix{\beta_0 \\ +(\beta_1w_{11})\delta_{11}+(\beta_1w_{12})\delta_{12}+...\\ ...\\ +(\beta_pw_{p1})\delta_{p1}+(\beta_pw_{p2})\delta_{p2}+... }$

$Score = (A - B\beta_0) -(B\beta_1w_{11})\delta_{11}-(B\beta_1w_{12})\delta_{12}-...\\ -...\\ -(B\beta_pw_{p1})\delta_{p1}-(B\beta_pw_{p2})\delta_{p2}$

知道线性表达式的两个参数A，B后就可以求每条记录（申请人）的分值。为了求得A，B，需要设定两个假设（分数的给定，很主观），求得A=521.8622；B=28.8539.

i. odds=1/15时，score=600分比率为1/15时，分值为600

ii. odds=1/15时，score=600分比率为1/15时，分值为600

计算出基础分接近600，而训练集的好坏比约为14，合理。

因为数据中“1”代表的是违约，直接建模预测，求的是“发生违约的概率”，log(odds)即为“坏好比”。为了符合常规理解，分数越高，信用越好，所有就调换“0”和“1”，使建模预测结果为“不发生违约的概率”，最后log(odds)即表示为“好坏比”。

trainWOE$y = train$y
glm.fit = glm(y~.,data = trainWOE,family = binomial(link = logit))
summary(glm.fit)
coe = glm.fit$coefficients
print(coe)

p <- 20/log(2)
q <- 600-20*log(15)/log(2)

$Score=q + p*(as.numeric(coe[1])+as.numeric(coe[2])*age.WOE\\ +as.numeric(coe[3])*NumberOfTime30.59DaysPastDueNotWorse.WOE\\ +p*as.numeric(coe[4])*MonthlyIncome.WOE\\ +p*as.numeric(coe[5])*NumberOfTime60.89DaysPastDueNotWorse.WOE\\ +p*as.numeric(coe[6])*NumberOfTimes90DaysLate.WOE\\ +p*as.numeric(coe[7])*NumberRealEstateLoansOrLine.WOE\\ +p*as.numeric(coe[8])*NumberOfDependents.WOE)$

个人总评分=基础分+各部分得分

基础分为:

base <- q + p*as.numeric(coe[1])
base

对各变量进行打分

age.score = round(p*age.WOE[[2]]*coe[2],0)
NumberOfTime30_59DaysPastDueNotWorse.score = round(p*NumberOfTime30_59DaysPastDueNotWorse.WOE[[2]]*coe[3],0)
MonthlyIncome.score = round(p*MonthlyIncome.WOE[[2]]*coe[4],0)
NumberOfTime60_89DaysPastDueNotWorse.score = round(p*NumberOfTime60_89DaysPastDueNotWorse.WOE[[2]]*coe[5],0)
NumberOfTimes90DaysLate.score = round(p*NumberOfTimes90DaysLate.WOE[[2]]*coe[6],0)
NumberRealEstateLoansOrLine.score = round(p*NumberRealEstateLoansOrLine.WOE[[2]]*coe[7],0)
NumberOfDependents.score = round(p*NumberOfDependents.WOE[[2]]*coe[8],0)

print(age.score)
print(NumberOfTime30_59DaysPastDueNotWorse.score)
print(MonthlyIncome.score)
print(NumberOfTime60_89DaysPastDueNotWorse.score)
print(NumberOfTimes90DaysLate.score)
print(NumberRealEstateLoansOrLine.score)
print(NumberOfDependents.score)