R语言广义线性模型中竖线含义 r语言广义线性混合模型

下面在此分享一下一次课程作业的答题思路及个人答题结果。如有错误欢迎指正。

*本文参考了很多文章。如有雷同纯属纯属巧合。

〇、题目

数据背景：数据1为用户申请手机分期的数据，共2000样本。数据2为某在线网贷平台数据，大家可以选择1个数据展开行业和意义写个报告。

数据集链接：R语言-基于广义线性模型的评分卡模型（数据集）资源

目的：建立基于广义线性模型的评分卡并评估其效果。大家也可以对比其他数据挖掘方法，并对这些方法同基于原始数据的逻辑回归，以及基于WOE的逻辑回归进行比较。

内容要求：
1. 介绍背景和想要解决的问题。
2. 数据介绍，Y怎么定义，X都有哪些，为什么选这些X。
3. 方法介绍。
4. 方法评估。
5. 总结讨论。

一、案例说明

1.1案例背景

近年来，随着科学技术的发展以及人们物质水平生活的升高，手机在我们生活中扮演着越来越重要的角色，因此，出门在外手机就成为了一项必需品。然而手机价格也水涨船高，于是越来越多人在买手机时，选择分期付款。目前信用卡分期付款最常见的期限有3期、6期和12期，期限越长，手续费越高。用户在借贷时，应尽量选择大平台，规避过高的利息以及其他风险。同时，也需及时还款，避免逾期的情况发生。逾期不仅会导致后续无法借款，还会影响个人征信，更严重可能会面临刑事问题。

1.2案例数据

本分析主要采用数据为2018年1月用于申请手机分期的放款客户的数据，共2001个样本。base_info表中有效指标共44个，以个人信息数据为主。

1.2.1好坏样本定义

坏样本即为逾期样本，好样本即为为逾期样本。

1.2.2观测点选择

基于样本数据，我们选择2018年9月为观测点。

1.2.3观察窗口的选择

由于我们要制作的是申请评分卡，并未获得信用局数据，故没有观察窗口的概念。

1.2.4表现窗口的选择

考虑到最短产品期数为6个月，故我们可以初步选择8个月为表现窗口。

1.3研究目的

基于此数据，我们需要结合各放款客户的相关特征，建立手机分期付款的申请评分卡模型，评估其还款能力。再据此判断当前用户出现违约的可能性，对于分数较高但未还款客户及时催收其还款，减少损失。

二、模型准备

2.0数据类型分析

通过对表的初步分析，两表数据集中的数据类型有：文本型、数值型以及逻辑型。考虑到目标变量为逻辑型变量，后续应采用逻辑回归的方法进行回归分析。 

2.1数据处理

2.1.1缺失值处理

对于base_info中数据，经过初步筛选，我们发现了在户籍类型中仅有6项不为空值，其余格中存在大量空值。因此该变量对研究意义不大，故删去户籍类型整列。

2.1.2异常值处理

经缺失值处理后，base_info中还有43个指标2001个样本。下面我们对10个数值类列变量中的值进行异常值检验。通过绘制箱线图观察，得到如下结果：

如图，我们可以发现在INCOME（个人收入）及LOAN_AMT（贷款金额）指标存在较大偏差。经检验，这两类数据并不服从正态分布，但是考虑到个人收入与贷款金额分布在人群中本身就具有随机性，故不做处理。

2.1.3数据划分

数据集进行划分，选择其中70%为训练集、30%为测试集，用训练集的数据进行模拟，测试集的数据进行验证。

2.2目标变量

在本案例中，我们选择STATUS列即截止至2018年9月逾期30天以上客户作为目标变量。0代表未逾期，1代表逾期。

2.3自变量选取

2.3.1相关性分析

由此除目标变量还有42列，这其余42个变量之间的相关关系，使用Spearman相关系数去表示相关关系的强弱情况并生成热力图如下：

根据热力图我们发现，有些变量具有强相关性，在后续回归时可以采用逐步回归的方法，减少变量间相互影响的作用。

【注：此处有很多变量为分类变量，按理说不应该按连续型变量进行相关系数分析，可以只做连续型变量的热力图】

2.3.2计算WOE及IV值

由于后续我们需要制作评分卡，故需要选择从原始数据生成WOE值，并用WOE值代替原始变量进行回归。为合理分段计算WOE值，在分箱时，需要先将各42列进行细分，再将相邻的且P(Y=1|X)相似的细分区合并在一起，避免维数过多引发“维数灾”。粗分后，计算他们的WOE值以及IV值，部分结果如下：

# install.packages('scorecard')
library(scorecard)

data <- read.csv('ibase_info.csv',header = TRUE, stringsAsFactors = F)
dt_s <- var_filter(data,y = 'STATUS')
set.seed(11)
bins <- woebin(dt_s, y="STATUS")

指标	类别	WOE	IV	IV_TOTAL
申请时间	[-Inf,20180103)	-0.13057	0.001024398	0.126118
	[20180103,20180119)	0.328706	0.056977843
	[20180119, Inf)	-0.41159	0.068115544
偿债比	[-Inf,1.2)	-1.24071	0.36215936	0.587857
	[1.2,1.26)	0.296959	0.009537048
	[1.26,1.34)	-0.70341	0.035896696
	[1.34, Inf)	0.579977	0.180264282
身份证命中同盾欺诈低级灰名单	[-Inf,2)	-0.12024	0.012194246	0.080823
	[2, Inf)	0.676708	0.068628769
…	…	…	…	…

其中，WOE越大，代表着该数据区间中因变量为1的比例越高，即客户违约的可能性越搞。变量的IV值代表着该变量对因变量的区分度，IV越大，该变量越有价值。但是由于直接凭IV值大小选取，可能避免不了多重共线性的影响，在这里我们先将IV值最小的10个变量进行剔除，后续再选择用逐步回归的方法选择变量。

三、模型建立

3.1变量选择

考虑到变量过多，可以采取逐步回归的方法，自动从可供选择的变量中选取最重要的几个变量，每次逐个引入自变量时，要保证其偏回归平方和经检验后是显著的。且每引入一个新变量后，要对旧的自变量逐个检验，剔除偏回归平方和不显著的自变量。由此建立逻辑回归分析的解释（预测）模型。

由逻辑回归的基本原理，我们将客户违约(Y=1)的概率表示为p，则正常(Y=0)的概率为1-p，因此可用得到：

则客户违约的概率可表示为：

逻辑回归计算模型如下：

其中，变量

是出现在最终模型中的自变量，即为入模指标。由于此时所有变量都用WOE转换进行了转换，将这些 表示为

，其中式中  为第i行第j个变量的WOE，为已知变量； 为二元变量，表示变量i是否取第j个值。

经逐步回归计算，回归分析结果如下：

dt_list <- split_df(dt_s)
train <- dt_list$train
test <- dt_list$test

woebin_plot(bins)

train_woe <-  woebin_ply(train, bins)
test_woe <-  woebin_ply(test, bins)
print(dim(train_woe))
print(dim(test_woe))
m1 <- glm( STATUS ~ ., family = binomial(), data = train_woe)
summary(m1)

# 逐步回归选择变量
m_step <- step(m1, direction="both", trace = FALSE)
m2 <- eval（m_step$call)
summary(m2)

library(car)
vif(m2)

变量	Estimate	p值	VIF
(Intercept)	-2.5034	< 2e-16 ***
申请时间_woe	0.7611	0.01339 *	1.030457
偿债比_woe	0.7404	2.82e-05 ***	1.177630
贷款金额_woe	0.5539	0.00623 **	1.229454
申请日年龄_woe	0.8282	0.02185 *	1.233900
手机号码使用年限_woe	0.8018	0.00472 **	1.118219
婚姻状况_woe	-0.7749	0.07464 .	1.361892
性别_woe	0.8250	0.01201 *	1.058399
单位性质_woe	1.4964	0.01709 *	1.033176
申请时间（秒）_woe	1.1525	0.00247 **	1.023388
7天内申请人在多个平台申请借款_woe	0.4213	0.04885 *	1.516915
12个月内申请人在多个平台申请借款_woe	-1.1472	0.05362 .	5.790108
18个月内申请人在多个平台申请借款_woe	1.1314	0.04162 *	5.238559
身份证命中同盾欺诈低级灰名单_woe	0.6715	0.04634 *	1.079546
注：***、**、*、.分别代表0.1%、1%、5%、10%的显著性水平。

由上表可以看出，所有指标均通过了显著性检验，解释能力良好。各指标VIF均较小，可以判断模型中多重共线性不显著。

采用“Boruta”法，获取自变量中对违约状态影响最显著的指标，并用箱图表示变量的重要性，如下图：

data <- read.csv('base_info.csv',header = TRUE, stringsAsFactors = F)
data <- as.data.frame(data)

# 定量指标筛选
# "Boruta"法
# install.packages('Boruta')
library(Boruta)
boruta_output<-Boruta(STATUS~.,data = na.omit(data),
                      doTrace=3)

boruta_signif<-names(boruta_output$finalDecision[
  boruta_output$finalDecision %in%c("Confirmed","Tentative",'Rejected')])
#获取自变量中确定的和实验性的指标
print(boruta_signif)
#Levels: Tentative Confirmed Rejected
#Confirmed坚定的；Tentative踌躇的；Rejected拒绝的
plot(boruta_output,cex.axis=.7,las=2,xlab="",main="Variable Importance")
#绘制变量显著性表示的箱图

与逐步回归选择的变量进行对比，具有较显著解释能力的解释变量基本都被选择为了入模指标，说明变量选择较为合理。

3.2建立模型

评分卡设定的分值刻度可以通过将分值表示为比率对数的线性表达式来定义，即可表示为下式：

其中，A,B为常数。式中的负号可以使得违约概率越低，得分越高。通常情况下，这是分值的理想变动方向，即高分值代表低风险，低分值代表高风险。

A、B的值可以通过两个已知或假设的分值带入计算得到。通常情况下，需要假定以下两个条件：

（1）给某个特定的比率设定特定的预期分值；

（2）确定比率翻番的分数（PDO）

根据以上的分析，我们首先假设比率为x的特定点的分值为P；则比率为2x的点的分值应该为P+PDO。代入式中，可以得到如下两个等式：

解之有：

评分卡刻度参数A和B确定以后，就可以计算比率和违约概率，以及对应的分值。通常将常数A称为补偿，常数B称为刻度。

评分卡的分值可表达为：

其中，模型参数 ， ，…， 可以用上面的建模参数进行拟合。

变量

是出现在最终模型中的自变量，即为入模指标。由于此时所有变量都用WOE转换进行了转换，将这些 表示为

，其中式中  为第i行第j个变量的WOE，为已知变量； 为二元变量，表示变量i是否取第j个值。

• 综上，带入拟合后的 ，最终评分卡公式为：

设基础分值等于

，变量 的第j行的分值取决于以下三个数值：

（1）刻度因子B；

（2）逻辑回归方程的参数 ；

（3）第i行的WOE值 。

给定比率x，假定基础分值为568，解出A、B的值并将其转化为标准的信用风险评分卡模型，如下表所示：

# 生成评分卡
card <- scorecard(bins, m2,)
card

指标	类别	分数
基础分值	-	568
申请时间	[-Inf,20180103)	7
	[20180103,20180119)	-18
	[20180119, Inf)	23
偿债比	[-Inf,1.2)	66
	[1.2,1.26)	-16
	[1.26,1.34)	38
	[1.34, Inf)	-31
贷款金额	[-Inf,1600)	35
	[1600,2800)	15
	[2800,3100)	-7
	[3100, Inf)	-30
申请日年龄	[-Inf,19)	6
	[19,20)	-41
	[20,22)	1
	[22,24)	-32
	[24,41)	10
	[41,44)	-4
	[44,48)	35
	48岁以上	-13
手机号码使用年限	[-Inf,10)	-37
	[10,15)	-1
	[15,25)	-24
	[25,85)	15
	[85, Inf)	32
婚姻状况	[-Inf,7002)	17
	[7002,7003)	-16
	[7003, Inf)	6
性别	19001	-13
	19002	31
单位性质	[-Inf,8005)	-57
	[8005,8007)	-1
	[8007,8013)	40
	[8013, Inf)	7
申请时间.精确到秒.	[-Inf,12)	32
	[12,15)	-6
	[15,16)	37
	[16,17)	-15
	[17,18)	-40
	[18, Inf)	7
X7天内申请人在多个平台申请借款	[-Inf,10)	14
	[10,15)	-7
	[15,20)	-33
	[20, Inf)	-26
X12个月内申请人在多个平台申请借款	-1	-58
	0	24
X18个月内申请人在多个平台申请借款	-1	65
	0	-23
身份证命中同盾欺诈低级灰名单	-1	6
	2	-33

*注：得分越高，说明客户信用越好，越不容易违约。得分越低，说明客户信用值不佳，越容易违约。

三、方法评估

3.1模型检验

3.1.1稳定性检验

对模型稳定性的检验，我们可以选择用PSI指标来衡量模型的预测值与实际值偏差大小，进而对模型进行评估。生成的PSI图如下：

train_score <- scorecard_ply(train, card, print_step = 0)
# 验证集评分
test_score <- scorecard_ply(test, card, print_step = 0)

print(train_score)
print(test_score)

# 模型的稳定性度量
# psi
psi_result <- perf_psi(
  score = list(train = train_score, test = test_score),
  label = list(train = train$STATUS, test = test$STATUS)
)

psi_result

由图知，PSI值为0.0132，远小于0.1，说明模型的稳定性很高。

3.1.2模型精确度检验

这里采用KS检验以及AUC检验，并绘制ROC图像。

# 模型性能验证 ks和roc
# 预测的概率
train_pred <- predict(m2, train_woe, type = 'response')
test_pred <- predict(m2, test_woe, type = 'response')

par(mfrow = c(2,1))
train_perf <- perf_eva(train_pred,train$STATUS, title = 'train',show_plot = c("ks", "roc") )

在训练集中，KS=0.5207>0.5,AUC=0.8015>0.75，说明模型有较好的区分度以及可信度，可以有效识别用户的“好坏”。

3.2与其他方法对比的拟合效果

将WOE代替原始变量的逻辑回归与基于原始数据的逻辑回归的进行比较，用测试集进行验证，两者的ROC曲线图如下：

test_perf <- perf_eva(test_pred,test$STATUS,title = 'test_woe',show_plot = c('roc'))

m <- glm( STATUS ~ ., family = binomial(), data = train)
m_step1 <- step(m, direction="both", trace = FALSE)
m22 <- eval（m_step1$call)
test_pred1 <- predict(m22, test, type = 'response')
test_perf1 <- perf_eva(test_pred1,test$STATUS,title = 'test',show_plot = c( "roc"))

可以发现，二者AUC值均大于0.5,说明拟合效果均较为良好。但用WOE拟合的AUC为0.7798，略大于用原始数据拟合的AUC，说明前者拟合效果更好。

3.3分数分布

生成建模样本和验证样本的评分分布如下：

#直方图
par(mfrow = c(1,1))
hist(train_score$score, col = rgb(1,0,0,0.2),freq = F)
lines(density(train_score$score), col = "red")

hist(test_score$score, col = rgb(0,1,0,0.2), freq = F, add = T)
lines(density(test_score$score),col = "green")

由图知，两者得分分布基本保持一致，说明模型建立较为良好。

3.4模型不足

①本文仅进行了一次样本抽样，在实操中，应多次抽样交叉检验，提升模型准确度。

②没有考虑hive表格中有效指标，具有局限性。

③数据分箱时只进行了一次粗分，应多次尝试选出粗分之后比细分时IV值下降最小的分箱方式。

四、研究总结

通过制作申请评分卡，我们可以根据收集上来的信息对用户“好坏”进行合理定义。在授信初期能准确识别借款人的信用水平，将显著减少损失，并减小后续存量用户管理的压力。在得到用户得分后，后续可以根据其他相关样本为其设置一个阈值分，若超过该得分，才可以申请借贷。

附  录

全部代码如下：

# install.packages('scorecard')
library(scorecard)

data <- read.csv('ibase_info.csv',header = TRUE, stringsAsFactors = F)


dt_s <- var_filter(data,y = 'STATUS')
set.seed(11)
bins <- woebin(dt_s, y="STATUS")
dt_list <- split_df(dt_s)
train <- dt_list$train
test <- dt_list$test


woebin_plot(bins)

train_woe <-  woebin_ply(train, bins)
test_woe <-  woebin_ply(test, bins)
print(dim(train_woe))
print(dim(test_woe))
m1 <- glm( STATUS ~ ., family = binomial(), data = train_woe)
summary(m1)

# 逐步回归选择变量
m_step <- step(m1, direction="both", trace = FALSE)
m2 <- eval（m_step$call)
summary(m2)

library(car)
vif(m2)

#===============================================================================
# 原始数据
m <- glm( STATUS ~ ., family = binomial(), data = train)
m_step1 <- step(m, direction="both", trace = FALSE)
m22 <- eval（m_step1$call)
test_pred1 <- predict(m22, test, type = 'response')
test_perf1 <- perf_eva(test_pred1,test$STATUS,title = 'test',show_plot = c( "roc"))

#===============================================================================

# 模型性能验证 ks和roc
# 预测的概率
train_pred <- predict(m2, train_woe, type = 'response')
test_pred <- predict(m2, test_woe, type = 'response')

par(mfrow = c(2,1))
train_perf <- perf_eva(train_pred,train$STATUS, title = 'train',show_plot = c("ks", "roc") )
test_perf <- perf_eva(test_pred,test$STATUS,title = 'test_woe',show_plot = c('roc'))


# 生成评分卡
card <- scorecard(bins, m2,)
card

# write.csv(card$身份证命中同盾欺诈低级灰名单, file = 'G:\\R语言用\\评分卡.csv')


train_score <- scorecard_ply(train, card, print_step = 0)
# 验证集评分
test_score <- scorecard_ply(test, card, print_step = 0)

print(train_score)
print(test_score)

# 模型的稳定性度量
# psi
psi_result <- perf_psi(
  score = list(train = train_score, test = test_score),
  label = list(train = train$STATUS, test = test$STATUS)
)
psi_result

#直方图
par(mfrow = c(1,1))
hist(train_score$score, col = rgb(1,0,0,0.2),freq = F)
lines(density(train_score$score), col = "red")

hist(test_score$score, col = rgb(0,1,0,0.2), freq = F, add = T)
lines(density(test_score$score),col = "green")

逐步回归所选变量的所有WOE、IV值如下：

指标	类别	woe	iv	iv_total
申请时间	[-Inf,20180103)	-0.13057	0.001024398	0.126118
	[20180103,20180119)	0.328706	0.056977843	0.126118
	[20180119, Inf)	-0.41159	0.068115544	0.126118
偿债比	[-Inf,1.2)	-1.24071	0.36215936	0.587857
	[1.2,1.26)	0.296959	0.009537048	0.587857
	[1.26,1.34)	-0.70341	0.035896696	0.587857
	[1.34, Inf)	0.579977	0.180264282	0.587857
贷款金额	[-Inf,1600)	-0.88088	0.13224792	0.346749
	[1600,2800)	-0.38099	0.041961336	0.346749
	[2800,3100)	0.168476	0.005877447	0.346749
	[3100, Inf)	0.739461	0.166662324	0.346749
申请日年龄	[-Inf,19)	-0.10264	0.00055474	0.097632
	[19,20)	0.693691	0.034703366	0.097632
	[20,22)	-0.0226	4.20E-05	0.097632
	[22,24)	0.533348	0.027558238	0.097632
	[24,41)	-0.16579	0.013826495	0.097632
	[41,44)	0.061662	0.000253506	0.097632
	[44,48)	-0.58563	0.017483859	0.097632
	48岁以上	0.218267	0.00320998	0.097632
手机号码使用年限	[-Inf,10)	0.638861	0.050896749	0.151473
	[10,15)	0.013864	2.74E-05	0.151473
	[15,25)	0.41989	0.036892913	0.151473
	[25,85)	-0.26637	0.027947456	0.151473
	[85, Inf)	-0.55284	0.035708752	0.151473
婚姻状况	[-Inf,7002)	0.300593	0.038532343	0.07819
	[7002,7003)	-0.28291	0.038738053	0.07819
	[7003, Inf)	0.105214	0.000919408	0.07819
性别	19001	0.216754	0.032999569	0.113423
	19002	-0.52826	0.080423885	0.113423
单位性质	[-Inf,8005)	0.530503	0.017745265	0.027923
	[8005,8007)	0.006657	2.95E-05	0.027923
	[8007,8013)	-0.37052	0.009347287	0.027923
	[8013, Inf)	-0.06329	0.00080127	0.027923
申请时间.精确到秒.	[-Inf,12)	-0.38639	0.017085748	0.074491
	[12,15)	0.067055	0.001359035	0.074491
	[15,16)	-0.44748	0.020519185	0.074491
	[16,17)	0.178526	0.004378311	0.074491
	[17,18)	0.480631	0.029804479	0.074491
	[18, Inf)	-0.08034	0.00134401	0.074491
X7天内申请人在多个平台申请借款	[-Inf,10)	-0.45605	0.116937579	0.33487
	[10,15)	0.232679	0.01026421	0.33487
	[15,20)	1.094336	0.122849061	0.33487
	[20, Inf)	0.847841	0.084818659	0.33487
X12个月内申请人在多个平台申请借款	[-Inf,0)	-0.70476	0.144243778	0.204724
	[0, Inf)	0.295498	0.060480015	0.204724
X18个月内申请人在多个平台申请借款	[-Inf,0)	-0.79796	0.16441761	0.223111
	[0, Inf)	0.284853	0.058692926	0.223111
身份证命中同盾欺诈低级灰名单	[-Inf,2)	-0.12024	0.012194246	0.080823
	[2, Inf)	0.676708	0.068628769	0.080823

其各变量分箱后的坏账率如下，基本满足单调性及U型趋势。