5 用parsnip拟合模型

5.1 AME数据集

parsnip包将R语言中不同的建模方法封装成统一的接口，使得用户可以快速、轻松地拟合不同类型的模型。我们将使用parsnip包拟合ames数据集中的回归模型，以房屋售价预测为例。

Note

tidymodels为不同的模型提供了统一的语法规则：

根据数学结构指定模型的类型，如线性回归、随机森林、KNN等。
模型的计算引擎(engine)，即指定建模使用的R包。这些软件包本身就是模型，parsnip使用这些软件包作为建模引擎，从而提供了一致的接口。
必要时指定模型的模式(mode)，即模型是回归还是分类。如果某个计算引擎只能处理其中一种类型，那在指定引擎的同时，也同时指定了模型的模式。

library(tidymodels)
library(modeldata)

# 载入数据
data(ames)
glimpse(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…

# 将因变量转换为常用的对数值
ames <- ames %>%
  mutate(Sale_Price = log10(Sale_Price))

# 划分训练集和测试集
set.seed(101)
ames_split <- initial_split(ames, prop = 0.80)
ames_train <- training(ames_split)
ames_test <- testing(ames_split)

5.1.1 线性回归模型示例

tidymodels_prefer() # 此函数用于设置 tidymodels 包的首选项，确保在多包环境中优先使用 tidymodels 的功能。
lm_model <-
  linear_reg() %>% # 指定模型类型为线性回归
  set_engine("lm") # 指定计算引擎为 lm

parsnip::translate(lm_model) # 将代码转换为软件语法包的详细信息

Linear Regression Model Specification (regression)

Computational engine: lm 

Model fit template:
stats::lm(formula = missing_arg(), data = missing_arg(), weights = missing_arg())

下面将数据、变量和模型联系起来，并在训练集上拟合模型，我们简单的探究房屋价格与经纬度的关系。

lm_form_fit <-
  lm_model %>%
  fit(Sale_Price ~ Longitude + Latitude, data = ames_train)
lm_form_fit

parsnip model object


Call:
stats::lm(formula = Sale_Price ~ Longitude + Latitude, data = data)

Coefficients:
(Intercept)    Longitude     Latitude  
   -301.337       -1.986        2.869

fit_xy()函数是fit()函数的变体版，适用于自变量为数据框，因变量为向量的情况。

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.869

5.1.2 随机森林模型

parsnip包对同一模型的不同计算引擎提供了统一的接口，使得用户可以轻松地切换计算引擎。

在随机森林模型中，我们需要设置每次随机选择的自变量个数、合并树的个数、每棵树的大小等参数，不同的随机森林扩展包(如randomForest、ranger，sparklyr等)对这些设置使用了不同的参数名。parsnip中的rand_forest()函数则使用了同一套参数名，使得用户可以轻松地切换计算引擎。如表 5.1 所示。

表 5.1: parsnip中随机森林参数

参数类型	parsnip
sampled predictors	mtry(每次分叉可以从多少个自变量中选择)
trees	trees(树的数目)
data points to split	min_n(每个叶节点包含的观测数下限，低于则不能继续拆分)

rf_model <-
  rand_forest(
    trees = 1000,
    min_n = 5
  ) %>%
  set_engine("ranger") %>%
  set_mode("regression")

parsnip::translate(rf_model)

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))

有一些参数是各个引擎特有的，这在parsnip中会自动设置一些合适的值，也可以在set_engine()函数中添加，如set_engine("ranger", verbose = TRUE)中的verbose = TRUE参数。

5.1.3 提取拟合结果

lm_form_fit %>%
  extract_fit_engine() %>%
  coef() # 回归系数

(Intercept)   Longitude    Latitude 
-301.337454   -1.986050    2.868594

lm_xy_fit %>%
  extract_fit_engine() %>%
  vcov() # 协方差矩阵

            (Intercept)     Longitude      Latitude
(Intercept)  203.734986  1.5807903348 -1.3252325337
Longitude      1.580790  0.0168109489 -0.0001562159
Latitude      -1.325233 -0.0001562159  0.0311793054

lm_xy_fit %>%
  extract_fit_engine() %>%
  summary() # 模型的整体信息


Call:
stats::lm(formula = ..y ~ ., 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

这些信息提取的函数都有各自用途，但返回的结果类型不统一，且不同引擎对信息提取函数的支持也不相同。我们需要将模型拟合的结果转换为统一的格式，以便于后续的分析和可视化，broom包中的tidy()函数正是为此而生。

lm_form_fit %>%
  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

5.1.4 预测

parsnip包中的predict()函数是以模型拟合的结果为输入，以new_data(即测试集)为输入的预测数据，返回预测值，输出的结果非常整洁统一：

输出结果总是一个tibble格式的数据框。
输出结果的行数与new_data的行数相同，如果输入数据有缺失值，则相应的结果为NA。
输出结果的变量名(列名)总是有明确的含义，具体规定如下：
- 因变量为数值型，要做预测，选type = "numeric"，结果变量为.pred。
- 因变量为数值型，要计算期望值的置信区间，选type = "conf_int"，结果变量为.pred_lower和.pred_upper。
- 因变量为数值型，要计算个别值的预测区间，选type = "pred_int"，结果变量为.pred_lower和.pred_upper。
- 因变量为类别型，要预测类别，选type = "class"，结果变量为.pred_class。
- 因变量为类别型，要计算每个类别的后验概率，选type = "prob"，结果变量为.pred_{类别水平名}。

lm_form_fit %>%
  predict(new_data = ames_test %>% slice(1:5))

# A tibble: 5 × 1
  .pred
  <dbl>
1  5.29
2  5.29
3  5.24
4  5.24
5  5.23

输出的结果为.pred列，表示预测的房屋价格，与输入数据的行一一对应。为增强输出结果的可读性，我们将输出的预测结果和预测数据进行合并。

ames_test_sm <- ames_test %>%
  slice(1:5)
ames_test_sm %>%
  select(Sale_Price) %>%
  bind_cols(
    predict(
      lm_form_fit,
      new_data = ames_test_sm
    )
  ) %>%
  bind_cols(
    predict(
      lm_form_fit,
      new_data = ames_test_sm,
      type = "conf_int"
    )
  )

# A tibble: 5 × 4
  Sale_Price .pred .pred_lower .pred_upper
       <dbl> <dbl>       <dbl>       <dbl>
1       5.28  5.29        5.28        5.30
2       5.29  5.29        5.28        5.30
3       5.26  5.24        5.23        5.26
4       4.94  5.24        5.23        5.25
5       5.11  5.23        5.23        5.24

这样统一的输出结果的好处很多，对于同一个问题，只要在指定模型时选择不同的引擎和模式即可。例如我们可以用决策树回归进行预测，代码与上述代码类似，只需将lm_model换成tree_model即可。

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,
    new_data = 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.15

5.1.5 支持的模型

tidymodels支持的模型和扩展包的详细列表见:

https://www.tidymodels.org/find/parsnip/

另外，与判别分析有关的扩展包支持放在了descrim扩展包中。

在RStudio中使用时， parsnip提供了一个RStudio addin，可以在Addins下拉菜单中找到，可以用“Generate parsnip model specification”功能生成模型设置代码。