R语言tidyverse数据处理建模案例

管道%>%

左连接left_join()

筛选行 filter(条件)

行排序arrange()

选择列select()

修改(计算)列mutate()

分组汇总group_by()、%>%、summarise() 计数：count()

数据处理案例：

企业上游业务量：

• 企业上游年均业务量:$X_1=12\times \frac{企业上游交易总次数}{企业数据的月份数}$

• 企业上游业务量年平均变化率：
  KaTeX parse error: Undefined control sequence: \notag at position 88: …数}{企业上游k月交易次数} \̲n̲o̲t̲a̲g̲ ̲

第$j$个企业上游业务量用$X_{1j}$表示，相应的企业上游业务量变化率用$X_{1j}^{’}(j=1,2,\cdots ,123)$表示.

企业下游业务量：

• 企业上游年均业务量:$X_2=12\times \frac{企业下游交易总次数}{企业数据的月份数}$

• 企业下游业务量年平均变化率：
  KaTeX parse error: Undefined control sequence: \notag at position 88: …数}{企业下游k月交易次数} \̲n̲o̲t̲a̲g̲ ̲

第$j$个企业下游业务量用$X_{2j}$表示，相应的企业上游业务量变化率用$X_{2j}^{’}(j=1,2,\cdots ,123)$表示.

企业毛利润：

• 企业年均毛利润$X_3=12\times \frac{企业总收益-企业直接总成本}{企业数据的月份数}$

• 企业的年均毛利润变化率:

KaTeX parse error: Undefined control sequence: \notag at position 90: …毛利润}{企业k月的毛利润} \̲n̲o̲t̲a̲g̲ ̲

第$j$个企业下游业务量用$X_{3j}$表示，相应的企业上游业务量变化率用$X_{3j}^{’}(j=1,2,\cdots ,123)$表示.

数据处理：

对上游业务量$X_{1j}$，下游业务量$X_{2j}$和毛利润率$X_{3j}$做标准化处理，即
$$X_{ij}^{\ast }=\frac{X_{ij}-{\overline{X}}_i}{{\sigma }_i},i=1,2,3;j=1,2,\cdots ,123$$
企业信贷风险的因素分析

(1)信誉评级量化(D级"一票否决")
$$\text{第 j }(j=1,2,3,\dots )\text{ 企业的信誉评级指标: }{X}_{0j}=\{\begin{array}{l}5，\text{ 信誉评级为 A 级，}\\ 3，\text{ 信誉评级为 B 级，}\\ 1，\text{ 信誉评级为 C 级，}\\ 0，\text{ 信誉评级为 D 级. }\end{array}$$
(2)是否有违约量化(违约"一票否决")
$$\text{第 j }(j=1,2,3,\dots )\text{ 企业的信誉评级指标: }{Y}_{0j}=\{\begin{array}{l}1，\text{ 企业无违约，}\\ 0，\text{ 企业有违约. }\end{array}$$
企业信贷风险的量化模型

(1)企业的实力指标

​ 综合实力指标为三项实力指标的动态加权求和：
$$S={\lambda }_1(X_1^{{}^{\prime }})X_1+{\lambda }_2(X_2^{{}^{\prime }})X_2+{\lambda }_3(X_3^{{}^{\prime }})X_3$$
​ 动态加权函数取偏大的S型分布：
$${\lambda }_i(X_i^{{}^{\prime }})=\{\begin{array}{ll}2-e^{-3|X_i^{{}^{\prime }}|}& ,X_i^{{}^{\prime }}⩾0\\ e^{-3|X_i^{{}^{\prime }}|}& ,X_i^{{}^{\prime }}<0\end{array}(i=1,2,3)$$
​ 各企业的综合实力指标值：
$$S_j={\lambda }_1(X_{1j}^{{}^{\prime }})X_{1j}+{\lambda }_2(X_{2j}^{{}^{\prime }})X_{2j}+{\lambda }_3(X_{3j}^{{}^{\prime }})X_{3j}(j=1,2,\dots ,123)$$
$X_1^{{}^{\prime }}$对应程序中的dfdX1，以此类推.

(2) 企业信誉度指标

​ 企业的信誉度指标由信誉评级和是否有违约构成，即
$$C_j=X_{0j}{Y}_{0j}(j=1,2,\dots ,123)$$

(3) 企业实力+信誉指标

​ 企业的信贷风险由企业的综合实力和信誉指标决定，而且银行对信誉评级为D和有违约记录的企业“一票否决”，则企业的实力+信誉指标：
$$S_{Cj}=C_j{S}_j(j=1,2,\cdots ,123)$$
(4) 企业的信贷风险指标

​ 企业的实力和信誉决定了信贷风险，信贷风险与实力+信贷指标不应该是线性关系，不难说明呈S型曲线的关系。

​ 利用生物学中常用的$Sigmoid$函数:
$$S_{Cj}{}^{\prime }=\{\begin{array}{ll}\frac{1}{1+e^{-S_{Cj}}}& ,当S_{Cj}\ne 0时,\\ 0& ,当S_{Cj}=0时\end{array}(j=1,2,,\cdots ,123).$$
则每个企业的信贷风险指标值：
$$R_j=1-S_{Cj}^{{}^{\prime }}\in [0,1](j=1,2,\cdots ,123)$$
​ 根据各企业信贷风险值R进行分类，不妨分为四类.事实上，企业的实力+信誉值越高，风险值越低，贷款违约的可能性越小，贷款利率就应越低.

各企业信贷风险的等级分类

信贷风险值	等级分级	企业数	各信誉评级分布
[0,0.2)	1级	63	A级26，B级20，C级20
[0.2,0.5]	2级	6	B级1，C级5
(0.5,1)	3级	27	A级1，B级16，C级10
1	4级	27	B级1，C级2，D级24

代码：

library(tidyverse)
library(readxl)
library(writexl)
library(tsibble) ##使用year、month函数

##读取数据
info <- read_xlsx("data/附件1：123家有信贷记录企业的相关数据.xlsx", sheet = 1)
input <- read_xlsx("data/附件1：123家有信贷记录企业的相关数据.xlsx", sheet = 2)
output <- read_xlsx("data/附件1：123家有信贷记录企业的相关数据.xlsx", sheet = 3)

##统计发票状态
input %>%
  count(发票状态)

output %>%
  count(发票状态)

## 开票日期列转化成日期型，删除作废发票
input <- input %>%
  mutate(开票日期 = yearmonth(开票日期)) %>%
  filter(发票状态 == "有效发票")
input

output <- output %>%
  mutate(开票日期 = yearmonth(开票日期)) %>%
  filter(发票状态 == "有效发票")
output

## 特征工程：构建新特征

# 企业上游业务量
dfX1 <- input %>%
  group_by(企业代号) %>%
  summarise(X1 = 12 * n() / n_distinct(开票日期))

# 企业上游业务量年均变化率 dx1,注意E104总共只有一笔交易

dfdX1 <- input %>%
  group_by(企业代号, 开票日期) %>%
  summarise(交易量 = n()) %>%
  mutate(月变化 = 交易量 - lag(交易量),
            月变化率 = 月变化 / lag(交易量)) %>%
  summarise(dX1 = 12 * mean(月变化率, na.rm = T))

# 企业下游业务量
dfX2 <- output %>%
  group_by(企业代号) %>%
  summarise(X2 = 12 * n() / n_distinct(开票日期))

# 企业上游业务量年均变化率 dx2

dfdX2 <- output %>%
  group_by(企业代号, 开票日期) %>%
  summarise(交易量 = n()) %>%
  mutate(月变化 = 交易量 - lag(交易量),
            月变化率 = 月变化 / lag(交易量)) %>%
  summarise(dX2 = 12 * mean(月变化率, na.rm = T))

# 企业年均毛利润X3
costs <- input %>%
  group_by(企业代号) %>%
  summarise(月份数 = n_distinct(开票日期), 总成本 = sum(金额))

sales <- output %>%
  group_by(企业代号) %>%
  summarise(总收益 = sum(金额))

dfX3 <- costs %>%
  left_join(sales, by = "企业代号") %>%
  mutate(X3 = 12 * (总收益 - 总成本) / 月份数) %>%
  select(-c(月份数, 总成本, 总收益))

# 企业年均毛利润变化率dX3
costs_mon <- input %>%
  group_by(企业代号, 开票日期) %>%
  summarise(月成本 = sum(金额))

sales_mon <- output %>%
  group_by(企业代号, 开票日期) %>%
  summarise(月收益 = sum(金额))

dfdX3 <- costs_mon %>%
  left_join(sales_mon, by = c("企业代号", "开票日期")) %>%
  mutate(月利润 = 月收益 - 月成本) %>%
  drop_na(月利润) %>%
  mutate(月变化率 = (月利润 - lag(月利润)) / lag(月利润)) %>%
  summarise(dX3 = 12 * mean(月变化率, na.rm = T))

## 将各个新特征合并到一个表中
df <- list(dfX1, dfdX1, dfX2, dfdX2, dfX3, dfdX3) %>%
  reduce(left_join, by = "企业代号")

write_xlsx(df, "data/企业信贷指标数据.xlsx")

## 用各列的均值插补缺失值
df <- df %>%
  mutate(across(where(is.numeric), ~ ifelse(is.nan(.x), NA, .x)),
         across(where(is.numeric), ~ naniar::impute_mean(.x)))

## 标准化处理
#注：用于数据相差很大时,可以避免归一化"大数吃小数"
df <- df %>%
  mutate(across(starts_with("X"), ~ (.x - mean(.x)) / sd(.x)))
df

## 企业实力指标：动态加权
# 注意，只有相互独立的指标才可以线性相加
SType = function(x) {
  ifelse(x >= 0, 2 - exp(-3 * abs(x)), exp(-3 * abs(x)))
}

ggplot() +
  xlim(-1, 1) +
  geom_function(fun = SType)

df <- df %>%
  mutate(S = SType(dX1) * X1 + SType(dX2) * X2 + SType(dX3) * dX3)
df

## 企业信誉度指标
rating <- info %>%
  mutate(X0 = case_when(信誉评级 == "A" ~ 5,
                        信誉评级 == "B" ~ 3,
                        信誉评级 == "C" ~ 1,
                        信誉评级 == "D" ~ 0,
                        ),
         Y0 = ifelse(是否违约 == "是", 0, 1),
         C = X0 * Y0) %>%
  select(企业代号, 信誉评级, C)
rating

## 企业实力 + 信誉综合指标
scores <- rating %>%
  left_join(df, by = "企业代号") %>%
  mutate(Sc = C * S)
scores
summary(scores$Sc)

##企业信贷风险指标
Sigmoid = function(x) {
  ifelse(x != 0, 1 / (1 + exp(-x)), 0)
}

ggplot() +
  xlim(-8, 8) +
  geom_function(fun = Sigmoid)

scores <- scores %>%
  mutate(R = 1 - Sigmoid(Sc))
scores

## 散点图
scores %>%
  mutate(企业代号 = parse_number(企业代号)) %>%
  ggplot(aes(企业代号, R)) +
  geom_point()
final <- scores %>%
  select(企业代号, 信誉评级, R) %>%
  mutate(等级分级 = case_when(R < 0.2 ~ "1级",
                              R <= 0.5 ~ "2级",
                              R < 1 ~ "3级",
                              R == 1 ~ "4级"))
## 统计分组频数
final %>%
  count(等级分级)
final %>%
  count(等级分级, 信誉评级) %>%
  pivot_wider(names_from = 信誉评级, values_from = n)

## 机器学习模型
df <- df %>%
  left_join(info, by = "企业代号") %>%
  rename(default = 是否违约, Rank = 信誉评级) %>%
  mutate(default = factor(default), Rank = factor(Rank))

library(mlr3verse) #此包不支持中文变量名,字符型变量需要变为因子

## 创建分类任务
task <- as_task_classif(df[-c(1, 8:9)], target = "Rank")

## 选择学习器,并设置两个超参数：最大深度，最小分支节点
learner <- lrn("classif.ranger", num.trees = 200, min.node.size = 20)

## 划分训练集,测试集
set.seed(123)
split <- partition(task, ratio = 0.8) #80%为训练集

## 训练模型
learner$train(task, row_ids = split$train)

## 模型预测
predictions <- learner$predict(task, row_ids = split$test)

## 模型评估
predictions$confusion #混淆矩阵
predictions$score(msr("classif.acc")) #准确率

#建议将信誉等级B和C合并变成一类

本文章内容为张敬信老师B站视频学习笔记