如何用R建模GBM？

本文使用梯度提升机(GBM)来找出是什么驱动了自行车租赁行为。

【需要加入数据人圈子，请加微信：luqin360】

与其他决策树算法相比，GBM是独一无二的，因为它按顺序构建模型，给予那些在先前模型中预测不佳的情况下的较高权重，从而逐步提高准确性，而不是简单地采用随机森林算法的所有模型的平均值。利用迭代方式减少错误以形成所谓的最终模型，GBM是用于分类和回归问题的有效且强大的算法。

R做GBM可以很好地获得探索图，包括参数贡献和部分依赖图，其提供模型中特征的值的效果的可视化。

加载包

完成这个分析所需要的包如下：

1library(rsample)      
2library(caret)        
3library(ggthemes)
4library(scales)
5library(wesanderson)
6library(tidyverse)
7library(gbm)
8library(Metrics)
9library(here)

运行先前的脚本以获得和处理数据集。

1base::source("code/01-import.R")
2base::source("code/02.3-wrangle.R")
3base::source("code/03.3-visualize.R") # we don't neeed to run 03-vizualize.R,
4# because it does not change the underlying data structure, but we wantto view 
5# the EDA outputs

第一部分所学习到的

探索性数据分析显示，在一定温度(~20摄氏度)下，自行车租金呈现下降趋势。

如图

而且节假日的租赁比非节假日要低。

因此，现在我们已经对BikeData中的变量有了一定的了解，我们可以使用广义增强回归模型(gbm)包对自行车租赁数据进行建模。

数据集

1BikeData %>% glimpse()

感兴趣的变量集

首先，我们想要为我们的模型构建一个数据框，其中包括我们的输出变量cnt(包括临时和注册自行车租赁总数)和我们想要帮助解释的特征。

 1BikeDataModel <- BikeData %>% dplyr::select(cnt,
 2                season_fct,
 3                yr_fct,
 4                month_fct,
 5                holiday_fct,
 6                weekday_fct,
 7                workingday_fct,
 8                weathersit_fct,
 9                temp,
10                atemp,
11                hum,
12                windspeed)
13BikeDataModel %>% dplyr::glimpse(78)

测试集和训练集划分

我们希望构建代表原始自行车数据集的训练和测试数据集。为此，我们将随机选择两个数据子集。我们还将指定采样过程，因此我们在BikeTest和BikeTrain数据集中获得相同比例的自行车租赁（cnt）（我们将数据随机分配到训练和测试集中，分成70％/ 30％）。

BikeTrain拥有我们将用于构建模型并展示其性能的数据。但是，因为我们的目标是尽可能使用“真实世界数据”进行泛化，我们将在我们的模型没有看到的数据（即BikeTest数据）上测试模型。

通过测试和训练数据集，我们可以1）构建并选择最佳模型，然后2）然后使用“新鲜”数据评估最终模型的性能。

1set.seed(123)
2BikeSplit <- initial_split(BikeDataModel, prop = .7)
3BikeTrain <- training(BikeSplit)
4BikeTest  <- testing(BikeSplit)

我们可以调用gbm函数并选择一些参数，包括交叉折叠验证（cv.folds = 10）。

交叉折叠验证随机将我们的训练数据划分为大小相对相等的k组。我们的模型将拟合排除一个折叠的所有组合数据集，然后使用保持即估计拟合的模型的误差。差。

每组用于测量模型误差，并在各组中计算平均值。使用R http://topepo.github.io/caret/index.html中 的caret包，可以调整shrinkage，interaction.depth，n.minobsinnode和n.trees的模型精度。

 1# model
 2set.seed(123)
 3bike_fit_1 <- gbm::gbm(cnt ~., 
 4             # the formula for the model (recall that the period means, "all 
 5             # variables in the data set")
 6             data = BikeTrain, 
 7             # data set
 8             verbose = TRUE, 
 9             # Logical indicating whether or not to print
10             #  out progress and performance indicators
11             shrinkage = 0.01, 
12             # a shrinkage parameter applied to each tree in the expansion. 
13             # Also known as the learning rate or step-size reduction; 0.001 
14             # to 0.1 usually work, but a smaller learning rate typically 
15             # requires more trees.
16             interaction.depth = 3, 
17             # Integer specifying the maximum depth of each tree (i.e., the 
18             # highest level of variable interactions allowed). A value of 1 
19             # implies an additive model, a value of 2 implies a model with up
20             #  to 2-way interactions
21             n.minobsinnode = 5,
22             # Integer specifying the minimum number of observations in the 
23             # terminal nodes of the trees. Note that this is the actual number 
24             # of observations, not the total weight.
25             n.trees = 5000, 
26             # Integer specifying the total number of trees to fit. This is 
27             # equivalent to the number of iterations and the number of basis 
28             # functions in the additive expansion.
29             cv.folds = 10
30             # Number of cross-validation folds to perform. If cv.folds>1 then
31             # gbm, in addition to the usual fit, will perform a 
32             # cross-validation, calculate an estimate of generalization error
33             #  returned in cv.error
34             )

1Distribution not specified, assuming gaussian ...
2Iter   TrainDeviance   ValidDeviance   StepSize   Improve
3     1  3887192.4522             nan     0.0100 53443.4349
4     2  3835212.0554             nan     0.0100 51509.5429
5     3  3781613.7580             nan     0.0100 46152.5706
6     4  3734439.6068             nan     0.0100 48076.6689
7     5  3684049.7222             nan     0.0100 47545.5212
8     6  3636533.0042             nan     0.0100 48434.4787
9     7  3589347.3456             nan     0.0100 49416.0562

注意:gbm::gbm()函数创建了大量输出(我们包括verbose = TRUE参数)。

存储和探索

现在我们可以检查模型对象(bike_fit_1)，从GBM产生的最优学习者数量和错误率开始。
我们可以使用gbm::gbm.perf()函数来查看每个学习者数量的错误率。

在下面的可视化中，我们可以看到蓝色的线表示使用交叉验证(cv)的最佳树数。GBM可能对过拟合很敏感，因此在我们的估计中使用= “cv”方法可以防止这种情况。
我们可以看到最优的树数是1794棵

做预测

现在我们可以使用带有测试集的predict（）函数和基于perf.gbm1估计的最佳树数预测我们的自行车租赁。

RMSE =均方根误差（RMSE）用于测量模型中的预测误差。顾名思义，这些错误通过平方来加权。RMSE的解释也非常简单，因为它与我们的结果变量具有相同的单位。

Metrics :: rmse计算两个数值向量之间的均方根误差，因此我们将使用它来比较我们的预测值和残差来计算模型的误差。

 1bike_prediction_1 <- stats::predict(
 2                           # the model from above
 3                          object = bike_fit_1, 
 4                          # the testing data
 5                          newdata = BikeTest,
 6                          # this is the number we calculated above
 7                          n.trees = perf_gbm1)
 8rmse_fit1 <- Metrics::rmse(actual = BikeTest$cnt, 
 9                           predicted = bike_prediction_1)
10print(rmse_fit1)

GBM提供部分依赖图，以探索模型中的特征与结果之间的相关性。例如，您可以在下图中看到环境温度与自行车租赁数量的增加相关，直到接近35度，此时骑车者不太可能租用自行车。

1gbm::plot.gbm(bike_fit_1, i.var = 9)

同样，我们可以看看自行车租赁的两个特征之间的交互作用。下面我们可以看到，尽管风速很大，但骑车者更有可能在周一之后租一辆自行车。

利用GBM提供的相对影响，我们可以直观地看到不同特征对自行车租赁预测的影响。首先，我们可以总结我们的模型，然后使用gbm::summary.gbm()将这些数据分配给tibble::as_tibble()函数。

1# summarize model
2BikeEffects <- tibble::as_tibble(gbm::summary.gbm(bike_fit_1, 
3                                         plotit = FALSE))
4BikeEffects %>% utils::head()

1 ## # A tibble: 6 x 2
2 ##   var            rel.inf
3 ##               
4 ## 1 atemp            35.3 
5 ## 2 yr_fct           22.1 
6 ## 3 month_fct         9.81
7 ## 4 hum               8.41
8 ## 5 season_fct        8.35
9 ## 6 weathersit_fct    3.91

这将创建一个新的数据集，其中包含var(一个包含模型中的变量的因素变量)和rel.inf(每个变量对模型预测的相对影响)。
然后，我们可以使用ggpplot和包含模型摘要(BikeEffects)的新数据框架来绘制10大特征。

 1# plot effects
 2BikeEffects %>% 
 3  # arrange descending to get the top influencers
 4  dplyr::arrange(desc(rel.inf)) %>%
 5  # sort to top 10
 6  dplyr::top_n(10) %>%
 7  # plot these data using columns
 8  ggplot(aes(x = forcats::fct_reorder(.f = var, 
 9                                      .x = rel.inf), 
10             y = rel.inf, 
11             fill = rel.inf)) +
12  geom_col() +
13  # flip
14  coord_flip() +
15  # format
16  scale_color_brewer(palette = "Dark2") +
17  theme_fivethirtyeight() +
18  theme(axis.title = element_text()) + 
19  xlab('Features') +
20  ylab('Relative Influence') +
21  ggtitle("Top 10 Drivers of Bike Rentals")

这让我们知道用于选择的变量：通过rel.inf选择

通过预测这些值并绘制差异，我们可以将我们预测的分布与实际自行车租赁进行比较。

 1# Predicted bike rentals
 2BikeTest$predicted <- base::as.integer(predict(bike_fit_1, 
 3                                         newdata = BikeTest, 
 4                                         n.trees = perf_gbm1))
 5
 6# plot predicted v actual
 7ggplot(BikeTest) +
 8  geom_point(aes(y = predicted, 
 9                 x = cnt, 
10                 color = predicted - cnt), alpha = 0.7) +
11  # add theme
12  theme_fivethirtyeight() +
13  # strip text
14  theme(axis.title = element_text()) + 
15  # add format to labels
16  scale_x_continuous(labels = comma) +
17  scale_y_continuous(labels = comma) +
18  # add axes/legend titles
19  scale_color_continuous(name = "Predicted - Actual", 
20                         labels = comma) +
21  ylab('Predicted Bike Rentals') +
22  xlab('Actual Bike Rentals') +
23  ggtitle('Predicted vs Actual Bike Rentals') 

我们可以看到我们的模型在预测自行车租赁方面做得相当不错。

我们学到了什么？

我们了解到环境温度是预测自行车租赁的最大影响因素，当温度达到~35度时，租赁数量会下降。我们还可以看到假期（或非假期）并不是预测自行车租赁的影响因素。

原文链接：
http://www.storybench.org/tidytuesday-bike-rentals-part-2-modeling-with-gradient-boosting-machine/

数据人才（ID：datarencai）

（一个帮助数据人才找工作的公众号，

也分享数据人才学习和生活的有趣事情。）

内容推荐