用tidymodels进行数学建模
卿成林
2025-01-17
1 概要
在这篇报告中，我先学习了使用tidymodels包进行建模包括数据准备、数据预处理、模型训练和模型评估；然后对一个较为复杂的例子进行建模并筛选出最佳模型进行预测。
2 介绍
经过查阅，我选择了tidyverse与tidymodels系列，tidyverse为R语言提供了一种高效统一的处理框架，而tidymodels基于tidyverse对机器学习和统计建模的各种R工具进行了更新和整合，并且将建模分成了如下四个步骤：
- 数据准备
- 数据预处理
- 模型的训练和调整
- 模型验证
相应的，tidymodels包括了以下的核心包：
- rsample：数据分类和重抽样
- recipes：用于特征工程的数据预处理
- parsnip：对一系列模型进行建模的统一接口
- tune：调整超参数，对预处理过程进行优化
- yardstick：使用性能指标度量模型的有效性
- workflows：将模型的预处理、建模、后处理整合在一起
于是我按照这个思路进行建模学习。
3 学习建模
3.1 parsnip建模
建模第一步应该是数据预处理，包括数据清洗和数据划分，数据清洗根据具体问题做出相应操作，一般包括删去或插补缺失数据、调整数据格式等，对于有监督学习，数据需划分为训练集与测试集，可以使用rsample包的initial_split函数；
parship包整合了绝大部分官方模型和第三方模型，将参数分为主参数和引擎参数，保持了参数的一致性，具体可以参看我下面的代码。
这个链接给出parship包支持的各种模型：https://www.tidymodels.org/find/parsnip/
常用的函数：
- fit函数进行拟合
- extract_*函数提取拟合结果
- predict函数对测试集进行预测
下面就是我们使用线性回归模型对Ames数据中房屋出售价格与经纬度关系进行训练与预测的程序
library(tidymodels)
tidymodels_prefer()
#载入数据ames
library(modeldata)
data(ames)
#对数化
ames <- ames |>
mutate(Sale_Price = log10(Sale_Price))
#按照出售价格分层划分数据
set.seed(123)
ames_split <- initial_split(ames, prop = 0.95,strata = Sale_Price)
ames_train <- training(ames_split)
ames_test <- testing(ames_split)
#设置模型
lm_model <-
linear_reg() |>
set_engine("lm")
#训练模型
lm_form_fit <-
lm_model |>
fit(Sale_Price ~ Longitude + Latitude,
data = ames_train)
#对测试集进行预测
ames_test_sm<-ames_test|>slice(1:10)
ames_test_sm |>
select(Sale_Price) |>
bind_cols(predict(
lm_form_fit,
ames_test_sm )) |>
bind_cols(predict(
lm_form_fit,
ames_test_sm,
type = "pred_int" ))## # A tibble: 10 × 4
## Sale_Price .pred .pred_lower .pred_upper
## <dbl> <dbl> <dbl> <dbl>
## 1 5.02 5.23 4.91 5.54
## 2 5.29 5.29 4.97 5.61
## 3 5.06 5.24 4.92 5.56
## 4 5.46 5.32 5.00 5.64
## 5 5.52 5.30 4.98 5.61
## 6 5.31 5.29 4.97 5.60
## 7 5.28 5.32 5.00 5.64
## 8 5.19 5.27 4.95 5.58
## 9 5.24 5.19 4.87 5.50
## 10 5.33 5.20 4.88 5.51可以看出代码十分简洁，结果也较为不错。
3.2 recipes特征工程
recipes包用来对数据进行预处理，比如数据的转换和编码、插补缺失值、对多维数据提取特征进行降维等。
常用步骤函数step_*()来进行处理。
下面这段代码对Ames数据进行了如下处理
- 对居住面积进行对数变换
- 将1%的房屋划分为other；类
- 将非数字格式列转换为数字格式
- 添加Gr_Liv_Area和Bldg_Type_的交叉作用
- 添加样条变换项
simple_ames <-
recipe(Sale_Price ~ Neighborhood + Gr_Liv_Area
+ Year_Built + Bldg_Type
+ Latitude + Longitude,
data = ames_train) |>
step_log(Gr_Liv_Area, base = 10) |>
step_other(Neighborhood, threshold = 0.01) |>
step_dummy(all_nominal_predictors()) |>
step_interact( ~ Gr_Liv_Area:starts_with("Bldg_Type_") ) |>
step_ns(Latitude, Longitude, deg_free = 20)3.3 workflow工作流
workflow将模型选择和数据预处理这两部分连接起来，形成一个workflow的对象；
常用的一些函数包括：
- add_formula 以公式形式添加预处理方法
- add_variables 以分类形式添加预处理方法
- add_recipes 将recipes处理后的数据添加进预处理方法
注意一个workflow对象只能有一个上述函数。。
前面对models的函数比如fit、extract_*、predict仍能对workflow类使用。
将前述预处理对象添加进workflow并训练、预测
#设置模型
lm_model <-
linear_reg() |>
set_engine("lm")
#将模型、recipes加入workflow
lm_wflow <- workflow() |>
add_model(lm_model) |>
add_recipe(simple_ames)
#只是设置数据，仍需拟合
lm_fit <- fit(lm_wflow, ames_train)
#对测试集进行预测
ames_test_sm<-ames_test|>slice(1:10)
ames_test_sm |>
select(Sale_Price) |>
bind_cols(predict(
lm_fit,
ames_test_sm )) |>
bind_cols(predict(
lm_fit,
ames_test_sm,
type = "pred_int" ))## # A tibble: 10 × 4
## Sale_Price .pred .pred_lower .pred_upper
## <dbl> <dbl> <dbl> <dbl>
## 1 5.02 5.07 4.92 5.22
## 2 5.29 5.28 5.13 5.44
## 3 5.06 5.06 4.90 5.21
## 4 5.46 5.45 5.30 5.60
## 5 5.52 5.54 5.39 5.69
## 6 5.31 5.33 5.18 5.48
## 7 5.28 5.30 5.15 5.45
## 8 5.19 5.23 5.08 5.38
## 9 5.24 5.20 5.05 5.35
## 10 5.33 5.30 5.14 5.45可以看出预处理后预测结果比之前好不少。
引入workflow最方便的一点的我觉得是可以同时放入多个模型，这个之后会用到。
3.4 yardstick评估模型效果
一旦我们训练了一个模型，我们就需要一种方法来衡量该模型对新数据的预测效果。yardstick包提供了很多评估模型效果的函数，并且这些函数十分有规律，第一个参数永远是数据集，第二个参数永远是真实结果，第三个参数永远是预测结果。
列出以下一些常用函数
通用函数
- metrics()只需提供数据集、真实结果、预测结果以及模型类别（回归、分类）
- metric_set()集合多个评估指标函数
- collect_metrics()收集参数搜索评估结果
回归函数 - rmse 均方根误差 - rsq R平方 - mae 平均绝对误差
分类函数 - recall 召回率 - accuracy 精度 - conf_mat 混淆矩阵 - roc_curve ROC曲线
我们对前面的模型进行评估
#预测
ames_test_res <- predict(
lm_fit,
new_data = ames_test |> select(-Sale_Price))
ames_test_res <- ames_test |>
select(Sale_Price) |>
bind_cols(ames_test_res)
#评估指标集合函数
ames_metrics <- metric_set(rmse, rsq, mae)
ames_metrics(
ames_test_res,
truth = Sale_Price,
estimate = .pred)## # A tibble: 3 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 rmse standard 0.0647
## 2 rsq standard 0.844
## 3 mae standard 0.0517下面我们使用randomforest法再对Ames数据集进行建模预测，以与线性模型对比
#设置模型
rf_model <-
rand_forest(trees = 1000) |>
set_engine("ranger") |>
set_mode("regression")
#建立workflow
rf_wflow <-
workflow() |>
add_model(rf_model) |>
add_formula(
Sale_Price ~ Neighborhood + Gr_Liv_Area
+ Year_Built + Bldg_Type
+ Latitude + Longitude)
#训练模型
rf_fit <- fit(rf_wflow, ames_train)
#预测
ames_test_res_rf <- predict(
rf_fit,
new_data = ames_test |> select(-Sale_Price))
ames_test_res_rf <- ames_test |>
select(Sale_Price) |>
bind_cols(ames_test_res_rf)
#评估模型
ames_metrics <- metric_set(rmse, rsq, mae)
ames_metrics(
ames_test_res_rf,
truth = Sale_Price,
estimate = .pred)## # A tibble: 3 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 rmse standard 0.0573
## 2 rsq standard 0.874
## 3 mae standard 0.0431可以看出randomforest模型比线性模型更好。
4 混凝土强度建模
至此我已经掌握使用tidymodels来建模并评估模型的基础方法，现在我们来按照上面的思路对一个具体的问题建立多个模型并筛选出较优模型。
考虑modeldata扩展包中的concrete数据。有1030个观测，因变量是混凝土压缩强度compressive_strength，有7个混凝土成分变量作为自变量，单位是每立方米千克数，另一自变量age是测试强度等待的天数。
4.1 划分数据
先观察数据
## tibble [1,030 × 9] (S3: tbl_df/tbl/data.frame)
## $ cement : num [1:1030] 540 540 332 332 199 ...
## $ blast_furnace_slag : num [1:1030] 0 0 142 142 132 ...
## $ fly_ash : num [1:1030] 0 0 0 0 0 0 0 0 0 0 ...
## $ water : num [1:1030] 162 162 228 228 192 228 228 228 228 228 ...
## $ superplasticizer : num [1:1030] 2.5 2.5 0 0 0 0 0 0 0 0 ...
## $ coarse_aggregate : num [1:1030] 1040 1055 932 932 978 ...
## $ fine_aggregate : num [1:1030] 676 676 594 594 826 ...
## $ age : int [1:1030] 28 28 270 365 360 90 365 28 28 28 ...
## $ compressive_strength: num [1:1030] 80 61.9 40.3 41 44.3 ...去掉重复的数据
concrete <-
concrete |>
group_by(across(-compressive_strength)) |>
summarize(
compressive_strength = mean(compressive_strength),
.groups = "drop")
nrow(concrete)## [1] 992能看出成功去除掉38个数据
分层抽样，对训练集产生重复5次的十折交叉验证集
4.2 预处理
有些模型要求自变量标准化，如神经网络、最近邻方法、支持向量机，另外一些模型则还需要使用正交多项式变换及其交互作用效应等自变量变换。分别产生两个recipes:
4.3 设置模型
先选择如下模型
library(rules)
library(baguette)
# 正则化的线性回归
linear_reg_spec <-
linear_reg(
penalty = tune(),
mixture = tune()) |>
set_engine('glmnet')
# 神经网络
nnet_spec <-
mlp(hidden_units = tune(), penalty = tune(), epochs = tune()) |>
set_engine("nnet") |>
set_mode("regression")
# 多元自适应回归样条
mars_spec <-
mars(
prod_degree = tune()) |>
set_engine('earth') |>
set_mode('regression')
# 支持向量机，使用径向基
svm_r_spec <-
svm_rbf(
cost = tune(),
rbf_sigma = tune() ) |>
set_engine('kernlab') |>
set_mode('regression')
# 支持向量机，使用多项式基
svm_p_spec <-
svm_poly(
cost = tune(),
degree = tune() ) |>
set_engine('kernlab') |>
set_mode('regression')
# 最近邻方法
knn_spec <-
nearest_neighbor(
neighbors = tune(),
weight_func = tune(),
dist_power = tune()) |>
set_engine('kknn') |>
set_mode('regression')
# 回归判别树
cart_spec <-
decision_tree(
min_n = tune(),
cost_complexity = tune()) |>
set_engine('rpart') |>
set_mode('regression')
# 装袋法
bag_cart_spec <-
bag_tree() |>
set_engine('rpart', times = 50L) |>
set_mode('regression')
# 随机森林
rf_spec <-
rand_forest(
mtry = tune(),
min_n = tune(),
trees = 1000) |>
set_engine('ranger') |>
set_mode('regression')
# XGBOOST
xgb_spec <-
boost_tree(
tree_depth = tune(),
trees = tune(),
learn_rate = tune(),
min_n = tune(),
loss_reduction = tune(),
sample_size = tune()) |>
set_engine('xgboost') |>
set_mode('regression')
# 基于规则的模型
cubist_spec <-
cubist_rules(
committees = tune(),
neighbors = tune()) |>
set_engine('Cubist')
#设置神经网络法的隐藏单元个数
nnet_param <-
nnet_spec |>
extract_parameter_set_dials() |>
update(hidden_units = hidden_units(c(1, 27)))4.4 建立workflow
用workflow集合将预处理和模型匹配
#建立需标准化自变量的workflow集合
normalized <-
workflow_set(
preproc = list(
normalized = normalized_rec),
models = list(
SVM_radial = svm_r_spec,
SVM_poly = svm_p_spec,
KNN = knn_spec,
neural_network = nnet_spec) )
#给nnet_spe加入调优参数
normalized <-
normalized |>
option_add(
param_info = nnet_param,
id = "normalized_neural_network")
#建立需自变量变换的workflow集合
with_features <-
workflow_set(
preproc = list(
full_quad = poly_recipe),
models = list(
linear_reg = linear_reg_spec,
KNN = knn_spec) )
#建立无需特殊处理，仅需指定因变量、自变量的workflow模型
model_vars <-
workflow_variables(
outcomes = compressive_strength,
predictors = everything())
no_pre_proc <-
workflow_set(
preproc = list(simple = model_vars),
models = list(
MARS = mars_spec,
CART = cart_spec,
CART_bagged = bag_cart_spec,
RF = rf_spec,
boosting = xgb_spec,
Cubist = cubist_spec) )
#将上面三个workflow集合合并成一个大的workflow集合
all_workflows <-
bind_rows(
no_pre_proc,
normalized,
with_features) |>
# 简化工作流程的id:
mutate(
wflow_id = gsub(
"(simple_)|(normalized_)", "", wflow_id))4.5 模型训练
我们采取一种赛跑的方式对模型进行调优，在调优过程中放弃表现很差的模型，以此加快调优的速度
library(doParallel)
library(finetune)
# 网格赛跑调优共用的选项
race_ctrl <- control_race(
save_pred = TRUE,
parallel_over = "everything",
save_workflow = TRUE )
# 产生一个集群对象并注册
cl <- makePSOCKcluster(6)
registerDoParallel(cl)
time0 <- proc.time()[3]
# 实际对每个工作流程进行调优
race_results <-
all_workflows |>
workflow_map(
"tune_race_anova",
seed = 1503,
resamples = concrete_folds,
grid = 25,
control = race_ctrl )
time1 <- proc.time()[3] - time0
cat("用时：", time1, "秒\n")## 用时： 556.21 秒## # A workflow set/tibble: 12 × 4
## wflow_id info option result
## <chr> <list> <list> <list>
## 1 MARS <tibble [1 × 4]> <opts[3]> <race[+]>
## 2 CART <tibble [1 × 4]> <opts[3]> <race[+]>
## 3 CART_bagged <tibble [1 × 4]> <opts[3]> <rsmp[+]>
## 4 RF <tibble [1 × 4]> <opts[3]> <race[+]>
## 5 boosting <tibble [1 × 4]> <opts[3]> <race[+]>
## 6 Cubist <tibble [1 × 4]> <opts[3]> <race[+]>
## 7 SVM_radial <tibble [1 × 4]> <opts[3]> <race[+]>
## 8 SVM_poly <tibble [1 × 4]> <opts[3]> <race[+]>
## 9 KNN <tibble [1 × 4]> <opts[3]> <race[+]>
## 10 neural_network <tibble [1 × 4]> <opts[4]> <race[+]>
## 11 full_quad_linear_reg <tibble [1 × 4]> <opts[3]> <race[+]>
## 12 full_quad_KNN <tibble [1 × 4]> <opts[3]> <race[+]>4.6 作图比较
autoplot(
race_results,
rank_metric = "rmse",
metric = "rmse",
select_best = TRUE) +
geom_text(aes(
y = mean - 1/2, label = wflow_id),
angle = 90, hjust = 1) +
lims(y = c(3.0, 9.5)) +
theme(legend.position = "none")
由图可以看出是XGBOOST模型最佳。
4.7 拟合
选择上面筛选出的最优模型XGBOOST，在训练集上拟合并在测试集上预测，绘制出预测值与真实值的对比图
best_results <-
race_results |>
# 提取一个工作流程的调优结果
extract_workflow_set_result("boosting") |>
# 从调优结果中选择最优的超参数组合
select_best(metric = "rmse")
#训练集上拟合并在测试集上预测
boosting_test_results <-
race_results |>
extract_workflow("boosting") |>
finalize_workflow(best_results) |>
last_fit(split = concrete_split)
#用collect_predictions()提取测试集上的预测值， 与真实值对比作图
boosting_test_results |>
collect_predictions() |>
ggplot(aes(
x = compressive_strength, y = .pred)) +
# 作y=x直线
geom_abline(color = "green", lty = 2) +
geom_point(alpha = 0.5) +
coord_obs_pred() +
labs(x = "observed", y = "predicted")