37 整洁建模
tidymodels包提供了针对统计建模用户简便易用的建模功能和使用规范。
R有许多建模用的扩展包，
不同的扩展包在解决类似的回归、判别问题时，
可能用了相近但不相同，
或者完全不同的界面接口。
tidymodels试图统一这些接口，
使得用户可以用相同的方法调用不同的计算包。
tidyr包的broom函数、nest、unnest函数提供了数据框分组建模的方便功能。
参见
- Max Kuhn and Julia Silge: Tidy Modeling with R https://www.tmwr.org/
37.1 建模一般步骤
在建模之前， 需要先获取数据。 统计建模用的数据一般都保存为数据框格式。
必须对数据有充分的了解， 否则在后续使用模型时， 就无法判断模型所对应的各种理论假定是否满足， 比如， 数据是否能够认为是独立的， 数据是否没有偏差地代表了所要研究的对象。
获得数据后， 一般需要花费时间进行“数据清理”。 数据中可能有许多错误、不一致、缺失内容， 都需要进行整理修改。 对缺失数据， 常见的办法是仅保留没有任何缺失内容的案例， 但是这常常会使得保留的数据量太小， 另外有可能删去的内容使得保留的数据带有偏差。 另一种对待缺失数据的方法是进行插补， 这又有简单插补， 如用均值代替， 与利用相关性插补等方法。
在建模之前， 应增加对数据的直观了解， 比如， 数据的收集方式， 每个属性（变量）的可取值范围， 变量分布情况， 变量之间的关系， 观测有没有明显分组聚集， 等等。 这些研究称为探索性数据分析（EDA）。
随后的工作是建模、评估、使用。 模型大致可以分为:
- 描述性模型， 如描述统计量，直方图，散点图，散点图局部曲线拟合。
- 推断性模型，如置信区间、假设检验、贝叶斯决策分析等。 推断性模型的结论一般是带有概率判断的。
- 预测性模型，如回归分析、判别分析等。主要关心预测的精度。 有些模型基于数学方程，如微分方程、回归模型， 有些模型所需的模型假设很少， 如最近邻回归。
建模的时候， 一般不仅仅考虑一个模型， 而是比较多个模型； 对同一个模型， 也有一些模型参数需要进行调整使得模型性能最优， 比如lasso回归中的光滑性系数的选择， 提升法(boosting)中学习速度和停止时间的选择。 模型比较和参数调优经常使用将数据交叉划分来评估性能的方法。
在研究模型表现时， 有可能会发现一些数据中的不足， 这时有可能需要对数据进行一些新的变换， 或进一步地去收集更多的数据。
在比较模型并筛选出少数较优模型后， 应对模型的性能进行更客观的判断， 如果是预报问题， 这通常需要使用建模时未使用过的数据来评估模型的表现。
37.1.1 训练集和测试集的划分
统计建模（统计学习）问题， 一般都需要将所有数据划分为两个子集：
- 训练集(training set)。 所有的探索性分析、模型试验、参数调整、模型比较都在训练集上进行， 从训练集产生一两个最成功的模型。
- 测试集(testing set)。 仅在训练完成后才将模型应用在测试集上， 评估模型的性能。 训练阶段不允许访问测试集， 一旦因为某种目的在训练阶段使用了测试集， 测试集就失效了。
划分训练集和测试集， 常用如下三种方法：
- 从\(n\)个观测中用随机无放回抽样抽取指定比例的训练集，
余集为测试集。
可以使用基本R函数
sample()获取训练集的下标集合。 - 对于因变量为分类变量且类别比例相差悬殊的情形， 可以按因变量类别进行分层抽样，组合成训练集， 然后余集为测试集。 这也适用于因变量为连续取值但分布重尾的情形。
- 对于沿时间记录的数据， 只能以前面一段时间为训练集， 后面剩余的时间为测试集， 不能选择中间的时间为测试集。
R扩展包rsample提供了对数据框划分训练集和测试集的函数initial_split。
为了用简单随机抽样划分，
用法如：
其中prop指定训练集比例，
da是原始数据框。
为了提取训练集和测试集，
用法如
为了按因变量y的不同类别进行分层抽样产生训练集，
用法如：
这也适用于y是连续取值的情形，
将自动按y的四分位数将y取值范围分成4份进行分层抽样。
对于数据框的各行是按时间顺序记录的情形， 划分训练集和测试集可用如：
有些数据， 如纵向数据(longitudinal data)， 每个受试者（观测单元）有多次观测， 这些观测不独立， 需要以受试者为单位进行随机抽取， 而不是以观测为单位进行抽取。
37.2 用parsnip建模
不同的R建模扩展包使用不同的访问接口， parsnip包则可以将这些扩展包作为计算引擎， 给用户提供一个统一的访问界面。
37.2.1 Ames数据例子
用modeldata扩展包的ames数据集为例。
## Rows: 2,930
## Columns: 74
## $ MS_SubClass <fct> One_Story_1946_and_Newer_All_Styles, One_Story_1946…
## $ MS_Zoning <fct> Residential_Low_Density, Residential_High_Density, …
## $ Lot_Frontage <dbl> 141, 80, 81, 93, 74, 78, 41, 43, 39, 60, 75, 0, 63,…
## $ Lot_Area <int> 31770, 11622, 14267, 11160, 13830, 9978, 4920, 5005…
## $ Street <fct> Pave, Pave, Pave, Pave, Pave, Pave, Pave, Pave, Pav…
## $ Alley <fct> No_Alley_Access, No_Alley_Access, No_Alley_Access, …
## $ Lot_Shape <fct> Slightly_Irregular, Regular, Slightly_Irregular, Re…
## $ Land_Contour <fct> Lvl, Lvl, Lvl, Lvl, Lvl, Lvl, Lvl, HLS, Lvl, Lvl, L…
## $ Utilities <fct> AllPub, AllPub, AllPub, AllPub, AllPub, AllPub, All…
## $ Lot_Config <fct> Corner, Inside, Corner, Corner, Inside, Inside, Ins…
## $ Land_Slope <fct> Gtl, Gtl, Gtl, Gtl, Gtl, Gtl, Gtl, Gtl, Gtl, Gtl, G…
## $ Neighborhood <fct> North_Ames, North_Ames, North_Ames, North_Ames, Gil…
## $ Condition_1 <fct> Norm, Feedr, Norm, Norm, Norm, Norm, Norm, Norm, No…
## $ Condition_2 <fct> Norm, Norm, Norm, Norm, Norm, Norm, Norm, Norm, Nor…
## $ Bldg_Type <fct> OneFam, OneFam, OneFam, OneFam, OneFam, OneFam, Twn…
## $ House_Style <fct> One_Story, One_Story, One_Story, One_Story, Two_Sto…
## $ Overall_Cond <fct> Average, Above_Average, Above_Average, Average, Ave…
## $ Year_Built <int> 1960, 1961, 1958, 1968, 1997, 1998, 2001, 1992, 199…
## $ Year_Remod_Add <int> 1960, 1961, 1958, 1968, 1998, 1998, 2001, 1992, 199…
## $ Roof_Style <fct> Hip, Gable, Hip, Hip, Gable, Gable, Gable, Gable, G…
## $ Roof_Matl <fct> CompShg, CompShg, CompShg, CompShg, CompShg, CompSh…
## $ Exterior_1st <fct> BrkFace, VinylSd, Wd Sdng, BrkFace, VinylSd, VinylS…
## $ Exterior_2nd <fct> Plywood, VinylSd, Wd Sdng, BrkFace, VinylSd, VinylS…
## $ Mas_Vnr_Type <fct> Stone, None, BrkFace, None, None, BrkFace, None, No…
## $ Mas_Vnr_Area <dbl> 112, 0, 108, 0, 0, 20, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6…
## $ Exter_Cond <fct> Typical, Typical, Typical, Typical, Typical, Typica…
## $ Foundation <fct> CBlock, CBlock, CBlock, CBlock, PConc, PConc, PConc…
## $ Bsmt_Cond <fct> Good, Typical, Typical, Typical, Typical, Typical, …
## $ Bsmt_Exposure <fct> Gd, No, No, No, No, No, Mn, No, No, No, No, No, No,…
## $ BsmtFin_Type_1 <fct> BLQ, Rec, ALQ, ALQ, GLQ, GLQ, GLQ, ALQ, GLQ, Unf, U…
## $ BsmtFin_SF_1 <dbl> 2, 6, 1, 1, 3, 3, 3, 1, 3, 7, 7, 1, 7, 3, 3, 1, 3, …
## $ BsmtFin_Type_2 <fct> Unf, LwQ, Unf, Unf, Unf, Unf, Unf, Unf, Unf, Unf, U…
## $ BsmtFin_SF_2 <dbl> 0, 144, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1120, 0…
## $ Bsmt_Unf_SF <dbl> 441, 270, 406, 1045, 137, 324, 722, 1017, 415, 994,…
## $ Total_Bsmt_SF <dbl> 1080, 882, 1329, 2110, 928, 926, 1338, 1280, 1595, …
## $ Heating <fct> GasA, GasA, GasA, GasA, GasA, GasA, GasA, GasA, Gas…
## $ Heating_QC <fct> Fair, Typical, Typical, Excellent, Good, Excellent,…
## $ Central_Air <fct> Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, Y, …
## $ Electrical <fct> SBrkr, SBrkr, SBrkr, SBrkr, SBrkr, SBrkr, SBrkr, SB…
## $ First_Flr_SF <int> 1656, 896, 1329, 2110, 928, 926, 1338, 1280, 1616, …
## $ Second_Flr_SF <int> 0, 0, 0, 0, 701, 678, 0, 0, 0, 776, 892, 0, 676, 0,…
## $ Gr_Liv_Area <int> 1656, 896, 1329, 2110, 1629, 1604, 1338, 1280, 1616…
## $ Bsmt_Full_Bath <dbl> 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0, …
## $ Bsmt_Half_Bath <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Full_Bath <int> 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 3, 2, …
## $ Half_Bath <int> 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, …
## $ Bedroom_AbvGr <int> 3, 2, 3, 3, 3, 3, 2, 2, 2, 3, 3, 3, 3, 2, 1, 4, 4, …
## $ Kitchen_AbvGr <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ TotRms_AbvGrd <int> 7, 5, 6, 8, 6, 7, 6, 5, 5, 7, 7, 6, 7, 5, 4, 12, 8,…
## $ Functional <fct> Typ, Typ, Typ, Typ, Typ, Typ, Typ, Typ, Typ, Typ, T…
## $ Fireplaces <int> 2, 0, 0, 2, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, …
## $ Garage_Type <fct> Attchd, Attchd, Attchd, Attchd, Attchd, Attchd, Att…
## $ Garage_Finish <fct> Fin, Unf, Unf, Fin, Fin, Fin, Fin, RFn, RFn, Fin, F…
## $ Garage_Cars <dbl> 2, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 2, …
## $ Garage_Area <dbl> 528, 730, 312, 522, 482, 470, 582, 506, 608, 442, 4…
## $ Garage_Cond <fct> Typical, Typical, Typical, Typical, Typical, Typica…
## $ Paved_Drive <fct> Partial_Pavement, Paved, Paved, Paved, Paved, Paved…
## $ Wood_Deck_SF <int> 210, 140, 393, 0, 212, 360, 0, 0, 237, 140, 157, 48…
## $ Open_Porch_SF <int> 62, 0, 36, 0, 34, 36, 0, 82, 152, 60, 84, 21, 75, 0…
## $ Enclosed_Porch <int> 0, 0, 0, 0, 0, 0, 170, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ Three_season_porch <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Screen_Porch <int> 0, 120, 0, 0, 0, 0, 0, 144, 0, 0, 0, 0, 0, 0, 140, …
## $ Pool_Area <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Pool_QC <fct> No_Pool, No_Pool, No_Pool, No_Pool, No_Pool, No_Poo…
## $ Fence <fct> No_Fence, Minimum_Privacy, No_Fence, No_Fence, Mini…
## $ Misc_Feature <fct> None, None, Gar2, None, None, None, None, None, Non…
## $ Misc_Val <int> 0, 0, 12500, 0, 0, 0, 0, 0, 0, 0, 0, 500, 0, 0, 0, …
## $ Mo_Sold <int> 5, 6, 6, 4, 3, 6, 4, 1, 3, 6, 4, 3, 5, 2, 6, 6, 6, …
## $ Year_Sold <int> 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 201…
## $ Sale_Type <fct> WD , WD , WD , WD , WD , WD , WD , WD , WD , WD , W…
## $ Sale_Condition <fct> Normal, Normal, Normal, Normal, Normal, Normal, Nor…
## $ Sale_Price <int> 215000, 105000, 172000, 244000, 189900, 195500, 213…
## $ Longitude <dbl> -93.61975, -93.61976, -93.61939, -93.61732, -93.638…
## $ Latitude <dbl> 42.05403, 42.05301, 42.05266, 42.05125, 42.06090, 4…后面将因变量变成房屋售价的常用对数值。
划分训练集和测试集：
set.seed(101)
ames_split <- initial_split(ames, prop = 0.80)
ames_train <- training(ames_split)
ames_test <- testing(ames_split)
dim(ames_train)## [1] 2344 7437.2.2 线性回归示例
parsnip包用类似的函数名和语法指定要使用的模型与计算引擎。 这会指定如下的选择：
- 模型的数学背景，如线性回归、随机森林、KNN等；
- 模型的计算引擎，这一般是规定要使用的R扩展包；
- 模型是回归还判别问题。 如果某个计算引擎只能处理其中一种问题， 则指定了计算引擎也选择了要处理的问题。
这些选择并不与数据集和变量绑定。
例如，
选择最小二乘的线性回归模型，
使用基本R的lm()函数作为计算引擎：
其中的tidymodels_prefer()用于载入的多个扩展包的名字冲突的解决，
可以用tidymodels默认的一些冲突解决方案。
下面将数据、变量与模型联系在一起， 在训练集上产生模型拟合结果：
## parsnip model object
##
##
## Call:
## stats::lm(formula = Sale_Price ~ Longitude + Latitude, data = data)
##
## Coefficients:
## (Intercept) Longitude Latitude
## -301.337 -1.986 2.869fit()函数的一个变种是fit_xy()，
适用于自变量保存为矩阵或数据框，
因变量保存为向量的情形，如：
lm_xy_fit <-
lm_model |>
fit_xy(
x = ames_train |> select(Longitude, Latitude),
y = ames_train |> pull(Sale_Price)
)
lm_xy_fit## parsnip model object
##
##
## Call:
## stats::lm(formula = ..y ~ ., data = data)
##
## Coefficients:
## (Intercept) Longitude Latitude
## -301.337 -1.986 2.869fit()和fit_xy()的区别是，
fit_xy()不对输入数据进行变换，
直接传递给计算引擎，
而fit()在计算引擎要求输入矩阵时，
会对输入数据进行变换，
比如将类别变量转换为哑变量。
37.2.3 随机森林示例
parsnip对于同一模型的不同计算引擎提供了模型参数的统一设置格式。
比如，
随机森林算法，
需要设置每次随机选择的自变量个数、合并的树的个数、每棵树的大小，
不同的随机森林扩展包如ranger, randomForest, sparklyr对这些设置使用不同的参数名称。
parsnip的rand_forest()函数则可以用统一的格式指定这些参数，
允许使用不同的引擎。
如：
这里min_n代表每个叶结点包含观测个数的下限，
低于此下限不能继续拆分，
trees代表树的数目，
mtry代表每次分叉可以从多少个自变量中随机选择。
可以用translate()函数查看parsnip的参数设置如何映射到计算引擎包本来的设置，
如：
## Random Forest Model Specification (regression)
##
## Main Arguments:
## trees = 1000
## min_n = 5
##
## Computational engine: ranger
##
## Model fit template:
## ranger::ranger(x = missing_arg(), y = missing_arg(), weights = missing_arg(),
## num.trees = 1000, min.node.size = min_rows(~5, x), num.threads = 1,
## verbose = FALSE, seed = sample.int(10^5, 1))可见在ranger包中trees变成了num.trees，
min_n变成了min.node.size。
有一些参数是各个引擎特有的，
这在parsnip中会自动设置一些合适的值，
也可以在set_engine()函数中添加，
如set_engine("ranger", verbose = TRUE)。
37.2.4 提取拟合结果
extract_fit_engine()从fit()和fit_xy()的结果中提取所使用引擎输出的拟合结果，
然后可以利用基本R本身的模型信息提取函数如coef(), vcov()等提取结果内容，如：
## (Intercept) Longitude Latitude
## -301.337454 -1.986050 2.868594## (Intercept) Longitude Latitude
## (Intercept) 203.734986 1.5807903348 -1.3252325337
## Longitude 1.580790 0.0168109489 -0.0001562159
## Latitude -1.325233 -0.0001562159 0.0311793054##
## Call:
## stats::lm(formula = Sale_Price ~ Longitude + Latitude, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.02459 -0.09740 -0.01793 0.10037 0.57726
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -301.3375 14.2736 -21.11 <2e-16 ***
## Longitude -1.9861 0.1297 -15.32 <2e-16 ***
## Latitude 2.8686 0.1766 16.25 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1597 on 2341 degrees of freedom
## Multiple R-squared: 0.1746, Adjusted R-squared: 0.1739
## F-statistic: 247.6 on 2 and 2341 DF, p-value: < 2.2e-16## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -301.337454 14.2735765 -21.11156 1.049258e-90
## Longitude -1.986050 0.1296570 -15.31772 1.514774e-50
## Latitude 2.868594 0.1765766 16.24561 2.569212e-56这些信息提取函数有各自的用处， 缺点是返回值类型很不统一， 而且不同引擎对信息提取函数的支持也不相同。 broom包的tidy函数可以将模型拟合结果统一转化为数据框形式， 方便了后续的利用。 如：
## # A tibble: 3 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -301. 14.3 -21.1 1.05e-90
## 2 Longitude -1.99 0.130 -15.3 1.51e-50
## 3 Latitude 2.87 0.177 16.2 2.57e-5637.2.5 预测
parsnip实现的predict()函数，
总是以模型拟合结果为输入，
用new_data输入要预测的数据，
输出结果十分规范：
- 结果总是一个tibble数据框；
- 结果的行数总是等于输入数据长度， 如果输入数据有缺失值， 则对应结果取缺失值；
- 数据框变量名总是有明显含义的。
这大大减轻了用户的学习成本， 因为许多机器学习扩展包提供的预测结果是非标准的， 需要用户自己想办法提取需要的预测内容。
输出结果数据框包含的变量规定如下：
- 若因变量为数值型，
要做点预测，
选
type = "numeric"， 结果变量为.pred； - 因变量为类别，
要预测类别，
选
type = "class"， 结果变量为.pred_class； - 若因变量为类别，
要计算每个类别的后验概率，
选
type = "prob"， 多个结果变量为.pred_{类别水平名}； - 若因变量为数值型，
要计算期望值的置信区间，
选
type = "conf_int"， 结果变量为.pred_lower和.pred_upper； - 若因变量为数值型，
要计算个别值的预测区间，
选
type = "pred_int"， 结果变量也是.pred_lower和.pred_upper。
比如， 计算线性回归模型在部分测试集上的预测值：
## # A tibble: 5 × 1
## .pred
## <dbl>
## 1 5.29
## 2 5.29
## 3 5.24
## 4 5.24
## 5 5.23结果是一个数据框，
其中用变量.pred保存因变量的点预测值，
各行与输入数据的各行次序一一对应。
可以用dplyr::bind_cols()将预测结果与输入的预测用数据进行横向合并，
如：
## # A tibble: 5 × 1
## .pred
## <dbl>
## 1 5.29
## 2 5.29
## 3 5.24
## 4 5.24
## 5 5.23ames_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: 5 × 4
## Sale_Price .pred .pred_lower .pred_upper
## <dbl> <dbl> <dbl> <dbl>
## 1 5.28 5.29 4.98 5.60
## 2 5.29 5.29 4.98 5.60
## 3 5.26 5.24 4.93 5.56
## 4 4.94 5.24 4.93 5.55
## 5 5.11 5.23 4.92 5.55这样统一的预测界面的好处是， 对同一个问题， 使用不同的统计方法和计算引擎， 输出结果格式是相同的， 只要在指定模型时选择不同的统计方法和计算引擎， 后续的预报代码不变。 例如， 用树回归作预测：
tree_model <-
decision_tree(min_n = 2) |>
set_engine("rpart") |>
set_mode("regression")
tree_fit <- tree_model |>
fit(
Sale_Price ~ Longitude + Latitude,
data = ames_train )
ames_test_sm |>
select(Sale_Price) |>
bind_cols(predict(
tree_fit, ames_test_sm ))## # A tibble: 5 × 2
## Sale_Price .pred
## <dbl> <dbl>
## 1 5.28 5.31
## 2 5.29 5.31
## 3 5.26 5.15
## 4 4.94 5.15
## 5 5.11 5.1537.3 用workflow管理建模流程
一个统计建模（统计学习）过程， 不仅仅是指定模型并估计参数的问题， 还包括自变量选择、缺失值插补、变量变换、判别阈值选择等问题， 称整个过程为一个工作流程(workflow)。 其它同类软件也有类似概念， 比如Python和Spark中有流水线(pipeline)的概念。 工作流程将基于训练集产生最优模型的过程完整地封装在一起管理， 也有利于管理研究项目。
37.3.1 工作流程基础
工作流程用来管理建模估计、建模前预处理、建模后的调整。
使用线性回归建模的简单例子如下。 工作流程至少需要指定一个模型。
然后可以产生初始的工作流程：
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: None
## Model: linear_reg()
##
## ── Model ───────────────────────────────────────────────────────────────────────
## Linear Regression Model Specification (regression)
##
## Computational engine: lm这里显示没有预处理步骤。 最简单的一种预处理步骤， 是指定一个回归公式：
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: linear_reg()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Sale_Price ~ Longitude + Latitude
##
## ── Model ───────────────────────────────────────────────────────────────────────
## Linear Regression Model Specification (regression)
##
## Computational engine: lm上面的lm_wflow是一个workflow对象，
这种对象支持用fit()函数拟合模型：
## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: linear_reg()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Sale_Price ~ Longitude + Latitude
##
## ── Model ───────────────────────────────────────────────────────────────────────
##
## Call:
## stats::lm(formula = ..y ~ ., data = data)
##
## Coefficients:
## (Intercept) Longitude Latitude
## -301.337 -1.986 2.869注意lm_fit仍是一个workflow对象，
但包含了拟合结果，
也可以当做拟合结果使用，
比如用predict()产生预测，
用broom::tidy()产生参数估计数据框：
## # A tibble: 3 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -301. 14.3 -21.1 1.05e-90
## 2 Longitude -1.99 0.130 -15.3 1.51e-50
## 3 Latitude 2.87 0.177 16.2 2.57e-56## # A tibble: 5 × 1
## .pred
## <dbl>
## 1 5.29
## 2 5.29
## 3 5.24
## 4 5.24
## 5 5.23工作流程对象， 包括已经拟合的工作流程对象， 其中的模型和预处理内容都可以删除或者更改，例如：
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: linear_reg()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Sale_Price ~ Longitude
##
## ── Model ───────────────────────────────────────────────────────────────────────
## Linear Regression Model Specification (regression)
##
## Computational engine: lm注意修改后的对象不包含拟合结果， 需要重新拟合。
37.3.2 指定变量
除了可以用add_formula()指定因变量和自变量，
还可以用add_variables()指定因变量和自变量，
指定变量集合的格式可以用tidyverse系统的tidyselect选择变量方式（参见23.8）。
如：
lm_wflow <- workflow() |>
add_model(lm_model) |>
add_variables(
outcome = Sale_Price,
predictors = c(Longtitude, Latitude) )选择自变量的例子还有：
如果指定的自变量中还包含了因变量， 因变量会自动忽略， 这样， 可以用如下方式指定所有其它变量作为自变量：
在自变量个数很多时，
tidyselect的指定变量方法可以更灵活地指定自变量集合，
而公式界面只能一一列举自变量，
或者选择使用所有除因变量以外的自变量。
但是，
add_variables()仅将选中的变量保存在数据框中传送给计算引擎，
不进行变换。
如果某些计算引擎需要用户自己将取属性值的自变量转换为哑变量，
则add_variables()不适用，
需要使用食谱(recipe)或者公式来指定因变量和自变量。
食谱是tidymodels提供的更强大的预处理功能。
37.3.3 工作流程中公式的作用
工作流程中用add_formula()添加公式，
可以用来指定因变量、自变量集合，
并在计算引擎需要时进行适当的变换，
比如将取属性值的自变量转换为哑变量，
将原始变量进行某些函数变换（如log(x)）或计算（如I((x-100)^2)），
或产生交叉项，等等。
但是，
不同的扩展包和统计学习方法对变换的要求不同，如：
- 基于树的方法对取属性值的自变量并不转换为哑变量；
- 某些特殊统计模型在公式中可能使用特殊的函数名说明变量用途，
比如生存分析模型可以用
strata(site)说明site变量用于分层， 而不是当做一个普通的自变量， 所以不会估计它的分类别的截距项； - 某些特殊统计模型推广了公式，
比如在线性混合效应模型中，
用
y ~ x | id指定id作为独立观测单位的标识， 同一id值的多个观测之间是相关的。
工作流程如何利用添加的公式呢？ 工作流程对公式的解释依赖于具体的模型也依赖于计算引擎（所利用的扩展包）， 所以工作流程力图模拟模型需要的数据变换， 如果无法实现， 就不对公式作任何改变。 下面举例说明。
37.3.3.1 树模型的公式转换
parsnip包对于树模型能够理解计算引擎对公式的不同理解， 比如， ranger和randomForest包要求不要将取属性值的变量转换为哑变量， 而xgboost包要求进行转换， 另一个C5.0引擎则不需要转换。 parsnip对于能支持的这些计算引擎（调用的扩展包）都可以将公式按照计算引擎的特殊要求进行适当处理。
37.3.3.2 公式内的特殊语法
有些特殊模型使用特殊的公式语法， 例如lme4包建立线性混合模型， 程序例如：
这时，
对每个受试者会拟合年龄项的不同截距项和不同斜率项。
基本R不支持这样的公式，比如model.matrix()函数不能基于这样的公式产生线性模型设计阵。
在工作流程中解决这个问题的办法在add_model()中同时指定对应的公式，
而用到的变量则使用add_variables()原样传递，
如：
## Warning: package 'multilevelmod' was built under R version 4.2.3data(Orthodont, package = "nlme")
multilevel_spec <- linear_reg() |>
set_engine("lmer")
multilevel_workflow <- workflow() |>
add_variables(
outcome = distance,
predictors = c(Sex, age, Subject) ) |>
add_model(
multilevel_spec,
formula = distance ~ Sex + (age | Subject) )
multilevel_fit <- fit(
multilevel_workflow,
data = Orthodont )
multilevel_fit## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Variables
## Model: linear_reg()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Outcomes: distance
## Predictors: c(Sex, age, Subject)
##
## ── Model ───────────────────────────────────────────────────────────────────────
## Linear mixed model fit by REML ['lmerMod']
## Formula: distance ~ Sex + (age | Subject)
## Data: data
## REML criterion at convergence: 471.1635
## Random effects:
## Groups Name Std.Dev. Corr
## Subject (Intercept) 7.3912
## age 0.6943 -0.97
## Residual 1.3100
## Number of obs: 108, groups: Subject, 27
## Fixed Effects:
## (Intercept) SexFemale
## 24.517 -2.14537.3.3.3 公式内的特殊函数
上述办法也能解决公式中有特定扩展包支持的特殊公式的问题。
比如，
在生存分析中，
可以用strata(site)表示site用来分层，
工作流程程序例如：
## Warning: package 'censored' was built under R version 4.2.3## Loading required package: survivalparametric_spec <- survival_reg()
parametric_workflow <- workflow() |>
add_variables(
outcome = c(fustat, futime),
predictors = c(age, rx) ) |>
add_model(
parametric_spec,
formula = Surv(futime, fustat) ~ age + strata(rx) )
parametric_fit <- fit(
parametric_workflow,
data = ovarian )
parametric_fit## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Variables
## Model: survival_reg()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Outcomes: c(fustat, futime)
## Predictors: c(age, rx)
##
## ── Model ───────────────────────────────────────────────────────────────────────
## Call:
## survival::survreg(formula = Surv(futime, fustat) ~ age + strata(rx),
## data = data, model = TRUE)
##
## Coefficients:
## (Intercept) age
## 12.8734120 -0.1033569
##
## Scale:
## rx=1 rx=2
## 0.7695509 0.4703602
##
## Loglik(model)= -89.4 Loglik(intercept only)= -97.1
## Chisq= 15.36 on 1 degrees of freedom, p= 8.88e-05
## n= 2637.3.4 批量建立工作流程
统计建模的过程， 一般都需要比较多个模型， 有时需要比较不同的变量子集。 workflowsets包可以指定预处理器列表和模型列表， 然后生成一批工作流程， 最后设法比较拟合结果。
以Ames数据为例， 考虑使用待售房屋方位的不同方式， 产生多个工作流程比较效果。 首先产生不同的使用房屋方位的列表：
location <- list(
longitude = Sale_Price ~ Longitude,
latitude = Sale_Price ~ Latitude,
coords = Sale_Price ~ Longitude + Latitude,
neighborhood = Sale_Price ~ Neighborhood )这些方式， 可以与不同模型交叉产生多个工作流程。 下面的例子仅使用了线性模型：
library(workflowsets)
location_models <- workflow_set(
preproc = location,
models = list(lm = lm_model) )
location_models## # A workflow set/tibble: 4 × 4
## wflow_id info option result
## <chr> <list> <list> <list>
## 1 longitude_lm <tibble [1 × 4]> <opts[0]> <list [0]>
## 2 latitude_lm <tibble [1 × 4]> <opts[0]> <list [0]>
## 3 coords_lm <tibble [1 × 4]> <opts[0]> <list [0]>
## 4 neighborhood_lm <tibble [1 × 4]> <opts[0]> <list [0]>可见，
workflow_sets()的结果是一个tibble数据框，
每一行是一个工作流程。
info列示例：
## # A tibble: 1 × 4
## workflow preproc model comment
## <list> <chr> <chr> <chr>
## 1 <workflow> formula linear_reg ""option列示例：
## an empty container for options提取某一行的工作流程示例：
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: linear_reg()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Sale_Price ~ Longitude + Latitude
##
## ── Model ───────────────────────────────────────────────────────────────────────
## Linear Regression Model Specification (regression)
##
## Computational engine: lm为了比较工作流程集合中的不同模型，
一般会使用交叉验证等重抽样方法，
重抽样设置与比较结果会保存在工作流程集合数据框的result列中。
这里先不进行重抽样比较，
先在训练数据上对每个模型进行拟合。
这可以用purrr包的map()函数，
以及dplyr包的mutate在工作流程集合数据框添加拟合结果列fit：
location_models <- location_models |>
mutate(
fit = map(
info,
~ fit(
.x$workflow[[1]],
ames_train)))
location_models## # A workflow set/tibble: 4 × 5
## wflow_id info option result fit
## <chr> <list> <list> <list> <list>
## 1 longitude_lm <tibble [1 × 4]> <opts[0]> <list [0]> <workflow>
## 2 latitude_lm <tibble [1 × 4]> <opts[0]> <list [0]> <workflow>
## 3 coords_lm <tibble [1 × 4]> <opts[0]> <list [0]> <workflow>
## 4 neighborhood_lm <tibble [1 × 4]> <opts[0]> <list [0]> <workflow>在上面的程序中，
mutate()输入location_models中的info列，
这是一个列表，
列表的每个元素是一个工作流程的信息，
每个工作流程信息是一个仅有一行的tibble数据框，
其中的workflow列是一个仅有一行的列表，
元素为工作流程，
但必须用.x$workflow[[1]]这样的格式取出，
而不是用.x$workflow取出。
.x是map()函数中用公式语法定义处理函数时的自变量名。
第一个工作流程的拟合结果：
## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: linear_reg()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Sale_Price ~ Longitude
##
## ── Model ───────────────────────────────────────────────────────────────────────
##
## Call:
## stats::lm(formula = ..y ~ ., data = data)
##
## Coefficients:
## (Intercept) Longitude
## -179.412 -1.972虽然使用mutate()和map()批量拟合模型很方便，
但是后面会给出更有用的批量拟合、比较模型的办法。
37.3.5 使用测试集
在一般的统计学习任务中， 测试集不应该用来比较不同模型， 而是在选择了最优模型后， 用来给出其预测性能的评判指标。 比较不同模型应该仅使用训练集。
现在假设已经有了最优模型，
比如是工作流程lm_wflow的结果，
则可以使用last_fit()在测试集上进行性能测试，
这时，
模型估计利用了整个训练集，
在测试集上测试其性能：
lm_model <- linear_reg() |>
set_engine("lm")
lm_wflow <- workflow() |>
add_model(lm_model) |>
add_formula(
Sale_Price ~ Longitude + Latitude )
final_lm_res <- last_fit(lm_wflow, ames_split)
final_lm_res## # Resampling results
## # Manual resampling
## # A tibble: 1 × 6
## splits id .metrics .notes .predictions .workflow
## <list> <chr> <list> <list> <list> <list>
## 1 <split [2344/586]> train/test split <tibble> <tibble> <tibble> <workflow>提取所使用的工作流程：
## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: linear_reg()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Sale_Price ~ Longitude + Latitude
##
## ── Model ───────────────────────────────────────────────────────────────────────
##
## Call:
## stats::lm(formula = ..y ~ ., data = data)
##
## Coefficients:
## (Intercept) Longitude Latitude
## -301.337 -1.986 2.869提取测试性能指标：
## # A tibble: 2 × 4
## .metric .estimator .estimate .config
## <chr> <chr> <dbl> <chr>
## 1 rmse standard 0.167 Preprocessor1_Model1
## 2 rsq standard 0.165 Preprocessor1_Model1提取在测试集上的预测结果， 仅显示前几项：
## # A tibble: 5 × 5
## id .pred .row Sale_Price .config
## <chr> <dbl> <int> <dbl> <chr>
## 1 train/test split 5.29 5 5.28 Preprocessor1_Model1
## 2 train/test split 5.29 6 5.29 Preprocessor1_Model1
## 3 train/test split 5.24 29 5.26 Preprocessor1_Model1
## 4 train/test split 5.24 32 4.94 Preprocessor1_Model1
## 5 train/test split 5.23 33 5.11 Preprocessor1_Model137.4 用recipe进行特征工程
训练数据集中的原始变量并不一定是最适合用来预报因变量的， 对原始变量进行一些变换、重新编码， 或者利用多个变量计算有意义的新变量， 都可能会改进模型的效果。 这样的做法称为“特征工程”(feature engineering)。 特征工程还可以筛选有作用的自变量， 进行降维， 插补缺失值， 将有偏分布转换成对称分布， 等等。 有些模型依赖于简单的欧式距离， 这时各个自变量应该进行标准化， 否则在计算距离时各个变量的影响会有很大的差异。 不同的模型对自变量形式要求也不同， 比如， 基于树的方法不要求分类变量转换为哑变量， 而回归方法则一般要把分类变量转换为哑变量表示（机器学习中称为独热(one-hot)编码）。
recipe包可以把多步特征工程和预处理操作整合为一个对象， 可以将其应用于多个数据集。 这是tidymodels扩展包组合的核心成员之一。
37.4.1 Ames数据集示例
已ames数据集为例， 进行简单的特征工程示例。 仅考虑如下几个自变量：
neighborhood, 分类变量，训练集中共有29个不同取值；Gr_Liv_Area，数值型变量，总的达标居住面积；Year_Built，建造年份；Bldg_Type，房屋类型，分类变量，可取：OneFam;TwoFmCon;Duplex;Twnhs;TwnhsE。
因变量Sale_Price已经变换成了常用对数值。
假设使用lm()函数计算线性回归:
该函数会自动将自变量提取出来组成为设计阵，
这会把分类变量(neighborhood和Bldg_Type)都变成哑变量形式，
变量Gr_Liv_Area会按公式要求用常用对数转换。
如果使用回归结果进行预测，
预测所用的新数据也会进行相同的转换。
一个食谱(recipe)用step_*函数指定类似的预处理步骤，
但暂不进行实际的计算。
例如，
下面的程序可以产生类似于上述lm()对数据进行的变量，
但仅规定变换操作，
暂不进行实际变换：
library(tidymodels)
tidymodels_prefer()
simple_ames <-
recipe(
Sale_Price ~ Neighborhood + Gr_Liv_Area
+ Year_Built + Bldg_Type,
data = ames_train) |>
step_log(Gr_Liv_Area, base = 10) |>
step_dummy(all_nominal_predictors())
simple_ames## Recipe
##
## Inputs:
##
## role #variables
## outcome 1
## predictor 4
##
## Operations:
##
## Log transformation on Gr_Liv_Area
## Dummy variables from all_nominal_predictors()上面的程序中，
用recipe添加了因变量和自变量集合（没有指定对自变量的变换），
以及训练数据集。
然后，用管道操作，
给recipe的结果用step_log对变量Gr_Liv_Area指定了作常用对数变换，
用step_dummy指定了将已选定的自变量集合中所有分类变量转换为哑变量形式。
在recipe环境下，
用all_nominal_predictors()表示所有分类自变量变量，
用all_numeric_predictors()表示所有数值型自变量，
用all_predictors()表示所有自变量，
用all_outcomes()表示所有因变量。
step_*函数中可以指定多个变量，
只要原样写变量名（不用撇号保护）并用逗号分隔。
在recipe()中用data =指定训练数据集，
这时仅从数据确定变量集合与变量类型，
但不进行实际的变量转换。
使用食谱处理变量变换， 而不是直接使用R的公式界面， 或者直接指定变量， 有如下优势：
- 如果使用了不同的模型，这些变换可以在不同模型之间共用；
- 食谱能做的变换比R的公式界面功能更强；
- 使用更方便，比如可以使用
all_nominal_predictors()表示所有分类自变量变量， 而公式中则没有这样的方法； - 对数据要进行的相同的预处理可以存放在一个单一的R对象中， 而不是在源程序中多次重复，甚至于在多个源文件中重复。
37.4.2 使用食谱作为预处理
数据预处理是工作流程的一部分， 在工作流程中加入食谱， 是更灵活的数据预处理指定方式。
工作流程只能有一种预处理方式，
不能在使用add_formula()或add_variables()后再使用食谱。
lm_model <- linear_reg() |>
set_engine("lm")
lm_wflow <- workflow() |>
add_model(lm_model) |>
add_recipe(simple_ames)
lm_wflow## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: linear_reg()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 2 Recipe Steps
##
## • step_log()
## • step_dummy()
##
## ── Model ───────────────────────────────────────────────────────────────────────
## Linear Regression Model Specification (regression)
##
## Computational engine: lm上面的食谱（预处理）仅进行了设置，
而没有进行实际的转换。
用fit()函数，
则可以将工作流程中的预处理进行实际计算，
并拟合模型：
在fit()函数中使用工作流程，
就会调用工作流程中的食谱，
进行实际的变量变换。
这些变换仅针对训练集。
用predict()对测试集的前三个观测计算预测值：
## Warning in predict.lm(object = object$fit, newdata = new_data, type =
## "response"): prediction from a rank-deficient fit may be misleading## # A tibble: 3 × 1
## .pred
## <dbl>
## 1 5.28
## 2 5.28
## 3 5.22在用predict()对测试集进行预测时，
也会在测试集上调用食谱进行预处理，
但这些预处理的参数还是使用从训练集得到的参数，
比如，
如果预处理是对变量作标准化，
即减去变量均值后除以变量标准差，
均值和标准差都是在训练集上计算，
用predict()对测试集进行预测时的变量标准化仍使用训练集上计算的变量均值和标准差。
可以用extract_recipe()加选项estimated = TRUE从拟合结果中提取已计算的预处理结果，如：
## Recipe
##
## Inputs:
##
## role #variables
## outcome 1
## predictor 4
##
## Training data contained 2344 data points and no missing data.
##
## Operations:
##
## Log transformation on Gr_Liv_Area [trained]
## Dummy variables from Neighborhood, Bldg_Type [trained]extract_fit_parsnip()从拟合结果中提取与模型和计算引擎有关的方面，
可以获得模型估计结果，
用tidy()函数可以提取整洁形式的模型参数估计：
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -0.6494257 | 0.2249441 | -2.887054 | 0.0039248 |
| Gr_Liv_Area | 0.6160497 | 0.0138745 | 44.401726 | 0.0000000 |
| Year_Built | 0.0019915 | 0.0001134 | 17.554954 | 0.0000000 |
| Neighborhood_College_Creek | 0.0159125 | 0.0081121 | 1.961565 | 0.0499331 |
| Neighborhood_Old_Town | -0.0290021 | 0.0082649 | -3.509056 | 0.0004583 |
37.4.3 典型的预处理
recipe包的step_*函数提供了各种常见的预处理能力，
这里举例说明。
37.4.3.1 将类别变量转码为数值型
step_dummy()可以将指定的类别变量（名义型、属性值）转换为哑变量编码。
对保存成R因子类型的类别变量，
也可以作一些水平的变化，
比如，
用step_unknown()将指定的类别变量的缺失值转换成一个因子水平值，
用step_novel()为指定的类别变量预留一个原来没有但将来可能出现的因子水平值，
step_other()可以统计各水平的频数，
将低频的水平值合并到一个other类别，
等等。
如：
simple_ames <-
recipe(Sale_Price ~ Neighborhood
+ Gr_Liv_Area + Year_Built + Bldg_Type,
data = ames_train) |>
step_log(Gr_Liv_Area, base = 10) |>
step_other(Neighborhood, threshold = 0.01) |>
step_dummy(all_nominal_predictors())其中的step_other()调用就将Neighborhood中总百分比小于\(1\%\)的低频取值合并成一个other类别。
注意：多数的模型需要将类别变量转换为数值型变量， 一般是转换为哑变量（参见32.10）， 但基于树的模型、基于规则的模型、朴素贝叶斯模型不需要转换。
step_dummy()进行标准的哑变量变换，
比如，
Neighborhood有5个类别，
就会转换成4个哑变量，
对应于5个类别中的后面4个类别，
命名格式是“变量名_类别值”。
每个哑变量都取值0或者1，
但每个观测仅有一个哑变量取1，
取值1的哑变量对应于相应的类别，
当4个哑变量全取0时对应于基线类别，
默认是因子的第一个水平。
可以用step_dummy()中的onehot = TRUE选项要求生成与类别个数相同的哑变量。
除了最常见的哑变量编码， recipes也支持其它的类别变量编码方法。
37.4.3.2 交互作用项
在R的公式中用a:b表示a和b的交互作用效应。
如果a和b都是数值型的，
则交互作用效应等于两者的乘积。
如果a是分类的，
b是数值型的，
则交互作用效应是a的每个类别分别计算b的系数（斜率）。
如果a和b都是分类的，
则对a, b的所有组合建立哑变量编码。
recipes中用step_interact()指定交互作用项，如：
simple_ames <-
recipe(Sale_Price ~ Neighborhood + Gr_Liv_Area
+ Year_Built + Bldg_Type,
data = ames_train) |>
step_log(Gr_Liv_Area, base = 10) |>
step_other(Neighborhood, threshold = 0.01) |>
step_dummy(all_nominal_predictors()) |>
# 注意下面的Gr_Liv_Area是已经做完对数变换的结果
step_interact( ~ Gr_Liv_Area:starts_with("Bldg_Type_") )上面的程序中step_interact()建立了数值型的Gr_Liv_Area的对数变换，
与已经做过哑变量编码的Bldg_Type之间的交互作用效应。
由变量a, b产生的交互作用效应变量命名为a_x_b的格式。
如果需要添加多个交互作用效应，
只要在公式的右侧用+号连接。
注意建立交互作用效应时， 如果有类别变量参与， 应该预先进行哑变量编码。 基本R的公式界面没有这样的要求。
37.4.3.3 样条函数项
在基本R的公式中可以用ns(x, df=3)这样的方式增加变量x的样条函数变换作为自变量。
这种变换将一列数值型变量x的值替换成若干个样条函数值的列，
提供灵活的非线性拟合能力。
自由度越大， 越容易产生过度拟合。 不同的自由度的拟合样例：
library(patchwork)
library(splines)
plot_smoother <- function(deg_free) {
ggplot(ames_train, aes(
x = Latitude, y = 10^Sale_Price)) +
geom_point(alpha = .2) +
scale_y_log10() +
geom_smooth(
method = lm,
formula = y ~ ns(x, df = deg_free),
color = "lightblue",
se = FALSE
) +
labs(
title = paste(deg_free, "Spline Terms"),
y = "Sale Price (USD)")
}
( plot_smoother(2) + plot_smoother(5) ) /
( plot_smoother(20) + plot_smoother(100) )
用recipes添加样条变换项的例子：
recipe(Sale_Price ~ Neighborhood + Gr_Liv_Area
+ Year_Built + Bldg_Type + Latitude,
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, deg_free = 20)37.4.3.4 特征提取
预处理经常需要根据旧有变量计算产生对预测因变量更有用的新变量。 比如， 主成分分析(PCA)就是一种消除旧有变量中的相关性并获得压缩结果的方法。 例如， ames数据中与面积大小有关的变量相关性强， 可以进行压缩， 程序如：
其中用了matches()指定变量名符合特定模式的变量子集。
如果各个变量量纲不同，
应该先用step_normalize()将变量转换到相近的取值范围。
除了主成分分析以外， recipes还支持其它一些变换方法。
37.4.3.5 对观测行进行抽样
在分类问题中， 如果两个类别的比例相差很大， 很容易忽略频率较少的类。 为此可以对训练样本的观测进行抽样， 产生新的训练样本：
- 可以减数抽样， 保持少数派观测不变， 对多数派样本仅抽取其中一部分；
- 可以增殖抽样， 保持多数派观测不变， 对少数派进行抽样增加， 增加方式可能是简单地重复， 也可能是产生与少数派表现相同但不同的样本。
- 可以对多数派抽取一部分， 对少数派增加样品。
themis包能够进行这些操作， 比如， 减数抽样可用：
这些对样本行的增删应该仅对训练样本进行，
测试样本和其它保留的样本不应该进行这些操作。
这些函数设置了一个skip参数，
默认为TRUE，
就是规定这样的原则。
还有一些操作样本行的函数，
如step_filter(),
step_sample(),
step_slice(),
step_arrange()。
37.4.4 对新数据应避免的预处理
如果需要对因变量进行变换， 不应该作为预处理步骤规定在食谱中， 这是因为在使用测试集时， 测试集中可以没有因变量存在。
如前所述，
如果对数据行进行了一些操作，
比如减数抽样，
应该仅对训练集进行，
但是测试集上不应该进行这些操作。
这样的行操作函数提供了skip = TRUE默认选项，
就是要求在predict()时忽略行操作预处理步骤。
当然，
fit()函数不会忽略这样的步骤。
如下的行操作函数仅对训练集有效：
step_adasyn(),
step_bsmote(),
step_downsample(),
step_filter(),
step_naomit(),
step_nearmiss(),
step_rose(),
step_sample(),
step_slice(),
step_smote(),
step_smotenc(),
step_tomek(),
step_upsample()。
37.4.5 食谱的整洁表示
可以用tidy()函数将食谱转换为数据框，
如：
ames_rec <-
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,
id = "nei_other") |>
step_dummy(all_nominal_predictors()) |>
step_interact(
~ Gr_Liv_Area:starts_with("Bldg_Type_") ) |>
step_ns(Latitude, Longitude, deg_free = 20)
tidy(ames_rec) |> knitr::kable()| number | operation | type | trained | skip | id |
|---|---|---|---|---|---|
| 1 | step | log | FALSE | FALSE | log_WPwWU |
| 2 | step | other | FALSE | FALSE | nei_other |
| 3 | step | dummy | FALSE | FALSE | dummy_Z4Mwx |
| 4 | step | interact | FALSE | FALSE | interact_pPAiA |
| 5 | step | ns | FALSE | FALSE | ns_LGEm4 |
数据集中的id是自动生成的，
可以在step_*函数中用id =认为指定，
如上面的step_other()调用。
下面利用上述预处理重新拟合：
lm_wflow <-
workflow() |>
add_model(lm_model) |>
add_recipe(ames_rec)
lm_fit <- fit(lm_wflow, ames_train)下面提取拟合过的step_other()预处理步的信息，
这时，
上面指定的id就可以用来从所有预处理步中定位：
estimated_recipe <-
lm_fit %>%
extract_recipe(estimated = TRUE)
tidy(estimated_recipe, id = "nei_other")## # A tibble: 20 × 3
## terms retained id
## <chr> <chr> <chr>
## 1 Neighborhood North_Ames nei_other
## 2 Neighborhood College_Creek nei_other
## 3 Neighborhood Old_Town nei_other
## 4 Neighborhood Edwards nei_other
## 5 Neighborhood Somerset nei_other
## 6 Neighborhood Northridge_Heights nei_other
## 7 Neighborhood Gilbert nei_other
## 8 Neighborhood Sawyer nei_other
## 9 Neighborhood Northwest_Ames nei_other
## 10 Neighborhood Sawyer_West nei_other
## 11 Neighborhood Mitchell nei_other
## 12 Neighborhood Brookside nei_other
## 13 Neighborhood Crawford nei_other
## 14 Neighborhood Iowa_DOT_and_Rail_Road nei_other
## 15 Neighborhood Timberland nei_other
## 16 Neighborhood Northridge nei_other
## 17 Neighborhood Stone_Brook nei_other
## 18 Neighborhood South_and_West_of_Iowa_State_University nei_other
## 19 Neighborhood Clear_Creek nei_other
## 20 Neighborhood Meadow_Village nei_other这个拟合过的预处理步显示结果包含了没有被合并的类别的列表。
tidy()函数除了用id定位某一个预处理步，
也可以用number =指定预处理步的序号。
对食谱调用tidy()函数，
可以选择某个预处理步，
显示与该预处理步有关的信息。
37.4.6 变量的角色
在公式中出现的变量， 基本的角色是自变量(predictor)和因变量(outcome)。 还可以有其它角色。
比如，
ames数据集的每一行可以有一个地址，
这虽然不是自变量也不是因变量，
但可以用来定位感兴趣的观测。
可以用add_role()、remove_role()、update_role()改变变量的角色。
如：
角色名称可以是任何字符串。
一个数据集变量也可以有多重角色，
可以用add_role()在已有角色时额外添加角色。
在数据集中保存非自变量、因变量角色的变量列， 好处是在进行行操作、随机抽样操作时， 可以同步地对这些列进行相同的操作， 如果使用数据集外部的一个变量保存这些信息（如地址）， 就不能自动跟随这些行操作。
37.4.7 更多预处理步骤
recipes和相关的包提供了上百个预处理方法， 参见https://www.tidymodels.org/find/。 还可以自己定义预处理步， 参见https://www.tidymodels.org/learn/develop/recipes/。
37.5 模型效果的评价指标计算
建模需要能够定量地评估模型效果， 用来比较不同的模型、对模型进行参数调优、更好地理解模型。 这里主要使用验证数据集来评估模型效果。 作为示例， 这一节将测试数据集当作验证数据集来使用， 但实际建模时应该使用单独的验证数据集， 或者使用随机抽样方法进行交叉验证。
选择什么样的定量度量来评估模型， 对最后选择的模型有关键的影响。 比如， 在回归分析中， 一般使用最小二乘估计， 这最小化了均方误差， 使得预测误差最小； 如果参数估计时最大化复相关系数平方\(R^2\)， 这实际上是最大化了因变量的真实值和拟合值的相关系数， 一般是不合理的做法。
yardstick是tidymodels扩展包集合中的核心组成包， 提供了各种模型评估度量。
37.5.1 模型评估与统计推断
本章对建模更强调模型的预测能力， 即使建模主要是为了进行统计推断， 其预测能力的评估也可以从另一侧面验证模型的可用性。
常见的推断性统计模型应用往往是不断在不同模型设定之间进行比较检验， 直到不能继续增加或减少解释项。 但是， 要注意这样还主要是一个拟合问题， 不管模型拟合的有多么好， 也不能说明这个模型就是客观规律的良好近似， 这是拟合问题的固有缺陷。 这时， 借助于交叉验证等方法进行预测能力的评估， 可以帮助我们认识模型的可信性。
37.5.2 回归问题的评估
在完成模型拟合后，
可以产生验证集或测试集上的预测结果，
将预测结果保存入数据集，
然后，
yardstick包的评估指标计算函数有统一的界面：
如果函数名为f，
保存预测结果的数据集为data，
保存因变量真实值的列为truth，
则计算评估指标的调用格式为
其中的...参数用来定位数据中的预测值。
下面用ames数据集举例说明，
先拟合模型：
ames_rec <-
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)
lm_model <- linear_reg() |> set_engine("lm")
lm_wflow <-
workflow() |>
add_model(lm_model) |>
add_recipe(ames_rec)
lm_fit <- fit(lm_wflow, ames_train)为了演示，
这里用ames_test当作验证数据集，
正常情况下应该单独生成验证数据集或者使用训练集的随机抽样交叉验证。
在验证数据集上进行预测，
产生因变量的预测值：
ames_test_res <- predict(
lm_fit,
new_data = ames_test |> select(-Sale_Price))
ames_test_res |> slice(1:5)## # A tibble: 5 × 1
## .pred
## <dbl>
## 1 5.28
## 2 5.27
## 3 5.22
## 4 5.08
## 5 5.14回归结果的点预测值，
在保存预测结果的数据框中命名为.pred。
可以将预测结果与验证集的原始数据中的因变量真实值一一对应起来保存：
ames_test_res <- ames_test |>
select(Sale_Price) |>
bind_cols(ames_test_res)
ames_test_res |> slice(1:5)## # A tibble: 5 × 2
## Sale_Price .pred
## <dbl> <dbl>
## 1 5.28 5.28
## 2 5.29 5.27
## 3 5.26 5.22
## 4 4.94 5.08
## 5 5.11 5.14可以浏览这个两列的数据集， 大致地感觉预测精度。 注意因变量真实值和预测值都是对数变换结果， 在评估模型效果时， 一般使用建模时的因变量形式， 而不是转换成原始形式再评估。
可以作散点图显示两者的关系：
ggplot(ames_test_res, aes(
x = Sale_Price, y = .pred)) +
geom_abline(lty = 2) + # 没有指定a, b则默认a=0, b=1
geom_point(alpha = 0.5) +
labs(y = "Predicted Sale Price (log10)",
x = "Sale Price (log10)") +
coord_obs_pred() # 令x、y轴使用相同长度单位
在有了保存因变量预测值和真值的数据集后，
可以用rmse()函数计算根均方误差：
\[
\text{RMSE} = \sqrt{\frac{1}{N} \sum_{i=1}^N
(y_i - \hat y_i)^2} .
\]
## # A tibble: 1 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 rmse standard 0.0828注意.pred不需要用撇号保护。
rmse()的如上调用方法是yardstick包中计算模型评估指标的函数的普遍坐法。
如果要同时计算多个模型评估指标，
因为每个评估指标对应于一个函数，
需要用metric_set()函数将若干个评估指标函数组合成一个评估指标集合函数，
调用这个函数就可以计算多个评估指标。
如：
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.0828
## 2 rsq standard 0.797
## 3 mae standard 0.0561这计算了根均方误差，\(R^2\)，平均绝对误差： \[ \text{MAE} = \frac{1}{N} \sum_{i=1}^N |y_i - \hat y_i| . \] RMSE和MAE与因变量量纲相同， 度量因变量和自变量的取值接近程度； 而\(R^2\)是一个取值在\([0, 1]\)区间的数值， \(R^2\)越大， 表示拟合越好。 \(R^2\)实际是因变量真实值与拟合值的相关系数的平方。
yardstick有意地取消了修正\(R^2\)的计算功能， 因为修正\(R^2\)主要是在同一个训练集上比较模型， 而更好的做法是在不同数据上评估模型预测效果。
下面用随机森林方法对ames数据集进行预测， 并与线性模型结果进行比较：
rf_model <-
rand_forest(trees = 1000) |>
set_engine("ranger") |>
set_mode("regression")
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.0749
## 2 rsq standard 0.840
## 3 mae standard 0.0480线性回归在测试集上的MAE为\(0.055\)， 随机森林的MAE为\(0.051\)， 随机森林略优。 但是这样的评估结果有随机性， 如果重新划分训练集、测试集， 评估结果的比较可能会改变。 更好的模型比较方法是利用再抽样技术， 见37.6。
37.5.3 二分类问题的评估
37.5.3.1 二分类问题的评估指标
二分类问题预测的评估有许多种指标， 各有侧重。
在医学检验中， 将有某种病的人正确检出的概率叫做敏感度(sensitivity)， 将无病的人判为无病的概率叫做特异度(specificity)。 将这样的术语对应到统计假设检验问题： \[ H_0: \text{无病} \longleftrightarrow H_a: \text{有病} \] 第一类错误是\(H_0\)成立但拒绝\(H_0\)的错误， 概率为\(\alpha\)； 第二类错误是\(H_a\)成立但承认\(H_0\)的错误， 概率为\(\beta\)。 则敏感度为\(1-\beta\)（检验功效）， 特异度为\(1-\alpha\)。
图37.1: 医学诊断四种情况表格
称\(\alpha\)为假阳性率(FPR)， 等于1减去特异度； 称敏感度为真阳性率(TPR)， 也称为召回率(recall)。
设图37.1中TN表示真阴性例数， 即无病判为无病个数； FN表示假阴性例数， 即有病判为无病例数； FP表示假阳性例数， 即无病判为有病例数； TP表示真阳性例数， 即有病判为有病例数。 这时，可以估计 \[\begin{aligned} \text{TPR(真阳性率)} =& 1 - \beta = \text{Sensitivity(敏感度)} = \text{Recall(查全率)} = \frac{TP}{TP + FN} ; \\ \text{FPR(假阳性率)} =& \alpha = \frac{FP}{FP + TN} ; \\ \text{TNR(真阴性率)} =& 1-\alpha = \text{Specificity(特异度)} = \frac{TN}{TN + FP}; \\ \text{FNR(假阴性率)} =& \beta = \frac{FN}{FN + TP} . \end{aligned}\]
另外， \[\begin{aligned} \text{Precision(查准率)} =& \text{Positive predictive value} = 1 - \text{FDR} = \frac{TP}{TP + FP}; \\ \text{FDR(错误发现率)} =& \frac{FP}{FP + TP} = 1 - \text{Precision}. \end{aligned}\]
F1指标是查准率与查全率的如下“调和平均”： \[ \frac{1}{\text{F1}} = \frac{1}{2} \left( \frac{1}{\text{Precision}} + \frac{1}{\text{Recall}} \right) . \] F1指标更重视将正例检出， 而忽视对反例的正确检出。
MCC(Mathew相关系数)度量\(2 \times 2\)列联表行、列变量的相关性， 与拟合优度卡方统计量等价，计算公式为： \[ \text{MCC} = \frac{\text{TN} \times \text{TP} - \text{FN} \times \text{FP} }{ \sqrt{ (\text{TP} + \text{FP}) (\text{TP} + \text{FN}) (\text{TN} + \text{FP}) (\text{TN} + \text{FN}) }} . \]
在两类的分类问题中， 将两类分别称为正例、反例， 或者阳性、阴性， 或者事件、非事件。 设自变量\(x\)得到的正例的后验概率为\(f(\boldsymbol x)\)， 任意取分界值\(t\)， 都可以按\(f(\boldsymbol x) > t\)和\(f(\boldsymbol x) < t\)将因变量值（标签）判为正例或者反例。 对每个\(t\)， 都可以在训练集上或者测试集上计算对应的假阳性率（健康人判入有病的比例）和真阳性率（病人判为有病的概率）。 \(t\)值取得越大， 相当于检验水平\(\alpha\)取值越小， 越不容易判为有病， 这时假阳性率低， 真阳性率也低； \(t\)值取得越小， 越容易判为有病， 相当于检验水平\(\alpha\)取值很大， 这时假阳性率很高，真阳性率也很高。 最好是对比较大的分界点\(t\)， 当假阳性率很低时， 真阳性率就很高。 以给定\(t\)后的假阳性率为横坐标， 真阳性率为纵坐标作曲线， 称为ROC(Receiver Operating Curve，受试者工作特征)曲线（例子见37.2）。 曲线越靠近左上角越好。 可以用ROC曲线下方的面积大小度量判别方法的优劣， 面积越大越好。 记这个面积为AUC(Area Under the Curve)。
37.5.3.2 心脏病诊断示例
为了演示二分类问题， 使用Heart数据。 这是心脏病诊断的数据， 因变量AHD为是否有心脏病， 试图用各个自变量预测（分类）。
读入Heart数据集，并去掉有缺失值的观测：
Heart <- read_csv("data/Heart.csv") |>
select(-1) |>
mutate(
ChestPain = factor(ChestPain, levels = c(
"asymptomatic", "nonanginal",
"nontypical", "typical")),
Thal = factor(Thal, levels = c(
"fixed", "normal", "reversable")),
AHD = factor(AHD, levels = c("Yes", "No")) ) |>
drop_na()## New names:
## Rows: 303 Columns: 15
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (3): ChestPain, Thal, AHD dbl (12): ...1, Age, Sex, RestBP, Chol, Fbs, RestECG,
## MaxHR, ExAng, Oldpeak,...
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`## tibble [297 × 14] (S3: tbl_df/tbl/data.frame)
## $ Age : num [1:297] 63 67 67 37 41 56 62 57 63 53 ...
## $ Sex : num [1:297] 1 1 1 1 0 1 0 0 1 1 ...
## $ ChestPain: Factor w/ 4 levels "asymptomatic",..: 4 1 1 2 3 3 1 1 1 1 ...
## $ RestBP : num [1:297] 145 160 120 130 130 120 140 120 130 140 ...
## $ Chol : num [1:297] 233 286 229 250 204 236 268 354 254 203 ...
## $ Fbs : num [1:297] 1 0 0 0 0 0 0 0 0 1 ...
## $ RestECG : num [1:297] 2 2 2 0 2 0 2 0 2 2 ...
## $ MaxHR : num [1:297] 150 108 129 187 172 178 160 163 147 155 ...
## $ ExAng : num [1:297] 0 1 1 0 0 0 0 1 0 1 ...
## $ Oldpeak : num [1:297] 2.3 1.5 2.6 3.5 1.4 0.8 3.6 0.6 1.4 3.1 ...
## $ Slope : num [1:297] 3 2 2 3 1 1 3 1 2 3 ...
## $ Ca : num [1:297] 0 3 2 0 0 0 2 0 1 0 ...
## $ Thal : Factor w/ 3 levels "fixed","normal",..: 1 2 3 2 2 2 2 2 3 3 ...
## $ AHD : Factor w/ 2 levels "Yes","No": 2 1 1 2 2 2 1 2 1 1 ...产生训练集和测试集：
set.seed(101)
split <- initial_split(Heart, prop = 0.60)
heart_train <- training(split)
heart_test <- testing(split)预处理步骤、建立工作流程、在训练集上拟合模型：
heart_rec <-
recipe(AHD ~ .,
data = heart_train) |>
step_dummy(all_nominal_predictors())
lgs_model <- logistic_reg() |> set_engine("glm")
lgs_wflow <-
workflow() |>
add_model(lgs_model) |>
add_recipe(heart_rec)
lgs_fit <- fit(lgs_wflow, heart_train)在测试集上预测， 生成预测类别、每类的后验概率， 并添加真实类别结果：
heart_test_res <- predict(lgs_fit,
new_data = heart_test |> select(-AHD)) |>
bind_cols(predict(lgs_fit,
new_data = heart_test |> select(-AHD),
type = "prob") ) |>
bind_cols(heart_test |> select(AHD))
heart_test_res |>
slice(1:5) |>
knitr::kable()| .pred_class | .pred_Yes | .pred_No | AHD |
|---|---|---|---|
| No | 0.1756546 | 0.8243454 | No |
| No | 0.4230829 | 0.5769171 | No |
| No | 0.0205135 | 0.9794865 | No |
| No | 0.1867585 | 0.8132415 | No |
| No | 0.0125266 | 0.9874734 | No |
显示混淆矩阵：
## Truth
## Prediction Yes No
## Yes 32 4
## No 15 68注意混淆矩阵有两种做法， 这里是机器学习中常用的做法， 真实类别为两列， 预测类别为两行。 相当于如下命令的结果：
## AHD
## .pred_class Yes No
## Yes 32 4
## No 15 68另一种做法是转置显示。
计算准确率， 即预测类别等于真实类别的样例的比例：
## # A tibble: 1 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.840Matthew相关系数：
## # A tibble: 1 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 mcc binary 0.665计算F1指标：
## # A tibble: 1 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 f_meas binary 0.771同时产生多个评估指标的示例：
classification_metrics <- metric_set(
accuracy, mcc, f_meas)
classification_metrics(heart_test_res,
truth = AHD,
estimate = .pred_class)## # A tibble: 3 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.840
## 2 mcc binary 0.665
## 3 f_meas binary 0.771准确率和MCC同时评估两类的正确预测，
而F1指标更强调能够将正例检出。
tidymodels默认第一个水平为正例，
可以用event_level选项指定一个水平作为正例（事件）。
在这一点上R语言的不同函数有不同的规定。
ROC曲线、计算AUC值，
则需要使用预测的后验概率而不是类别。
用roc_curve()计算ROC曲线上点的坐标，如：
## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## ℹ The deprecated feature was likely used in the yardstick package.
## Please report the issue at <]8;;https://github.com/tidymodels/yardstick/issueshttps://github.com/tidymodels/yardstick/issues]8;;>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.| .threshold | specificity | sensitivity |
|---|---|---|
| -Inf | 0.0000000 | 1 |
| 0.0033245 | 0.0000000 | 1 |
| 0.0072132 | 0.0138889 | 1 |
| 0.0081136 | 0.0277778 | 1 |
| 0.0113795 | 0.0416667 | 1 |
用autoplot()从上面的数据集作AOC图：
图37.2: 心脏病诊断建模ROC曲线
计算AUC值：
## # A tibble: 1 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 roc_auc binary 0.90037.5.4 多类别分类问题的评估
用鸢尾花数据作为示例。 这里不拆分训练集和测试集， 用交叉验证方法进行模型评估。 关于交叉验证方法， 见37.6.2.1。
用随机森林作判别， 建立工作流程并拟合：
data(iris)
iris_rec <-
recipe(Species ~ .,
data = iris)
rf_model <- rand_forest(
trees = 100,
min_n = 5 ) |>
set_engine("ranger") |>
set_mode("classification")
rf_wflow <-
workflow() |>
add_model(rf_model) |>
add_recipe(iris_rec)产生交叉验证分组， 用交叉验证样本拟合并评估：
set.seed(101)
iris_folds <- vfold_cv(iris, v = 10)
keep_pred <- control_resamples(
save_pred = TRUE, save_workflow = TRUE)
iris_cv <- rf_wflow |>
fit_resamples(resamples = iris_folds,
control = keep_pred)## Warning: package 'ranger' was built under R version 4.2.3从交叉验证结果中提取预测值、后验概率：
## tibble [150 × 8] (S3: tbl_df/tbl/data.frame)
## $ id : chr [1:150] "Fold01" "Fold01" "Fold01" "Fold01" ...
## $ .pred_setosa : num [1:150] 1 1 1 0 0 0.002 0 0.002 0 0 ...
## $ .pred_versicolor: num [1:150] 0 0 0 0.851 1 ...
## $ .pred_virginica : num [1:150] 0 0 0 0.149 0 ...
## $ .row : int [1:150] 5 9 46 53 83 90 93 99 116 124 ...
## $ .pred_class : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 2 2 2 2 2 3 3 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 2 2 2 2 2 3 3 ...
## $ .config : chr [1:150] "Preprocessor1_Model1" "Preprocessor1_Model1" "Preprocessor1_Model1" "Preprocessor1_Model1" ...显示混淆矩阵：
## Truth
## Prediction setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 47 3
## virginica 0 3 47准确率：
## # A tibble: 1 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy multiclass 0.96MCC指标也可以推广到多类别分类，如：
## # A tibble: 1 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 mcc multiclass 0.94原来适用于二分类的评价指标都可以设法推广到多类别分类问题， 推广方法有（设有\(k\)个类别）：
- 宏平均法，将每个类别与“其它”类别作二分类计算， 将获得的\(k\)个评估指标计算平均值。
- 宏加权法，类似宏平均法， 但是计算平均值时按每类的样例数加权。
- 微平均法，对每个类分别计算对评估指标的贡献， 将这些贡献汇总计算成一个单一指标。
用estimator选项可以指定推广方法。
例如，
sensitivity()在二分类问题中计算敏感度（将正例识别为正例的概率），
用宏平均法推广到多分类：
## # A tibble: 1 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 sensitivity macro 0.96在上例中， 每一个类别都轮流作为正例， 然后将所得的3个敏感度平均。
宏加权法用estimator = "macro_weighted"，
微平均法用estimator = "micro"。
计算AUC的推广指标，使用Hand and Till(2001)文章的方法：
## # A tibble: 1 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 roc_auc hand_till 0.996文献：
- Hand, D, and R Till. 2001. “A Simple Generalisation of the Area Under the ROC Curve for Multiple Class Classification Problems.” Machine Learning 45 (August): 171–86.
37.6 用再抽样技术计算模型评价指标
37.6.1 拟合效果
再抽样技术是指从同一个样本数据集中多次抽样产生多个样本， 这些样本当然不能做到相互独立， 但是也可以用来增大样本容量， boostrap是典型的再抽样技术， 而统计建模、机器学习中的交叉验证(cross validation)也是一种再抽样技术。
将全体样本划分为训练集(training set)和测试集(test set)后， 测试集不能用来训练模型、比较模型， 只能在已经选定模型后计算模型性能评估指标。 所有的模型参数估计、超参数调优、模型比较都只能使用训练集。 如果完全用训练集上的预测效果评估模型性能， 则更复杂的模型会占有绝对优势， 比如， 平面上横坐标不相等的\(10\)个点， 可以用一个9阶多项式完美拟合， 但是这样凑出来的模型对新的数据完全没有意义。 这样的问题称为“过度拟合”问题。 见32.6的示例。 所以， 除了少量模型复杂度很低的模型以外， 不应该使用在训练集上拟合的效果来评估模型性能。
在训练集上， 为了评估模型， 必须使用建模时没有用到的数据计算评估指标。 可以将训练集继续随机拆分， 拆分的两部分称为训练集和“验证集”（validation set）， 在训练集上估计模型， 在验证集上计算模型性能评估指标。 但是， 当数据量比较小的时候， 计算得到的评估指标严重依赖于拆分的随机性， 变化很大。 各种再抽样方法可以克服这个问题。
37.6.2 再抽样方法
再抽样方法一般重复地从训练样本生成多组分析集(analysis set)和评估集(assessment set)对， 在生成的每个分析集上估计模型参数， 然后在相应的评估集上计算模型评估指标， 最后汇总这些评估指标为一个数值。 再抽样仅在训练集上进行， 不允许对测试集作再抽样。 分析集和评估集的作用类似于训练集和测试集， 但分析集和评估集是从训练集中获得， 而且可以用随机抽样方法获得多个分析集、评估集对。 也有人将抽取出的分析集和评估集称为训练集和验证集。
分析集、评估集、测试集两两不相交， 但是不同的分析集、评估集对之间可以有交集。 因为估计模型在分析集上进行， 评估模型在评估集上进行， 评估集与分析集没有共同数据， 所以这样的评估方法比用训练集上的拟合效果评估更为客观。 但是， 因为不同的分析集、评估集对之间有交集， 所以这种方法的评估结果可靠性还是比不上使用完全独立的测试集。
37.6.2.1 交叉验证方法
设训练集有\(n=10000\)个样例。 将每个样例随机地分配一个\(1:10\)的号码， 则可以将训练集分为10“折”（fold）， 就是10个子集， 每个子集大约有1000个样例。 用分配的号码作为样例所属子集（折）编号， 这可以产生10对分析集、评估集：
- 第一对，分析集由编号\(2:10\)的折组成， 评估集为编号\(1\)的折；
- 第二对，分析集由编号\(1, 3:10\)的折组成， 评估集为编号\(2\)的折；
- …………
这种方法称为十折交叉验证(ten-fold cross validation)。
划分折时， 最简单的是简单随机抽样， 给每个观测随机分配一个\(1:k\)的号码（\(k\)是折数）。 也可以使用分层抽样方法， 在某个分类变量的每一类别内部进行均匀分折。
在最后计算评估指标时， 取较多的折数会减小评估指标的估计偏差，但增大标准误差； 取较少的折数会增大评估指标的估计偏差，但能减小标准误差。 常用十折和五折。
用tidymodels产生这样的再抽样方案的程序示例：
## # 10-fold cross-validation
## # A tibble: 10 × 2
## splits id
## <list> <chr>
## 1 <split [2109/235]> Fold01
## 2 <split [2109/235]> Fold02
## 3 <split [2109/235]> Fold03
## 4 <split [2109/235]> Fold04
## 5 <split [2110/234]> Fold05
## 6 <split [2110/234]> Fold06
## 7 <split [2110/234]> Fold07
## 8 <split [2110/234]> Fold08
## 9 <split [2110/234]> Fold09
## 10 <split [2110/234]> Fold10结果是一个rset类型的数据集，
保存为一个两列的tibble数据框，
id列用于区分这些分析集、验证集对，
而splits列是元素为列表的列，
内容是产生的多个分析集、验证集对，
如第一折相应的分析集、验证集对：
## <Analysis/Assess/Total>
## <2109/235/2344>用analysis()取出某一对的分析集，
用validation()取出测试集，如：
## [1] 2109 74## [1] 235 74这说明第一折对应的分析集由2109个样例， 评估集有235个样例。
许多使用交叉验证的方法会自动进行再抽样并计算评估结果， 所以往往不需要人为生成再抽样的多个分析集、评估集对。
除了上述的交叉验证方法以外， 还有一些交叉验证方法的变种。
重复交叉验证方法：
如果仅取一个分析集和一个评估集， 则样本量较小的情况下评估结果会随着随机抽样的不同而变化很大。 即使用了\(k\)折的交叉验证， 评估结果仍可能有较大的随机性。 为此， 可以将这样的操作重复\(r\)次， 则性能评估结果的标准误差可以降低到原来的\(\frac{1}{\sqrt{r}}\)。
程序是在vfolds()函数中添加一个repeats =参数，
如：
留一法：
将每个样例都作为评估集一次， 对应的分析集则是剩余的\(n-1\)个样例， 这样一共有\(n\)个分析、评估集对。 除了样本量特别小的情况以外， 这种方法从计算效率和统计性质上都是不好的， 但是这种方法提示了交叉验证方法的产生。
重复随机模拟生成分析集、评估集：
如果将训练集随机划分为分析集、评估集， 这样计算出来的评估指标有较大的随机性； 但是， 如果将这样的步骤重复多次并汇总评估结果， 则评估结果比较稳定。
与通常的交叉验证的区别是， 交叉验证的不同折对应的分析集是互不相交的， 但这样重复随机模拟生成的评估集之间是有交集的。
生成重复随机模拟分析、评估集的程序示例：
## # Monte Carlo cross-validation (0.9/0.1) with 10 resamples
## # A tibble: 10 × 2
## splits id
## <list> <chr>
## 1 <split [2109/235]> Resample01
## 2 <split [2109/235]> Resample02
## 3 <split [2109/235]> Resample03
## 4 <split [2109/235]> Resample04
## 5 <split [2109/235]> Resample05
## 6 <split [2109/235]> Resample06
## 7 <split [2109/235]> Resample07
## 8 <split [2109/235]> Resample08
## 9 <split [2109/235]> Resample09
## 10 <split [2109/235]> Resample1037.6.2.2 仅取一个验证集
当数据量较大时， 可以将原始数据随机划分为训练集、验证集、测试集三个部分， 在训练集上估计模型参数， 在验证集上比较模型、寻找最优超参数， 用测试集给出最终的模型性能指标。
用validation_split()可以从原来的训练集划分出训练集和验证集，
这个结果与vfold_cv()结果类型相同，
但仅有一个观测行，如：
## # Validation Set Split (0.75/0.25)
## # A tibble: 1 × 2
## splits id
## <list> <chr>
## 1 <split [1758/586]> validation划分的训练集和验证集的大小：
## [1] 1758 74## [1] 586 7437.6.2.3 自助法
自助法(bootstraping)在统计中用来评估统计量分布。 训练集的一个bootstrap样本是从训练集的\(n\)个样例中等概有放回抽取\(n\)次形成的一个数据集。 显然， 这样的数据集中有许多样例是重复的。 每个样例被选入bootstrap样本的概率是\(63.2\%\)， 所以与每个bootstrap样本相对应， 有大约\(36.8\%\)的样例没有被选入该bootstrap样本， 称这些落选的样例组成的数据集为袋外样本(out-of-bag sample, OOB)。 这样可以用bootstrap样本作为训练集（分析集）， 用没有被选入的样例构成的数据集作为验证集（评估集）， 这样的抽样可用重复\(B\)次， 得到\(B\)对分析、评估集。
用rsample包产生多个bootstrap分析、评估集对的程序示例：
## # Bootstrap sampling
## # A tibble: 5 × 2
## splits id
## <list> <chr>
## 1 <split [2344/875]> Bootstrap1
## 2 <split [2344/865]> Bootstrap2
## 3 <split [2344/860]> Bootstrap3
## 4 <split [2344/870]> Bootstrap4
## 5 <split [2344/871]> Bootstrap5这种方法产生的模型评估指标有很低的标准误差， 但对模型精度的估计会偏向于低估。
37.6.2.4 时间序列的预测验证
如果原始数据集的各个样例行是按时间次序记录的， 则不能随机划分训练集、测试集、分析集、评估集。 测试集必须整体位于比训练集更晚的时间； 每一对分析集、评估集， 评估集对应的时间必须晚于分析集。
假设有\(n = 10000\)个时间点的数据， 则可以将前\(n_1 = 8000\)个样例行作为训练集， 将后\(n_2 = 2000\)个样例行作为测试集。 为了在训练集上评估不同模型的性能， 可以滚动地用前面\(t\)个作为分析集估计模型参数， 用后面\(n_1-t\)个或者其中一部分作为评估集， 评估模型性能。 构成分析集有两种方法：
- 以第\(1\)到第\(t\)行为分析集， 估计模型参数， 然后预测第\(t+1\)行（单步预测问题）， 或者预测第\(t+1\)到第\(t+m\)行（多步预测问题）。 对\(t = n_3\)到\(n_4\)（如\(t\)从4000到7999）的每一个时间点进行如上的估计、预测， 最后汇总\(n_4 - n_3 + 1\)组预测的评估结果。 这可以称为可变窗宽的滚动窗口预测评估。
- 上面的步骤中， 也可以在使用\(1\tilde t\)行作为分析集后， 下一步使用\(1 \tilde t+k\)行作为分析集， 即分析集每次增加\(k\)个时间点， 而不是仅增加一个时间点。
- 以第\(t - n_5 + 1\)到第\(t\)行为分析集， 比如以第\(t-999\)到第\(t\)行的1000行为分析集， 预测第\(t+1\)或者第\(t+1\)到第\(t+m\)行。 这可以称为固定窗宽的滚动窗口预测评估。
下面的样例程序从1到\(t = 6 \times 30\)处开始作为分析集， 评估集使用分析集后面的30个点， 而且下一个分析集跳过29个时间点（即下一分析集使用1到\(t+30\)时间点）。
time_slices <-
tibble(x = 1:365) %>%
rolling_origin(
initial = 6 * 30, # 初始的分析集大小
assess = 30, # 分析集大小
skip = 29, # 从一个分析集到下一个分析集跳过的点，缺省为0
cumulative = FALSE) # 分析集保持固定大小180
data_range <- function(x) {
summarize(x, first = min(x), last = max(x), n = n())
}下面的程序计算每一组分析集、评估集的时间起始、结束位置：
bind_cols(
map_dfr(time_slices$splits,
~ analysis(.x) %>% data_range()) |>
rename(ana_first = first, ana_last = last, ana_n = n),
map_dfr(time_slices$splits,
~ assessment(.x) %>% data_range()) |>
rename(ass_first = first, ass_last = last, ass_n = n) )## # A tibble: 6 × 6
## ana_first ana_last ana_n ass_first ass_last ass_n
## <int> <int> <int> <int> <int> <int>
## 1 1 180 180 181 210 30
## 2 31 210 180 211 240 30
## 3 61 240 180 241 270 30
## 4 91 270 180 271 300 30
## 5 121 300 180 301 330 30
## 6 151 330 180 331 360 30可见分析集始终保持了180个时间点大小， 评估集始终保持了30个时间点大小， 最后的5个时间点不够作为评估集所以弃用， 分析集的末端每次前进\(29+1=30\)个时间点。
37.6.3 性能评估
可以用多组分析集、评估集来评估模型性能。 对每一组分析集、评估集， 应使用如下步骤：
- 在分析集上计算预处理需要参数并进行预处理， 然后在分析集上拟合模型， 获得模型参数估计；
- 用分析集上得到的预处理参数对评估集进行预处理， 用分析集上估计的模型参数对评估集进行预测， 然后计算模型性能指标。
- 最后，将所有分析集、评估集上获得的性能指标进行汇总， 得到单一指标。 tidymodels的指标汇总方法是简单地计算平均值。
对上一小节的各种再抽样方法产生的rset集合，
可以用fit_resamples(workflow, resamples)来对每一对分析、评估集进行估计、评估计算，
可以用选项metrics指定要计算的一个或多个评估指标，
用选项control规定一些输出选项，
比如要求输出每个评估样例的预测值。
例如（重复前面的例子）：
data(iris)
iris_rec <-
recipe(Species ~ .,
data = iris)
rf_model <- rand_forest(
trees = 100,
min_n = 5 ) |>
set_engine("ranger") |>
set_mode("classification")
rf_wflow <-
workflow() |>
add_model(rf_model) |>
add_recipe(iris_rec)
set.seed(101)
iris_folds <- vfold_cv(iris, v = 10)
keep_pred <- control_resamples(
save_pred = TRUE, save_workflow = TRUE)
iris_cv <- rf_wflow |>
fit_resamples(resamples = iris_folds,
control = keep_pred)
iris_cv## # Resampling results
## # 10-fold cross-validation
## # A tibble: 10 × 5
## splits id .metrics .notes .predictions
## <list> <chr> <list> <list> <list>
## 1 <split [135/15]> Fold01 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [15 × 7]>
## 2 <split [135/15]> Fold02 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [15 × 7]>
## 3 <split [135/15]> Fold03 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [15 × 7]>
## 4 <split [135/15]> Fold04 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [15 × 7]>
## 5 <split [135/15]> Fold05 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [15 × 7]>
## 6 <split [135/15]> Fold06 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [15 × 7]>
## 7 <split [135/15]> Fold07 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [15 × 7]>
## 8 <split [135/15]> Fold08 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [15 × 7]>
## 9 <split [135/15]> Fold09 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [15 × 7]>
## 10 <split [135/15]> Fold10 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [15 × 7]>fit_resamples的结果是一个tibble数据框，
包含了每一组分析、评估集上计算的性能指标，
以及评估集每个样例的预测结果。
汇总计算性能指标的方法（使用结果中的mean作为汇总结果）：
## # A tibble: 2 × 6
## .metric .estimator mean n std_err .config
## <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 accuracy multiclass 0.96 10 0.0147 Preprocessor1_Model1
## 2 roc_auc hand_till 0.996 10 0.00201 Preprocessor1_Model1用交叉验证计算的准确率为\(96.00\%\)， AUC(多类别情况下的Hand-Till推广)为\(0.9962\)。
如果在collect_metrics()中加选项summarize = FALSE，
则仅提取每一对分析、评估集的性能指标，
而不计算汇总统计量。
下面的程序将所有评估集的样例的预测都提取出来：
## # A tibble: 150 × 8
## id .pred_setosa .pred_versicolor .pred_virginica .row .pred_class Species
## <chr> <dbl> <dbl> <dbl> <int> <fct> <fct>
## 1 Fold… 1 0 0 5 setosa setosa
## 2 Fold… 1 0 0 9 setosa setosa
## 3 Fold… 1 0 0 46 setosa setosa
## 4 Fold… 0 0.851 0.149 53 versicolor versic…
## 5 Fold… 0 1 0 83 versicolor versic…
## 6 Fold… 0.002 0.996 0.002 90 versicolor versic…
## 7 Fold… 0 1 0 93 versicolor versic…
## 8 Fold… 0.002 0.974 0.024 99 versicolor versic…
## 9 Fold… 0 0 1 116 virginica virgin…
## 10 Fold… 0 0.184 0.816 124 virginica virgin…
## # ℹ 140 more rows
## # ℹ 1 more variable: .config <chr>其中的预测结果列的命名与parsnip包的predict()函数的结果数据框预测结果列命名规则相同。
其中的变量.row可以定位到评估集样例所在的原始数据行。
如果仅使用了一个验证集，
仍可以使用fit_resamples()函数进行估计、评估，
用collect_metrics()提取评估指标值。
37.6.4 利用并行计算功能
tune包的fit_resamples()等函数可以利用并行计算功能，
而且因为在不同的分析、评估集对上的计算（除了时间序列以外）都是互不相关的，
所以可以很容易地将其分解到不同的CPU核心，或操作系统线程，或计算机集群上完成计算。
在单一计算机上， 可以用如下的框架进行并行计算：
37.6.5 保存再抽样建模、预测结果
再抽样方法一般用于超参数调优和模型比较， 所以其中的参数估计结果和在评估集上的预测值默认是不保存的。 在选定最优模型后， 可以使用整个训练集重新估计模型参数。
fit_resamples()函数的control参数可以提供一个用control_resamples()生成的设置，
在control_resamples()中设置save_pred = TRUE可以保存在评估集上的预测结果,
设置save_workflow = TRUE可以保存在分析集上拟合的工作流程结果。
还可以在control_resamples()中设置extract =，
这需要指定一个函数f(x)，
其中x是在分析集上拟合的工作流程结果，
用f(x)返回要从拟合结果中提取的信息。
这种做法比较麻烦，但是更为灵活。
如：
get_model <- function(x) {
extract_fit_parsnip(x) |> tidy()
}
ctrl <- control_resamples(extract = get_model)
lm_wflow <-
workflow() |>
add_recipe(ames_rec) |>
add_model(linear_reg() |> set_engine("lm"))
lm_res <- lm_wflow |>
fit_resamples(resamples = ames_folds, control = ctrl)
lm_res## # Resampling results
## # 10-fold cross-validation
## # A tibble: 10 × 5
## splits id .metrics .notes .extracts
## <list> <chr> <list> <list> <list>
## 1 <split [2109/235]> Fold01 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [1 × 2]>
## 2 <split [2109/235]> Fold02 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [1 × 2]>
## 3 <split [2109/235]> Fold03 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [1 × 2]>
## 4 <split [2109/235]> Fold04 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [1 × 2]>
## 5 <split [2110/234]> Fold05 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [1 × 2]>
## 6 <split [2110/234]> Fold06 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [1 × 2]>
## 7 <split [2110/234]> Fold07 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [1 × 2]>
## 8 <split [2110/234]> Fold08 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [1 × 2]>
## 9 <split [2110/234]> Fold09 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [1 × 2]>
## 10 <split [2110/234]> Fold10 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [1 × 2]>查看extract选项产生的结果对应于第一折的部分：
## # A tibble: 1 × 2
## .extracts .config
## <list> <chr>
## 1 <tibble [71 × 5]> Preprocessor1_Model1## [[1]]
## # A tibble: 71 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -0.0804 0.292 -0.275 7.83e- 1
## 2 Gr_Liv_Area 0.625 0.0161 38.8 1.27e-246
## 3 Year_Built 0.00176 0.000137 12.9 1.76e- 36
## 4 Neighborhood_College_Creek -0.0501 0.0347 -1.44 1.49e- 1
## 5 Neighborhood_Old_Town -0.0421 0.0129 -3.25 1.16e- 3
## 6 Neighborhood_Edwards -0.123 0.0279 -4.39 1.18e- 5
## 7 Neighborhood_Somerset 0.0392 0.0191 2.05 4.03e- 2
## 8 Neighborhood_Northridge_Heights 0.0995 0.0280 3.55 3.90e- 4
## 9 Neighborhood_Gilbert -0.00391 0.0221 -0.177 8.59e- 1
## 10 Neighborhood_Sawyer -0.144 0.0266 -5.40 7.37e- 8
## # ℹ 61 more rows37.7 用再抽样技术比较模型
模型比较， 可能是同一类模型的不同自变量集合、不同预处理方法的比较， 也可能是不同类型模型的比较。 一般是用再抽样方法计算一系列的模型性能指标(如RMSE)， 然后在模型之间比较这些性能指标。
因为涉及到多个模型， tidymodels提供了“工作流程集合”的类型， 可以同时规定多种模型的建模步骤。 需要更详细地了解不同的性能指标， 以便在比较时选用。 还有一些统计推断方法可以用来比较模型。
37.7.1 工作流程集合
用ames数据， 使用不同的预处理步骤和不同模型， 可以生成包含多个工作流程的工作流程集合。
下面的程序使用了不同的工作流程， 仅使用线性回归：
basic_rec <-
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())
interaction_rec <-
basic_rec |>
step_interact( ~ Gr_Liv_Area:starts_with("Bldg_Type_") )
spline_rec <-
interaction_rec |>
step_ns(Latitude, Longitude, deg_free = 50)
preproc <- list(
basic = basic_rec,
interact = interaction_rec,
splines = spline_rec )
lm_models <- workflow_set(
preproc,
list(lm = linear_reg()))
lm_models## # A workflow set/tibble: 3 × 4
## wflow_id info option result
## <chr> <list> <list> <list>
## 1 basic_lm <tibble [1 × 4]> <opts[0]> <list [0]>
## 2 interact_lm <tibble [1 × 4]> <opts[0]> <list [0]>
## 3 splines_lm <tibble [1 × 4]> <opts[0]> <list [0]>对工作流程集合中每个工作流程，
可以使用workflow_map()，
将fit_resamples()函数应用到每个工作流程上。
lm_models <-
lm_models |>
workflow_map(
"fit_resamples",
# 下面是`workflow_map()`的选项:
seed = 1101, verbose = TRUE,
# 下面是`fit_resamples()`的选项:
resamples = ames_folds,
control = keep_pred)## i 1 of 3 resampling: basic_lm## ✔ 1 of 3 resampling: basic_lm (1.6s)## i 2 of 3 resampling: interact_lm## ✔ 2 of 3 resampling: interact_lm (1.9s)## i 3 of 3 resampling: splines_lm## ✔ 3 of 3 resampling: splines_lm (2.7s)## # A workflow set/tibble: 3 × 4
## wflow_id info option result
## <chr> <list> <list> <list>
## 1 basic_lm <tibble [1 × 4]> <opts[2]> <rsmp[+]>
## 2 interact_lm <tibble [1 × 4]> <opts[2]> <rsmp[+]>
## 3 splines_lm <tibble [1 × 4]> <opts[2]> <rsmp[+]>结果中的result包含fit_resamples()函数计算出的再抽样性能指标，
option包含了再抽样的分析、评估集对集合，
与fit_resamples()所用的control选项值。
可以用collect_metrics()从拟合过的工作流程集合对象提取所有的性能指标值，
然后可以筛选比较，如：
## # A tibble: 3 × 9
## wflow_id .config preproc model .metric .estimator mean n std_err
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <int> <dbl>
## 1 basic_lm Preprocesso… recipe line… rmse standard 0.0783 10 0.00211
## 2 interact_lm Preprocesso… recipe line… rmse standard 0.0778 10 0.00212
## 3 splines_lm Preprocesso… recipe line… rmse standard 0.0765 10 0.00186从mean列的结果比较，
最后一个食谱（包含交互作用和样条项）较优。
对于随机森林方法， 可以生成这个方法的工作流程集合， 然后与线性回归方法的工作流程合并， 最后比较：
rf_model <-
rand_forest(trees = 1000) |>
set_engine("ranger") |>
set_mode("regression")
rf_wflow <-
workflow() |>
add_model(rf_model) |>
add_formula(
Sale_Price ~ Neighborhood + Gr_Liv_Area
+ Year_Built + Bldg_Type
+ Latitude + Longitude)
rf_res <- fit_resamples(
rf_wflow,
resamples = ames_folds,
control = keep_pred)
four_models <-
as_workflow_set(
random_forest = rf_res) |>
bind_rows(lm_models)
four_models## # A workflow set/tibble: 4 × 4
## wflow_id info option result
## <chr> <list> <list> <list>
## 1 random_forest <tibble [1 × 4]> <opts[0]> <rsmp[+]>
## 2 basic_lm <tibble [1 × 4]> <opts[2]> <rsmp[+]>
## 3 interact_lm <tibble [1 × 4]> <opts[2]> <rsmp[+]>
## 4 splines_lm <tibble [1 × 4]> <opts[2]> <rsmp[+]>比较4个模型的RMSE：
## # A tibble: 4 × 9
## wflow_id .config preproc model .metric .estimator mean n std_err
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <int> <dbl>
## 1 random_forest Preproces… formula rand… rmse standard 0.0714 10 0.00220
## 2 basic_lm Preproces… recipe line… rmse standard 0.0783 10 0.00211
## 3 interact_lm Preproces… recipe line… rmse standard 0.0778 10 0.00212
## 4 splines_lm Preproces… recipe line… rmse standard 0.0765 10 0.00186随机森林模型结果较优。
下面的图形比较了不同工作流程的\(R^2\)统计量：
library(ggrepel)
autoplot(
four_models,
metric = "rsq") +
geom_text_repel(
aes(label = wflow_id),
nudge_x = 1/8,
nudge_y = 1/100) +
theme(legend.position = "none")
因为每个工作流程用多组分析、评估集合计算了多个\(R^2\)结果， 所以可以做成置信区间的图形。 这里的\(R^2\)定义为因变量和预测值的相关系数的平方。
37.7.2 用再抽样技术计算的模型性能指标的选择
使用不同的性能指标比较各个模型， 最后选出的最优模型可能会不一致。 这就需要我们了解这些性能指标， 然后根据需要选取。
因为不同的模型使用了相同的交叉验证分层数据， 所以这些模型计算得到的同一性能指标的相同交叉验证层之间的值有很强的正相关性。 这样，在用统计检验比较两个模型时， 这种正相关性会干扰模型比较。 比如， 不能用独立两样本检验比较两个模型的\(RMSE\)指标是否相等， 强正相关性会使得两个模型的指标之间的差异变小， 使得检验不出显著差异。
作为显著性检验的替代， 可以预先确定一个“实际效应”阈值， 当差别超过这个阈值时就认为两个模型有差异。 虽然这种做法比较主观， 但是在选择模型时很有用。
37.7.4 用贝叶斯方法比较模型
仍可以用性能指标作为因变量， 建立方差分析模型， 但是假设模型系数、误差标准差有先验分布， 然后计算后验分布。 可以将不同的交叉验证层作为一个随机效应截距。 tidyposterior包可以用来进行这样的贝叶斯计算。
比较4个模型的\(R^2\)的例子：
## Warning: package 'tidyposterior' was built under R version 4.2.3## Warning: package 'rstanarm' was built under R version 4.2.3## Loading required package: Rcpp##
## Attaching package: 'Rcpp'## The following object is masked from 'package:rsample':
##
## populate## This is rstanarm version 2.21.4## - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!## - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.## - For execution on a local, multicore CPU with excess RAM we recommend calling## options(mc.cores = parallel::detectCores())rsq_anova <-
perf_mod(
four_models,
metric = "rsq",
prior_intercept = rstanarm::student_t(df = 1),
chains = 4,
iter = 5000,
seed = 1102
)用\(R^2\)的后验分布的随机样本代表后验分布：
## Rows: 40,000
## Columns: 2
## $ model <chr> "random_forest", "basic_lm", "interact_lm", "splines_lm", "r…
## $ posterior <dbl> 0.8308552, 0.7984321, 0.8001101, 0.8058113, 0.8340578, 0.798…四个模型的\(R^2\)的后验分布图形：
model_post |>
mutate(model = forcats::fct_inorder(model)) |>
ggplot(aes(x = posterior)) +
geom_histogram(
bins = 50,
color = "white",
fill = "blue",
alpha = 0.4) +
facet_wrap(~ model, ncol = 1)
各个模型的\(R^2\)置信区间图形：
library(ggrepel)
autoplot(rsq_anova) +
geom_text_repel(
aes(label = workflow),
nudge_x = 1/8, nudge_y = 1/100) +
theme(legend.position = "none")
不同模型的后验分布随机样本是按行一一对应的， 所以很容易生成比较用的样本。 例如， 比较基本的线性模型和加入交互项、样条项的线性模型：
rqs_diff <-
contrast_models(
rsq_anova,
list_1 = "splines_lm",
list_2 = "basic_lm",
seed = 1104)
rqs_diff |>
as_tibble() |>
ggplot(aes(x = difference)) +
geom_vline(xintercept = 0, lty = 2) +
geom_histogram(
bins = 50, color = "white",
fill = "red", alpha = 0.4)
这是`list_1 - list_2\(的\)R^2$后验样本的直方图。
计算差值的可信区间(credible interval)：
## # A tibble: 1 × 6
## contrast probability mean lower upper size
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 splines_lm vs basic_lm 1 0.00899 0.00543 0.0126 0其中的概率是差值大于0的概率。 差值的均值仅有\(0.0085\)， 如果要求差值显著超过\(0.02\)才算是有效的差距， 可以如下计算：
## # A tibble: 1 × 4
## contrast pract_neg pract_equiv pract_pos
## <chr> <dbl> <dbl> <dbl>
## 1 splines_lm vs basic_lm 0 1 0其中pract_equiv列是差值的后验样本中位于[-size, size]区间内的比例，
可以认为以\(0.02\)为有效差距阈值，
则两个方法之间没有有效差距。
可以用图形表示多组之间是否有有效差距：
autoplot(
rsq_anova, type = "ROPE", size = 0.02) +
geom_text_repel(aes(label = workflow)) +
theme(legend.position = "none")
图形说明，以\(0.02\)为有效差值的阈值， 则三个线性模型都和随机森林模型有有效差距， 而这三个之间没有有效差距。
为了提高贝叶斯推断的精度， 可以使用重复交叉验证方法生成更多的交叉验证样本， 这可以缩短贝叶斯可信区间长度。
37.8 模型参数调优和过度拟合问题
37.8.1 概念
在建模时， 有些模型参数是直接从训练数据估计的； 而其它一些参数， 称为调节参数(tuning parameters)或者超参数(hypterparmeters)， 不能直接估计而是需要预先指定。
考虑线性回归模型 \[ y_i = a + b x_i + e_i, \ i=1,2,\dots,n . \] 其中的参数\(a, b\)可以用最小二乘估计法直接从训练集得到估计结果。
有一些参数则不能直接从训练集估计，称为“调节参数”或者“超参数”。 比如：
- 随机森林算法在训练集上训练之前， 需要指定树的深度。
- GBM算法需要指定树的棵数、树的深度、学习率。
- lasso回归需要指定\(L^1\)惩罚项的系数， 称为调节参数， 不能直接估计。
- 在多层神经网络模型中， 网络的层数、每一层的节点数、选用的激活函数， 都需要预先确定， 不能从训练集估计。
- 深度神经网络中梯度下降算法的改进版本， 需要确定学习率，动量系数， 运行的轮数(epochs)等， 也需要调优，但不能从训练集估计。
有时需要对预处理步骤调优，如：
- 在主成分分析、偏最小二乘中， 可以对原始变量作变换， 突出其线性相关性； 可以选择输出的降维变量个数。
- 在对缺失值进行插补时， 可以对某一变量用最近邻方法插补， 这时需要确定使用的近邻个数。
有一些统计模型有结构上的不同设置， 这可以看成是模型比较问题， 也可以广义地看成是模型调优问题：
- 对于二值因变量建立广义线性模型， 联系函数一般使用logit函数， 但也可以使用probit或者余重对数(complementary log-log)。
- 在线性混合模型中， 随机效应的方差阵可以指定为常数方差对角阵、 对角阵或者自回归的、Toeplitz的。
不能作为参数调优问题的一个典型例子是贝叶斯分析中的先验分布选择。 先验分布代表了研究者对参数的先验知识， 这不能从训练集中调优， 否则就违背了贝叶斯分析的原则。
另外， 在随机森林算法和bagging算法中， 树的棵数也不需要调优， 只要取得比较大（对随机森林，达到几千，对bagging，50到100即可）， 性能就没有进一步的改进， 也不会造成过度拟合问题。
37.8.2 参数调优的标准
为了参数调优应该选用什么样的性能指标？ 这依赖于考虑的模型以及建模的目的。 如果主要目的是进行统计推断， 而模型选择或比较可以通过统计推断方法进行， 可以利用统计假设检验、置信区间来进行模型比较， 一般这不需要使用再抽样技术，只用训练集上的检验结果即可； 但是， 如果建模目的是预报， 则预报精度是更重要的考虑， 应使用再抽样技术进行比较， 统计推断得到的模型比较结论与按预报精度得到结论有可能不一致。
37.8.3 超参数选择的问题
超参数在许多模型中控制的是模型的复杂度。 模型过于简单， 则模型对训练集和测试集的预报效果都不好； 模型越复杂， 越能够抓住数据中的复杂规律， 但是也可能产生过度拟合， 即抓住了训练集中一些随机发生的非规律性的模式， 在测试集上表现很差。 见32.6的示例。
在选择超参数时必须使用再抽样技术， 把估计普通参数的数据与选择超参数的数据分开， 否则选择超参数就会倾向于完美拟合训练集， 产生过度拟合， 在测试集上表现会很差。
37.8.4 超参数优化的两种常用策略
超参数优化经常使用网格搜索法(grid search)和迭代法(iterative search)。
网格法有些像是方差分析中的全面试验设计， 每个要优化的超参数选取若干个最可能的值， 然后将这些超参数进行完全组合。 比如， 超参数A取值为\(0.1, 0.2\)， 超参数B取值\(2, 5\)， 超参数C取值\(10, 100, 1000\)， 则网格设计表如下：
| A | B | c |
|---|---|---|
| 0.1 | 2 | 10 |
| 0.1 | 2 | 100 |
| 0.1 | 2 | 1000 |
| 0.1 | 5 | 10 |
| 0.1 | 5 | 100 |
| 0.1 | 5 | 1000 |
| 0.2 | 2 | 10 |
| 0.2 | 2 | 100 |
| 0.2 | 2 | 1000 |
| 0.2 | 5 | 10 |
| 0.2 | 5 | 100 |
| 0.2 | 5 | 1000 |
共有\(2 \times 2 \times 3 = 12\)个网格点， 需要对12种超参数组合用再抽样技术进行比较。
除了使用以上的完全组合作为网格外， 还可以使用不规则的网格， 类似于拟随机数设计， 使得网格点在空间中既分布比较均匀， 又互相远离， 称为“空间填充设计”。
迭代搜索法则从某个初始的参数组合的出发， 针对指定的目标函数（即利用再抽样方法计算的模型性能指标）进行优化迭代， 原则上可以使用任意的优化迭代方法， 但是要注意到目标函数利用了再抽样方法， 所以目标函数值有随机误差， 不能要求太高的优化精度， 一般可以规定搜索使用的时间阈值或者迭代步数阈值。
一种迭代方法是结合了网格法的简单， 与迭代搜索法的灵活性， 即指定一个参数组合的较大网格， 然后在这些网格点上随机地跳跃， 直到达到规定的运行时间或者搜索步数， 使用已搜索到的最优超参数组合作为输出。
还可以用网格法搜索迭代初始值， 然后用网格法得到的参数组合作为初值进行迭代搜索。
37.8.5 调优超参数管理
tidymodels系统用dials扩展包管理超参数组合。 因为tidymodels对同一问题支持多种模型， 对同一模型支持多种计算引擎， 所以将超参数分为两个级别：
- 主要参数， 这是同一模型的每个计算引擎都要设置的， 只不过使用的参数名称不一定相同， parsnip对这些参数进行了统一化， 只要使用parsnip规定的参数名称即可。
- 引擎特有的参数，
这些参数一般不需要进行优化。
如需指定，
可以在
set_engine()函数中与计算引擎一同指定。
例如：
rand_forest(trees = 2000, min_n = 10) # trees和min_n是主要参数
set_engine("ranger",
regularization.factor = 0.5) # regularization.factor是引擎特有参数在指定模型的函数(如上面的rand_forest())中，
可以指定超参数的值，
也可以将超参数说明为= tune()，
则是规定该超参数通过超参数调优获得。
如：
tune()函数并不实际运行，
只起到说明作用，
后面要讲到的模型调优函数则能够识别工作流程中的这样的说明，
进行实际的参数调优。
这些参数调优函数能够识别这些要调优的参数，
自动获得其取值范围和常用值。
用extract_parameter_set_dials()可以查看哪些超参数需要调优：
## Collection of 2 parameters for tuning
##
## identifier type object
## mtry mtry nparam[?]
## min_n min_n nparam[+]
##
## Model parameters needing finalization:
## # Randomly Selected Predictors ('mtry')
##
## See `?dials::finalize` or `?dials::update.parameters` for more information.nparam[?]表示取值范围待定的数值型参数，
它对应的mtry参数是随机森林算法每一次分叉时允许使用的自变量个数，
这个数值依赖于训练数据集中的自变量个数。
params[+]表示要调优的min_n参数是一个数值型参数且范围可确定。
结果中的identifier在有多于一处需要对同一参数调优时可以区分不同的调优参数，
比如，
有两个自变量都需要用ns()函数转换，
则需要有两个自由度参数，
程序如：
ames_rec <-
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 = tune()) |>
step_dummy(all_nominal_predictors()) |>
step_interact( ~ Gr_Liv_Area:starts_with("Bldg_Type_") ) |>
step_ns(Longitude, deg_free = tune("longitude df")) |>
step_ns(Latitude, deg_free = tune("latitude df"))
recipes_param <- extract_parameter_set_dials(ames_rec)
recipes_param## Collection of 3 parameters for tuning
##
## identifier type object
## threshold threshold nparam[+]
## longitude df deg_free nparam[+]
## latitude df deg_free nparam[+]其中tune()函数增加了一个用来作为调节参数标识符的输入项。
预处理步骤和模型中都可以有超参数， 将二者结合到工作流程中以后， 这些超参数都可以看到， 如：
wflow_param <-
workflow() |>
add_recipe(ames_rec) |>
add_model(rf_mod) |>
extract_parameter_set_dials()
wflow_param## Collection of 5 parameters for tuning
##
## identifier type object
## mtry mtry nparam[?]
## min_n min_n nparam[+]
## threshold threshold nparam[+]
## longitude df deg_free nparam[+]
## latitude df deg_free nparam[+]
##
## Model parameters needing finalization:
## # Randomly Selected Predictors ('mtry')
##
## See `?dials::finalize` or `?dials::update.parameters` for more information.wflow_param这样的对象被超参数调优函数使用。
查看某一个超参数的设置：
## Threshold (quantitative)
## Range: [0, 0.1]可以修改wflow_param这样的对象，
如：
wflow_param <- wflow_param |>
update(threshold = threshold(c(0.8, 1.0)))
wflow_param |>
extract_parameter_dials("threshold")## Threshold (quantitative)
## Range: [0.8, 1]又如，
人为确定mtry参数的取值范围：
wflow_param <- wflow_param |>
update(mtry = mtry(c(1, 70)))
wflow_param |>
extract_parameter_dials("mtry")## # Randomly Selected Predictors (quantitative)
## Range: [1, 70]注意ames数据集共有74个变量，
最多可以有73个自变量。
但是，如果预处理步骤或者算法可能会增删自变量，
这样的范围指定有可能不适合，
可以对调优参数对象使用finalize()函数，
指定一个训练集获得实际可用的参数范围，如：
pca_rec <-
recipe(Sale_Price ~ ., data = ames_train) |>
# 将与面积有关的自变量用主成分分析法压缩:
step_normalize(contains("SF")) |>
# 用方差累计贡献率选择主成分个数
step_pca(contains("SF"), threshold = .95)
updated_param <-
workflow() |>
add_model(rf_mod) |>
add_recipe(pca_rec) |>
extract_parameter_set_dials() |>
# 这一步可以确定mtry的范围:
finalize(ames_train |> select(-Sale_Price))
updated_param## Collection of 2 parameters for tuning
##
## identifier type object
## mtry mtry nparam[+]
## min_n min_n nparam[+]## # Randomly Selected Predictors (quantitative)
## Range: [1, 73]引擎特有的超参数也可以用extract_parameter_set_dials()管理。
有些超参数使用变换后的值更方便，
需要了解这些变换，
人为指定范围时也要使用变换后的范围，
如：
## Amount of Regularization (quantitative)
## Transformer: log-10 [1e-100, Inf]
## Range (transformed scale): [-10, 0]这个参数使用时用了实际参数值的常用对数变换，
指定范围时需要在变换后指定。
比如，
指定penalty在\([0.1, 1]\)范围，
需要作常用对数变换，
指定范围为\([-1, 0]\)，
如：
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1003 0.1771 0.3200 0.3919 0.5754 0.9922上例指定了penalty在\([0.1, 1]\)范围，
并随机生成了1000个值，
计算这些值的简单统计量（在没有对数变换的基础上计算）。
可以用trans = NULL指定在未变换的基础上指定范围，
如：
## Amount of Regularization (quantitative)
## Range: [0.1, 1]37.9 超参数调优的网格搜索法
网格搜索法预先确定要试验的超参数组合， 并对这些组合逐一使用再抽样技术计算模型性能参数， 从而找到最优超参数组合。
37.9.1 规则网格与不规则网格
规则网格法产生的超参数网格类似于方差分析中的完全试验设计， 由每个超参数选择少量水平值作完全组合产生网格。 非规则网格则是在超参数组合空间内直接布点。
以单隐藏层神经网络为例。 要调节的超参数包括隐藏层单元数、拟合迭代的轮数(epochs)、权重衰减系数。
指定模型、引擎，
在模型函数mlp()中指定主要超参数，
三个主要超参数都需要调优：
mlp_spec <-
mlp(
hidden_units = tune(),
penalty = tune(),
epochs = tune()) |>
set_engine("nnet", trace = 0) |>
set_mode("classification")用extract_parameter_set_dials()提取超参数集合：
用extract_parameter_dials()显示各个超参数设置：
## # Hidden Units (quantitative)
## Range: [1, 10]## Amount of Regularization (quantitative)
## Transformer: log-10 [1e-100, Inf]
## Range (transformed scale): [-10, 0]## # Epochs (quantitative)
## Range: [10, 1000]37.9.1.1 使用规则网格
前面的mlp_param仅是列出了要调优的参数以及每个参数的可取范围，
但是没有给出具体的参数组合。
规则网格法则对每个超参数指定若干个典型取值，
然后完全组合。
tidyr包的crossing()函数可以生成这样的完全组合，
如：
## # A tibble: 12 × 3
## hidden_units penalty epochs
## <int> <dbl> <dbl>
## 1 1 0 100
## 2 1 0 200
## 3 1 0.1 100
## 4 1 0.1 200
## 5 2 0 100
## 6 2 0 200
## 7 2 0.1 100
## 8 2 0.1 200
## 9 3 0 100
## 10 3 0 200
## 11 3 0.1 100
## 12 3 0.1 200实际上， dials包的调优函数给出的超参数设置可以自动产生规则网格， 所以仅在用户有经验时才人为指定规则网格， 如：
## # A tibble: 8 × 3
## hidden_units penalty epochs
## <int> <dbl> <int>
## 1 1 0.0000000001 10
## 2 10 0.0000000001 10
## 3 1 1 10
## 4 10 1 10
## 5 1 0.0000000001 1000
## 6 10 0.0000000001 1000
## 7 1 1 1000
## 8 10 1 1000其中的levels输入指定每个超参数要取几个典型水平值。
也可以每个超参数分别指定水平值个数，如：
## # A tibble: 12 × 3
## hidden_units penalty epochs
## <int> <dbl> <int>
## 1 1 0.0000000001 10
## 2 5 0.0000000001 10
## 3 10 0.0000000001 10
## 4 1 1 10
## 5 5 1 10
## 6 10 1 10
## 7 1 0.0000000001 1000
## 8 5 0.0000000001 1000
## 9 10 0.0000000001 1000
## 10 1 1 1000
## 11 5 1 1000
## 12 10 1 1000也可以不进行完全试验， 而是在上述的完全组合中选一部分， 比如正交设计法。
当超参数个数较多时， 对大多数模型来说规则网格法耗时较多， 也有一些模型例外。 规则网格法优点是每个超参数的调优作用比较容易理解， 需要的话也可以研究一些交互作用效应。
37.9.1.2 非规则网格法
最简单的非规则网格方法是在超参数取值范围内进行随机抽样。
这样理论上可以比较均匀，
但是结果具有随机性，
实际上可能有些点距离较近，
而有些参数空间位置附近没有网格点。
用grid_random()产生这样的网格。
如：
## hidden_units penalty epochs
## Min. : 1.000 Min. :0.0000000 Min. : 10.0
## 1st Qu.: 3.000 1st Qu.:0.0000000 1st Qu.:265.8
## Median : 5.000 Median :0.0000061 Median :497.0
## Mean : 5.381 Mean :0.0437435 Mean :509.5
## 3rd Qu.: 8.000 3rd Qu.:0.0026854 3rd Qu.:761.0
## Max. :10.000 Max. :0.9814405 Max. :999.0其中的size是非规则网格点个数。
如果某些超参数进行了变换，
随机抽样是在变换后的范围内均匀抽样的。
更合理的做法是用“空间填充设计”， 其目标是尽可能覆盖整个超参数空间范围， 又避免相距过近的点。 具体方法有拉丁方设计， 最大熵设计， 最大投影设计， 等等。
dials包能使用拉丁方设计与最大熵设计。 使用拉丁方设计如：
library(ggforce)
set.seed(1303)
mlp_param |>
grid_latin_hypercube(
size = 20, original = FALSE) |>
ggplot(aes(x = .panel_x, y = .panel_y)) +
geom_point() +
geom_blank() +
facet_matrix(
vars(hidden_units, penalty, epochs),
layer.diag = 2) +
labs(title = "Latin Hypercube design with 20 candidates")
tune扩展包默认使用最大熵设计产生非规则网格，
函数为grid_max_entropy()
这种方法覆盖超参数空间效果好，
能更快搜索最优超参数组合。
37.9.2 评估不同超参数组合的性能
使用网格搜索法寻找最优超参数组合， 本质上是对每一超参数组合利用重抽样方法计算模型性能指标。
以modeldata包的cells数据集作为示例。
数据集中有2019例肺癌病人的医学影像数据，
因变量class取两个类别值，
56个自变量则是细胞不同部位的大小、形状、密度测量值，
自变量存在较高相关性，
自变量分布有的不对称。
建模目的是预测分类。
上面删去的变量case用来区分每个病例，
既不是因变量也不是自变量。
本节仅演示参数调优问题， 所以不进行训练集、测试集拆分， 使用所有观测作为训练集。 使用十折交叉验证进行参数调优：
因为自变量的自相关性，
预处理增加一个用主成分分析降维的步骤，
保留主成分个数也作为预处理步骤的待调优超参数。
因为有些自变量分布不对称，
这使得基于方差的主成分方法受到离群值的影响较大，
所以利用了step_YeoJohnson()将自变量分布变得更加对称。
因为获得的主成分的方差是递减的，
所以在产生主成分以后再次用step_normalize()进行变量标准化。
程序如下：
mlp_rec <-
recipe(class ~ ., data = cells) |>
step_YeoJohnson(all_numeric_predictors()) |>
step_normalize(all_numeric_predictors()) |>
step_pca(
all_numeric_predictors(),
num_comp = tune()) |>
step_normalize(all_numeric_predictors())
mlp_spec <-
mlp(
hidden_units = tune(),
penalty = tune(),
epochs = tune()) |>
set_engine("nnet", trace = 0) |>
set_mode("classification")
mlp_wflow <-
workflow() |>
add_model(mlp_spec) |>
add_recipe(mlp_rec)下面提取用tune()指定的待调优超参数集合，
并人为修改部分超参数的可选范围：
mlp_param <-
mlp_wflow |>
extract_parameter_set_dials() |>
update(
epochs = epochs(c(50, 200)),
num_comp = num_comp(c(0, 40))
)其中的num_comp是PCA预处理步的保留主成分个数。
在step_pca()中，
如果主成分个数选取为0个，
就作为使用原始变量而不进行维数压缩的处理。
下面利用tune_grid()函数在已经定义的超参数网格上进行超参数调优。
这与fit_resamples()本质上相同，
但对多组超参数组合进行计算。
tune_grid()以工作流程和一个再抽样对象（如交叉验证对象）为输入。
可以指定要计算的性能指标集合。
最好利用并行计算功能。 程序如：
library(doParallel)
# 产生一个集群对象并注册
cl <- makePSOCKcluster(8)
registerDoParallel(cl)
roc_res <- metric_set(roc_auc) # 性能指标集合
set.seed(1305)
mlp_reg_tune <-
mlp_wflow |>
tune_grid(
cell_folds,
grid = mlp_param |> grid_regular(levels = 3),
metrics = roc_res )
# 取消集群对象
stopCluster(cl)
# 卸载并行库
detach("package:doParallel")
mlp_reg_tune## # Tuning results
## # 10-fold cross-validation
## # A tibble: 10 × 4
## splits id .metrics .notes
## <list> <chr> <list> <list>
## 1 <split [1817/202]> Fold01 <tibble [81 × 8]> <tibble [0 × 3]>
## 2 <split [1817/202]> Fold02 <tibble [81 × 8]> <tibble [0 × 3]>
## 3 <split [1817/202]> Fold03 <tibble [81 × 8]> <tibble [0 × 3]>
## 4 <split [1817/202]> Fold04 <tibble [81 × 8]> <tibble [0 × 3]>
## 5 <split [1817/202]> Fold05 <tibble [81 × 8]> <tibble [0 × 3]>
## 6 <split [1817/202]> Fold06 <tibble [81 × 8]> <tibble [0 × 3]>
## 7 <split [1817/202]> Fold07 <tibble [81 × 8]> <tibble [0 × 3]>
## 8 <split [1817/202]> Fold08 <tibble [81 × 8]> <tibble [0 × 3]>
## 9 <split [1817/202]> Fold09 <tibble [81 × 8]> <tibble [0 × 3]>
## 10 <split [1818/201]> Fold10 <tibble [81 × 8]> <tibble [0 × 3]>共有4个调优超参数， 对于规则网格， 试验设计类似完全试验设计， 可以作图显示不同参数的性能比较：
从图形明显看出， 超参数“正则性惩罚的量”以最高值1为最优。 epochs的差异不明显。 主成分个数好像以0个（不作主成分变换）为优。 隐藏单元的个数的影响与其它超参数有交互作用， 且仅在正则性惩罚非最优时影响较大。
用show_best()可以从tune_grid()的输出结果中显示若干个性能相近的最优模型：
| hidden_units | penalty | epochs | num_comp | .metric | mean | n | std_err | .config |
|---|---|---|---|---|---|---|---|---|
| 5 | 1 | 50 | 0 | roc_auc | 0.8973350 | 10 | 0.0085720 | Preprocessor1_Model08 |
| 10 | 1 | 125 | 0 | roc_auc | 0.8945163 | 10 | 0.0089850 | Preprocessor1_Model18 |
| 10 | 1 | 50 | 0 | roc_auc | 0.8944248 | 10 | 0.0096021 | Preprocessor1_Model09 |
| 5 | 1 | 200 | 0 | roc_auc | 0.8938910 | 10 | 0.0078434 | Preprocessor1_Model26 |
| 5 | 1 | 125 | 0 | roc_auc | 0.8924622 | 10 | 0.0082247 | Preprocessor1_Model17 |
仍使用结果表格中的mean列作为不同超参数组合的评价指标。
从结果看，
有必要增加penalty超参数的值进行网格搜索。
为了使用空间填充设计的非规则网格搜索法，
可以在tune_grid()中指定grid =一个网格点个数，
或者取为grid_max_entropy()这样的网格生成函数结果。
如：
set.seed(1306)
mlp_sfd_tune <-
mlp_wflow |>
tune_grid(
cell_folds,
grid = 20,
param_info = mlp_param, # 指定调整过范围的超参数集合
metrics = roc_res
)
mlp_sfd_tune## # Tuning results
## # 10-fold cross-validation
## # A tibble: 10 × 4
## splits id .metrics .notes
## <list> <chr> <list> <list>
## 1 <split [1817/202]> Fold01 <tibble [20 × 8]> <tibble [0 × 3]>
## 2 <split [1817/202]> Fold02 <tibble [20 × 8]> <tibble [0 × 3]>
## 3 <split [1817/202]> Fold03 <tibble [20 × 8]> <tibble [0 × 3]>
## 4 <split [1817/202]> Fold04 <tibble [20 × 8]> <tibble [0 × 3]>
## 5 <split [1817/202]> Fold05 <tibble [20 × 8]> <tibble [0 × 3]>
## 6 <split [1817/202]> Fold06 <tibble [20 × 8]> <tibble [0 × 3]>
## 7 <split [1817/202]> Fold07 <tibble [20 × 8]> <tibble [0 × 3]>
## 8 <split [1817/202]> Fold08 <tibble [20 × 8]> <tibble [0 × 3]>
## 9 <split [1817/202]> Fold09 <tibble [20 × 8]> <tibble [0 × 3]>
## 10 <split [1818/201]> Fold10 <tibble [20 × 8]> <tibble [0 × 3]>对非规则网络的调优结果，
用autoplot()作图，
类似于假设了一个可加模型，
作因变量（性能指标）对各个自变量的散点图：
因为实际上超参数之间可能有交互作用， 所以上述的图形可能不完全反映超参数组合的影响。
用show_best()显示最优的几个网格点：
## # A tibble: 5 × 9
## hidden_units penalty epochs num_comp .metric mean n std_err .config
## <int> <dbl> <int> <int> <chr> <dbl> <int> <dbl> <chr>
## 1 8 0.594 97 22 roc_auc 0.880 10 0.00998 Prepro…
## 2 3 0.00000000649 135 8 roc_auc 0.878 10 0.00955 Prepro…
## 3 9 0.141 177 11 roc_auc 0.873 10 0.0104 Prepro…
## 4 8 0.0000000103 74 9 roc_auc 0.869 10 0.00761 Prepro…
## 5 6 0.00581 129 15 roc_auc 0.865 10 0.00658 Prepro…其中最优的AUC值\(0.88\)， 比规则网格达到的\(0.89\)略差。 规则网格计算了81个网格点， 这里仅使用了20个网格点。
在用网格法进行参数调优时， 还应该考虑多个性能指标， 而不是完全按照同一个性能指标。 另外， 还可以不选择搜索到的性能指标最优的组合， 而是选用次优但是更简单的组合。 是否“更简单”可以用代表模型复杂度的参数来衡量。
参数调优仅为了搜索最优的超参数组合，
所以在tune_grid()运行完以后一般不需要保存评估集上的因变量预测值和各个分析、评估集上计算的性能指标，
如果需要可以用control参数输入函数control_grid()的设置结果。
在control_grid()中设置extract = TRUE可以保留预处理、拟合参数结果，
设置save_pred = TRUE可以保留预测值，
在调优后用collect_predictions()提取预测值。
37.9.3 确定最优超参数
tune_grid()函数仅在网格点上计算模型性能指标，
show_best()函数可以显示最优的若干组超参数组合，
但都用来不在训练集上估计模型，
所以为了确定某个参数组合用来估计模型，
可以人为记录超参数组合然后直接在建模函数值指定超参数值，
还可以用select_*()函数从tune_grid()结果中挑选超参数组合。
例如，select_best()取出指定的性能指标最优的超参数组合：
## # A tibble: 1 × 5
## hidden_units penalty epochs num_comp .config
## <int> <dbl> <int> <int> <chr>
## 1 5 1 50 0 Preprocessor1_Model08这用一个单行数据框形式返回了AUC参数最优的超参数组合。
也可以人为输入超参数组合到一个单行的数据框。
得到单行数据框形式的超参数组合后，
用finalize_workflow()函数进行模型估计：
logistic_param <- tibble(
num_comp = 0, # 主成分个数
epochs = 125, # 训练轮数
hidden_units = 1, # 隐藏单元个数
penalty = 1 # 正则性惩罚项因子
)
final_mlp_wflow <-
mlp_wflow |>
finalize_workflow(logistic_param)
final_mlp_wflow## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: mlp()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 4 Recipe Steps
##
## • step_YeoJohnson()
## • step_normalize()
## • step_pca()
## • step_normalize()
##
## ── Model ───────────────────────────────────────────────────────────────────────
## Single Layer Neural Network Model Specification (classification)
##
## Main Arguments:
## hidden_units = 1
## penalty = 1
## epochs = 125
##
## Engine-Specific Arguments:
## trace = 0
##
## Computational engine: nnet注意上面是超参数固定的工作流程， 还没有进行模型拟合。 在训练集（现在是整个数据集）拟合模型：
如果没有使用工作流程，
而是使用了模型，
则可以使用finalize_model()在模型中固定超参数值。
如果使用了食谱，
可以用finalize_recipe()在食谱中固定预处理超参数值。
37.9.4 调优工作的辅助设置参数工具
usemodels包可以用数据框和模型公式作为输入， 产生模型调优需要的R代码， 所以是一种元编程工具。
例如， 对ames数据集产生xgboost建模用的代码：
## Warning: package 'usemodels' was built under R version 4.2.3use_xgboost(
Sale_Price ~ Neighborhood + Gr_Liv_Area
+ Year_Built + Bldg_Type
+ Latitude + Longitude,
data = ames_train,
# verbose选项指示增加说明生成的代码的注释
verbose = TRUE)## xgboost_recipe <-
## recipe(formula = Sale_Price ~ Neighborhood + Gr_Liv_Area + Year_Built + Bldg_Type +
## Latitude + Longitude, data = ames_train) %>%
## step_novel(all_nominal_predictors()) %>%
## ## This model requires the predictors to be numeric. The most common
## ## method to convert qualitative predictors to numeric is to create
## ## binary indicator variables (aka dummy variables) from these
## ## predictors. However, for this model, binary indicator variables can be
## ## made for each of the levels of the factors (known as 'one-hot
## ## encoding').
## step_dummy(all_nominal_predictors(), one_hot = TRUE) %>%
## step_zv(all_predictors())
##
## xgboost_spec <-
## boost_tree(trees = tune(), min_n = tune(), tree_depth = tune(), learn_rate = tune(),
## loss_reduction = tune(), sample_size = tune()) %>%
## set_mode("regression") %>%
## set_engine("xgboost")
##
## xgboost_workflow <-
## workflow() %>%
## add_recipe(xgboost_recipe) %>%
## add_model(xgboost_spec)
##
## set.seed(60916)
## xgboost_tune <-
## tune_grid(xgboost_workflow, resamples = stop("add your rsample object"), grid = stop("add number of candidate points"))在生成的代码的食谱部分，
step_novel()对每个分类自变量指定一个训练集中未出现的新水平。
step_dummy()将所有分类自变量转换为哑变量，
且哑变量个数等于类水平数。
step_zv()删去取常数值的自变量。
最后的tune_grid()调用是一个模板，
不能直接运行，
需要将resamples参数设置为某个再抽样分析、评估集对集合，
将grid参数设置为某个规则或非规则超参数网格。
如果在use_xgboost()调用中设置tune = FALSE，
则会生成模型拟合用的R代码，
而不是超参数调优用的R代码。
37.9.5 网格搜索加速工具
参数调优总是比较耗时的， 这里介绍一些加速网格搜索法参数调优的辅助工具。
37.9.5.1 优化子模型
有些模型可以一次拟合就获得多个调优参数对应的性能指标， lasso回归的调节参数就可以这样调优。
偏最小二乘也有这样的特点。 偏最小二乘与主成分分析类似， 但有因变量存在， 在最大化自变量线性组合的方差的同时， 使得因变量与自变量线性组合的相关性达到最大。 一个超参数是要保留的自变量线性组合个数， 一次拟合可以产生保留个数为任意值时的预测结果。
许多常用模型有这样的性质，如：
- 提升方法，可以对不同迭代次数一次性输出对应的预测值。
- 正则化方法，如glmnet包的lasso等功能， 可以对不同惩罚系数值一次性输出对应的预测值。
- 多元适应回归样条(MARS, multivariate adaptive regression splines)， 调节参数是要保留的项数， 这可以一次拟合对不同的项数同时输出相应的预测值。
tune扩展包在参数调优时， 会自动利用这些技术。
比如， 使用C5.0作为基础算法的提升法， 将迭代次数作为唯一调优参数时， 可以一次拟合就得到所有迭代次数对应的预测值：
c5_spec <-
boost_tree(trees = tune()) |>
set_engine("C5.0") |>
set_mode("classification")
set.seed(1307)
c5_spec |>
tune_grid(
class ~ .,
resamples = cell_folds,
grid = data.frame(trees = 1:100),
metrics = roc_res )## Warning: package 'C50' was built under R version 4.2.3## # Tuning results
## # 10-fold cross-validation
## # A tibble: 10 × 4
## splits id .metrics .notes
## <list> <chr> <list> <list>
## 1 <split [1817/202]> Fold01 <tibble [100 × 5]> <tibble [0 × 3]>
## 2 <split [1817/202]> Fold02 <tibble [100 × 5]> <tibble [0 × 3]>
## 3 <split [1817/202]> Fold03 <tibble [100 × 5]> <tibble [0 × 3]>
## 4 <split [1817/202]> Fold04 <tibble [100 × 5]> <tibble [0 × 3]>
## 5 <split [1817/202]> Fold05 <tibble [100 × 5]> <tibble [0 × 3]>
## 6 <split [1817/202]> Fold06 <tibble [100 × 5]> <tibble [0 × 3]>
## 7 <split [1817/202]> Fold07 <tibble [100 × 5]> <tibble [0 × 3]>
## 8 <split [1817/202]> Fold08 <tibble [100 × 5]> <tibble [0 × 3]>
## 9 <split [1817/202]> Fold09 <tibble [100 × 5]> <tibble [0 × 3]>
## 10 <split [1818/201]> Fold10 <tibble [100 × 5]> <tibble [0 × 3]>37.9.5.2 利用并行计算功能
tidymodels不自动利用并行计算功能， 而是需要在程序中指定并注册计算集群， 见37.6.4。
进行参数调优过程可以看成如下的两层循环：
# 外层循环，针对不同的分析、评估集对
for(rs in resamples){
# 提取分析、评估集对
# 进行预处理
# 内层循环，针对不同的超参数组合网格点
for(mod in configarations) {
# 在分析集上用mod超参数组合估计模型
# 在评估集上预测，并计算性能指标
}
}注意针对超参数组合网格点的循环是在内层循环而不是外层循环。
tune_grid()在并行计算时，
默认仅对外层循环并行化；
可以在tune_grid()中设置参数parallel_over = "everything"，
则可以对外层、内层循环的组合进行并行计算。
在预处理计算量不大时，
选用此选项可以充分利用并行计算节点数目，
比仅在外层循环并行化略有改进。
并行计算时要避免使用全局变量，
因为工人节点可能无法访问主控节点的全局变量值。
rlang扩展包的!!运算符可以在设置工作流程参数时将全局变量的值插入到函数调用中。
37.9.5.3 赛跑方法
网格点需要将所有的网格点和所有的分析、评估集对都计算完成， 才能评估不同超参数组合的性能。 如果能在中间就排除某些一定会表现不佳的模型， 就可以加速搜索过程。 在机器学习中这种方法称为“赛跑方法”。 可以不是对所有的分析、评估集对都进行计算， 而是先选一小部分分析、评估集对， 在这样的子集上比较所有的参数组合， 预先将表现特别差的组合删去， 然后再完成其它的分析、评估集对的计算。
程序示例：
library(finetune)
set.seed(1308)
mlp_sfd_race <-
mlp_wflow |>
tune_race_anova(
cell_folds,
grid = 20,
param_info = mlp_param,
metrics = roc_res,
control = control_race(
verbose_elim = TRUE)
)
show_best(mlp_sfd_race, n = 10)
#> # A tibble: 2 × 10
#> hidden_units penalty epochs num_comp .metric .estimator mean n std_err
#> <int> <dbl> <int> <int> <chr> <chr> <dbl> <int> <dbl>
#> 1 8 0.814 177 15 roc_auc binary 0.887 10 0.0103
#> 2 3 0.0402 151 10 roc_auc binary 0.885 10 0.00810
#> # ℹ 1 more variable: .config <chr>调优结果仅能显示没有被中途删除的最优组合。
37.10 超参数调优的迭代搜索法
当超参数个数较多时， 网格法需要大量的格子点， 否则在超参数空间中的布点就过于稀疏， 无法抓住最优点。 这时， 可以考虑迭代搜索法。 迭代搜索法的缺点是不一定找到全局近似最优。
各种优化算法都可以使用， 这一节仅介绍贝叶斯优化和模拟退火优化。 以支持向量机建模为例。
37.10.1 支持向量机
使用cells数据集， 使用径向基函数， 要调优的超参数是损失值和径向基的尺度参数\(\sigma\)。
径向基使用点积（向量内积）， 所以应该将自变量标准化。 还可以使用PCA方法压缩自变量维数， 这里为了将超参数个数控制为2个就没有使用PCA。
针对cells数据集建立食谱、SVM模型、工作流程：
svm_rec <-
recipe(class ~ ., data = cells) |>
# 将自变量分布变得较为对称
step_YeoJohnson(all_numeric_predictors()) |>
step_normalize(all_numeric_predictors())
svm_spec <-
# SVM，基于径向基(radial basis function)
svm_rbf(cost = tune(), rbf_sigma = tune()) |>
set_engine("kernlab") |>
set_mode("classification")
svm_wflow <-
workflow() |>
add_model(svm_spec) |>
add_recipe(svm_rec)生成交叉验证集：
查看两个超参数的取值范围：
## Cost (quantitative)
## Transformer: log-2 [1e-100, Inf]
## Range (transformed scale): [-10, 5]## Radial Basis Function sigma (quantitative)
## Transformer: log-10 [1e-100, Inf]
## Range (transformed scale): [-10, 0]微调\(\sigma\)的范围：
svm_param <-
svm_wflow |>
extract_parameter_set_dials() |>
update(
rbf_sigma = rbf_sigma(c(-7, -1)))后面的两个算法需要一些搜索点的性能指标已知， 所以先构建一个小的网格搜索：
roc_res <- metric_set(roc_auc) # 性能指标集合
set.seed(1401)
start_grid <-
svm_param |>
update(
cost = cost(c(-6, 1)),
rbf_sigma = rbf_sigma(c(-6, -4))
) |>
grid_regular(levels = 2)
set.seed(1402)
svm_initial <-
svm_wflow |>
tune_grid(
resamples = cell_folds,
grid = start_grid,
metrics = roc_res)
collect_metrics(svm_initial)## # A tibble: 4 × 8
## cost rbf_sigma .metric .estimator mean n std_err .config
## <dbl> <dbl> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 0.0156 0.000001 roc_auc binary 0.864 10 0.00864 Preprocessor1_Model1
## 2 2 0.000001 roc_auc binary 0.863 10 0.00867 Preprocessor1_Model2
## 3 0.0156 0.0001 roc_auc binary 0.863 10 0.00862 Preprocessor1_Model3
## 4 2 0.0001 roc_auc binary 0.866 10 0.00855 Preprocessor1_Model4网格点上的AUC值（mean列）差别不大，
将作为迭代法的初值。
37.10.2 贝叶斯优化
贝叶斯优化基于已经计算过的超参数组合上的目标函数值（性能指标）， 计算出下一个要尝试的超参数组合并用再抽样技术计算其性能指标， 如此进行直到达到最大迭代次数或者不再有改进为止。
贝叶斯优化方法依赖于目标函数的模型与下一个尝试点的推荐算法。
37.10.2.1 高斯过程模型
高斯过程假设所有的观测值服从多元正态分布。 可以用均值函数和协方差函数描述联合分布。 协方差函数的常用模型为 \[ \text{Cov}(y_i, y_j) = \exp\left( -\frac{1}{2} | \boldsymbol x_i - \boldsymbol x_j |^2 \right) + \sigma_{ij}^2 . \] 这个模型意味着， 当超参数(\(\boldsymbol x_i\))组合距离增加时， 对应的性能指标的协方差以负指数速度衰减。 \(\sigma_{ii}=0\)。
使用这样一个概率模型的好处是， 对未来的\(\boldsymbol x_{n+1}\)上的\(y_{n+1}\)的预测， 可以是预测其高斯分布， 用均值和方差限定， 而不是仅有点预测或者区间预测。
对于本例， 已经有\(n=4\)个初始观测， 于是可以拟合高斯过程模型， 预测给出第5个点， 计算其性能指标， 重新拟合高斯过程模型， 预测给出第6个点， 如此重复。
37.10.2.2 购置函数
在估计了高斯过程模型以后， 如何产生推荐的下一个超参数组合点？ 可以在一个较大网格超参数上预测对应的性能指标的均值和方差， 并按某种方法在均值大和方差小之间权衡， 给出下一个点的推荐值。 用来权衡的函数称为购置函数。 在均值和方差之间权衡， 可以看成是偏向探索新区域和在旧有熟悉区域择优的权衡。 已知点较稀疏的新区域， 其预测分布的方差较大， 需要布点查看是否较优； 已知点较密集的熟悉区域， 其预测分布的方差较小， 可以用均值来择优。
一种常用方法是用改进量期望值作为选择依据， 仅考虑候选点比当前最优值可能的改进量， 而忽略其低于当前最优值的分布， 这样方差大的候选点有可能因此获胜。 这是tidymodels的默认设置。
37.10.2.3 用tune_bayes()函数进行迭代搜索
tune_bayes()函数基于贝叶斯优化方法进行迭代搜索，
它的主要输入仍为工作流程、再抽样对象，
除此以外还有如下选项：
iter规定最大迭代次数；initial可以是网格搜索的结果用作初值， 也可以指定一个整数， 指定用空间填充设计进行非规则网格搜索的点数， 还可以指定一个赛跑函数。objective指定在均值和方差之间权衡的购置函数， 默认为改进量期望值exp_improve()。param_info()，给出待优化超参数集合， 包括范围、变换。 可用用缺省值。
使用control参数输入一些迭代搜索设置，
用control_bayes()生成这些设置。
对cells数据集， SVM径向基模型， 用贝叶斯方法对两个超参数调优：
library(doParallel)
# 产生一个集群对象并注册
cl <- makePSOCKcluster(4)
registerDoParallel(cl)
ctrl <- control_bayes(verbose = TRUE)
set.seed(1403)
svm_bo <-
svm_wflow |>
tune_bayes(
# 再抽样对象，十折交叉验证数据集
resamples = cell_folds,
# AUC作为性能指标
metrics = roc_res,
# 初值是用网格搜索产生的4个点
initial = svm_initial,
# 限制了范围的两个超参数
param_info = svm_param,
# 最大迭代次数
iter = 25,
control = ctrl
)
# 取消集群对象
stopCluster(cl)注意， 迭代法的两层循环和网格法的两层循环不同。 网格法的外层循环是不同的再抽样折， 内层循环是不同的超参数组合； 迭代法的外层循环是顺序变化的超参数组合， 内层循环是不同的再抽样折。
上面的程序用了verbose选项，
输出内容过于详细，
没有显示这些输出。
从输出可以看出，
初始的AUC值为\(0.8659\)，
有些迭代步改进了指标值，
有些步并没有改进指标值。
查看调优结果：
# A tibble: 5 × 9
cost rbf_sigma .metric .estimator mean n std_err .config .iter
<dbl> <dbl> <chr> <chr> <dbl> <int> <dbl> <chr> <int>
1 16.9 0.00267 roc_auc binary 0.899 10 0.00804 Iter13 13
2 15.2 0.00283 roc_auc binary 0.899 10 0.00810 Iter16 16
3 13.5 0.00262 roc_auc binary 0.898 10 0.00780 Iter15 15
4 16.1 0.00227 roc_auc binary 0.898 10 0.00787 Iter17 17
5 13.6 0.00316 roc_auc binary 0.898 10 0.00822 Iter12 12最优值为\(0.8985\)。
按迭代次序对不同折的AUC指标作图：
可以看出在前10次迭代就已经接近了最优指标值。
显示各次迭代的超参数点变化，这是每个超参数的变化序列：
37.10.3 模拟退火优化
模拟退火算法是一种全局优化算法， 它避免落入局部最优的做法是有一定概率跳到目标函数变差的位置。
算法基于超参数空间内的随机游动。 用随机游动在当前点局域内产生下一个候选点， 如果目标函数（性能指标）改善， 则使用下一候选点； 如果目标函数值变差， 则以一定正概率使用此候选点， 结果越差， 接受此候选点的概率也越差； 随着迭代步数的增加， 接受较差候选点的概率也随之下降。
cooling_coef控制接受较差点的概率，
此参数越小则越容易接受较差的点。
算法中如果连续出现较差的点， 则可以退回到已发现的最优点重新开始。
用finetune包的tune_sim_anneal()函数实现超参数调优的模拟退火方法。
对cells数据集，
进行参数调优：
library(finetune)
library(doParallel)
# 产生一个集群对象并注册
cl <- makePSOCKcluster(4)
registerDoParallel(cl)
ctrl_sa <- control_sim_anneal(
verbose = TRUE,
# 如果10步都是变差，就从已有的最优点重新开始
no_improve = 10L)
set.seed(1404)
svm_sa <-
svm_wflow |>
tune_sim_anneal(
# 指定再抽样对象，即cells的十折交叉验证数据集
resamples = cell_folds,
# 指定优化的性能指标，是AUC
metrics = roc_res,
# 用一个网格优化结果作为初值
initial = svm_initial,
# 指定要优化的参数、范围、变换
param_info = svm_param,
# 超过最大迭代步数就停止
iter = 50,
# 算法控制
control = ctrl_sa
)
# 取消集群对象
stopCluster(cl)程序中用了verbose选项详细显示迭代过程，
这些显示略去。
最后得到的AUC最优值是\(0.8985\)，
与贝叶斯搜索方法AUC结果相同。