基于R的Score Card开发流程记录

最新推荐文章于 2021-03-05 12:09:23 发布

abbyatlarge

最新推荐文章于 2021-03-05 12:09:23 发布

阅读量1.1k

点赞数

文章标签： data analysis R credit analysis

本文链接：https://blog.csdn.net/abbyatlarge/article/details/82414963

版权

本文记录了从数据清洗到输出各变量得分的评分卡建立过程，建模数据来自数据挖掘竞赛网站Kaggle上由捷信消费金融提供的数据集，数据来源Data Source。
捷信提供了不同属性的多个数据集，为使建模过程更加清晰，本次建模数据仅采用其中的application_train.csv数据集，同时需要下载数据字典文件HomeCredit_columns_description.csv备查。

本次建模用到以下Package，数据读入后首先检查数据概况：

options(scipen=10, digits=5)
data_path="data_path"
library(data.table)
library(dplyr)
library(corrplot)
library(mice)
library(caret)
library(ggplot2)
library(InformationValue)
library(smbinning)
library(pROC)
data_raw<-fread(data_path,header=TRUE,sep=",",skip=0,fill=TRUE)
summary(data_raw)

共计122各变量，307511条记录，其中第一个变量SK_ID_CURR是编号，第二个变量Target是标签，不涉及数据清洗的过程。剩余变量仍有120个，数量较多，可以首先查看变量特征。:View(arrange(cln_info(data_raw),type,-incmplt))。执行该命令后，按照变量类型，缺失比例降序排列显示变量信息，下表是部分结果：
这里写图片描述
120个自变量中，字符型变量有16个，数值型有104个，首先处理字符型变量中相关系数过高的变量。
执行命令data_cln<-cln_prcss1(data_raw,0.3,'TARGET')的作用是剔除data_raw中自变量相关性高于0.3的数值型自变量，并保证因变量不会被剔除。

自变量取值波动过小会出现以下警告，不影响建模

Warning message:

In cor(temp, use = “pairwise.complete.obs”) : 标准差为零

执行函数后，剔除相关系数高于0.3的数值型变量，data_cln数据集包含63个变量，即剔除了59个数值型变量，下面考虑剩余变量的缺失问题。
执行命令data_cln<-cln_prcss2(data_cln),其作用是处理缺失数据，具体规则为数值型缺失数据用中位数替换，字符型用众数，这一规则可以避免异常值干扰，异常值在分bin中统一处理。
执行View(arrange(cln_info(cln_prcss2(data_cln)),-unq))命令查看按照变量取值的unique值降序排序的变量信息，下表是部分结果：
这里写图片描述
下面开始切分数据集，由于样本不均衡——Target为1的样本约占全部样本的8%——建模过程采用采用caret包中createDataPartition函数均衡抽样得到训练集和测试集。

# TARGET样本不均衡
table(data_cln$TARGET)/length(data_cln$TARGET)
# 采用caret包中createDataPartition函数分割数据
set.seed(1) 
spliter<-createDataPartition(data_cln$SK_ID_CURR,time=1,p=0.7,list=FALSE) 
train<-data_cln[spliter,] # 训练集
test<-data_cln[-spliter,] # 测试集

抽取所有数据的70%作为训练集，剩余作为测试集。开始基于数据集train选取变量，讲数据集切分为train_num，和train_chr，前者包含所有数值型变量，后者包含所有字符型变量。首先考虑数值型变量：

# 选取定量指标
# 方法一，自变量的逐步回归法
base.mod<-lm(TARGET~1,data = train_num)
# 获取线性回归模型的截距
all.mod<-lm(TARGET~.,data = train_num)
# 获取完整的线性回归模型
step.mod<-step(base.mod,scope = list(lower=base.mod,upper=all.mod),
               direction = "both",trace = 0,steps = 1000)
# 采用双向逐步回归法，筛选变量
shortlistedVars<-names(unlist(step.mod[[1]]))
# 获取逐步回归得到的变量列表
shortlistedVars<-shortlistedVars[!shortlistedVars %in% "(Intercept)"]
# 删除逐步回归的截距
print(shortlistedVars)
# 输出逐步回归后得到的变量
rm(all.mod,base.mod,step.mod)

# 方法二，计算变量间的相对重要性，回归法
library(relaimpo)
lm.mod<-lm(TARGET~.,data = train_num)  # 线性回归
relImportance<-calc.relimp(lm.mod,type = "lmg",rela = TRUE)
# 计算自变量间的相对重要性，通过所有可能的自变量顺序组合计算序列平方和的均值
print(sort(relImportance$lmg,decreasing = TRUE))
# 排序并输出自变量间的相对重要性
rm(lm.mod)

# 方法三，自变量间的广义交叉验证法
library(earth)
mars.mod<-earth(TARGET~.,data = train_num)
ev<-evimp(mars.mod)
print(ev)
# 经过自变量间的广义交叉验证后，获取自变量的重要性
rm(mars.mod)

# # 方法四，随机森林法
# # train_num全量数据跑不了，100000条在本地是极限了
# library(party)
# cf1<-cforest(TARGET~.,data = train_num,
#              controls = cforest_unbiased(mtry=2,ntree=50))
# varimp(cf1)
# varimp(cf1,conditional = TRUE)
# #经过变量间的相关系数调整后，获取自变量的重要性
# varimpAUC(cf1)

综合以上方法，选取以下数值变量参与建模：

chosen_num<-c('EXT_SOURCE_3',
              'EXT_SOURCE_2',
              'EXT_SOURCE_1',
              'DAYS_ID_PUBLISH',
              'DAYS_LAST_PHONE_CHANGE',
              'DAYS_REGISTRATION',
              'AMT_GOODS_PRICE')

之后考虑字符型变量的选择，根据各变量的IV值决定是否剔除：

IV	预测能力
<0.03	无预测能力
0.03~0.09	低
0.1~0.29	中
0.3~0.49	高
>=0.5	极高

# 选取定性指标
factor_vars<-names(train_chr[,2:ncol(train_chr)])
# 获取所有名义变量
all_iv<-data.table(vars=rep(NA,ncol(train_chr)-1),
                   IV=rep(0,ncol(train_chr)-1),
                   strength=rep(NA,ncol(train_chr)-1))
# 初始化待输出的数据框
k=1
for(factor_var in factor_vars)
{
  temp<-InformationValue::IV(X=as.factor(train_chr[[factor_var]]),Y=as.factor(train_chr$TARGET))
  all_iv$vars[k]<-factor_var
  all_iv$IV[k]<-temp[1]
  # 计算每个指标的IV值
  all_iv$strength[k]<-attr(temp,"howgood")
  # 提取每个IV指标的描述
  k=k+1
}
all_iv<-all_iv[order(-all_iv$IV),]
print(all_iv)

选择IV>0.03的字符型变量参与建模：

chosen_chr<-c('OCCUPATION_TYPE',
              'ORGANIZATION_TYPE',
              'NAME_INCOME_TYPE',
              'NAME_EDUCATION_TYPE',
              'CODE_GENDER',
              'FLAG_DOCUMENT_3')
train_cln<-train[,c('TARGET',chosen_num,chosen_chr),with=FALSE]

建模数据的清洗及筛选工作就做完了，Train_cln将在下面的步骤中进行分bin和woe化，之后就可以用于建立Logistic回归模型。

以变量EXT_SOURCE_2为例，展示一个数值型变量的分bin及woe化过程：

# 初始化一个list用于存储分箱的结果
bin_result<-list()

# EXT_SOURCE_2
# 观察最优分箱
smbin<-smbin_info(df=train_cln,y='TARGET',x='EXT_SOURCE_2',p=0.05)
# 自变函数smbin_info()输出分bin信息，包括cutpoints，Odds，woe等，并画图展示

这里写图片描述

# 画图看分布
p<-ggplot(data=train_cln,mapping=aes(x =EXT_SOURCE_2,fill=factor(TARGET)))
p+geom_density(alpha=0.3)+scale_x_continuous(limits=c(min(data=train_cln$EXT_SOURCE_2), max(data=train_cln$EXT_SOURCE_2)), breaks=seq(min(data=train_cln$EXT_SOURCE_2), max(data=train_cln$EXT_SOURCE_2),(max(data=train_cln$EXT_SOURCE_2)-min(data=train_cln$EXT_SOURCE_2))/10))+geom_vline(xintercept=smbin$smbin_cuts)+theme(axis.text.x=element_text(angle=90))# 好坏样本分布密度
# p+geom_histogram(binwidth=0.02)+xlim(0,1)# 样本分布柱状图
# 分箱woe单调，记录分箱结果
bin_result$EXT_SOURCE_2_cuts<-smbin$smbin_cuts
bin_result$EXT_SOURCE_2_tags<-smbin$woe_info$Cutpoint[1:(length(smbin$smbin_cuts)+1)]
bin_result$EXT_SOURCE_2_woes<-smbin$woe_info$WoE[1:(length(smbin$smbin_cuts)+1)]

这里写图片描述

# 分箱woe单调，记录分箱结果
bin_result$EXT_SOURCE_2_cuts<-smbin$smbin_cuts
bin_result$EXT_SOURCE_2_tags<-smbin$woe_info$Cutpoint[1:(length(smbin$smbin_cuts)+1)]
bin_result$EXT_SOURCE_2_woes<-smbin$woe_info$WoE[1:(length(smbin$smbin_cuts)+1)]
# 原始数据分箱
train_cln$EXT_SOURCE_2<-binning(df=train_cln$EXT_SOURCE_2,
                                cuts=bin_result$EXT_SOURCE_2_cuts,
                                tags=bin_result$EXT_SOURCE_2_tags)
# 建模数据woe化
train_model$EXT_SOURCE_2<-woelizing(df=train_cln$EXT_SOURCE_2,
                                  tags=bin_result$EXT_SOURCE_2_tags,
                                  woes=bin_result$EXT_SOURCE_2_woes)

以变量OCCUPATION_TYPE为例，展示一个字符型变量的分bin及woe化过程：

# 定性变量降维并woe化
# OCCUPATION_TYPE
# 观察分类特征
temp<-arrange(WOETable(X=as.factor(train_cln$OCCUPATION_TYPE), Y=train_cln$TARGET),-WOE)
print(temp)
p<-ggplot(data=train_cln,mapping=aes(x = reorder(OCCUPATION_TYPE,TARGET,length),fill=factor(TARGET))) 
p+geom_bar(stat="count")+coord_flip()#
p<-ggplot(data=temp,mapping=aes(x=reorder(CAT,-WOE),y=WOE))
p+geom_bar(stat='identity')+theme(axis.text.x=element_text(angle=90))

这里写图片描述

# 记录降维结果
bin_result$OCCUPATION_TYPE_cuts<-c('Low-skill Laborers|Drivers|Waiters/barmen staff|Security staff|Cooking staff|Laborers|Sales staff|Cleaning staff',
                                   'Realty agents|Secretaries|IT staff|Private service staff|Medicine staff|High skill tech staff|Core staff|Managers|HR staff|Accountants',
                                   '')
bin_result$OCCUPATION_TYPE_tags<-c('OCCUPATION_1',
                                   'OCCUPATION_2',
                                   'OCCUPATION_blank')
# 降维
temp<-train_cln$OCCUPATION_TYPE
temp<-gsub(bin_result$OCCUPATION_TYPE_cuts[1],bin_result$OCCUPATION_TYPE_tags[1],temp)
temp<-gsub(bin_result$OCCUPATION_TYPE_cuts[2],bin_result$OCCUPATION_TYPE_tags[2],temp)
temp[temp=='']<-bin_result$OCCUPATION_TYPE_tags[3] # 空白占比太多了，单独成为一类
train_cln$OCCUPATION_TYPE<-temp
temp<-arrange(WOETable(X=as.factor(train_cln$OCCUPATION_TYPE), Y=train_cln$TARGET),-WOE)
p<-ggplot(data=temp,mapping=aes(x=reorder(CAT,-WOE),y=WOE))
p+geom_bar(stat='identity')+theme(axis.text.x=element_text(angle=90))
bin_result$OCCUPATION_TYPE_woes<-temp$WOE
# 建模数据woe化
train_model$OCCUPATION_TYPE<-woelizing(df=train_cln$OCCUPATION_TYPE,
                                     tags=bin_result$OCCUPATION_TYPE_tags,
                                     woes=bin_result$OCCUPATION_TYPE_woes)

所有变量woe化之后查看变量相关性：

View(cln_info(train_model))
head(train_model)
cor1<-cor(train_model,use='pairwise.complete.obs')
corrplot(cor1,method ='number')

这里写图片描述
ORGANIZATION_TYPEOCCUPATION_TYPE及NAME_INCOME_TYPE相关系数分别为0.46和0.65，建模数据剔除后两者；
ORGANIZATION_TYPE和DAYS_ID_PUBLISH相关系数为0.21，暂不处理。

chosen_chr<-chosen_chr[!chosen_chr %in% c('OCCUPATION_TYPE','NAME_INCOME_TYPE')]
chosen_var<-c('TARGET',chosen_num,chosen_chr)
train_model<-train_model[,chosen_var,with=FALSE]
bin_result<-bin_result[!grepl('OCCUPATION_TYPE|NAME_INCOME_TYPE',names(bin_result))]

截至这一步，建模数据train_model准备好了，同时bin_result记录了所有分bin及woe化的过程。
下面可以开始使用train_model数据集进行逻辑回归。

# 逻辑回归
model_LR<-glm(TARGET~.,data=train_model,family = binomial())
summary(model_LR)
# 画训练集的ROC和KS
pre_act<-data.table(prediction=exp(predict(model_LR,train_model))/(1+exp(predict(model_LR,train_model))),actual=train_model$TARGET)
# 画ROC曲线，训练集上AUC=0.741
model_LR_roc <- roc(pre_act$actual,pre_act$prediction)
plot.roc(model_LR_roc,print.auc=TRUE,
         axes=TRUE,xlab="1 - Specificity",
         grid=TRUE,grid.col='dark gray',
         print.thres=TRUE)

这里写图片描述

# 画KS曲线
pre_act_KS<-data.table(x=seq(0,1,0.001),
                       tpr=rep(0,1001),
                       fpr=rep(0,1001),
                       diff=rep(0,1001))
for (i in 1:1000){
  pre_act_KS$tpr[i+1]=nrow(pre_act[prediction<=i*0.001&actual==0])/nrow(pre_act[actual==0])
  pre_act_KS$fpr[i+1]=nrow(pre_act[prediction<=i*0.001&actual==1])/nrow(pre_act[actual==1])
  pre_act_KS$diff[i+1]=pre_act_KS$tpr[i+1]-pre_act_KS$fpr[i+1]
}

p<-ggplot(data=pre_act_KS)
p+geom_line(mapping=aes(x=x,y=tpr),col='green')+geom_line(mapping=aes(x=x,y=fpr),col='red')+geom_line(mapping=aes(x=x,y=diff),col='blue')+scale_x_continuous(limits=c(0,1), breaks=seq(0,1,0.1))+geom_point(mapping=aes(x=(which(pre_act_KS$diff==max(pre_act_KS$diff))-1)/1000,y=max(pre_act_KS$diff)),col='blue')+geom_text(aes(x=(which(pre_act_KS$diff==max(pre_act_KS$diff))-1)/1000+0.1,y=max(pre_act_KS$diff),label=paste('KS: ',round(max(pre_act_KS$diff),3),sep='')), size=4, vjust=0)
# 模型的KS在训练集上是0.359

这里写图片描述
下面开始交叉验证，将前面准备好的验证数据集test按照处理训练集类似的方法分bin并woe化后画出模型在验证集上的ROC曲线和KS曲线：

模型的KS在训练集上是0.359，在验证集上是0.360，KS在验证集上升高很少见，感觉可能是均衡抽样的功劳？
下面开始求评分卡，假设评分卡的分的公式为score=alpha-beta*log(odds)，其中odds=p/(1-p)即model_LR的预测值。
考虑odds=1/100分数为700，每下降50分odds翻一倍，得到方程组

700=alpha-beta*log(1/100)

50=beta*log(2)

解上述方程得到：

beta=50/log(2)

alpha=700+beta*log(1/100)

以EXT_SOURCE_3_scores为例，计算各变量各bin得分

# EXT_SOURCE_3_scores的各分箱得分
bin_result$EXT_SOURCE_3_scores<-bin_result$EXT_SOURCE_3_woes*(-beta)*model_LR$coefficients['EXT_SOURCE_3']

# score=alpha-beta*(Intercept)+各bin得分=560.414+各bin得分
score_min<-560.414
score_max<-560.414
for (i in names(bin_result[grepl('scores',names(bin_result))])){
  score_min<-score_min+min(unlist(bin_result[i]))
  print(min(unlist(bin_result[i])))
  score_max<-score_max+max(unlist(bin_result[i]))
  print(max(unlist(bin_result[i])))
}

得到各变量各bin得分，保存在bin_result$EXT_SOURCE_3_scores变量中，同时求得模型的理论最高分=826.3624，理论低分=316.2758。
之后给所有数据打分，并在得分上绘制ROC和KS图：
这里写图片描述

由于目前没有业务数据可以进一步验证模型，评分卡建模到这里暂时结束。以下是各分数段样本的Odds和违约率：

score_band	odds_predict	p_predict	odds_calculate	p_calculate
500	0.16	0.138	0.371	0.271
525	0.113	0.102	0.17	0.145
550	0.08	0.074	0.12	0.108
575	0.057	0.054	0.068	0.064
600	0.04	0.038	0.056	0.053
625	0.028	0.028	0.034	0.033
650	0.02	0.02	0.031	0.03
675	0.014	0.014	0.013	0.013
700	0.01	0.01	0.007	0.007

参考资料：
信用风险评级模型的开发
 信用卡评分模型（R语言）
数据挖掘模型中的IV和WOE详解
 真假正负例、混淆矩阵、ROC曲线、召回率、准确率、F值、AP
ROC曲线-阈值评价标准

abbyatlarge

关注

0
点赞
踩
8

收藏

觉得还不错? 一键收藏
3
评论
基于R的Score Card开发流程记录

本文记录了从数据清洗到输出各变量得分的评分卡建立过程，建模数据来自数据挖掘竞赛网站Kaggle上由捷信消费金融（Home Credit Home）提供的数据集，数据来源Data Source。捷信提供了不同属性的多个数据集，为使建模过程更加清晰，本次建模数据仅采用其中的application_train.csv数据集，同时需要下载数据字典文件HomeCredit_columns_descrip...
复制链接

扫一扫