# How to model with gradient boosting machine in R

The tutorial is part 2 of our #tidytuesday post from last week, which explored bike rental data from Washington, D.C. Check it out here. The following tutorial will use a gradient boosting machine (GBM) to figure out what drives bike rental behavior.

GBM is unique compared to other decision tree algorithms because it builds models sequentially with higher weights given to those cases that were poorly predicted in previous models, thus improving accuracy incrementally instead of simply taking an average of all models like a random forest algorithm would. By reducing the error iteratively to produce what will become the final model, GBM is an efficient and powerful algorithm for classification and regression problems.

Implementing GBM in R allows for a nice selection of exploratory plots including parameter contribution, and partial dependence plots which provide a visual representation of the effect across values of a feature in the model.

# Load the packages

The packages below are needed to complete this analysis.

# packages -------------------------------------------------------------- library(rsample) library(caret) library(ggthemes) library(scales) library(wesanderson) library(tidyverse) library(gbm) library(Metrics) library(here)

Run previous scripts to import and wrangle data.

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

## What did we learn in part 1?

The exploratory data analysis showed us that bike rentals seem to drop off at a certain temperature (~20˚C).

ggRentalsByTempAnd that rentals were lower on holidays compared to non-holidays.

So now that we have some understanding of the variables in `BikeData`

, we can use the Generalized Boosted Regression Models (`gbm`

) package to model the bike rental data.

## The data

BikeData %>% glimpse()

## Observations: 731 ## Variables: 30 ## $ weekday6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2… ## $ weekday_chr "Saturday", "Sunday", "Monday", "Tuesday", "Wednesda… ## $ weekday_fct Saturday, Sunday, Monday, Tuesday, Wednesday, Thursd… ## $ holiday 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0… ## $ holiday_chr "Non-Holiday", "Non-Holiday", "Non-Holiday", "Non-Ho… ## $ holiday_fct Non-Holiday, Non-Holiday, Non-Holiday, Non-Holiday, … ## $ season 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1… ## $ season_chr "Spring", "Spring", "Spring", "Spring", "Spring", "S… ## $ season_fct Spring, Spring, Spring, Spring, Spring, Spring, Spri… ## $ workingday 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1… ## $ workingday_chr "Non-Working Day", "Non-Working Day", "Working Day",… ## $ workingday_fct Non-Working Day, Non-Working Day, Working Day, Worki… ## $ month_ord Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Jan, Ja… ## $ month_fct January, January, January, January, January, January… ## $ yr 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0… ## $ yr_chr "2011", "2011", "2011", "2011", "2011", "2011", "201… ## $ yr_fct 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011… ## $ weathersit 2, 2, 1, 1, 1, 1, 2, 2, 1, 1, 2, 1, 1, 1, 2, 1, 2, 2… ## $ weathersit_chr "Clouds/Mist", "Clouds/Mist", "Good", "Good", "Good"… ## $ weathersit_fct Clouds/Mist, Clouds/Mist, Good, Good, Good, Good, Cl… ## $ instant 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1… ## $ dteday 2011-01-01, 2011-01-02, 2011-01-03, 2011-01-04, 201… ## $ mnth 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1… ## $ temp 8, 9, 1, 1, 2, 1, 1, 0, -1, 0, 0, 0, 0, 0, 2, 2, 0, … ## $ atemp 28.36325, 28.02713, 22.43977, 23.21215, 23.79518, 23… ## $ hum 80, 69, 43, 59, 43, 51, 49, 53, 43, 48, 68, 59, 47, … ## $ windspeed 10, 16, 16, 10, 12, 6, 11, 17, 24, 14, 8, 20, 20, 8,… ## $ casual 331, 131, 120, 108, 82, 88, 148, 68, 54, 41, 43, 25,… ## $ registered 654, 670, 1229, 1454, 1518, 1518, 1362, 891, 768, 12… ## $ cnt 985, 801, 1349, 1562, 1600, 1606, 1510, 959, 822, 13…

## Variables of interest

First we want build a data frame for our model, which includes our outcome variable `cnt`

(the ‘*count of total rental bikes including both casual and registered*’) and the features we want to help explain:

`season_fct`

= Factor w/ 4 levels, “Spring”, “Summer”, “Fall”, “Winter”

`yr_fct`

= Factor w/ 2 levels, “2011”, “2012”

`month_fct`

= Factor w/ 12 levels, “January”, “February”, “March”, “April”, “May”, “June”, “July”, “August”, “September”, “October”, “November”, “December”

`holiday_fct`

= Factor w/ 2 levels “Non-Holiday”, “Holiday”

`weekday_fct`

= Factor w/ 7 levels, “Sunday”, “Monday”, “Tuesday”, “Wednesday”, “Thursday”, “Friday”,”Saturday”

`workingday_fct`

= Factor w/ 2 levels, “Non-Working Day”, “Working Day”

`weathersit_fct`

= Factor w/ 3 levels, “Good”, “Clouds/Mist”, “Rain/Snow/Storm”

`temp`

= int [1:731] 8 9 1 1 2 1 1 0 -1 0

`atemp`

= num [1:731] 28.4 28 22.4 23.2 23.8 ..

`hum`

= int [1:731] 80 69 43 59 43 51 49 53 43 48 …

`windspeed`

= int [1:731] 10 16 16 10 12 6 11 17 24 14 …

BikeDataModel <- BikeData %>% dplyr::select(cnt, season_fct, yr_fct, month_fct, holiday_fct, weekday_fct, workingday_fct, weathersit_fct, temp, atemp, hum, windspeed) BikeDataModel %>% dplyr::glimpse(78)

## Observations: 731 ## Variables: 12 ## $ cnt985, 801, 1349, 1562, 1600, 1606, 1510, 959, 822, 13… ## $ season_fct Spring, Spring, Spring, Spring, Spring, Spring, Spri… ## $ yr_fct 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011… ## $ month_fct January, January, January, January, January, January… ## $ holiday_fct Non-Holiday, Non-Holiday, Non-Holiday, Non-Holiday, … ## $ weekday_fct Saturday, Sunday, Monday, Tuesday, Wednesday, Thursd… ## $ workingday_fct Non-Working Day, Non-Working Day, Working Day, Worki… ## $ weathersit_fct Clouds/Mist, Clouds/Mist, Good, Good, Good, Good, Cl… ## $ temp 8, 9, 1, 1, 2, 1, 1, 0, -1, 0, 0, 0, 0, 0, 2, 2, 0, … ## $ atemp 28.36325, 28.02713, 22.43977, 23.21215, 23.79518, 23… ## $ hum 80, 69, 43, 59, 43, 51, 49, 53, 43, 48, 68, 59, 47, … ## $ windspeed 10, 16, 16, 10, 12, 6, 11, 17, 24, 14, 8, 20, 20, 8,…

## The testing and training split

We want to build training and testing data sets that are representative of the original `bike`

data set. To achieve this, we will randomly select observations for two subsets of data. We’ll also specify the sampling process, so we get an equal proportion of bike rentals (`cnt`

) in our `BikeTest`

and `BikeTrain`

data sets (we’ll randomize our data into training and test sets with a 70% / 30% split).

The `BikeTrain`

has the data we will use to build a model and demonstrate it’s performance. However, because our goal is to generalize as much as possible using ‘real world data’, we’ll be testing the model on data our model *hasn’t* seen (i.e. the `BikeTest`

data).

Having testing and training data sets allows us to 1) build and select the best model, and 2) then assess the final model’s performance using ‘fresh’ data.

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

We can call the `gbm`

function and select a number of parameters, including cross-fold validation (`cv.folds = 10`

).

Cross-fold validation randomly divides our training data into `k`

sets that are relatively equal in size. Our model will be fit using all the sets with the exclusion of the first fold. The model error of the fit is estimated with the hold-out sets.

Each set is used to measure the model error and an average is calculated across the various sets. `shrinkage`

, `interaction.depth`

, `n.minobsinnode`

and `n.trees`

can be adjusted for model accuracy using the caret package in R http://topepo.github.io/caret/index.html.

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

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

**NOTE:** The `gbm::gbm()`

function creates a lot of output (we included the `verbose = TRUE`

argument).

## Store and explore

Now we can inspect the model object (`bike_fit_1`

) beginning with the optimal number of learners that GBM has produced and the error rate.

We can use the `gbm::gbm.perf()`

function to see the the error rate at each number of learners.

# model performance perf_gbm1 = gbm.perf(bike_fit_1, method = "cv")

In the visualization below we can see that the blue line represents the optimal number of trees with our cross validation (`cv`

). GBM can be sensitive to over-fitting, so using the `method = "cv"`

in our estimate protects against this.

We can see that the optimal number of trees is `1794`

## Make predictions

Now we can predict our bike rentals using the `predict()`

function with our test set and the optimal number of trees based on our `perf.gbm1`

estimate.

`RMSE`

= The root mean squared error (`RMSE`

) is used to measure the prediction error in our model(s). As the name suggests, these errors are weighted by means of squaring them. The RMSE is also pretty straightforward to interpret, because it's in the same units as our outcome variable. Additional attractive qualities include the fact equally captures the overestimates and underestimates, and the misses are penalized according to their relative size.

`Metrics::rmse`

computes the root mean squared error between two numeric vectors, so we will use it to compare our predicted values with the residuals to calculate the error of our model.

bike_prediction_1 <- stats::predict( # the model from above object = bike_fit_1, # the testing data newdata = BikeTest, # this is the number we calculated above n.trees = perf_gbm1) rmse_fit1 <- Metrics::rmse(actual = BikeTest$cnt, predicted = bike_prediction_1) print(rmse_fit1)

GBM offers partial dependency plots to explore the correlations between a feature in our model and our outcome. For example, you can see in the graph below that ambient temperature is associated with increased numbers of bike rentals until close to 35 degrees when riders tend to be less likely to rent a bike.

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

Similarly we can look at the interaction of two features on bike rentals. Below we can see that riders are more likely to rent a bike after Monday, despite wind speed.

gbm::plot.gbm(bike_fit_1, i.var = c(5, 11))

We can visualize the impact of different features on predicting bike rentals using the relative influence provided by GBM. First we can summarize our model then assign these data to a `tibble`

using `gbm::summary.gbm()`

and passing this to the `tibble::as_tibble()`

function.

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

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

This creates a new data set with `var`

, a factor variable with the variables in our model, and `rel.inf`

, the relative influence each variable had on our model predictions.

We can then plot the top ten features by impact using `ggpplot`

and our new data frame containing our model summary (`BikeEffects`

).

# plot effects BikeEffects %>% # arrange descending to get the top influencers dplyr::arrange(desc(rel.inf)) %>% # sort to top 10 dplyr::top_n(10) %>% # plot these data using columns ggplot(aes(x = forcats::fct_reorder(.f = var, .x = rel.inf), y = rel.inf, fill = rel.inf)) + geom_col() + # flip coord_flip() + # format scale_color_brewer(palette = "Dark2") + theme_fivethirtyeight() + theme(axis.title = element_text()) + xlab('Features') + ylab('Relative Influence') + ggtitle("Top 10 Drivers of Bike Rentals")

This let's us know what is being used for selection: `Selecting by rel.inf`

We can visualize the distribution of our predicted compared with actual bike rentals by predicting these values and plotting the difference.

# Predicted bike rentals BikeTest$predicted <- base::as.integer(predict(bike_fit_1, newdata = BikeTest, n.trees = perf_gbm1)) # plot predicted v actual ggplot(BikeTest) + geom_point(aes(y = predicted, x = cnt, color = predicted - cnt), alpha = 0.7) + # add theme theme_fivethirtyeight() + # strip text theme(axis.title = element_text()) + # add format to labels scale_x_continuous(labels = comma) + scale_y_continuous(labels = comma) + # add axes/legend titles scale_color_continuous(name = "Predicted - Actual", labels = comma) + ylab('Predicted Bike Rentals') + xlab('Actual Bike Rentals') + ggtitle('Predicted vs Actual Bike Rentals')

We can see that our model did a fairly good job predicting bike rentals.

## What did we learn?

We learned the ambient temperature is the largest influencer for predicting bike rentals and that rental numbers go down when the temperature reaches ~35 degrees. We can also see holiday (or non-holiday) was not much of an influencer for predicting bike rentals.

- How to model with gradient boosting machine in R - April 9, 2019
- Exploring bike rental behavior using R - April 2, 2019

This is a nice post for learning. Could you please show, how can we see/save the training period simulation data. Then we have the opportunity to calculate different matrix for training period also.

Thank you for this article. Bootstrapping is commonly used with big datasets to help remove the chance of a single training/testing sample skewing the overall result. How would we apply a loop to the model so as to generate the average result after n boostrapping repititions?

Hi Peter

You mention visualizations and graphs, but I only see the code. Is it just me or didn’t you include the graphics?