<jupyter_start><jupyter_text>Gradient Boosting Machines  * * *![Alt text](./images/shrubs.jpg "GBMs")* * *Image Source: brucecompany.com  Introduction[Gradient boosting](https://en.wikipedia.org/wiki/Gradient_boosting) is a machine learning technique for regression and classification problems, which produces a prediction model in the form of an ensemble of weak prediction models, typically decision trees. It builds the model in an iterative fashion like other boosting methods do, and it generalizes them by allowing optimization of an arbitrary differentiable loss function.It is recommended that you read through the accompanying [Classification and Regression Trees Tutorial](decision-trees.ipynb) for an overview of decision trees.  HistoryBoosting is one of the most powerful learning ideas introduced in the last twenty years. It was originally designed for classification problems, but it can be extended to regression as well. The motivation for boosting was a procedure that combines the outputs of many "weak" classifiers to produce a powerful "committee."  A weak classifier (e.g. decision tree) is one whose error rate is only slightly better than random guessing. [AdaBoost](https://en.wikipedia.org/wiki/AdaBoost) short for "Adaptive Boosting", is a machine learning meta-algorithm formulated by [<NAME>](https://en.wikipedia.org/wiki/Yoav_Freund) and [<NAME>](https://en.wikipedia.org/wiki/Robert_Schapire) in 1996, which is now considered to be a special case of Gradient Boosting.  There are [some differences](http://stats.stackexchange.com/questions/164233/intuitive-explanations-of-differences-between-gradient-boosting-trees-gbm-ad) between the AdaBoost algorithm and modern Gradient Boosting.  In the AdaBoost algorithm, the "shortcomings" of existing weak learners are identified by high-weight data points, however in Gradient Boosting, the shortcomings are identified by gradients.The idea of gradient boosting originated in the observation by Leo Breiman that boosting can be interpreted as an optimization algorithm on a suitable cost function. Explicit regression gradient boosting algorithms were subsequently developed by <NAME>, simultaneously with the more general functional gradient boosting perspective of <NAME>, <NAME>, <NAME> and <NAME>. The latter two papers introduced the abstract view of boosting algorithms as iterative functional gradient descent algorithms. That is, algorithms that optimize a cost function over function space by iteratively choosing a function (weak hypothesis) that points in the negative gradient direction. This functional gradient view of boosting has led to the development of boosting algorithms in many areas of machine learning and statistics beyond regression and classification.In general, in terms of model performance, we have the following heirarchy:$$Boosting > Random \: Forest > Bagging > Single \: Tree$$  Boosting[Gradient boosting](https://en.wikipedia.org/wiki/Gradient_boosting) is a machine learning technique for regression and classification problems, which produces a prediction model in the form of an ensemble of weak prediction models, typically decision trees. It builds the model in a stage-wise fashion like other boosting methods do, and it generalizes them by allowing optimization of an arbitrary differentiable loss function.The purpose of boosting is to sequentially apply the weak classification algorithm to repeatedly modified versions of the data, thereby producing a sequence of weak classifiers $G_m(x)$, $m = 1, 2, ... , M$.  Stagewise Additive ModelingBoosting builds an additive model:$$F(x) = \sum_{m=1}^M \beta_m b(x; \gamma_m)$$where $b(x; \gamma_m)$ is a tree and $\gamma_m$ parameterizes the splits.  With boosting, the parameters, $(\beta_m, \gamma_m)$ are fit in a *stagewise* fashion.  This slows the process down, and overfits less quickly.  AdaBoost- AdaBoost builds an additive logistic regression model by stagewise fitting.- AdaBoost uses an exponential loss function of the form, $L(y, F(x)) = exp(-yF(x))$, similar to the negative binomial log-likelihood loss.- The principal attraction of the exponential loss in the context of additive modeling is computational; it leads to the simple modular reweighting- Instead of fitting trees to residuals, the special form of the exponential loss function in AdaBoost leads to fitting trees to weighted versions of the original data.![Alt text](./images/boosting_diagram.png "Adaboost M1")Source: Elements of Statistical Learning![Alt text](./images/adaboost_m1.png "Adaboost M1")Source: Elements of Statistical Learning  Gradient Boosting Algorithm- The idea of gradient boosting originated in the observation by Leo Breiman that boosting can be interpreted as an optimization algorithm on a suitable cost function.- Explicit regression gradient boosting algorithms were subsequently developed by <NAME> simultaneously with the more general functional gradient boosting perspective of <NAME>, <NAME>, <NAME> and <NAME>. - The latter two papers introduced the abstract view of boosting algorithms as iterative functional gradient descent algorithms. That is, algorithms that optimize a cost function over function space by iteratively choosing a function (weak hypothesis) that points in the negative gradient direction. - This functional gradient view of boosting has led to the development of boosting algorithms in many areas of machine learning and statistics beyond regression and classification.Friedman's Gradient Boosting Algorithm for a generic loss function, $L(y_i, \gamma)$:![Alt text](./images/friedman_gbm.png "GBM Algorithm")Source: Elements of Statistical Learning Loss Functions and Gradients![Alt text](./images/gbm_gradients_loss.png "GBM Algorithm")Source: Elements of Statistical LearningThe optimal number of iterations, T, and the learning rate, λ, depend on each other.  Stochastic GBM[Stochastic Gradient Boosting](https://statweb.stanford.edu/~jhf/ftp/stobst.pdf) (Friedman, 2002) proposed the stochastic gradient boosting algorithm that simply samples uniformly without replacement from the dataset before estimating the next gradient step. He found that this additional step greatly improved performance. Some implemenations of Stochastic GBM have both column and row sampling (per split and per tree) for better generalization.  XGBoost and H2O are two implemenations that have a per-tree column and row sampling.  This is reason that these implementations are popular among competitive data mining competitors (e.g. Kaggle).  Practical Tips- It's more common to grow shorter trees ("shrubs" or "stumps") in GBM than you do in Random Forest.- It's useful to try a variety of column sample (and column sample per tree) rates.  - Don't assume that the set of optimal tuning parameters for one implementation of GBM will carry over and also be optimal in a different GBM implementation.  Resources - [<NAME> - Gradient Boosting & Random Forests at H2O World 2014](https://www.youtube.com/watch?v=wPqtzj5VZus&index=16&list=PLNtMya54qvOFQhSZ4IKKXRbMkyLMn0caa) (YouTube) - [<NAME> - Data Science of GBM (2013)](http://www.slideshare.net/0xdata/gbm-27891077) (slides) - [<NAME> - Gradient Boosting Method and Random Forest at H2O World 2015](https://www.youtube.com/watch?v=9wn1f-30_ZY) (YouTube) - [<NAME> - Gradient Boosted Regression Trees in scikit-learn at PyData London 2014](https://www.youtube.com/watch?v=IXZKgIsZRm0) (YouTube) - [<NAME> and <NAME> - Gradient boosting machines, a tutorial](http://journal.frontiersin.org/article/10.3389/fnbot.2013.00021/full) (blog post) ***  GBM Software in RThis is not a comprehensive list of GBM software in R, however, we detail a few of the most popular implementations below: [gbm](https://cran.r-project.org/web/packages/gbm/index.html), [xgboost](https://cran.r-project.org/web/packages/xgboost/index.html) and [h2o](https://cran.r-project.org/web/packages/gamboostLSS/index.html).  The [CRAN Machine Learning Task View](https://cran.r-project.org/web/views/MachineLearning.html) lists the following projects as well.  The Hinge-loss is optimized by the boosting implementation in package [bst](https://cran.r-project.org/web/packages/bst/index.html). Package [GAMBoost](https://cran.r-project.org/web/packages/GAMBoost/index.html) can be used to fit generalized additive models by a boosting algorithm. An extensible boosting framework for generalized linear, additive and nonparametric models is available in package [mboost](https://cran.r-project.org/web/packages/mboost/index.html). Likelihood-based boosting for Cox models is implemented in [CoxBoost](https://cran.r-project.org/web/packages/CoxBoost/index.html) and for mixed models in [GMMBoost](https://cran.r-project.org/web/packages/GMMBoost/index.html). GAMLSS models can be fitted using boosting by [gamboostLSS](https://cran.r-project.org/web/packages/gamboostLSS/index.html).   gbmAuthors: Originally written by <NAME>, added to by various authors, currently maintained by <NAME>Backend: C++The [gbm](https://github.com/gbm-developers/gbm) R package is an implementation of extensions to Freund and Schapire's AdaBoost algorithm and Friedman's gradient boosting machine.  This is the original R implementation of GBM. A presentation is available [here](http://www.slideshare.net/mark_landry/gbm-package-in-r) by <NAME>. Features:- Stochastic GBM.- Supports up to 1024 factor levels.- Supports Classification and regression trees.- Includes regression methods for:   - least squares  - absolute loss  - t-distribution loss   - quantile regression  - logistic   - multinomial logistic  - Poisson  - Cox proportional hazards partial likelihood  - AdaBoost exponential loss  - Huberized hinge loss  - Learning to Rank measures ([LambdaMart](https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/tr-2008-109.pdf))- Out-of-bag estimator for the optimal number of iterations is provided.- Easy to overfit since early stopping functionality is not automated in this package.- If internal cross-validation is used, this can be parallelized to all cores on the machine.- Currently undergoing a major refactoring & rewrite (and has been for some time).- GPL-2/3 License.<jupyter_code>#install.packages("gbm")
#install.packages("cvAUC")
library(gbm)
library(cvAUC)
# Load 2-class HIGGS dataset
train <- read.csv("data/higgs_train_10k.csv")
test <- read.csv("data/higgs_test_5k.csv")
set.seed(1)
model <- gbm(formula = response ~ ., 
             distribution = "bernoulli",
             data = train,
             n.trees = 70,
             interaction.depth = 5,
             shrinkage = 0.3,
             bag.fraction = 0.5,
             train.fraction = 1.0,
             n.cores = NULL)  #will use all cores by default
print(model)
# Generate predictions on test dataset
preds <- predict(model, newdata = test, n.trees = 70)
labels <- test[,"response"]

# Compute AUC on the test set
cvAUC::AUC(predictions = preds, labels = labels)<jupyter_output><empty_output><jupyter_text>xgboostAuthors: <NAME>, <NAME>, <NAME>Backend: C++The [xgboost](https://cran.r-project.org/web/packages/xgboost/index.html) R package provides an R API to "Extreme Gradient Boosting", which is an efficient implementation of gradient boosting framework. [Parameter tuning guide](http://www.analyticsvidhya.com/blog/2016/03/complete-guide-parameter-tuning-xgboost-with-codes-python/) and more resources [here](https://github.com/dmlc/xgboost/tree/master/demo).  The xgboost package is quite popular on [Kaggle](http://blog.kaggle.com/tag/xgboost/) for data mining competitions.Features:- Stochastic GBM with  column and row sampling (per split and per tree) for better generalization.- Includes efficient linear model solver and tree learning algorithms. - Parallel computation on a single machine. - Supports various objective functions, including regression, classification and ranking. - The package is made to be extensible, so that users are also allowed to define their own objectives easily.- Apache 2.0 License.<jupyter_code>#install.packages("xgboost")
#install.packages("cvAUC")
library(xgboost)
library(Matrix)
library(cvAUC)
# Load 2-class HIGGS dataset
train <- read.csv("data/higgs_train_10k.csv")
test <- read.csv("data/higgs_test_5k.csv")
# Set seed because we column-sample
set.seed(1)

y <- "response"
train.mx <- sparse.model.matrix(response ~ ., train)
test.mx <- sparse.model.matrix(response ~ ., test)
dtrain <- xgb.DMatrix(train.mx, label = train[,y])
dtest <- xgb.DMatrix(test.mx, label = test[,y])

train.gdbt <- xgb.train(params = list(objective = "binary:logistic",
                                      #num_class = 2,
                                      #eval_metric = "mlogloss",
                                      eta = 0.3,
                                      max_depth = 5,
                                      subsample = 1,
                                      colsample_bytree = 0.5), 
                                      data = dtrain, 
                                      nrounds = 70, 
                                      watchlist = list(train = dtrain, test = dtest))
# Generate predictions on test dataset
preds <- predict(train.gdbt, newdata = dtest)
labels <- test[,y]

# Compute AUC on the test set
cvAUC::AUC(predictions = preds, labels = labels)
#Advanced functionality of xgboost
#install.packages("Ckmeans.1d.dp")
library(Ckmeans.1d.dp)

# Compute feature importance matrix
names <- dimnames(data.matrix(train[,-1]))[[2]]
importance_matrix <- xgb.importance(names, model = train.gdbt)

# Plot feature importance
xgb.plot.importance(importance_matrix[1:10,])<jupyter_output>Warning message:
In set(allTrees, i = which(allTrees[, Feature] != "Leaf"), j = "Yes.Feature", : Supplied 1665 items to be assigned to 1584 items of column 'Yes.Feature' (81 unused)Warning message:
In set(allTrees, i = which(allTrees[, Feature] != "Leaf"), j = "Yes.Cover", : Supplied 1665 items to be assigned to 1584 items of column 'Yes.Cover' (81 unused)Warning message:
In set(allTrees, i = which(allTrees[, Feature] != "Leaf"), j = "Yes.Quality", : Supplied 1665 items to be assigned to 1584 items of column 'Yes.Quality' (81 unused)Warning message:
In set(allTrees, i = which(allTrees[, Feature] != "Leaf"), j = "No.Feature", : Supplied 1665 items to be assigned to 1584 items of column 'No.Feature' (81 unused)Warning message:
In set(allTrees, i = which(allTrees[, Feature] != "Leaf"), j = "No.Cover", : Supplied 1665 items to be assigned to 1584 items of column 'No.Cover' (81 unused)Warning message:
In set(allTrees, i = which(allTrees[, Feature] != "Leaf"), j = "No.Quality", : Supplied 16[...]<jupyter_text>h2oAuthors: [<NAME>](https://www.linkedin.com/in/candel), [<NAME>](http://www.cliffc.org/blog/), H2O.ai contributorsBackend: Java[H2O GBM Tuning guide by <NAME>](https://github.com/h2oai/h2o-3/blob/master/h2o-docs/src/product/tutorials/gbm/gbmTuning.Rmd) and [H2O GBM Vignette](http://docs.h2o.ai/h2o/latest-stable/h2o-docs/booklets/GBMBooklet.pdf).Features:- Distributed and parallelized computation on either a single node or a multi-node cluster.- Automatic early stopping based on convergence of user-specied metrics to user-specied relative tolerance.- Stochastic GBM with column and  row sampling  (per split and per tree) for better generalization.- Support  for  exponential  families  (Poisson,  Gamma,  Tweedie)  and  loss functions in addition to binomial (Bernoulli), Gaussian and multinomial distributions, such as Quantile regression (including Laplace).- Grid search for hyperparameter optimization and model selection.- Data-distributed, which means the entire dataset does not need to fit into memory on a single node, hence scales to any size training set.- Uses histogram approximations of continuous variables for speedup.- Uses dynamic binning - bin limits are reset at each tree level based on the split bins' min and max values discovered during the last pass.- Uses squared error to determine optimal splits.- Distributed implementation details outlined in a [blog post](http://blog.h2o.ai/2013/10/building-distributed-gbm-h2o/) by <NAME>.- Unlimited factor levels.- Multiclass trees (one for each class) built in parallel with each other.- Apache 2.0 Licensed.- Model export in plain Java code for deployment in production environments.- GUI for training & model eval/viz (H2O Flow).<jupyter_code>#install.packages("h2o")
library(h2o)
#h2o.shutdown(prompt = FALSE)  #if required
h2o.init(nthreads = -1)  #Start a local H2O cluster using nthreads = num available cores
# Load 10-class MNIST dataset
train <- h2o.importFile("data/higgs_train_10k.csv")
test <- h2o.importFile("data/higgs_test_5k.csv")
print(dim(train))
print(dim(test))
# Identity the response column
y <- "response"

# Identify the predictor columns
x <- setdiff(names(train), y)

# Convert response to factor
train[,y] <- as.factor(train[,y])
test[,y] <- as.factor(test[,y])
# Train an H2O GBM model
model <- h2o.gbm(x = x,
                 y = y,
                 training_frame = train,
                 ntrees = 70,
                 learn_rate = 0.3,
                 sample_rate = 1.0,
                 max_depth = 5,
                 col_sample_rate_per_tree = 0.5,
                 seed = 1)
# Get model performance on a test set
perf <- h2o.performance(model, test)
print(perf)
# To retreive individual metrics
h2o.auc(perf)
# Print confusion matrix
h2o.confusionMatrix(perf)
# Plot scoring history over time
plot(model)
# Retreive feature importance
vi <- h2o.varimp(model)
vi[1:10,]
# Plot feature importance
barplot(vi$scaled_importance,
        names.arg = vi$variable,
        space = 1,
        las = 2,
        main = "Variable Importance: H2O GBM")<jupyter_output><empty_output><jupyter_text>Note that all models, data and model metrics can be viewed via the [H2O Flow GUI](http://127.0.0.1:54321/flow/index.html), which should already be running since you started the H2O cluster with `h2o.init()`.<jupyter_code># Early stopping example
# Keep in mind that when you use early stopping, you should pass a validation set
# Since the validation set is used to detmine the stopping point, a separate test set should be used for model eval

#fit <- h2o.gbm(x = x,
#               y = y,
#               training_frame = train,
#               model_id = "gbm_fit3",
#               validation_frame = valid,  #only used if stopping_rounds > 0
#               ntrees = 500,
#               score_tree_interval = 5,      #used for early stopping
#               stopping_rounds = 3,          #used for early stopping
#               stopping_metric = "misclassification", #used for early stopping
#               stopping_tolerance = 0.0005,  #used for early stopping
#               seed = 1)
# GBM hyperparamters
gbm_params <- list(learn_rate = seq(0.01, 0.1, 0.01),
                   max_depth = seq(2, 10, 1),
                   sample_rate = seq(0.5, 1.0, 0.1),
                   col_sample_rate = seq(0.1, 1.0, 0.1))
search_criteria <- list(strategy = "RandomDiscrete", 
                         max_models = 20)

# Train and validate a grid of GBMs
gbm_grid <- h2o.grid("gbm", x = x, y = y,
                      grid_id = "gbm_grid",
                      training_frame = train,
                      validation_frame = test,  #test frame will only be used to calculate metrics
                      ntrees = 70,
                      seed = 1,
                      hyper_params = gbm_params,
                      search_criteria = search_criteria)

gbm_gridperf <- h2o.getGrid(grid_id = "gbm_grid", 
                            sort_by = "auc", 
                            decreasing = TRUE)
print(gbm_gridperf)<jupyter_output>|======================================================================| 100%
H2O Grid Details
================

Grid ID: gbm_grid 
Used hyper parameters: 
  -  sample_rate 
  -  max_depth 
  -  learn_rate 
  -  col_sample_rate 
Number of models: 20 
Number of failed models: 0 

Hyper-Parameter Search Summary: ordered by decreasing auc
   sample_rate max_depth learn_rate col_sample_rate         model_ids
1          0.8         6       0.06               1  gbm_grid_model_9
2          0.8         7       0.05             0.8  gbm_grid_model_7
3          0.9         7       0.08             0.6 gbm_grid_model_16
4          0.9         9       0.04             0.4 gbm_grid_model_14
5          0.7         5       0.07             0.5 gbm_grid_model_11
6          0.5         7       0.05             0.3 gbm_grid_model_17
7          0.5         6       0.03             0.8  gbm_grid_model_3
8          0.5         9       0.01             0.9  gbm_grid_model_6
9          0.6         9      [...]
{'fraction_non_alphanumeric': 0.09414266383143244, 'fraction_numerical': 0.01942779818287261, 'mean_word_length': 3.5389339767492873, 'pattern_counts': {'":': 0, '<': 74, '<?xml version=': 0, '>': 47, 'https://': 25, 'lorem ipsum': 0, 'www.': 9, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '320115', 'license': 'Apache-2.0', 'n_tokens_mistral': 6192, 'n_tokens_neox': 5601, 'n_words': 2162, 'path': 'gradient-boosting-machines.ipynb', 'repo_name': 'purushothamgowthu/useR-machine-learning-tutorial', 'size': '20886'}