利用R语言对贷款客户作风险评估（下）——零膨胀回归分析

利用R语言对贷款客户作风险评估（下）——零膨胀回归分析

前言

上一篇的分类预测是决定好坏客户的初步判断, 不足以直接决策, 因此还需要进一步分析. 通过随机森林, 对影响好坏客户的解释变量的重要性进行排序. 由结果可以得出, 六个月内的还款状态是决定客户是非为好客户的最为重要的影响因素.因此, 本部分将由六个月的还款状态产生新的新的因变量, 即逾期还款次数, 其他变量作为解释变量, 从而研究影响逾期还款次数的因素.

初步分析

将数据中的6个月的还款状态做合并处理得到新的因变量逾期还款次数, 并删除还款状态6个解释变量；再将因子变量性别、教育程度和婚姻转化为哑变量. 由此得到20个解释变量, 1个计数因变量的数据. 因变量。逾期还款次数的计数分布如图， 计数分布不均衡, 零计数所占的比例为66.44%.

先使用传统的计数分布模型对数据进行拟合数据, 然后与简单线性回归做比较.为了比较模型, 将全部数据按因变量分层抽样抽取5份, 使得每份样本中因变量的分布均匀, 然后采用5折交叉验证得到模型的测试集和训练集的拟合结果. 模型的拟合结果如表：
结果表明, 虽然线性回归模型理论上并不适用, 其结果优于泊松回归和负二项回归, 表明传统的计数分布模型不如简单的线性回归有效。零点观测的比例较大, 远远超过了poisson分布中零值的概率分布所允许的范围. 当事件发生数中含有大量的零值, 我们把这种现象称为零膨胀现象, 并且这一类特殊的计数数据广泛的存在于自然、社会科学的实际研究中.

零膨胀数据

由于零事件所占比例过高, 传统模型的预测能力不足以解决, 将会导致有偏的统计推断结果, 所以需要找到适合的模型用来拟合这一类计数数据。
零膨胀模型, 由取值为零的计数数据和取值服从某特殊分布的计数数据各占一定比例组成的混合分布, 然后分别建立logit回归模型和一般的计数回归模型. 其中,logit 回归模型拟合的是解释变量与事件发生与否的关系, 而计数回归模型拟合的是解释变量与事件发生次数的关系。

实现的R代码如下：

library(pscl)
library(MASS)
card=read.csv("card.csv")
card[,-c(2,3,16:21)]<-scale(card[,-c(2,3,16:21)])#标准化
library(ggplot2)
tx<-ggplot(data=card,mapping=aes(x=yqcs))+geom_bar(stat='count',fill ='pink',colour = 'grey')

fit1<-glm(yqcs~.,family="poisson",data=card,control=list(maxit=1000))
summary(fit1)

fit2=lm(yqcs~.,data=card)
summary(fit2)
#library(arm) 
#bayesglm(yqcs~.,family="poisson",data=card)

#负二项回归
fit3<-glm(yqcs~., family = negative.binomial(theta=1),data=card) 

#过度离势检验
mean(card$yqcs) #0.8342
var(card$yqcs)  #2.415858
library(qcc)
qcc.overdispersion.test(card$yqcs,type="poisson")

#ZIP回归
library(pscl)
fit4<-zeroinfl(yqcs~.|.,card)
fit5<-zeroinfl(yqcs~.|.,card,EM=TRUE)
summary(fit4)


#zip结果有两部分：计数模型--泊松回归；零膨胀模型
#检验两个模型的参数都具有统计意义，即含有预测变量的模型优于null
#模型，可以用lmtest()检验(似然比检验)。
library(zoo)
library(lmtest)
lrtest(fit4)#卡方值为6357.6，自由度为40，对应的P值为2.2e-16<0.0001
#说明模型具有统计显著性

#但是以上的检验没有给出与传统的泊松回归模型相比，零膨胀模型是否有改进
#检验zip优于poisson
vuong(fit1,fit4)

fit0<-hurdle(yqcs~.|., data=card,dist="poisson")
summary(fit0)
lglk0<-fit0$loglik#df=42

#####分割训练集和测试集 ,用NMSE比较模型
##(1)留出法：按照Y进行分层抽样，可以保证训练集和测试集中Y的比例是一致的
library(caret)
set.seed(123)
table(card$yqcs)
table(card$yqcs)/nrow(card)#应变量各个水平的比例
index=createFolds(card$yqcs,k=5,list=FALSE)
mse0<-rep(0,5)
mse<-rep(0,5)

#泊松回归
for(i in 1:5){
  train<-card[-which(index==i),]
  test<-card[which(index==i),]
  a=glm(yqcs~.,family="poisson",data=train)
  y0=predict(a,train)#对训练集预测
  y1=predict(a,test)#对测试集预测
  mse0[i]=mean((train[,21]-y0)^2)/mean((train[,21]-mean(train[,21]))^2)
  mse[i]=mean((test[,21]-y1)^2)/mean((test[,21]-mean(test[,21]))^2)
}
nmse01<-mean(mse0)#2.712316
nmse1<-mean(mse)#2.745967

#负二项分布
for(i in 1:5){
  train<-card[-which(index==i),]
  test<-card[which(index==i),]
  a<-glm(yqcs~., family = negative.binomial(theta=1),data=train)
  y0=predict(a,train)#对训练集预测
  y1=predict(a,test)#对测试集预测
  mse0[i]=mean((train[,21]-y0)^2)/mean((train[,21]-mean(train[,21]))^2)
  mse[i]=mean((test[,21]-y1)^2)/mean((test[,21]-mean(test[,21]))^2)
}
nmse02<-mean(mse0)#2.172339
nmse2<-mean(mse)#2.193163

#线性回归
for(i in 1:5){
  train<-card[-which(index==i),]
  test<-card[which(index==i),]
  a=lm(yqcs~.,data=train)
  y0=predict(a,train)#对训练集预测
  y1=predict(a,test)#对测试集预测
  mse0[i]=mean((train[,21]-y0)^2)/mean((train[,21]-mean(train[,21]))^2)
  mse[i]=mean((test[,21]-y1)^2)/mean((test[,21]-mean(test[,21]))^2)
}
nmse03<-mean(mse0)#0.9061491
nmse3<-mean(mse)#0.9093126
#a<-glm(yqcs~., family = negative.binomial(theta=1),data=card)


#ZIP回归
for(i in 1:5){
  train<-card[-which(index==i),]
  test<-card[which(index==i),]
  a<-zeroinfl(yqcs~.|.,train)
  y0=predict(a,train)#对训练集预测
  y1=predict(a,test)#对测试集预测
  mse0[i]=mean((train[,21]-y0)^2)/mean((train[,21]-mean(train[,21]))^2)
  mse[i]=mean((test[,21]-y1)^2)/mean((test[,21]-mean(test[,21]))^2)
}
nmse05<-mean(mse0)#0.8282108
nmse5<-mean(mse)#0.8308125

#ZINB回归
for(i in 1:5){
  train<-card[-which(index==i),]
  test<-card[which(index==i),]
  a<-zeroinfl(yqcs~.|.,train,dist = "negbin")
  y0=predict(a,train)#对训练集预测
  y1=predict(a,test)#对测试集预测
  mse0[i]=mean((train[,21]-y0)^2)/mean((train[,21]-mean(train[,21]))^2)
  mse[i]=mean((test[,21]-y1)^2)/mean((test[,21]-mean(test[,21]))^2)
}
nmse05<-mean(mse0)#0.8354855
nmse5<-mean(mse)#0.836005