基于R语言的用户贷款风险预测

作者介绍：皮吉斯，热爱R语言，知乎号：皮吉斯
已授权转载

本次的分析数据来自kaggle数据竞赛平台的 “give me some credit” 竞赛项目。任务是提高模型精度AUC。
本次分析用到了多种算法，分别有：逻辑回归，cart决策树，神经网络，xgboost，随机森林。通过多种模型相互对比，最终根据auc选出最好的模型。

分析步骤：

1. 导入数据
2. 数据清洗及准备
3. 模型建立
4. 模型评估
5. 多模型对比

导入数据

cs_training <- read.csv("cs_training.csv")
names(cs_training) 
a <- cs_training 
a1 <- cs_training 
colnames(a)<-c('id','response','x1','age','x2','debtratio','monthlyincome','x3','x4','x5','x6','x7')

导入数据，并对列名重命名，方便分析。
通过summary了解数据的整体情况，可以看到monthliincome和x7变量有缺失值

数据清洗

#age变量
age<-a$age
var_age<-c(var="age",
mean=mean(age,na.rm=TRUE),#na.rm=TRUE去除NA的影响
median=median(age,na.rm=TRUE),
quantile(age,c(0,0.01,0.1,0.25,0.5,0.75,0.9,0.99,1),na.rm=TRUE),  
max=max(age,na.rm=TRUE),
missing=sum(is.na(age)) )
View(var_age)

var	mean	median	0	0.01	0.1	0.25	0.5	0.75	0.9	0.99	1	max	missing
age	52.2952066666667	52	0	24	33	41	52	63	72	87	109	109	0

这样可以看到各个变量的数据分布情况

# 查看异常值
boxplot(age~response,data = a,horizontal = T ,frame = F ,        col = "lightgray",main = "Distribution")

上图可以看到，数据存在异常值。
处理异常值：通常采用盖帽法，即用数据分布在1%的数据覆盖在1%以下的数据，用在99%的数据覆盖99%以上的数据。

block<-function(x,lower=T,upper=T){
if(lower){
q1<-quantile(x,0.01)
x[x<=q1]<-q1
}  
if(upper){
q99<-quantile(x,0.99)
x[x>q99]<-q99
}
return(x)
}
boxplot(age~response,data = a,horizontal = T ,frame = F ,        col = "lightgray",main = "Distribution")

经过处理，异常值大量减少

# x1
xa1<-a$x1
var_x1<-c(var='xa1',
mean=mean(xa1,na.rm = T),
median=median(xa1,na.rm = T),
quantile(xa1,c(0,0.01,0.1,0.5,0.75,0.9,0.99,1),na.rm = T), max=max(xa1,na.rm = T), 
miss=sum(is.na(xa1)) )
boxplot(x1~response,data=a,horizontal=T, frame=F, col="lightgray",main="Distribution-x1") 
#对X1变量进行处理 
a$x1<-block(a$x1)

# x2
a$x2 
summary(a$x2) 
xa2<-a$x2 
var_x2<-c(var='xa2',
mean=mean(xa2,na.rm = T), 
median=median(xa2,na.rm = T), 
quantile(xa2,c(0,0.01,0.1,0.5,0.75,0.9,0.99,1),na.rm = T),  max=max(xa2,na.rm = T), 
miss=sum(is.na(xa2)) ) 
boxplot(x2~response,data=a,horizontal=T, frame=F,  col="lightgray",main="Distribution") 
# 对x2变量进行处理 
a$x2<-block(a$x2) 
a$x1<-round(a$x1,2)
# debtratio
summary(a$debtratio) 
ratio<-a$debtratio 
var_ratio<-c(var='debtratio', 
mean=mean(ratio,na.rm = T), 
median=median(ratio,na.rm = T), 
quantile(ratio,c(0,0.01,0.1,0.5,0.75,0.9,0.99,1),na.rm = T),  max=max(ratio,na.rm = T), 
miss=sum(is.na(ratio)) ) 
hist(a$debtratio)
# 因为debtratio是百分比，对异常值的处理要结合实际
a$debtratio<-ifelse(a$debtratio>1,1,a$debtratio)
# monthlyincome
income<-a$monthlyincome 
var_income<-c(var='income',
mean=mean(income,na.rm = T), 
median=median(income,na.rm = T), 
quantile(income,c(0,0.01,0.1,0.5,0.75,0.9,0.99,1),na.rm = T),
max=max(income,na.rm = T), 
miss=sum(is.na(income)) ) 
hist(income) 
boxplot(monthlyincome~response,data=a,horizontal=T, frame=F,  col="lightgray",main="Distribution-x2")
# 对缺失值处理
a$monthlyincome<-ifelse(is.na(a$monthlyincome)==T,6670.2,a$monthlyincome)
# 对异常值处理
a$monthlyincome<-block(a$monthlyincome)

# x3
summary(a$x3) 
x3<-a$x3 
var_x3<-c(var='x3',mean=mean(x3,na.rm = T),  
median=median(x3,na.rm = T),          
quantile(x3,c(0,0.01,0.1,0.5,0.75,0.9,0.99,1),na.rm = T),   max=max(x3,na.rm = T),  
miss=sum(is.na(x3)) ) 
names(a) 
# 对x3进行处理 
a$x3<-block(a$x3)

# x4-延迟90天的次数
summary(a$x4) 
x4<-a$x4 
var_x4<-c(var='x4',mean=mean(x4,na.rm = T),
median=median(x4,na.rm = T),          
quantile(x4,c(0,0.01,0.1,0.5,0.75,0.9,0.99,1),na.rm = T),  max=max(x4,na.rm = T),  
miss=sum(is.na(x4)) ) 
# 对x4进行处理 
a$x4<-block(a$x4)

# x5-贷款额度
summary(a$x5) 
x5<-a$x5 
var_x5<-c(var='x5',mean=mean(x5,na.rm = T),   
median=median(x5,na.rm = T),          
quantile(x5,c(0,0.01,0.1,0.5,0.75,0.9,0.99,1),na.rm = T),  max=max(x5,na.rm = T),  
miss=sum(is.na(x5)) ) 
# 对x5异常值进行采用盖帽法处理 
a$x5<-block(a$x5) 
boxplot(a$x5) 
boxplot(x5~response , data = a , horizontal = T, main="Distribution-x5")

# x6
summary(a$x6) 
x6<-a$x6 
var_x6<-c(var='x6',mean=mean(x6,na.rm = T), 
median=median(x6,na.rm = T),  
quantile(x6,c(0,0.01,0.1,0.5,0.75,0.9,0.99,1),na.rm = T),  max=max(x6,na.rm = T),  
miss=sum(is.na(x6)) ) 
# 对x6进行盖帽法 
a$x6<-block(a$x6)

# x7
summary(a$x7) 
x7<-a$x7 
var_x7<-c(var='x7',mean=mean(x7,na.rm = T),  
median=median(x7,na.rm = T),  
quantile(x7,c(0,0.01,0.1,0.5,0.75,0.9,0.99,1),na.rm = T),  max=max(x7,na.rm = T),  
miss=sum(is.na(x7)) ) 
hist(a$x7,freq=F)
# 对缺失值处理
a$x7<-ifelse(is.na(a$x7)==T,0.76,a$x7)
# 对异常值处理
a$x7<-block(a$x7) 
summary(a)

#response变量
# 为了方便理解，将1作为违约，0表示不违约
a$response<-as.numeric(!as.logical(a$response))

建模

数据分组

table(a$response)

0	1
10026	139974

数据正负比例不平衡，我们要对数据进行smote处理，smote算法的思想是合成新的少数类样本，合成的策略是对每个少数类样本a，从它的最近邻中随机选一个样本b，然后在a、b之间的连线上随机选一点作为新合成的少数类样本。

library(lattice) 
library(grid) 
library(DMwR) 
a$response <- factor(a$response) 
a$response <- as.factor(ifelse(a$response == 0 , "yes" , "no")) 
newData <- SMOTE(response ~ ., a) 
table(newData$response) 
newData$response <- ifelse(newData$response == "yes" , 0 , 1) 
newData$response <- as.numeric(newData$response) 
str(newData)

# 接下来进行分组
index <- sample(1:nrow(newData) , nrow(newData)*0.7) 
train <- newData[index,] 
test <- newData[-index,] 
table(train$response) 
table(test$response)

table(train$response)

0	1
21213	27914

table(test$response)

0	1
8865	12190

逻辑回归

建立模型

glm_model <- glm(response~. , data = train  , family = binomial()) 
summary(glm_model)
#p值全部显著

使用逐步法剔除变量

step(glm_model , direction = "both")

VIF多重共线性检验

library(carData) 
library(car) 
vif(glm_model) #一般认为VIF值大于2的话，表明变量间存在共线性。此时没有大于2的值,各个变量间相互独立

预测

train_pred <- predict(glm_model , newdata = train , type = "response") 
test_pred <- predict(glm_model , newdata = test , type = "response")

模型评估（auc=0.842466）

library(gplots)
library(ROCR)
test_prob <- prediction(test_pred,  test$response) 
test_perf <- performance(test_prob , "tpr" , "fpr") 
test_auc <- performance(test_prob ,"auc")
plot(test_perf) 
abline(a=0,b=1)

cart决策树模型

模型建立

library(rpart) 
library(rpart.plot) 
library(RWeka) 
cart_dt<-rpart.control(minsplit = 50,maxdepth = 4,xval=10,cp=0) 
cart_model <- rpart(response~.,  data = train , method = "anova" , control = cart_dt) 
rpart.plot(cart_model,digits=3)

预测

test_pred <- predict(cart_model  , newdata = test)

模型评估（auc=0.3670627）

test_prob <- prediction(test_pred,  test$response) 
test_perf <- performance(test_prob , "tpr" , "fpr") 
test_auc <- performance(test_prob ,"auc")
plot(test_perf) 
abline(a=0,b=1)

nnet包的单隐层BP网络

对数据进行标准化

#对数据进行标准化 
names(train) 
train[,2:11] <- scale(train[,2:11]) 
test[,2:11] <- scale(test[,2:11])

模型建立

library(nnet) 
train_nnet<-nnet(response~., linout =F,size=10, decay=0.0076, maxit=200,data = train) 
# size隐节点数,decay权值衰减率

预测并将概率调整为0和1

test_pred<-predict(train_nnet, test) 
test_pred[test_pred<0.5]=0 
test_pred[test_pred>=0.5]=1 
table(test_pred)

test_pred

0	1
14523	6532

计算准确率

test_pred_rate<-sum(test_pred==test$response)/length(test$response)

模型评估（auc=0.7848754）

library(gplots) 
library(ROCR) 
test_prob <- prediction(test_pred,  test$response) 
test_perf <- performance(test_prob , "tpr" , "fpr") 
test_auc <- performance(test_prob ,"auc")
plot(test_perf) 
abline(a=0,b=1)

模型调参

##输入数据的预测变量必须是二分类的。且所有变量只包含模型输入输出变量。
##调参时如果变量过多size不宜过大。
##构建调参函数network()。 构建调参函数network()
##构建调参函数network()。 

network<-function(formula,data,size,adjust,decay=0,maxit=200,scale=TRUE,samplerate=0.7,seed=1,linout=FALSE,ifplot=TRUE){
library(nnet) 
##规范输出变量为0,1   
yvar<-colnames(data)==(all.vars(formula)[1])  
levels(data[,yvar])<-c(0,1)  ##抽样建立训练集和测试集  
set.seed(seed)  
select<-sample(1:nrow(data),nrow(data)*samplerate)  train=data[select,]  
test=data[-select,]  ##根据给定判断进行标准化  
if(scale==T){    
xvar<-colnames(data)!=(all.vars(formula)[1])    train[,xvar]=scale(train[,xvar])    test[,xvar]=scale(test[,xvar])  
}  ##循环使用nnet训练调参  
obj<-eval(parse(text = adjust))  
auc<-data.frame()  
for(i in obj){    
if(adjust=="size"){       
mynnet<-nnet(formula,size=i,linout=linout,decay=decay,                   
maxit=maxit,trace=FALSE,data=train)    
}    
else if(adjust=="decay"){     
mynnet<-nnet(formula,size=size,linout=linout,decay=i,                   
maxit=maxit,trace=FALSE,data=train)} ##调用之前的ROC()得到对应参数的AUC值    
objcolname<-all.vars(formula)[1]    
auc0<-ROC(model=mynnet,train=train,test=test,              
objcolname=objcolname,ifplot=F)    ##输出指定参数不同值对应的数据框    
out<-data.frame(i,auc0)    
auc<-rbind(auc,out)  }    
names(auc)<-c(adjust,"Train_auc","Test_auc")  
if(ifplot==T){    
library(plotrix)    
twoord.plot(auc[,1],auc$Train_auc,auc[,1],auc$Test_auc,lcol=4,rcol=2,xlab=adjust,ylab="Train_auc",  rylab="Test_auc",type=c("l","b"),lab=c(15,5,10))  }   
return(auc) } 
auc<-network(response~.,data=train,size=4:16,adjust="size",             
decay=0.0001,maxit=200,scale=T) #发现当隐藏节点数为6的时候，auc最高 
auc<-network(response~.,data=train,size=11,adjust="decay",             
decay=c(0,seq(0.0001,0.01,0.0003)),maxit=200) 
#发现decay等于0.0025的时候， auc最高

建模2.0

#调参后的模型 #size =11 , decay = 0.0013 , maxit = 200 
library(nnet) 
train_nnet <- nnet(response~. , linout = F , size = 4 , decay = 0.0031 , maxit = 200 , data = train)

模型评估（auc = 0.7923195）

# 预测为0，1 
test_pred<-predict(train_nnet, test) test_pred[test_pred<0.5]=0 
test_pred[test_pred>=0.5]=1 
table(test_pred) #计算准确率 
test_pred_rate<-sum(test_pred==test$response)/length(test$response) # 模型评估 
library(gplots) 
library(ROCR) 
test_prob <- prediction(test_pred,  test$response) 
test_perf <- performance(test_prob , "tpr" , "fpr") 
test_auc <- performance(test_prob ,"auc")
plot(test_perf)
abline(a=0,b=1)

xgboost算法

模型建立

library(xgboost) 
xgb_model <- xgboost(data = data.matrix(train[,-1]) ,                              
label=data.matrix(train$response) ,                                        
max_depth = 6, eta=  0.3 , nthread = 2 , nrounds = 15 ,                             
objective = "binary:logistic" , seed = 123)

预测

test_prob <- predict(xgb_model , data.matrix(test[,-1]))

模型评估（auc = 0.860617）

test_pred <- prediction(test_prob , test$response) 
test_perf <- performance(test_pred , "tpr" , "fpr") 
test_auc <- performance(test_pred , "auc")
plot(test_perf)
abline(a=0,b=1)

随机森林

模型建立

library(randomForest) 
random_model <- randomForest(response~. , data = train) 
random_model

模型预测

test_prob <- predict(random_model , test , type= "response")

模型评估（auc = 0.8531173）

test_pred <- prediction(test_prob , test$response) 
test_perf <- performance(test_pred , "tpr" , "fpr") 
test_auc <- performance(test_pred , "auc")
plot(test_perf)
abline(a=0,b=1)

综上，通过多种模型对比可以看到xgboost算法的模型精确度是最高的达到0.860617。