一、项目简述
(一)项目目标
基于WOE转换及逻辑回归模型制作一张反映贷款人信用风险水平的信用评分卡。
(二)数据来源
数据来自Kaggle的cs-training.csv样本，该样本有15万条的个人消费类贷款数据。 数据下载地址
(三)数据说明
变量名(原始)为文件中变量名，由于原始变量名较长，所以就转换成较短的y、x1、x2…代替。
| 变量名(原始) | 变量名 | 变量类型 | 变量描述 |
|---|---|---|---|
| id | integer | 序号 | |
| SeriousDlqin2yrs | y | Y/N | 超过90天或更糟的逾期拖欠 |
| RevolvingUtilizationOfUnsecuredLines | x1 | percentage | 无担保放款的循环利用：除了不动产和像车贷那样除以信用额度总和的无分期付款债务的信用卡和个人信用额度总额 |
| age | x2 | integer | 借款人当时的年龄 |
| NumberOfTime30-59DaysPastDueNotWorse | x3 | integer | 35-59天逾期但不糟糕次数 |
| DebtRatio | x4 | percentage | 负债比率 |
| MonthlyIncome | x5 | real | 月收入 |
| NumberOfOpenCreditLinesAndLoans | x6 | integer | 开放式信贷和贷款数量，开放式贷款（分期付款如汽车贷款或抵押贷款）和信贷（如信用卡）的数量 |
| NumberOfTimes90DaysLate | x7 | integer | 90天逾期次数：借款者有90天或更高逾期的次数 |
| NumberRealEstateLoansOrLines | x8 | integer | 不动产贷款或额度数量：抵押贷款和不动产放款包括房屋净值信贷额度 |
| NumberOfTime60-89DaysPastDueNotWorse | x9 | integer | 60-89天逾期但不糟糕次数：借款人在在过去两年内有60-89天逾期还款但不糟糕的次数 |
| NumberOfDependents | x10 | integer | 家属数量：不包括本人在内的家属数量 |
二、建模过程
(一)数据读入及重命名
#读入训练数据集
traindata<-read.csv("cs-training.csv",header = TRUE)
#原始变量名太长,重命名
names(traindata)=c("id","y","x1","x2","x3","x4","x5","x6","x7","x8","x9","x10")
#查看前6行数据,检查变量名是否更新完成
head(traindata)得到下面的结果，可以看到，变量名已经更新完成
> head(traindata)
id y x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
1 1 1 0.7661266 45 2 0.80298213 9120 13 0 6 0 2
2 2 0 0.9571510 40 0 0.12187620 2600 4 0 0 0 1
3 3 0 0.6581801 38 1 0.08511338 3042 2 1 0 0 0
4 4 0 0.2338098 30 0 0.03604968 3300 5 0 0 0 0
5 5 0 0.9072394 49 1 0.02492570 63588 7 0 1 0 0
6 6 0 0.2131787 74 0 0.37560697 3500 3 0 1 0 1(二)缺失值处理
library(VIM)
matrixplot(traindata) VIM包中的matrixplot可以把缺失值的分布展示出来，下图中浅色表示值比较小，深色表示值比较大，而红色表示缺失值。
library(mice)
md.pattern(traindata)mice包中的md.pattern可以按变量统计缺失值个数，最后一行数字表示所在列变量的缺失值个数，可以看到变量x5和x10包含缺失值。(x5:月收入；x10:家属数量)
> md.pattern(traindata)
id y x1 x2 x3 x4 x6 x7 x8 x9 x10 x5
120269 1 1 1 1 1 1 1 1 1 1 1 1 0
25807 1 1 1 1 1 1 1 1 1 1 1 0 1
3924 1 1 1 1 1 1 1 1 1 1 0 0 2
0 0 0 0 0 0 0 0 0 0 3924 29731 33655上述x5和x10变量缺失值占比较高，直接移除会损失较多数据，所以这里我们采用多重插补的方法对缺失值进行填补。
imp<-mice(traindata,met="cart",m=1)#该方法只对数值变量进行插补，分类变量的缺失值保留
traindata_imp<-complete(imp)此外也可以用DMwR包中的knnImputation()进行缺失值填补，不过速度较慢，代码如下
traindata<-knnImputation(traindata,k=5,meth = "median")填补完缺失值后再次调用md.pattern,可以看到已经没有缺失值了。
> md.pattern(traindata_imp)
id y x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
[1,] 1 1 1 1 1 1 1 1 1 1 1 1 0
[2,] 0 0 0 0 0 0 0 0 0 0 0 0 0(三)异常值处理
分别对x2～x10进行异常值检验，一般来说异常值需要结合业务知识来将不符合业务尝试的数据剔除掉，这就要求建模人员熟悉每个变量的数据生产口径，同时也要掌握基本的异常值检测技术。
x2变量为‘借款人当时的年龄’，通过下面的语句发现年龄中包含0岁的借款人，明显是异常值，应该剔除掉。
＃画出x2分布的箱线图
boxplot(traindata_imp$x2)
#列举x2的数据
> unique(traindata_imp$x2)
[1] 45 40 38 30 49 74 57 39 27 51 46 76 64 78 53 43 25 32 58 50 69 24 28 62 42 75
[27] 26 52 41 81 31 68 70 73 29 55 35 72 60 67 36 56 37 66 83 34 44 48 61 80 47 59
[53] 77 63 54 33 79 65 86 92 23 87 71 22 90 97 84 82 91 89 85 88 21 93 96 99 94 95
[79] 101 98 103 102 107 105 0 109
剔除样本中x2等于0的异常值
traindata_imp=traindata_imp[-which(traindata_imp$x2==0),]此外也可以用boxplot.stats()函数剔除异常值,该函数适用于整体分布较为对称的数据，长尾比较严重的数据不适合，其中默认coef＝1.5,可以适当调整coef值，代码如下
traindata_imp=traindata_imp[-which(traindata_imp$x2 %in% boxplot.stats(traindata_imp$x2,coef = 1.5)$out),]剔除完x2的异常值后，同样观察x3、x7、x9的箱线图，发现三个变量靠近100的区域有两个异常值
boxplot(traindata_imp$x3,traindata_imp$x7,traindata_imp$x9,names = c('x3','x7','x9'))
剔除x3、x7、x9的异常值
traindata_imp=traindata_imp[-which(traindata_imp$x3>95|traindata_imp$x7>95|traindata_imp$x9>95),]x4的原始变量名为DebtRatio(负债比率，即全部负债与全部资产的比率)，会计学有个公式：资产＝负债＋所有者权益，而负债比率=负债／资产＊100％，所有者权益>=0,所以负债比例应该不大于1,剔除掉x4中大于1的异常值
traindata_imp=traindata_imp[-which(traindata_imp$x4>1),]观察剩余变量x5、x6、x8、x10，无明确异常值故不做处理。
(四)变量分析
逻辑回归模型需要检验各变量之间是否存在多重共线性问题，通过corrplot()可以打印出各变量两两之间的相关性，可以看到各变量之间的相关性都比较小，之后建模时还可以用VIF来检验多重共线性的问题。
cor1<-cor(traindata_imp[,2:12])
library(corrplot)
corrplot(cor1,method = "number")
(五)变量筛选
使用createDataPartition将样本分隔成70%的训练集和30%的测试集，通过p值控制分割的比例。
set.seed(100)
splitIndex<-createDataPartition(traindata_imp$y,time=1,p=0.7,list=FALSE)
train<-traindata_imp[splitIndex,]
test<-traindata_imp[-splitIndex,] 查看原样本、训练集、测试集中目标变量y的分布，可以看到好坏分布基本一致。
> prop.table(table(traindata_imp$y))
0 1
0.93359856 0.06640144
> prop.table(table(train$y))
0 1
0.93393158 0.06606842
> prop.table(table(test$y))
0 1
0.9328215 0.0671785 建立逻辑回归模型
fit=glm(y~.,data = train[,-1],family = binomial(link = 'logit'))
#查看模型结果
summary(fit)建立逻辑回归模型，查看模型结果，结果如下,可以知道x1的p值>0.05,未通过检验：
> summary(fit)
Call:
glm(formula = y ~ ., family = binomial(link = "logit"), data = train[,
-1])
Deviance Residuals:
Min 1Q Median 3Q Max
-4.6971 -0.3370 -0.2770 -0.2216 3.6352
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.191e+00 6.663e-02 -32.883 < 2e-16 ***
x1 -4.556e-05 9.275e-05 -0.491 0.623281
x2 -2.423e-02 1.254e-03 -19.313 < 2e-16 ***
x3 5.440e-01 1.506e-02 36.134 < 2e-16 ***
x4 1.039e+00 8.401e-02 12.369 < 2e-16 ***
x5 -1.062e-05 4.111e-06 -2.584 0.009774 **
x6 -8.819e-03 3.822e-03 -2.308 0.021022 *
x7 8.815e-01 2.364e-02 37.296 < 2e-16 ***
x8 -7.468e-02 2.061e-02 -3.623 0.000291 ***
x9 7.092e-01 3.221e-02 22.020 < 2e-16 ***
x10 6.217e-02 1.314e-02 4.732 2.23e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 39050 on 80234 degrees of freedom
Residual deviance: 31787 on 80224 degrees of freedom
AIC: 31809
Number of Fisher Scoring iterations: 6去除x1后重现建模，公式如下
fit2=glm(y ~ x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10,data = train[,-1],family = binomial(link = 'logit'))查看新建模型，结果如下，所有变量的p值都通过检验
> summary(fit2)
Call:
glm(formula = y ~ x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10,
family = binomial(link = "logit"), data = train[, -1])
Deviance Residuals:
Min 1Q Median 3Q Max
-4.6973 -0.3371 -0.2770 -0.2217 3.6384
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.191e+00 6.663e-02 -32.887 < 2e-16 ***
x2 -2.422e-02 1.254e-03 -19.309 < 2e-16 ***
x3 5.441e-01 1.506e-02 36.137 < 2e-16 ***
x4 1.039e+00 8.401e-02 12.364 < 2e-16 ***
x5 -1.067e-05 4.112e-06 -2.595 0.009470 **
x6 -8.788e-03 3.821e-03 -2.300 0.021467 *
x7 8.816e-01 2.363e-02 37.299 < 2e-16 ***
x8 -7.467e-02 2.061e-02 -3.622 0.000292 ***
x9 7.093e-01 3.221e-02 22.021 < 2e-16 ***
x10 6.220e-02 1.314e-02 4.735 2.2e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 39050 on 80234 degrees of freedom
Residual deviance: 31788 on 80225 degrees of freedom
AIC: 31808
Number of Fisher Scoring iterations: 6对各变量进行VIF多重共线性检验(当0<VIF<10，不存在多重共线性；当10≤VIF<100，存在较强的多重共线性；当VIF≥100，存在严重多重共线性),可以看到各变量之间不存在多重共线性。
> vif(fit2)
x2 x3 x4 x5 x6 x7 x8 x9 x10
1.109821 1.103973 1.630034 1.430896 1.457423 1.116286 1.841333 1.104896 1.057597 通过测试集检验逻辑回归模型的效果，评估模型效果一般用ROC曲线，ROC横轴为特异率，纵轴为敏感度，代码如下
#用模型对测试集进行预测
test$lg_p <-predict(fit2,test)
test$p <-1/(1+exp(-1*test$lg_p))
#画出模型的ROC曲线
library(pROC)
modelroc <- roc(test$y,test$p)
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)AUC（Area Under Curve）被定义为ROC曲线下的面积，显然这个面积的数值不会大于1，又由于ROC曲线一般都处于y=x这条直线的上方，所以AUC的取值范围在0.5和1之间，AUC越接近于1模型预测效果越好。
可以看到，本模型的AUC为0.808，说明该模型预测效果较为不错。
经过上述过程，最后经过筛选后最后进入信用评分卡模型的变量为 x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 。
(六)对变量进行WOE转换
WOE称为证据权重(Weight of Evidence),计算公式为：ln(违约件占比／正常件占比)，违约件占比高于正常件时，WOE为负数。WOE的绝对值越高，说明该组别好坏客户的区分能力越强，各组之间的WOE值差距应该尽量拉开。
对变量进行分箱
R语言中woe()可以对变量自动计算woe值，不过该函数自动分组功能并不好用，所以我们先手动对变量进行分组。
x2、x4、x5、x6为连续变量，用cut()将每个变量分为等深的10组，x2分组代码如下，其它连续变量代码类似。
#x2分成10组
#查看x2变量分布
plot(density(traindata_imp$x2))
#提取x2的十分位数
x2_q=quantile(traindata_imp$x2,probs = seq(0, 1, 0.1),names=FALSE)
#用无穷小和无穷大分别替换最小值和最大值
x2_q[1]=-Inf
x2_q[length(x2_q)]=Inf
#生成分组组号
traindata_imp$x2_rank=cut(traindata_imp$x2,x2_q,labels = F)
#查看分组样本数
table(traindata_imp$x2_rank)x3、x7、x8、x9、x10为非连续变量，手动将每个变量分为5组(x9分成4组)，x3分组代码如下，其它非连续变量代码类似。
#x3分成5组
#查看x3变量分布
table(traindata_imp$x3)
#x3非连续变量，将其分为5组，变量值>=4的样本归为一组
traindata_imp$x3_rank=traindata_imp$x3+1
traindata_imp$x3_rank[which(traindata_imp$x3_rank>5)]=5
#查看分组样本数
table(traindata_imp$x3_rank)计算各分组WOE值
之前说过本文采用的WOE的计算公式为：ln(违约件占比／正常件占比)，而R语言woe()中WOE值的计算公式为ln(违约件占比／正常件占比)＊100，所以要在woe()的计算结果上在乘以0.01进行转换。
#计算各分组的WOE值
woe_2=woe(Data=traindata_imp,"x2_rank",FALSE,"y",5,Bad=0,Good=1)
woe_3=woe(Data=traindata_imp,"x3_rank",FALSE,"y",5,Bad=0,Good=1)
woe_4=woe(Data=traindata_imp,"x4_rank",FALSE,"y",5,Bad=0,Good=1)
woe_5=woe(Data=traindata_imp,"x5_rank",FALSE,"y",5,Bad=0,Good=1)
woe_6=woe(Data=traindata_imp,"x6_rank",FALSE,"y",5,Bad=0,Good=1)
woe_7=woe(Data=traindata_imp,"x7_rank",FALSE,"y",5,Bad=0,Good=1)
woe_8=woe(Data=traindata_imp,"x8_rank",FALSE,"y",5,Bad=0,Good=1)
woe_9=woe(Data=traindata_imp,"x9_rank",FALSE,"y",5,Bad=0,Good=1)
woe_10=woe(Data=traindata_imp,"x10_rank",FALSE,"y",5,Bad=0,Good=1)以woe_2为例，结果中WOE所在的列就是分组对应自动计算出来的WOE值。
> woe_2
BIN BAD GOOD TOTAL BAD% GOOD% TOTAL% WOE IV BAD_SPLIT GOOD_SPLIT
1 1 10795 1287 12082 0.101 0.169 0.105 51.5 0.035 0.893 0.107
2 2 10433 1095 11528 0.097 0.144 0.101 39.5 0.019 0.905 0.095
3 3 11339 1114 12453 0.106 0.146 0.109 32.0 0.013 0.911 0.089
4 4 10419 871 11290 0.097 0.114 0.098 16.1 0.003 0.923 0.077
5 5 10900 885 11785 0.102 0.116 0.103 12.9 0.002 0.925 0.075
6 6 10219 740 10959 0.095 0.097 0.096 2.1 0.000 0.932 0.068
7 7 9735 540 10275 0.091 0.071 0.090 -24.8 0.005 0.947 0.053
8 8 11881 516 12397 0.111 0.068 0.108 -49.0 0.021 0.958 0.042
9 9 10457 304 10761 0.098 0.040 0.094 -89.6 0.052 0.972 0.028
10 10 10832 259 11091 0.101 0.034 0.097 -108.9 0.073 0.977 0.023提取并转换各分组的WOE值
woe2=woe_2$WOE*0.01
woe3=woe_3$WOE*0.01
woe4=woe_4$WOE*0.01
woe5=woe_5$WOE*0.01
woe6=woe_6$WOE*0.01
woe7=woe_7$WOE*0.01
woe8=woe_8$WOE*0.01
woe9=woe_9$WOE*0.01
woe10=woe_10$WOE*0.01对样本进行WOE赋值
traindata_imp$x2_woe=woe2[traindata_imp$x2_rank]
traindata_imp$x3_woe=woe3[traindata_imp$x3_rank]
traindata_imp$x4_woe=woe4[traindata_imp$x4_rank]
traindata_imp$x5_woe=woe5[traindata_imp$x5_rank]
traindata_imp$x6_woe=woe6[traindata_imp$x6_rank]
traindata_imp$x7_woe=woe7[traindata_imp$x7_rank]
traindata_imp$x8_woe=woe8[traindata_imp$x8_rank]
traindata_imp$x9_woe=woe9[traindata_imp$x9_rank]
traindata_imp$x10_woe=woe10[traindata_imp$x10_rank] 以上就是对各变量进行WOE转换的全部过程。
(七)评分卡设定
标准评分卡采用的格式是评分卡中的每一个变量都遵循一系列IF-THEN法则，变量的值决定了该变量所分配的分值，总分就是各变量分值的和。
信用评分公式为：Score=A-Bln(odds)
A为基础分
B为权重值
odds=p/(1-p),即为坏好比，坏好比越高，score越低
ln(odds)=(b0+b1·x1＋b2·x2+……) b0、b1…即为我们要算出的逻辑回归模型的系数为了求出A和B，需要下面两个设定：
- 当odds＝1:1时，设定score＝500分；
- 当odds＝1:2时，设定score＝530分；(每提高30分，坏好比降低一半)
500分和30分坏好比降一半这两个设定可以根据自身的要求自己调整。这两个设定就给了我们两个公式：
- 500=A-Bln(1/1)
- 530=A-Bln(1/2)
可以计算得到A=500, B=43.28085
(八)计算各变量分数
计算WOE转化后的各变量的系数
#生成训练集和测试集
train<-traindata_imp[splitIndex,]
test<-traindata_imp[-splitIndex,]
#建立逻辑回归模型
fit3=glm(y ~ x2_woe + x3_woe + x4_woe + x5_woe + x6_woe + x7_woe + x8_woe + x9_woe + x10_woe,data = train[,-1],family = binomial(link = 'logit'))
summary(fit3)
#将变量参数赋给coe
coe = fit3$coefficients计算各变量的分数
1、计算基础分数：
A=500
B=30/log(2)
base_score=A-B*coe[1]2、计算训练集各变量的信用评分：
train$x2_score=(-1)*B*train$x2_woe*coe[2]
train$x3_score=(-1)*B*train$x3_woe*coe[3]
train$x4_score=(-1)*B*train$x4_woe*coe[4]
train$x5_score=(-1)*B*train$x5_woe*coe[5]
train$x6_score=(-1)*B*train$x6_woe*coe[6]
train$x7_score=(-1)*B*train$x7_woe*coe[7]
train$x8_score=(-1)*-B*train$x8_woe*coe[8]
train$x9_score=(-1)*B*train$x9_woe*coe[9]
train$x10_score=(-1)*B*train$x10_woe*coe[10]3、计算训练集各观测样本的总分，总评分＝基础分＋各变量分数和。
train$score=base_score+train$x2_score+train$x3_score+train$x4_score+train$x5_score+train$x6_score+train$x7_score+train$x8_score+train$x9_score+train$x10_score4、同样的，计算测试集各变量的信用评分：
test$x2_score=(-1)*B*test$x2_woe*coe[2]
test$x3_score=(-1)*B*test$x3_woe*coe[3]
test$x4_score=(-1)*B*test$x4_woe*coe[4]
test$x5_score=(-1)*B*test$x5_woe*coe[5]
test$x6_score=(-1)*B*test$x6_woe*coe[6]
test$x7_score=(-1)*B*test$x7_woe*coe[7]
test$x8_score=(-1)*-B*test$x8_woe*coe[8]
test$x9_score=(-1)*B*test$x9_woe*coe[9]
test$x10_score=(-1)*B*test$x10_woe*coe[10]5、计算测试集各观测样本的总分：
test$score=base_score+test$x2_score+test$x3_score+test$x4_score+test$x5_score+test$x6_score+test$x7_score+test$x8_score+test$x9_score+test$x10_score
summary(test$score)(九)打印评分卡
需要先将各个变量的实际变量变化范围和信用评分组合在一起，代码如下，其它变量代码一致：
x2_s=tapply(train$x2_score,train$x2_rank,max)
x2_t=tapply(train$x2,train$x2_rank,range)
names(x2_s)=x2_t最后将各变量评分规则合并成一个list：
score_list=list(x2=x2_s,x3=x3_s,x4=x4_s,x5=x5_s,x6=x6_s,x7=x7_s,x8=x8_s,x9=x9_s,x10=x10_s)最后把评分卡打印出来，每个变量的区间下面对应着这个区间的信用评分：
> score_list
$x2
c(21, 32) c(33, 38) c(39, 43) c(44, 47) c(48, 51) c(52, 55) c(56, 59) c(60, 64) c(65, 71) c(72, 102)
-13.0443064 -10.0048563 -8.1052001 -4.0779288 -3.2674088 -0.5319038 6.2815301 12.4110876 22.6945602 27.5830090
$x3
c(0, 0) c(1, 1) c(2, 2) c(3, 3) c(4, 13)
14.84904 -24.81542 -45.98320 -55.89214 -71.37305
$x4
c(0, 0.013594562) c(0.013597281, 0.090909091) c(0.09092792, 0.159954299) c(0.159956689, 0.219178082)
20.8327762 -0.4528864 1.8568344 5.2081940
c(0.219180026, 0.275072557) c(0.27507641, 0.333288895) c(0.333296177, 0.399334443) c(0.399336476, 0.483461963)
10.5522540 10.0087903 2.8078959 -2.2191435
c(0.483471394, 0.621075785) c(0.621126291, 1)
-12.8166862 -22.8707651
$x5
c(0, 2400) c(2401, 3208) c(3209, 4000) c(4001, 4800) c(4801, 5563) c(5564, 6500) c(6501, 7700)
-5.0877876 -5.9124137 -3.1429147 -2.8628530 -0.3267387 0.9490980 3.2829455
c(7701, 9282) c(9283, 11880) c(11881, 1794060)
5.2122595 6.7681578 6.7681578
$x6
c(0, 3) c(4, 4) 5:6 c(7, 7) c(8, 8) c(9, 9) c(10, 10) 11:12 c(13, 15) c(16, 57)
-12.8976555 -0.3313627 1.3254508 3.9763523 8.6409194 3.1352008 3.6704790 3.0587325 2.2175811 1.0960458
$x7
c(0, 0) c(1, 1) c(2, 2) c(3, 3) c(4, 17)
10.68024 -55.38219 -73.26872 -85.26963 -105.91234
$x8
c(0, 0) c(1, 1) c(2, 2) c(3, 3) c(4, 29)
6.644972 -5.523458 -6.364594 -2.803786 10.009515
$x9
c(0, 0) c(1, 1) c(2, 2) c(3, 11)
5.197519 -34.669594 -51.021340 -64.861925
$x10
c(0, 0) c(1, 1) c(2, 2) c(3, 3) c(4, 20)
3.582490 -1.529603 -3.985017 -5.937273 -9.278247 三、评分卡模型评估指标
一般信用评分卡模型评估指标大致可以分成两类：
- 一类是预测能力指标，表示模型对违约事件的预测能力，如AUC、K-S指标、GINI系数和分离度；
- 一类是稳定性指标，表示模型在训练样本和测试样本中预测能力的一致性，如roc.test和PSI指标。
(一)ROC曲线
ROC曲线上面已经介绍过了，主要通过AUC来衡量模型的预测能力，AUC越接近1模型预测效果越好。
r1=roc(train$y1,train$score, plot=TRUE, print.thres=F, print.auc=F,col="3",lwd=5)
r2=roc(test$y1,test$score, plot=TRUE, print.thres=TRUE, print.auc=TRUE,add=TRUE,col="1")
legend("bottomright", legend=c("train", "test"), col=c("3", "1"), lwd=2)
roc.test可以用来比较训练样本和测试样本的roc是否一致，p-value>0.05说明两者一致。
> roc.test(r1,r2)
DeLong's test for two ROC curves
data: r1 and r2
D = 0.85911, df = 65072, p-value = 0.3903
alternative hypothesis: true difference in AUC is not equal to 0
sample estimates:
AUC of roc1 AUC of roc2
0.8123469 0.8071067(二)K-S指标
K-S指标主要是测量违约和正常用户分布的最大差距，按分数顺序分别画出违约用户和正常用户的累积百分比线图，两者之间的最大差距值就是K-S指标。计算公式为max（正常用户累计－违约用户累计），K-S指标越大表示模型的区分能力越强。
pred_Tr<-prediction(train$score,train$y1)
tpr<-performance(pred_Tr,measure='tpr')@y.values[[1]]
fpr<-performance(pred_Tr,measure='fpr')@y.values[[1]]
ks<-(tpr-fpr)
depth<-performance(pred_Tr,measure='rpp')@y.values[[1]]
plot(depth,ks,type='l',main='train k-s',ylab='ks',xlab='depth',ylim=c(0,1))
lines(depth,tpr,col='2')
lines(depth,fpr,col='3')
kslable<-c(paste("KS:",round(max(ks),3),sep=""),paste("tpr:",round(tpr[which(ks==max(ks))],3),sep=""),paste("fpr:",round(fpr[which(ks==max(ks))],3),sep=""))
legend(0.0,1.0,c(kslable),bty="n",ncol = 1)
legend(0.4,1.0,legend=c("TPR", "FPR","TPR-FPR"), col=c("2", "3","1"),cex=0.6,lwd=2)
(三)基尼系数
GINI系数是正常用户累计百分比与违约用户累计百分比通过模型选择与随机选择时的差值，对角线下面的那条红色曲线也叫洛伦兹曲线，对角线表示完全没有好坏区分能力的随记选择结果，两者之间的区域除以对角线以下的面积即为GINI系数，GINI系数越大表示模型区分能力越强。
pred_Tr<-prediction(train$score,train$y1)
fnr<-performance(pred_Tr,measure='fnr')@y.values[[1]] #FN/(FN+TP)=FN/P
tnr<-performance(pred_Tr,measure='tnr')@y.values[[1]] #TN/(TN+FP)=TN/N
plot(tnr,fnr,type='l',main='train_Gini',xlab='x-tnr',ylab='y-fnr',yaxs="i",xaxs="i",col=2)
abline(a=0,b=1,col=1)
kslable1<-paste("Gini:",round(Gini(fnr),3),sep="")
legend(0.5,0.5,c(kslable1),bty="n",ncol = 1)
legend("bottomright", legend=c("RandomSelection", "ModelPick"), col=c("1", "2"), cex=0.75,lwd=2)
(四)分离度指标
分离度主要是基于评分的好坏用户的整体分布的差别，通常分离度越高，评分模型的排序能力越强，计算公式如下图所示：
train_s0=train$score[train$y==0]
train_s1=train$score[train$y==1]
plot(density(train_s0),main ="train_divergence",col=1)
lines(density(train_s1),col=2)
divergence=((mean(train_s0)-mean(train_s1))^2)/(0.5*(var(train_s0)+var(train_s1)))
kslable<-paste("divergence:",round(divergence,3),sep="")
legend(400,0.013,c(kslable),bty="n",ncol = 1)
legend(350,0.017, legend=c("Goods", "Bads"), col=c("1", "2"), cex=0.75,lwd=2)
(五)稳定度指标PSI
PSI的计算公式为：sum(（实际占比-预期占比）* ln(实际占比/预期占比)),实际占比就是各分数组在测试集的占比，预期占比就是各分数组在训练集的占比。一般认为psi小于0.1时候模型稳定性很高，0.1-0.25一般，大于0.25模型稳定性差，建议重做。
#提取train$score的十分位数
score_q=quantile(train$score,probs = seq(0, 1, 0.1),names=FALSE)
#用无穷小和无穷大分别替换最小值和最大值
score_q[1]=-Inf
score_q[length(score_q)]=Inf
#生成分组组号
train$score_rank=cut(train$score,score_q,labels = F)
test$score_rank=cut(test$score,score_q,labels = F)
#计算PSI
a=b=c=0
for(i in 1:10){
a=sum(train$score_rank==i)/length(train$score);
b=sum(test$score_rank==i)/length(test$score);
c=c+(b-a)*log(b/a);
}
#画出分布图
plot(density(train$score),main='PSI',col=2,lwd=5)
lines(density(test$score),col=1,lwd=2)
kslable<-paste("PSI:",round(c,5),sep="")
legend(400,0.013,c(kslable),bty="n",ncol = 1)
legend(350,0.017, legend=c("train", "test"), col=c("2", "1"), cex=0.75,lwd=2)
附录：参考文献
- 信用卡评分模型(R语言)_程松
- 《互联网金融时代消费信贷评分建模与应用》_单良,茆小林
- 《SAS开发经典案例解析》_杨池然
- Scorecard 评分卡模型_textboy
版权声明：本文为大枫原创文章，转载请注明来源。
