0.原起 🔗
本新手打算整理一下建模型的基本方法,做点系统性的梳理认识下自身的不足,查漏补缺。不过并不打算一篇笔记写完全部,先初步写个非常基础的做二分类的建模流程,随后慢慢对模型可解释性、特征工程、效果验证等几方面内容单独写写。另外,虽然对 caret、mlr3 早有耳闻,但本新手都没接触过,有很多很优秀的方法和包都听说过但一直没来得及学习,希望未来的日子可以好好学学。
这篇笔记使用的 R 包如下:
- data.table:数据处理。
- DataExplorer:数据探索,主要用来批量绘图,由于此包依赖 ggplot2,数据量大了出图慢。
- echarts4r/DT:由于平时工作用的是 Rstudio Server,服务器上没有装 DataExplorer 这个包,实际上更多用基础的图表包来探索数据。
- scorecard:计算IV值,WOE转换,绘制 ROC曲线、K-S 曲线图。
- xgboost:训练模型。
1.数据集介绍 🔗
这是一个典型的用于做二分类预测的数据集,是在 Kaggle 上随意找的1,地址是https://www.kaggle.com/datasets/ifteshanajnin/carinsuranceclaimprediction-classification。
建模场景是:根据车险保单的车辆信息预测未来6月内是否会出险理赔。训练集有58592条数据,测试集有39063条数据,其中训练集有结果标签,而测试集无。
-
数据介绍
- 没有任何时间相关字段
- 没有缺失值,也没有系统误差造成的异常值
- 仅有车辆本身相关信息,没有任何与投保人自身属性或行为相关信息
- 也没有保单相关信息,如保费、保额、险类等。车险保费定价和有无出险理赔过相关,保费保额也是能反映出一些车辆的信息的。
-
字段介绍
- is_claim:结果标签,投保人未来6月内是否提出索赔
- policy_id:保单投保人唯一标识
- policy_tenure:保单期限
- age_of_car:车龄(按年,标准化处理后)
- age_of_policyholder:投保人年龄(按年,标准化处理后)
- area_cluster:投保人所在区域(重新编码处理后)
- population_density:城市人口密度(投保人城市)
其他字段
- make:汽车制造商的编码
- segment:汽车车型(A/ B1/ B2/ C1/ C2)(C1指手动挡小型汽车,C2指自动挡小型汽车,B1指中型客车,B2指大型货车,A1指大型客车)
- model:汽车的某种类型编码(从M1到M9分别限定了车辆载客座位数、载货总质量)
- fuel_type:汽车使用的燃料类型
- max_torque:汽车产生的最大扭矩
- max_power:汽车产生的最大功率
- engine_type:汽车发动机类型
- airbags:车内安全气囊数量
- is_esc:是否存在电子稳定控制系统(ESC)
- is_adjustable_steering:汽车方向盘是否可调整
- is_tpms:是否存在轮胎压力监测系统(TPMS)
- is_parking_sensors:是否存在停车传感器
- is_parking_camera:是否存在停车摄像头
- rear_brakes_type:汽车制动器类型
- displacement:汽车发送机排量
- cylinder:汽车发送机气缸数
- transmission_type:汽车变速器类型
- gear_box:车内齿轮数
- steering_type:汽车动力转向类型
- turning_radius:车辆转弯所需米数
- length:汽车长度
- width:汽车宽度
- height:汽车高度
- gross_weight:汽车最大载重
- is_front_fog_lights:是否有前雾灯
- is_rear_window_wiper:是否有后窗雨刷器
- is_rear_window_washer:是否有后窗清洗器
- is_rear_window_defogger:是否有后窗除雾器
- is_brake_assist:是否有制动辅助功能
- is_power_door_locks:是否有电动门锁
- is_central_locking:是否能中央控锁
- is_power_steering:是否能动力转向
- is_driver_seat_height_adjustable:是否可调整驾驶座高度
- is_day_night_rear_view_mirror:是否有后视镜
- is_ecw:是否有发动机检查警告(ECW)
- is_speed_alert:是否有速度报警
- ncap_rating:安全评级
保单期限(policy_tenure)这个字段有点奇怪,一般车险的保险保障期限都是一年,而数据集中做过标准化处理的保单期限分布在0.002735至1.396641之间。猜测这个字段的含义可能是一张保单未来剩余的保障期限,不然的话也可能是数据中混合了多种不同的车险险类。
2.探索性数据分析 🔗
在正式开始探索数据之前,需要先了解一下数据的基本概况。比如数据中好坏样本的比例,比例越失衡,建模型难度越高。要是碰到坏样本占比在千分之一或万分之一的情况,可能得先考虑清楚这份数据到底适不适合建立一个二分类的模型,如果不做二分类的话,是不是做异常检测更合适。另外,在建模型之前,确定一个清晰明确的验证目标,以及在上线后始终保持逻辑一致的标签定义,是更重要的事。比如先厘清这份数据的预测目标的标签定义为何是“未来6月内是否提出理赔”,如果只是预测短期出险的话为何是“6月内”而不是“3月内”或“1年内”,这里“理赔”的含义是否包含保险金欺诈。
-
查看训练集标签
0:1 = 54844:3748
,整体理赔比例是6.39%,属于非平衡数据集。 -
一般需要看看按年、按月等不同时间维度的数据量变化趋势,比对一下各年数据差异。由于此数据集中没有时间字段,所以暂略。
-
看看有无拆分子模型的必要。一般越平衡的数据,建模越容易。如果数据拆分后,好坏样本比例差异很大,不妨一试。
-
看看数据缺失情况,若有,可删、可填补。有没有异常值,若有,可删、可修正。
library(data.table)
library(echarts4r)
library(DT)
library(scorecard)
library(xgboost)
library(DataExplorer)
train <- fread('baseline1/train.csv')
test <- fread('baseline1/test.csv')
train$is_claim <- as.numeric(train$is_claim)
train$type <- "train"
test$type <- "test"
test$is_claim <- NA
data <- rbind(train, test)
#最大扭矩、最大功率分别拆成两个字段
data$max_power_1 <-
sapply(strsplit(as.character(data$max_power), 'bhp@'), "[", 1)
data$max_power_2 <-
gsub("rpm", "", gsub("^.+bhp@", "", data$max_power))
data$max_torque_1 <-
sapply(strsplit(data$max_torque, 'Nm@'), "[", 1)
data$max_torque_2 <-
gsub("rpm", "", gsub("^.+Nm@", "", data$max_torque))
# 根据字段含义,部分字段可统一转换为因子型
char <- c(
'area_cluster','make','segment','model','fuel_type','max_torque','max_power','engine_type','is_esc','is_adjustable_steering',
'is_tpms','is_parking_sensors','is_parking_camera','rear_brakes_type','transmission_type','gear_box','steering_type','is_front_fog_lights',
'is_rear_window_wiper','is_rear_window_washer','is_rear_window_defogger','is_brake_assist','is_power_door_locks','is_central_locking',
'is_power_steering','is_driver_seat_height_adjustable','is_day_night_rear_view_mirror','is_ecw','is_speed_alert','ncap_rating','type')
# 根据字段含义,部分字段可统一转换为数值型
num <- c(
'policy_tenure','age_of_car','age_of_policyholder','population_density','airbags','displacement','cylinder','turning_radius',
'length','width','height','gross_weight','max_power_1','max_power_2','max_torque_1','max_torque_2')
# head(data[,..num])
# head(data[,..char])
data <-
data[, (char) := lapply(.SD, as.factor), .SDcols = char][, (num) := lapply(.SD, as.numeric), .SDcols = num]
# sapply(data, function(x) sum(is.na(x))) # 计算缺失值
# summary(data)
2.1.批量查看数据分布 🔗
用 DataExplorer 包来批量查看数据的一些基本情况。一番比较查看后,很容易发现,这份数据的训练集和测试集的数据分布相当一致,绝大多数字段的理赔比例都分布得很均匀,构造特征的空间不大。
- 数值型数据查看直方图或箱线图
- 类别型数据查看柱状图、条形图、饼图
- 比对一下训练集、测试集数据差异
plot_bar(data[, ..char])
plot_histogram(data[, ..num])
# 比较训练集、测试集数据分布有没有差异
# plot_bar(data[type == 'train', ..char])
# plot_bar(data[type == 'test', ..char])
# plot_histogram(data[type == 'train', ..num])
# plot_histogram(data[type == 'test', ..num])
# 比较训练集、测试集在数据取值比例上有无差异
plot_bar(data[, ..char], by = "type")
plot_boxplot(data, by = "type")
# 查看训练集中各字段在标签上的比例有无差异
train$is_claim <- as.factor(train$is_claim)
plot_bar(train, by = "is_claim")
plot_bar(train, by = "is_claim", by_position = "dodge")
去年看到一种叫做对抗验证的思路,用于确认训练集测试集的数据分布是否一致,大意是将数据来源于训练集还是测试集作为标签来训练模型,如果模型的 AUC 值接近0.5,那么训练集测试集的数据分布便算是一致的。
2.2.进一步探索 🔗
如果是分析熟悉的业务数据,可以直接根据数据表现来验证猜想。但是笔者对车辆性能一点都不熟,两眼一摸瞎,干脆先算一下所有特征的 IV 值(Information Value,信息值),看看有没有哪些特征对结果影响大点。
iv <- scorecard::iv(data[type=='train', -c('policy_id', 'type')], y = 'is_claim')
datatable(iv)
只有四个 IV 值超过0.02的字段,其中投保人所在区域(area_cluster)和城市人口密度(population_density)的 IV 值居然完全相等,两个字段的取值可能是一一对应的。还有不少 IV 值小于0.02的字段,其 IV 值也完全一致,那些字段的具体取值大概都是一一对应的。
table(data$population_density, data$area_cluster)
- age_of_car:数值型,车龄(按年,标准化处理后)。怎么车龄为0的客户理赔比例最高呢?
折叠代码块
tmp <- data[type == 'train', by = .(age_of_car), .(
cnt = uniqueN(policy_id),
cnt_0 = uniqueN(policy_id[is_claim == '0']),
cnt_1 = uniqueN(policy_id[is_claim == '1'])
)][, ':='(bili = round(cnt_1 / cnt, 4))]
tmp |>
e_charts(age_of_car) |>
e_bar(cnt_0,
name = '理赔样本数',
stack = 'group',
y_index = 0) |>
e_bar(cnt_1,
name = '未理赔样本数',
stack = 'group',
y_index = 0) |>
e_line(bili, name = '比例', y_index = 1) |>
e_tooltip(trigger = 'axis') |>
e_title(text = 'age_of_car') |>
e_x_axis(max = 1) |>
e_mark_line("比例", data = list(yAxis = 0.0639))
- age_of_policyholder:数值型,投保人年龄(按年,标准化处理后)。随着投保人年龄增长,理赔比例在波动……
折叠代码块
tmp <- data[type == 'train', by = .(age_of_policyholder), .(
cnt = uniqueN(policy_id),
cnt_0 = uniqueN(policy_id[is_claim == '0']),
cnt_1 = uniqueN(policy_id[is_claim == '1'])
)][, ':='(bili = round(cnt_1 / cnt, 4))]
tmp |>
e_charts(age_of_policyholder) |>
e_bar(cnt_0,
name = '理赔样本数',
stack = 'group',
y_index = 0) |>
e_bar(cnt_1,
name = '未理赔样本数',
stack = 'group',
y_index = 0) |>
e_line(bili, name = '比例', y_index = 1) |>
e_tooltip(trigger = 'axis') |>
e_title(text = 'age_of_policyholder') |>
e_x_axis(max = 1) |>
e_mark_line("比例", data = list(yAxis = 0.0639))
- area_cluster:字符型,投保人所在区域(重新编码处理后)。部分区域样本数远超其他。
折叠代码块
tmp <-
data[type == 'train', by = .(area_cluster), .(
cnt = uniqueN(policy_id),
cnt_0 = uniqueN(policy_id[is_claim == '0']),
cnt_1 = uniqueN(policy_id[is_claim == '1'])
)][, ':='(bili = round(cnt_1 / cnt, 4))]
tmp |>
e_charts(area_cluster) |>
e_bar(cnt_0,
name = '理赔样本数',
stack = 'group',
y_index = 0) |>
e_bar(cnt_1,
name = '未理赔样本数',
stack = 'group',
y_index = 0) |>
e_line(bili, name = '比例', y_index = 1) |>
e_tooltip(trigger = 'axis')|>
e_title(text = 'area_cluster')|>
e_mark_line("比例", data = list(yAxis = 0.0639)) |>
e_x_axis(axisLabel = list(interval = 0))
- 维度交叉看看,比如投保人所在区域和投保人年龄或车龄的分布,暂不展开。这里就是看看整体数据中有没有哪一部分的理赔比例高出6.39%。如果能找出这样的一个组成部分,且其背后代表的规律能始终保持一致,那么大概率就是找出来一个很有用的特征。笔者的个人习惯是很喜欢分析数据的组成成分和变化趋势……
3.特征筛选 🔗
在筛选特征之前,先得在原数据基础上做充分的探索,在分析的基础上构造新特征。一般情况下,需要上线应用的模型,在构造新特征有几点需要特别注意:
-
不能泄露标签。当前取用历史数据建模预测未来时,所使用的数据应能还原历史。
-
谨慎地验证特征的稳定性。在探索数据时,应谨慎观察分析得到的数据规律是否能在多数情况下保持一致。建模型实际上就是探寻数据中隐藏的规律,积累大量的偶然性规律后虽不能直接得到必然性结论,但可以发现或者说预测出提高偶然规律的概率。不稳定的特征如果只是偶然看起来管用,实则就是无用。
-
构造的特征太多,后面筛选特征时也麻烦。
-
在训练数据和验证数据中要保持各方面一致。举个奇葩例子,如果某些特征刚好只在保留的时间外验证数据中效果极好,可能让验证效果看起来很好。
3.0.构造特征 🔗
由于这份数据构造特征的空间很有限,加上笔者水平也很有限,本小节写得十分简略。
-
方法一:针对类别型特征,将标签替换为数字后两两相乘,或者直接把各字段中的字符串拼接在一起,看能否得到 IV 值超过0.02的新特征。
-
方法二:针对数值型特征,计算不同维度下或者不同维度交叉组合下的件数(count)、金额(sum),或计算最大值、最小值、平均值、中位数等。
-
方法三:结合对数据的探索和认知来构造新特征。比如仔细看看特征两两组合以后,哪些组合情况下理赔比例高于总体样本,且占比也不低(约大于0.1)的,然后单独构造新特征。依此方法共构建了 f1/f2/f3/f4/f5等5个新特征。
折叠代码块
# 方法三 举个栗子
# tmp <-
# data[type == 'train', by = .(area_cluster, make), .(
# cnt = uniqueN(policy_id),
# cnt_0 = uniqueN(policy_id[is_claim == '0']),
# cnt_1 = uniqueN(policy_id[is_claim == '1'])
# )][, ':='(bili1 = round(cnt_1 / cnt, 4),
# bili2 = round(cnt / 58592, 4))]
#
# datatable(tmp,
# filter = list(
# position = 'top',
# # 可选参数值有"none", "bottom", "top"
# clear = TRUE,
# # 是否展示在筛选框中输入过滤条件后的清除按钮
# plain = FALSE
# )) |>
# formatStyle(columns = c(6), color = styleInterval(c(0.0639), c('black', 'red'))) |>
# formatStyle(columns = c(7), color = styleInterval(c(0.1), c('black', 'green')))
data$f1 <- ifelse(data$age_of_car == 0, 1, 0)
data$f2 <-
ifelse(data$area_cluster %in% c('C2', 'C3', 'C4', 'C8', 'C14', 'C18', 'C19', 'C21', 'C22'), 1, 0)
data$f3 <-
ifelse((
data$area_cluster %in% c('C8', 'C2', 'C14') &
data$segment %in% c('B2', 'C2')) |
(data$area_cluster == 'C3' &
data$segment %in% c('A', 'B2')) |
(data$area_cluster == 'C8' &
data$segment == 'C1') |
(data$area_cluster == 'C12' & data$segment == 'B2'), 1, 0)
data$f4 <-
ifelse((
data$make == '1' &
data$area_cluster %in% c('C8', 'C3', 'C2', 'C14', 'C19', 'C4')) |
(data$make == '3' & data$area_cluster %in% c('C8', 'C2', 'C14')), 1, 0)
data$f5 <-
ifelse((
data$area_cluster %in% c('C8', 'C2', 'C3', 'C14') &
data$fuel_type %in% c('Petrol', 'Diesel')) | (data$area_cluster %in% c('C3', 'C14') &
data$fuel_type == 'CNG'), 1, 0)
3.1.数据分箱 🔗
使用 scorecard 包做数据分箱,就是将特征进行离散化处理,比如将连续型特征处理成离散型;或者将类别较多的特征使其部分类别合并,从而处理成类别更少的。这样做的好处有几点:一是特征更稳定,模型也更稳定;二是离散化的特征更好解释;三是可以提升计算速度,比如说如果手动用笔计算两个矩阵相乘,肯定是更稀疏或者更简单的矩阵算起来更快的。
####
#1.进行WOE转换时,需要先去掉一些id类、文本类、时间类的字段,
#2.最好将类别特别多的特征跟普通特征分开执行,不然的话全部一起执行会很慢
####
df <- data[type == 'train', -c('policy_id','type')]
# WOEBIN
bins = woebin(df,
y = "is_claim",
method = "tree",
print_info = TRUE)
# 将结果从列表转换为数据框,方便后续用 DT 表格查看结果
new <-
unlist(lapply(bins, function(x)
if (is.data.frame(x))
list(x)
else
x), recursive = FALSE)
newbin <- do.call(rbind.data.frame, new)
# 一般仅保留IV值大于等于0.02的特征
datatable(
newbin[total_iv >= 0.01, ],
extensions = 'RowGroup',
options = list(
rowGroup = list(dataSrc = 1),
columnDefs = list(list(
targets = '_all',
className = 'dt-center'
))
),
selection = 'none') |>
formatRound(columns = c(4, 7:10), digits = 2) |>
formatStyle(columns = c(10), color = 'red')
# unique(newbin[total_iv >= 0.01, ]$variable)
feature <-
c(
"is_claim",
"policy_tenure",
"age_of_car",
"age_of_policyholder",
"area_cluster",
"population_density",
"f1", "f2", "f3", "f4", "f5")
# 分箱后有些特征的类别仍然较多的话,可以将类别的取值范围接近且 IV 值较小的合并
# 手动调整分箱
break_adj = list(
age_of_car = c(0.01, 0.02, 0.15),
population_density = c(15000, 20000, 30000, 40000))
bins_adj = woebin(df[,..feature],
y = "is_claim",
breaks_list = break_adj,
print_info = FALSE)
# 比对一下调整前后的差异
woebin_plot(bins$age_of_car)
woebin_plot(bins_adj$age_of_car)
# WOE转换,得到分箱后的特征
bins_df = data.table::rbindlist(bins_adj)
# 这里的bins_df的格式是数据框,也可以用前面筛选了 IV>0.01 的 newbin
df.woe = woebin_ply(df[,..feature], bins = bins_df, to = 'woe')
3.2.相关系数 🔗
查看特征之间的相关系数,是为了找出强相关(大于0.85)的特征。笔者的理解是,在高等代数里面,如果一个矩阵里面有完全相等的两列的话,矩阵的行列式等于0,行列式非满秩的矩阵是奇异矩阵,要解方程求解的时候可能会有无穷解或者无解。如果一个矩阵里面存在强相关的特征,也会导致类似情况的发生。为了避免无穷解或无解,强相关的特征要进行删减或转换。
- 在特征数量较少的情况下,可以直接用 corrplot 包的 corrplot 函数来绘制相关系数的图。
library(corrplot)
matrix <- cor(df.woe)
corrplot(matrix, method = "number")
- 但是在特征数量极多的情况下,也可以灵活运用 DT 表格来查看数据。
data.matrix <- as.data.table(matrix)
data.matrix$variable <- colnames(data.matrix)
iv <- iv(df.woe, y = 'is_claim')
data.matrix <- iv[data.matrix, on = 'variable']
datatable(
data.matrix,
rownames = FALSE,
extensions = 'FixedColumns',
options = list(
pageLength = 16,
dom = 't',
scrollX = TRUE,
fixedColumns = list(leftColumns = 2, rightColumns = 1),
columnDefs = list(list(
targets = '_all',
className = 'dt-center')))) |>
formatRound(columns = c(2:13), digits = 2) |>
formatStyle(columns = c(2),backgroundColor = 'pink')|>
formatStyle(columns = c(3:13), backgroundColor = styleInterval(c(-0.85, 0.85), c('green', 'white', 'red')))
3.3.逐步回归 🔗
用 bigstep 包的 stepwise 函数做逐步回归来筛选特征。这个包有个限制,一次运行的数据量不能超过100万。
library(bigstep)
library(dplyr)
y <- df.woe%>% pull(is_claim)
x <- df.woe[, -c('is_claim')]
datal <- prepare_data(y, x, type = "logistic")
step_lr <-
datal %>% reduce_matrix() %>% fast_forward() %>% stepwise()
summary(step_lr)
##
## Call:
## glm(formula = y ~ ., family = binomial(), data = as.data.frame(Xm[,
## -1]))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.6408 -0.4112 -0.3515 -0.2835 3.0102
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.68618 0.01743 -154.115 < 2e-16 ***
## policy_tenure_woe 1.07944 0.05643 19.130 < 2e-16 ***
## age_of_car_woe 1.25917 0.07028 17.917 < 2e-16 ***
## f2_woe 0.81153 0.11520 7.044 1.86e-12 ***
## age_of_policyholder_woe 0.64019 0.16008 3.999 6.36e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27860 on 58591 degrees of freedom
## Residual deviance: 27085 on 58587 degrees of freedom
## AIC: 27095
##
## Number of Fisher Scoring iterations: 6
3.4.K折交叉验证 🔗
用 glmnet 包的 cv.glmnet 函数做交叉验证来筛选特征。这个函数默认设置nfolds = 10
,即10折交叉验证。
library(Matrix)
library(glmnet)
#生成矩阵和标签
x <-df.woe[,-"is_claim"]
x <- data.matrix(x)
y <- df.woe$is_claim
# 设定随机种子
set.seed(2023)
fit.cv <- cv.glmnet(x, y, family = "binomial", type.measure = "auc")
#查看系数
#plot(fit.cv)
fit.cv$lambda.1se
#去掉系数没有具体数值的特征
coef(fit.cv, s = "lambda.1se")
下面结果中系数有数值的特征,就是筛选出来可以保留的。
## 11 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) -2.66135153
## policy_tenure_woe 0.71153340
## age_of_car_woe 0.65452866
## age_of_policyholder_woe .
## area_cluster_woe .
## population_density_woe .
## f1_woe 0.19768138
## f2_woe 0.19161653
## f3_woe 0.03489065
## f4_woe .
## f5_woe .
3.5.逻辑回归 🔗
查看用逻辑回归训练的模型基础效果,看看筛选出来的特征是否充分。AUC 值只有0.62,这说明特征是相当不够的。
#正式训练,拆分训练集和验证集
set.seed(2023)
ind <- sample(2, nrow(df.woe), replace = TRUE, prob = c(0.7, 0.3))
train <- df.woe[ind == 1,]
valid <- df.woe[ind == 2,]
#根据逐步回归结果去掉
#根据十折交叉验证结果去掉
#去掉系数显著性检验不通过
#去掉IV值超过0.5的特征
#去掉IV值小于0.02的特征
# 训练
model1 <- glmnet::glm(
#is_claim ~ policy_tenure_woe + age_of_car_woe + f2_woe + age_of_policyholder_woe
is_claim ~ policy_tenure_woe + age_of_car_woe + f1_woe + f2_woe + f3_woe
, family = binomial(link = 'logit'),
data = train)
summary(model1)
#预测
pred.train <- predict(model1, type = "response", newdata = train)
pred.valid <- predict(model1, type = "response", newdata = valid)
# 训练集
perf.train = scorecard::perf_eva(
pred = pred.train,
label = train$is_claim,
show_plot = c('ks', 'roc'))
# 验证集
perf.test = scorecard::perf_eva(
pred = pred.valid,
label = valid$is_claim,
show_plot = c('ks', 'roc'))
4.训练模型 🔗
训练模型时会遇到模型过拟合和欠拟合的问题,一般是先减少欠拟合,再减少过拟合。比如,先调整参数使 ROC 曲线下的 AUC 值不断增大,再调整参数使训练集、验证集 KS 值相差不超过0.01。
4.1.调参 🔗
library(xgboost)
##拆分数据集
set.seed(2022)
ind <- sample(2, nrow(df.woe), replace = TRUE, prob = c(0.7, 0.3))
factor <- c('is_claim',
'policy_tenure_woe',
'age_of_car_woe',
'f1_woe',
'f2_woe',
'f3_woe')
train <-
df.woe[ind == 1, ..factor]
valid <-
df.woe[ind == 2, ..factor]
#拆出训练集和验证集
trainlabel <- train[, 'is_claim']
validlabel <- valid[, 'is_claim']
traindata <- train[,-'is_claim']
validdata <- valid[,-'is_claim']
#若要使用xgb.train,需要先将原来的数据转化为矩阵格式
traindata <- as.matrix(traindata)
trainlabel <- as.matrix(trainlabel)
dtrain <- xgb.DMatrix(data = traindata, label = trainlabel)
validdata <- as.matrix(validdata)
validlabel <- as.matrix(validlabel)
dvalid <- xgb.DMatrix(data = validdata, label = validlabel)
#寻参得到的参数
set.seed(2152)
param1 = list(
booster = "gbtree",
objective = "binary:logistic",
eval_metric = "auc",
eta = 0.02,
max_depth = 6,
min_child_weight = 20,
gamma = 0.125,
subsample = 0.89,
colsample_bytree = 0.51,
n_estimators = 1000,
alpha = 0.1,
lambda = 1,
verbose = 1)
watchlist <- list(train = dtrain, test = dvalid)
bst1 <-
xgb.train(
data = dtrain,
params = param1,
nrounds = 1200,
# nthread = 16, # 默认1,值越大,执行越快
verbose = 1,
watchlist = watchlist,
early_stopping_rounds = 20
)
#预测结果
pred_train = predict(bst1, newdata = as.matrix(traindata),
ntree_limit = 1000)
pred_valid = predict(bst1, newdata = as.matrix(validdata),
ntree_limit = 1000)
#训练集
perf.test = scorecard::perf_eva(
pred = pred_train,
label = trainlabel,
show_plot = c('ks', 'roc', 'gain','lift'))
#验证集
perf.test = scorecard::perf_eva(
pred = pred_valid,
label = validlabel,
show_plot = c('ks', 'roc', 'gain','lift'))
4.2.网格寻参 🔗
下面这段寻参代码不记得是在哪里看到,然后笔者只是略微改了一点就拿来用了。寻参过程总是很慢,笔者也许应该去了解一下 foreach、parallel 那些实施并行计算的包,然后来改进这部分的执行效率。
# 寻参时,可不必拆分数据集
train <- df.woe
trainlabel <- train[, 'is_claim']
traindata <- train[, -is_claim]
traindata <- as.matrix(traindata)
trainlabel <- as.matrix(trainlabel)
dtrain <- xgb.DMatrix(data = traindata, label = trainlabel)
best_param = list()
best_seednumber = 1234
# best_logloss = Inf
# best_logloss_index = 0
best_auc = Inf
best_auc_index = 0
for (iter in 1:100) {
param <- list(
objective = "binary:logistic",
# eval_metric = c("logloss"),
eval_metric = "auc",
max_depth = sample(5:8, 1),
eta = runif(1, .01, .3),
gamma = runif(1, 0.0, 0.2),
subsample = runif(1, .6, .9),
colsample_bytree = runif(1, .5, .8),
min_child_weight = sample(1:40, 1)
)
cv.nround = 200 #运行轮次
cv.nfold = 5
seed.number = sample.int(10000, 1)[[1]]
set.seed(seed.number)
mdcv <- xgb.cv(
data = dtrain,
params = param,
nthread = 16, # 默认1,值越大,执行越快
nfold = cv.nfold,
nrounds = cv.nround,
verbose = T,
early_stop_round = 8,
maximize = FALSE
)
# min_logloss = min(mdcv[, test_logloss_mean])
# min_logloss_index = which.min(mdcv[, test_logloss_mean])
#
# if (min_logloss < best_logloss) {
# best_logloss = min_logloss
# best_logloss_index = min_logloss_index
# best_seednumber = seed.number
# best_param = param
# }
min_auc = min(mdcv$evaluation_log[, test_auc_mean])
min_auc_index = which.min(mdcv$evaluation_log[, test_auc_mean])
if (min_auc < best_auc) {
best_auc = min_auc
best_auc_index = min_auc_index
best_seednumber = seed.number
best_param = param
}
}
best_param
best_seednumber
# nround=best_auc_index
# set.seed(best_seednumber)
# md <-
# xgb.train(
# data = dtrain,
# params = best_param,
# nrounds = nround,
# nthread = 16)
4.3.特征重要性 🔗
可以用 xgb.importance 函数查看特征重要性。讲真,笔者也看不懂下面这些数值,只知道 Gain 值自动从高到低排序,代表模型中这些特征重要性逐渐下降。
importance <- xgb.importance(feature_names = colnames(dtrain), model = bst1)
head(importance)
## Feature Gain Cover Frequency
## 1: age_of_car_woe 0.36945417 0.44952249 0.45714286
## 2: policy_tenure_woe 0.31074011 0.18814007 0.21428571
## 3: f1_woe 0.13229080 0.09504380 0.08571429
## 4: f3_woe 0.12491070 0.18622903 0.17142857
## 5: f2_woe 0.06260423 0.08106461 0.07142857
刚翻了下xgboost 包最新的文档,可以查看 SHAP 值,还能画图。试了试,报错,不明所以中……笔者对 SHAP 很陌生,连同上网翻的资料一、资料二一并先记在这里,回头再琢磨琢磨单独写一篇。
4.4.lime 🔗
用 lime 包可以查看单个样本具体受到哪些特征影响,是正相关还是负相关,影响程度多大等等。
- 蓝色表示正相关,红色表示负相关。
- case 表示样本的序号,Label 是指是否在6月内索赔的结果(1表示是,0表示否)。
- 条形的长度表示特征重要性的大小(即权重,下图中限制仅展示重要性最高的前5个特征)。
需要注意的是,由于 lime 包不一定支持 R 中所有的模型,所以要事先用 model_type 函数检查一下,看需要解释的模型是否被支持。否则,强行解释会报下面这个错。
Error: The class of model must have a model_type method. See ?model_type to get an overview of models supported out of the box.
library(lime)
# 检查前面 xgboost 模型 bst1 是否可以支持
model_type(bst1)
# lime 不支持 xgboost 包的特殊数据格式,得转换成普通数据框
traindata <- as.data.frame(traindata)
explainer <- lime(traindata, bst1)
#查看单个样本
validdata<-as.data.frame(validdata)
explanation <-
lime::explain(validdata[1,],
explainer,
n_labels = 1,
n_features = 5)
plot_features(explanation)
5.验证 🔗
在建好模型以后,在正式预测之前,需要做一番验证:
-
数据分布是否保持一致。
-
特征是否稳定。假如一个特征在训练集IV值较低,在验证集IV值较高,会导致验证集效果比训练集更好。但是一个特征IV值变化很大是很有风险的,因为不敢确定该特征在未来的测试集上是否依然表现良好。
PSI特征稳定性 🔗
library(scorecard)
pred_list2 = list(pred_train, pred_valid)
names(pred_list2) <- c("train", "valid")
label_list2 = list(train$label, valid)
names(label_list2) <- c("train", "valid")
psi2 = perf_psi(score = pred_list2, label = label_list2)
psi2$psi # psi data frame
# 分数分布
library(Cairo)
Cairo(
file = "Cairo_PNG_psi_dpi.png",
type = "png",
width = 1200,
height = 800
)
psi2$pic
dev.off()
IV特征重要性 🔗
iv_train <- iv(train, y = "label")
iv_valid <- iv(valid, y = "label")
colnames(iv_train)[2]<-"iv_train"
colnames(iv_valid)[2]<-"iv_valid"
iv.data <- iv_train[iv_valid, on = "variable"]
iv.data$iv_train_l <-
ifelse(iv.data$iv_train < 0.02, '无预测力',
ifelse(iv.data$iv_train < 0.1,'弱',
ifelse(iv.data$iv_train < 0.3, '中等', '强')))
iv.data$iv_valid_l <-
ifelse(iv.data$iv_valid < 0.02, '无预测力',
ifelse(iv.data$iv_valid < 0.1, '弱',
ifelse(iv.data$iv_valid < 0.3, '中等', '强')))
iv.data
6.预测 🔗
上节建模时设定的目标函数是“binary:logistic”,评价指标是“auc”,预测时输出的结果是概率,即0至1之间的数值。通常在预测概率的基础上再设定阈值,判定好与坏的结果标签。实际业务中设定阈值通常与验证结果所花费的成本相关。若想预测结果直接是0或1,须将目标函数改为“binary:hinge”,评价指标也一并修改。
# 用分箱的结果对测试数据集做 WOE 转换
dtest <- data[type == 'test', ..feature]
dtest.woe = scorecard::woebin_ply(dtest, bins = bins_df, to = 'woe')
# 用训练好的 xgboost 模型 bst1 进行预测
testdata <- as.matrix(dtest.woe[, c('policy_tenure_woe',
'age_of_car_woe',
'f1_woe',
'f2_woe',
'f3_woe')])
pred_test = predict(bst1, newdata = testdata, ntree_limit = 1000)
-
这份数据最奇葩的一点就是,它真的一点都不奇葩。笔者几次三番想换份数据,又无法干脆利落地做到轻易放弃,唉。 ↩︎