In this lab, we will cover bagging, random forest, gradient boosting and extreme boosting for regression problems. We use the same Boston Housing data. In the next lab, we discuss classification problem with these state-of-the art machine linear algorithms.

# load Boston data
library(MASS)
library(tidyverse)
data(Boston)
index <- sample(nrow(Boston),nrow(Boston)*0.70)
boston.train <- Boston[index,]
boston.test <- Boston[-index,]

Bagging

Bagging stands for Bootstrap and Aggregating. It employs the idea of bootstrap but the purpose is not to study bias and standard errors of estimates. Instead, the goal of Bagging is to improve prediction accuracy. It fits a tree for each bootsrap sample, and then aggregate the predicted values from all these different trees. For more details, you may look at Wikepedia, or you can find the original paper Leo Breiman (1996).

An available R package, ipred, provides functions to perform Bagging. You need to install this package if you didn’t do it before.

library(ipred)

Fit tree with bagging on Boston training data, and calculate MSE on testing sample.

boston.bag<- bagging(medv~., data = boston.train, nbagg=100)
boston.bag
## 
## Bagging regression trees with 100 bootstrap replications 
## 
## Call: bagging.data.frame(formula = medv ~ ., data = boston.train, nbagg = 100)

Prediction on testing sample.

boston.bag.pred<- predict(boston.bag, newdata = boston.test)
mean((boston.test$medv-boston.bag.pred)^2)
## [1] 9.427629

Comparing with a single tree.

library(rpart)
boston.tree<- rpart(medv~., data = boston.train)
boston.tree.pred<- predict(boston.tree, newdata = boston.test)
mean((boston.test$medv-boston.tree.pred)^2)
## [1] 13.17775

How many trees are good?

ntree<- c(seq(10, 200, 10))
MSE.test<- rep(0, length(ntree))
for(i in 1:length(ntree)){
  boston.bag1<- bagging(medv~., data = boston.train, nbagg=ntree[i])
  boston.bag.pred1<- predict(boston.bag1, newdata = boston.test)
  MSE.test[i]<- mean((boston.test$medv-boston.bag.pred1)^2)
}
plot(ntree, MSE.test, type = 'l', col=2, lwd=2)

By fitting the Bagging multiple times and predicting the testing sample, we can draw the following boxplot to show the variance of the prediction error at different number of trees.

ntree<- c(1, 3, 5, seq(10, 200, 10))
MSE.test<- matrix(0, length(ntree), 50)
for(k in 1:50){
  for(i in 1:length(ntree)){
    boston.bag1<- bagging(medv~., data = boston.train, nbagg=ntree[i])
    boston.bag.pred1<- predict(boston.bag1, newdata = boston.test)
    MSE.test[i,k]<- mean((boston.test$medv-boston.bag.pred1)^2)
  }
}
boxplot(t(MSE.test), names=ntree, xlab="Number of Tree", ylab="Test MSE")
lines(apply(MSE.test, 1, mean), col="red", lty=2, lwd=2)

Out-of-bag (OOB) prediction

The out-of-bag prediction is similar to LOOCV. We use full sample. In every bootstrap, the unused sample serves as testing sample, and testing error is calculated. In the end, OOB error, root mean squared error by default, is obtained

boston.bag.oob<- bagging(medv~., data = boston.train, coob=T, nbagg=100)
boston.bag.oob
## 
## Bagging regression trees with 100 bootstrap replications 
## 
## Call: bagging.data.frame(formula = medv ~ ., data = boston.train, coob = T, 
##     nbagg = 100)
## 
## Out-of-bag estimate of root mean squared error:  4.3826

Random Forests

Random forest is an extension of Bagging, but it makes significant improvement in terms of prediction. The idea of random forests is to randomly select \(m\) out of \(p\) predictors as candidate variables for each split in each tree. Commonly, \(m=\sqrt{p}\). The reason of doing this is that it can decorrelates the trees such that it reduces variance when we aggregate the trees. You may refer Wikipedia and the tutorial on the author’s website.

Fitting a Random Forest model

We start with Boston Housing data.

library(randomForest)
boston.rf<- randomForest(medv~., data = boston.train, importance=TRUE)
boston.rf
## 
## Call:
##  randomForest(formula = medv ~ ., data = boston.train, importance = TRUE) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 4
## 
##           Mean of squared residuals: 11.94228
##                     % Var explained: 86.86
# How "Mean of squared residuals" (MSE) is computed
mean((boston.train$medv-boston.rf$predicted)^2)
## [1] 11.94228
# How "% Var explained" (Rsq) is computed
1-sum((boston.train$medv-boston.rf$predicted)^2)/sum((boston.train$medv-mean(boston.train$medv))^2)
## [1] 0.8686456
# $mse shows the mse of each ntree
head(boston.rf$mse)
## [1] 46.00002 45.89441 37.08099 30.73221 31.78830 27.54704
# $rsq shows the rsq of each ntree
head(boston.rf$rsq)
## [1] 0.4940408 0.4952024 0.5921421 0.6619731 0.6503570 0.6970071
# $oob.times shows how many times each observation was in oob sample
head(boston.rf$oob.times)
## [1] 168 181 205 196 154 178

By default, \(m=p/3\) for regression tree, and \(m=\sqrt{p}\) for classification problem. You can change it by specifying mtry=. You can also specify number of trees by ntree=. The default is 500.

Variable importance

The argument importance=TRUE allows us to see the variable importance.

boston.rf$importance
##            %IncMSE IncNodePurity
## crim    11.4352905     2108.0747
## zn       1.0484739      304.5956
## indus    6.9150658     1907.8394
## chas     0.6587992      200.5951
## nox      9.6245681     1820.1729
## rm      32.0872548     8464.1884
## age      4.2687199      815.8375
## dis      7.8986617     1975.7752
## rad      1.6944793      300.3245
## tax      4.4648937     1070.5355
## ptratio  7.1949002     1928.5114
## black    1.6563274      677.6526
## lstat   71.0559749     9723.5300

Here “%IncMSE” is the average percentage increase in MSE over all trees. Specifically, For each variable, the algorithm does: 1) Take the out-of-bag (OOB) samples for each tree. 2) Randomly permute the values of that variable only. 3) Recompute the prediction error (MSE) on those OOB samples. 4) Compare this MSE with the original OOB MSE.

“IncNodePurity” is the total reduction in impurity across all trees attributable to that variable. Bigger values mean the variable is frequently used in splits that substantially reduce residual variance.

This variable importance plot can be visualized.

varImpPlot(boston.rf)

OOB Error

The fitted randomForest actually saves all OOB errors for each ntree value from 1 to 500. We can make a plot to see how the OOB error changes with different ntree.

plot(boston.rf$mse, type='l', col=2, lwd=2, xlab = "ntree", ylab = "OOB Error")

Prediction

Prediction on the testing sample.

boston.rf.pred<- predict(boston.rf, boston.test)
mean((boston.test$medv-boston.rf.pred)^2)
## [1] 7.191286

Try different \(m\)

As we mentioned before, the number of candidate predictors in each split is \(m\approx \sqrt{p}=\sqrt{13}\approx 4\). We can also specify \(m\) with argument mtry. Now let’s see how the OOB error changes with mtry.

oob.err<- rep(0, 13)
for(i in 1:13){
  fit<- randomForest(medv~., data = boston.train, mtry=i)
  oob.err[i]<- fit$mse[500]
  cat(i, " ")
}
## 1  2  3  4  5  6  7  8  9  10  11  12  13
plot(oob.err, pch=15, col = "blue", type = "b", ylab = "OOB MSE", xlab = "mtry")
abline(v=which.min(oob.err), lty=2, col="red")

Exercise: Create a plot displaying the OOB error across ntree=1, …, 500, and mtry= 1, …, 13. (You can draw 13 lines in different color representing each \(m\)).

Boosting

Boosting builds a number of small trees, and each time, the response is the residual from last tree. It is a sequential procedure. We use gbm package to build boosted trees.

library(gbm)
?gbm
boston.boost<- gbm(medv~., data = boston.train, 
                   distribution = "gaussian", 
                   n.trees = 10000, 
                   shrinkage = 0.01, 
                   interaction.depth = 3,
                   cv.folds = 5, n.cores = 5)

Note that we need to specify distribution = "gaussian" if we are working on regression tree. The default is Bernoulli distribution for binary classification problem. n.trees is the number of small trees we fit. We need to choose this parameter carefully because it may results in overfitting if the number is too large. shrinkage is another tuning parameter that controls how much contribution each tree makes. interaction.depth is how many splits of each tree we want. All those tuning parameters can be chosen from cross-validation. The idea is that we don’t want overfitting.

Find the optimal number of trees

When we fit the model, we specified cv, which conducts an internal cross-validation. This is used to find optimal number of trees. The function gbm.perf shows it.

best_iter <- gbm.perf(boston.boost, method = "cv")

best_iter
## [1] 2978

There are two curves in the figure. One is the training error and the other is testing error. You should know which one is which.

Variable importance

summary(boston.boost)

##             var    rel.inf
## lstat     lstat 39.9746010
## rm           rm 25.8778093
## dis         dis  9.3671112
## crim       crim  5.5863927
## nox         nox  4.5535221
## age         age  3.7936479
## black     black  2.9831058
## ptratio ptratio  2.9032183
## tax         tax  2.0003707
## indus     indus  1.1962384
## rad         rad  0.8181993
## chas       chas  0.8179984
## zn           zn  0.1277847

Partial dependence / effect plots

The fitted boosted tree also gives the relation between response and each predictor.

par(mfrow=c(1,2))
plot(boston.boost, i="lstat", n.trees = best_iter)

plot(boston.boost, i="rm", n.trees = best_iter)

Prediction

Prediction on testing sample.

boston.boost.pred.test<- predict(boston.boost, boston.test, n.trees = best_iter)
mean((boston.test$medv-boston.boost.pred.test)^2)
## [1] 8.337923

XGBoost

XGBoost (eXtreme Gradient Boosting) is an efficient implementation of gradient boosted trees.
Conceptually it is the same idea as gradient boosting:

  1. Fit many small regression trees in sequence.
  2. Each new tree tries to correct the residuals of the previous trees.
  3. Average all trees to form the final prediction.

What makes XGBoost attractive in practice is that it adds:

  • Regularization on tree complexity (to reduce overfitting),
  • Efficient handling of sparse / missing values,
  • Parallelization and many engineering tricks that make it much faster and more scalable than a “vanilla” gradient boosting implementation.
library(xgboost)
# Prepare matrices for XGBoost (it works on numeric matrices, not data frames)
x_train <- as.matrix(boston.train[, -14])
y_train <- boston.train[, 14]
x_test <- as.matrix(boston.test[, -14])
y_test <- boston.test[, 14]

# Put data in DMatrix objects (XGBoost's internal format)
dtrain <- xgb.DMatrix(data = x_train, label = y_train)
dtest <- xgb.DMatrix(data = x_test, label = y_test)

Fitting a basic XGBoost model

We start with a simple model and a fixed set of hyperparameters.
Here we control:

  • max_depth: complexity of each tree
  • eta: learning rate (how much each tree changes the prediction)
  • nrounds: number of boosting iterations (number of trees)
params <- list(
  objective = "reg:squarederror", # squared-error loss for regression
  eta       = 0.1,                # learning rate
  max_depth = 4,                  # tree depth
  subsample = 0.8,                # row subsampling
  colsample_bytree = 0.8          # column subsampling
)

fit.xgboost.reg <- xgboost(
  data = dtrain,
  params = params,
  nrounds = 100,
  verbose = 0
)

Cross-validation and choosing the number of trees

We do not want to guess how many trees (nrounds) to use.
XGBoost has built-in cross-validation. We also use early stopping: stop when the validation error stops improving.

set.seed(2025)

cv.nrounds <- 500   # upper bound; early stopping will usually stop earlier

xgboost.cv <- xgb.cv(
  params  = params,
  data    = dtrain,
  nrounds = cv.nrounds,
  nfold   = 10,
  metrics = "rmse",
  early_stopping_rounds = 20,   # stop if no improvement for 20 rounds
  verbose = 0
)

# Best number of trees according to CV
best_nrounds <- xgboost.cv$best_iteration
best_nrounds
## [1] 141
# Plot test RMSE over boosting iterations
plot(xgboost.cv$evaluation_log$iter,
     xgboost.cv$evaluation_log$test_rmse_mean,
     type = "l",
     xlab = "Boosting iteration",
     ylab = "CV test RMSE")
abline(v = best_nrounds, lty = 2)

The performance stays the same after 100 iterations. We choose nrounds=100 and use xgboost() to fit the model. You may also grid search learning rate eta using the similar way.

Prediction

pred.xgboost<- predict(fit.xgboost.reg, newdata = as.matrix(boston.test[,-14]))
mean((boston.test$medv-pred.xgboost)^2)
## [1] 8.484655

Interpreting the fitted XGBoost model

XGBoost is still a tree-based model, so we can look at variable importance and partial dependence to understand what it is doing.

# Variable importance
xgb_imp <- xgb.importance(model = fit.xgboost.reg)
xgb_imp
##     Feature         Gain       Cover  Frequency
##      <char>        <num>       <num>      <num>
##  1:   lstat 0.4990388214 0.148194383 0.10902896
##  2:      rm 0.2693138987 0.162294449 0.14565588
##  3:     dis 0.0620014488 0.153996066 0.13287905
##  4:    crim 0.0538924486 0.082175541 0.16439523
##  5:     nox 0.0365863983 0.071560078 0.08347530
##  6: ptratio 0.0271247454 0.051074569 0.04088586
##  7:     tax 0.0146260763 0.060495568 0.04684838
##  8:   black 0.0114994338 0.084672241 0.09284497
##  9:     age 0.0105208206 0.094093240 0.09625213
## 10:   indus 0.0079069660 0.051047626 0.04429302
## 11:     rad 0.0047835357 0.030858487 0.02385009
## 12:    chas 0.0020417851 0.006987166 0.00681431
## 13:      zn 0.0006636212 0.002550585 0.01277683
# Plot top variables
xgb.plot.importance(xgb_imp, top_n = 10)