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

本文记录了从数据清洗到输出各变量得分的评分卡建立过程,建模数据来自数据挖掘竞赛网站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_bandodds_predictp_predictodds_calculatep_calculate
5000.160.1380.3710.271
5250.1130.1020.170.145
5500.080.0740.120.108
5750.0570.0540.0680.064
6000.040.0380.0560.053
6250.0280.0280.0340.033
6500.020.020.0310.03
6750.0140.0140.0130.013
7000.010.010.0070.007

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

  • 0
    点赞
  • 8
    收藏
    觉得还不错? 一键收藏
  • 3
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值