用R做评分卡模型
評分卡模型流程
參考書：《信用風險評分卡研究》
數據讀取和EDA
Kaggle數據cs-training，共15萬條
https://www.kaggle.com/c/GiveMeSomeCredit/data
rawdata<-read.csv(file.choose(),header = T)#讀取原始數據
data<-rawdata[,-1]#去除序列號列
#剔除缺失值超過20%或單一值超過80%的變量
l=c()
for(j in 2:ncol(data))
{
if(length(which(is.na(data[,j])))/nrow(data)>=0.2|max(table(data[,j]))>=0.8*nrow(data))
? ? {l=c(l,j)}
}
#剔除有缺失值的行（個體）
data=data[,-l]
l=c()
for(i in 1:nrow(data))
{
? if(length(which(is.na(data[i,])))>0)
? ? {l=c(l,i)}
}
data=data[-l,]
注：數據EDA分析和處理要根據具體數據和業務去做，這里主要介紹評分卡模型，數據處理簡化處理
數據準備（計算woe、分箱）
連續變量的數據準備
定義qua函數：返回變量的K等分位點
#fwd[1]-1是為了方便統一后續的區間編寫（前開后閉），輸出的結果會改寫回來
qua<-function(var,K)
{
? fwd=round(unique(quantile(var,0:K/K)),4)
? fwd[1]=fwd[1]-1
? return(fwd)
}
這里數據較少，大家可以觀察每個變量的分布，自行分箱。
定義函數binning：對連續變量分箱（K等分），返回變量指標var_index(list格式):WOE(每個個體在該變量上的woe值），index（每個箱的最小值，最大值，壞的數量，壞的比率，好的數量，好的比率，總數量，woe值，iv值），IV（該變量的iv值），quantitl（分位點）
#usqua=T:使用默認K等分分箱，否則自行定義分位點fwd（注：包括前后端點）
#若某一箱全為好或全為壞，則woe為0，iv為0
binning<-function(var,K,fwd,usequa=T)
{
? if(usequa)
? {
? ? fwd=qua(var,K)
? }
? numallgood=length(which(data[,1]==0))
? numallbad=length(which(data[,1]==1))
? WOE=data.frame(matrix(0,nrow(data),1));var_index=list()
? num_bad=0;num_good=0;rate_bad=0;rate_good=0;woe=0;iv=0;IV=0
? index=data.frame(matrix(0,length(fwd)-1,9))
? colnames(index)=c("min","max","num_bad","rate_bad","num_good","rate_good","count","woe","iv")
? for(k in 1:(length(fwd)-1))
? {
? ? num_bad=length(which(data[,1][which(var>fwd[k]&var<=fwd[k+1])]==1))
? ? num_good=length(which(data[,1][which(var>fwd[k]&var<=fwd[k+1])]==0))
? ? rate_bad=num_bad/numallbad
? ? rate_good=num_good/numallgood
? ? if(num_bad!=0&num_good!=0)
? ? {
? ? ? woe=log(rate_bad/rate_good)
? ? ? iv=(rate_bad-rate_good)*woe
? ? }
? ? if(num_bad==0|num_good==0)
? ? {
? ? ? woe=0;iv=0
? ? }
? ? WOE[which(var>fwd[k]&var<=fwd[k+1]),1]=woe
? ? index$min[k]=fwd[k];index$max[k]=fwd[k+1]
? ? index$num_bad[k]=num_bad;index$rate_bad[k]=rate_bad
index$num_good[k]=num_good;index$rate_good[k]=rate_good
? ? index$count[k]=num_bad+num_good
? ? index$woe[k]=woe;index$iv[k]=iv
? ? IV=IV+iv
? }
? index$min[1]=index$min[1]+1
? var_index[[1]]=WOE
? var_index[[2]]=index
? var_index[[3]]=IV
? var_index[[4]]=fwd
? names(var_index)=c("WOE","index","IV","quantitl")
? return(var_index)
}
運行代碼
age_index=binning(data$age,K=10,usequa=T)
age_index
得到結果：
定義合箱函數closing：根據woe相鄰最近的兩箱進行合并
var 是某一變量，var_index是binning函數對var返回的結果
closing<-function(var,var_index)
{
? fwd=var_index$quantitl
? woe=var_index$index$woe
? diff_woe=abs(diff(woe))
? num_min=which.min(diff_woe)
? fwd=fwd[-(num_min+1)]
? var_index_byclosing=binning(var,fwd=fwd,usequa = F)
? return(var_index_byclosing)
}
運行代碼
age_index_byclosing=closing(data$age,age_index)
age_index_byclosing
得到結果：由10箱合并為9箱
定義按單調性函數clobymono：根據變量各箱woe值的單調性合箱，非單調且箱數大于K_min時循環合箱
（注：有些變量U型，或倒U型，這里暫不考慮）
clobymono<-function(var,K,K_min)
{
? var_index=binning(var,K,usequa = T)
? woe=var_index$index$woe
? while(length(which(woe!=sort(woe)))>0&
? ? ? ? length(which(woe!=sort(woe,decreasing = T)))>0&
? ? ? ? length(var_index$quantitl)>=K_min+1)
? {
? ? var_index=closing(var,var_index)
? ? woe=var_index$index$woe
? ? if(length(var_index$quantitl)==(K_min+1))
? ? {break}
? }
? return(var_index)
}
運行代碼
age_index=clobymono(data$age,10,5)
age_index
得到結果
這里結果與binning(data$age,10,ueage=T)相同，是因為woe單調，無需合箱。我們帶入另一個變量再觀察：
運行代碼
MonthlyIncome_index=clobymono(data$MonthlyIncome,10,5)
MonthlyIncome_index
得到結果
變量MonthlyIncome一開始等分為10箱，直到合為6箱woe才單調
定義連續變量指標函數convarindex
#all continuous variable's index?
convarindex<-function(data,K,K_min)
{
? convar_index=list()
? WOEmatrix=as.data.frame(data[,1],norw=nrow(data),ncol=1)
? colnames(WOEmatrix)[1]=colnames(data)[1]
? for(j in 2:ncol(data))
? {
? ? var_index=clobymono(data[,j],K,K_min)
? ? WOEmatrix=cbind(WOEmatrix,var_index$WOE)
? ? colnames(WOEmatrix)[j]=colnames(data)[j]
? ? convar_index[[j-1]]=var_index[-1]
? ? names(convar_index)[j-1]=colnames(data)[j]
? }
? convar_index[[length(convar_index)+1]]=WOEmatrix
? names(convar_index)[length(convar_index)]="WOEmatrix"
? return(convar_index)
}
這里取4個連續變量演示，data的第一列為目標值
運行代碼：
convar_index=convarindex(data[,1:5],10,5)
convar_index
得到結果：所有連續變量指標和連續變量woe矩陣
離散變量的數據準備
定義分類函數classifing：返回離散變量指標
classifing<-function(var)
{
? classify<-unique(var)
? numallgood=length(which(data[,1]==0))
? numallbad=length(which(data[,1]==1))
? WOE=data.frame(matrix(0,nrow(data),1));var_index=list()
? num_bad=0;num_good=0;rate_bad=0;rate_good=0;woe=0;iv=0; ?IV=0
? index=data.frame(matrix(0,length(classify),8))
? colnames(index)=c("classify","num_bad","rate_bad","num_good","rate_good","count","woe","iv")
? for(k in 1:length(classify))
? {
? ? num_bad=length(which(data[,1][which(var==classify[k])]==1))
? ? num_good=length(which(data[,1][which(var==classify[k])]==0))
? ? rate_bad=num_bad/numallbad
? ? rate_good=num_good/numallgood
? ? if(num_bad!=0&num_good!=0)
? ? {
? ? ? woe=log(rate_bad/rate_good)
? ? ? iv=(rate_bad-rate_good)*woe
? ? }
? ? if(num_bad==0|num_good==0)
? ? {
? ? ? woe=0;iv=0
? ? }
? ? WOE[which(var==classify[k]),1]=woe
? ? index$classify[k]=classify[k]
? ? index$num_bad[k]=num_bad;index$rate_bad[k]=rate_bad
? ? index$num_good[k]=num_good;index$rate_good[k]=rate_good
? ? index$count[k]=num_bad+num_good
? ? index$woe[k]=woe;index$iv[k]=iv
? ? IV=IV+iv
? }
? var_index[[1]]=WOE
? var_index[[2]]=index
? var_index[[3]]=IV
? var_index[[4]]=classify
? names(var_index)=c("WOE","index","IV","classify")
? return(var_index)
}
運行代碼：
numofdep_index=classifing(data$NumberOfDependents)
numofdep_index
得到結果：
共13類，與連續變量相同，若某一類全為好或全為壞，則woe為0，iv為0（NumberOfDependents其實也是離散變量，這里只是作為演示）
定義所有離散變量指標函數clavarindex，與連續變量一致，返回所有變量指標和woe矩陣（合為list）
clavarindex<-function(data)
{
? cladata_index=list()
? WOEmatrix=data.frame(data[,1])
? for(j in 2:ncol(data))
? {
? ? classifyvar_index=classifing(data[,j])
? ? cladata_index[[j-1]]=classifyvar_index[-1]
? ? WOEmatrix=cbind(WOEmatrix,classifyvar_index[[1]])
? ? colnames(WOEmatrix)[j]=colnames(data)[j]
? }
? cladata_index[[length(cladata_index)+1]]=WOEmatrix
? names(cladata_index)=c(colnames(data)[-1],"WOEmatrix")
? return(cladata_index)
}
運行代碼：
cladata_index=clavarindex(data[,c(1,6,7)])
cladata_index
這里，假設第6和第7列為離散變量，第一列為目標值，結果與連續變量一致，這里不再給出圖形結果。
以上出去數據處理以外，從qua函數開始的所有函數都是為了下面allvarindex函數調取所用：
定義所有變量指標函數allvarindex:返回所有變量指標和woe矩陣
（注：concol為連續變量的列，不包括第一列目標值列）
#combining continuous and classify indexes
allvarindex<-function(data,K,K_min,concol)
{
? conindex=convarindex(data[,c(1,concol)],K,K_min)
? claindex=clavarindex(data[,-concol])
? all_index=append(conindex[-length(conindex)],claindex[-length(claindex)])
? WOEmatrix=cbind(conindex$WOEmatrix,claindex$WOEmatrix[,-1])
? colnames(WOEmatrix)=c(colnames(conindex$WOEmatrix),colnames(claindex$WOEmatrix)[-1])
? all_index[[length(all_index)+1]]=WOEmatrix
? names(all_index)[length(all_index)]="WOEmatrix"
? return(all_index)
}
運行代碼：
all_index=allvarindex(data,10,5,c(2:7))
all_index
給出all_index的內容：
得出所有變量的指標和WOEmatrix（woe矩陣）
可以看出變量NumberRealEstateLoansOrLines的值為0、1、2的較多
變量選擇
下面對woe矩陣建立逐步logistic回歸模型：
運行代碼：
library(caret)
set.seed(1234)
woeMatrix<-all_index$WOEmatrix
splitIndex<-createDataPartition(woeMatrix[,1],time=1,p=0.5,list=FALSE)?
train<-woeMatrix[splitIndex,];traindata=data[splitIndex,]?
test<-woeMatrix[-splitIndex,];testdata=data[splitIndex,]
glm_fit=glm(woeMatrix[,1]~.,family = binomial,woeMatrix[,-1])
step_glm_fit=step(glm_fit)
pre=predict(step_glm_fit,test,type="response")
注：train、test分別為woe矩陣的訓練集和測試集；traindata、testdata分別為數據data的訓練集和測試集。
運行代碼
summary(step_glm_fit)
1
結果如下：
逐步回歸后不通過的變量不再出現在模型中，在轉化為評分卡時，只調用留下來的變量的woe。
模型評估
ROC曲線
計算ROC畫出圖形：
library(pROC)
modelroc <- roc(test[,1],pre)
plot(modelroc, print.auc=T, auc.polygon=T, grid=c(0.2, 0.1),
? ? ?grid.col=c("red", "yellow"), max.auc.polygon=T,
? ? ?auc.polygon.col="seagreen2", print.thres=T)
結果如下：
KS曲線
計算KS
定義KS函數ksindex：這里pre是對test的預測，所以y應該是test[,1]
ksindex<-function(y,pre,K=10)
{
? per_bad=c();per_good=c()
? y2=y[order(rank(pre),decreasing = T)]
? numallgood=length(which(y2==0))
? numallbad=length(which(y2==1))
? for(i in 1:K)
? {
? ? fwd=floor(quantile(0:length(y2),0:K/K))
? ? per_bad[i]=length(which(y2[(fwd[i]+1):fwd[i+1]]==1))*100/numallbad
? ? per_good[i]=length(which(y2[(fwd[i]+1):fwd[i+1]]==0))*100/numallgood
? }
? acc_per_bad=c();acc_per_good=c()
? for(i in 1:K)
? {
? ? acc_per_bad[i]=sum(per_bad[1:i])
? ? acc_per_good[i]=sum(per_good[1:i])
? }
? per_diff=acc_per_bad-acc_per_good
? ksindex=data.frame(acc_per_bad,acc_per_good,per_diff)
? colnames(ksindex)=c("acc_per_bad","acc_per_good","per_diff")
? return(ksindex)
}
運行代碼：
ks_index=ksindex(test[,1],pre)
ks_index
有如下結果：
定義ks作圖函數ksplot：
ksplot<-function(ksindex,goodcol="blue",badcol="darkgreen",kscol="red")
{
? plot(1:nrow(ksindex),ksindex$acc_per_good,type="l",col=goodcol,lwd=2,ylab ="acc_per ?%",
? ? ? ?xlim=c(1,10),ylim=c(0,100),xlab="Decili",main="KS curve")
? points(1:nrow(ksindex),ksindex$acc_per_good,cex=0.3)
? text(floor(nrow(ksindex)*0.2)+1,ksindex$acc_per_good[floor(nrow(ksindex)*0.2)]-10,"goodrate",col = goodcol)
? lines(1:nrow(ksindex),ksindex$acc_per_bad,"l",col=badcol,lwd=2)
? points(1:nrow(ksindex),ksindex$acc_per_bad,cex=0.3)
? text(floor(nrow(ksindex)*0.8)-1,ksindex$acc_per_bad[floor(nrow(ksindex)*0.8)]-10,"badrate",col=badcol)
? lines(1:nrow(ksindex),ksindex$per_diff,col=kscol,lwd=2)
? points(1:nrow(ksindex),ksindex$per_diff,cex=0.3)
? points(which.max(ksindex$per_diff),ksindex$per_diff[which.max(ksindex$per_diff)],cex=1,pch=16)
? text(which.max(ksindex$per_diff),ksindex$per_diff[which.max(ksindex$per_diff)]+5,
? ? ? ?paste("KS=",round(ksindex$per_diff[which.max(ksindex$per_diff)],2),"%",sep = ""),col=kscol)
}
運行代碼：
ksplot(ks_index)
1
輸出圖形
評分卡刻度
定義計算變量woe函數（個體某一變量的值對應的woe值，若離散變量的類不在訓練集中，則輸出woe為0）get_varwoe:
get_varwoe<-function(colvar,var_value)
{
? if(names(colvar)[3]=="quantitl")
? {
? ? k=min(which(var_value<=colvar$quantitl))
? ? if(k==1)
? ? {k=2}
? ? if(length(k)==0)
? ? {k=length(colvar$quantitl)}
? ? woe_var=colvar$index$woe[k-1]
? }
? if(names(colvar)[3]=="classify")
? {
? ? if(var_value%in%colvar$classify)
? ? {
? ? ? woe_var=colvar$index$woe[var_value==colvar$classify]
? ? }
? else
? ? {woe_var=0;print("new classify")}
? }
? return(woe_var)
}
例如：得出第一個測試樣本的age值所對應的woe值；
運行代碼：
get_varwoe(all_index$age,test$age[1])
得出結果：
定義計算變量得分函數var_score：colvar_index為allvarindex函數返回值的某一元素變量的index元素
var_score<-function(colvar_index)
{
? df<-data.frame()
? if(colnames(colvar_index)[1]=="min")
? {
? ? df[1,1]=paste("[",as.character(round(colvar_index$min[1],4)),",",
? ? ? ? ? ? ? ? ? as.character(round(colvar_index$max[1],4)),"]",sep = "")
? ? df[2:nrow(colvar_index),1]=paste("(",as.character(round(colvar_index$min[2:nrow(colvar_index)],4)),
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?",",
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?as.character(round(colvar_index$max[2:nrow(colvar_index)],4)),"]",
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?sep = "")
? }
? if(colnames(colvar_index)[1]=="classify")
? {
? ? df=colvar_index$classify
? }
? df=cbind(df,colvar_index$woe)
? rownames(df)=as.character(c(1:nrow(df)))
? return(df)
}
運行代碼：
var_score(all_index$age$index)
輸出如下：
定義評分卡函數score_card:
注：若theta取1，表明theta為默認值：樣本中的壞好比
model_fit為模型擬合對象，allvar_index為函數allvarindex返回的結果，若list取真，則輸出list格式，若list為假輸出data.frame格式。
scorecard<-function(y,theta=1,PDO,p0,model_fit,allvar_index,list=T)
{
? if (theta==1)
? {
? ? theta=length(which(y==1))/length(which(y==0))
? }
? B=PDO/log(2);A=p0+B*log(theta)
? beta<-model_fit$coefficients
? namesofmodelcoef=names(beta)
? allvar_index=allvar_index[names(allvar_index)%in%namesofmodelcoef]
? score_card=list()
? varscore=data.frame()
? for(i in 1:length(allvar_index))
? {
? ? varscore=var_score(allvar_index[[i]]$index)
? ? varscore[,2]=floor(-varscore[,2]*B*beta[i+1]+(A-B*beta[1])/length(allvar_index))
? ? varscore=cbind(as.data.frame(rep(names(allvar_index)[i],nrow(varscore))),varscore)
? ? colnames(varscore)=c("varname","interval","score")
? ? score_card[[i]]=varscore
? ? names(score_card)[i]=names(allvar_index)[i]
? }
? if(!list)
? {
? ? sc=score_card[[1]]
? ? for(i in 2:length(score_card))
? ? {
? ? ? sc=rbind(sc,score_card[[i]])
? ? }
? ? score_card=sc
? ? rownames(score_card)=as.character(1:nrow(score_card))
? }
? return(score_card)
}
運行代碼：
score_card=scorecard(data[,1],1,50,580,step_glm_fit,all_index,T)
score_card
輸出結果：
運行代碼：
score_card=scorecard(data[,1],1,50,580,step_glm_fit,all_index,F)
score_card
輸出結果：
代碼流程總結
數據預處理完成后，按以下代碼即可得到評分卡：
#計算woe等相關指標
all_index=allvarindex(data,K=10,K_min=5,concol=c(2:7))
#建立logistic回歸模型
library(caret)
set.seed(1234)
woeMatrix<-all_index$WOEmatrix
splitIndex<-createDataPartition(woeMatrix[,1],time=1,p=0.5,list=FALSE)?
train<-woeMatrix[splitIndex,];traindata=data[splitIndex,]?
test<-woeMatrix[-splitIndex,];testdata=data[splitIndex,]
glm_fit=glm(woeMatrix[,1]~.,family = binomial,woeMatrix[,-1])
step_glm_fit=step(glm_fit)
pre=predict(step_glm_fit,test,type="response")
#觀察ROC曲線和KS曲線
library(pROC)
modelroc <- roc(test[,1],pre)
plot(modelroc, print.auc=T, auc.polygon=T, grid=c(0.2, 0.1),
? ? ?grid.col=c("red", "yellow"), max.auc.polygon=T,
? ? ?auc.polygon.col="seagreen2", print.thres=T)
ks_index=ksindex(test[,1],pre)
ksplot(ks_index)
#woe轉化為評分卡
score_card=scorecard(data[,1],1,50,580,step_glm_fit,all_index,F)
score_card
當然，這只是個不精細的通用代碼，在數據預處理、變量分析與篩選的過程中往往要精細化，比如觀察每個變量的分布，對異常值的分箱考慮等。Kaggle的cs-training.csv的變量都是連續變量，但家屬數量“類別”較少，這里就當作離散變量演示。
?
總結
- 上一篇： 评分卡模型剖析之一（woe、IV、ROC
- 下一篇： 大数据风控-提高授信审查效率，做好这7点