Advanced Tree Models (for Classification) – Bagging, Random Forests, and Boosting

In this lab, we show how to do classification using the tree-based ensemble methods.

library(tidyverse)
# load credit card data
credit.data <- read_csv("https://www.dropbox.com/s/tnoo06n8m842uit/credit_card_default.csv?dl=1")
credit.data<- rename(credit.data, default=`default payment next month`)
# convert categorical data to factor
credit.data$SEX<- as.factor(credit.data$SEX)
credit.data$EDUCATION<- as.factor(credit.data$EDUCATION)
credit.data$MARRIAGE<- as.factor(credit.data$MARRIAGE)
# random splitting
index <- sample(nrow(credit.data),nrow(credit.data)*0.80)
credit.train = credit.data[index,]
credit.test = credit.data[-index,]

# load MNIST data
digit<- data.matrix(read_csv("https://www.dropbox.com/s/ulujvi2a4ykfzju/train.csv?dl=1"))
# random splitting
index<- sample(1:nrow(digit), 0.6*nrow(digit))
train<- digit[index,]
test<- digit[-index,]
# standardize X
train.x <- train[,-1] #remove 'label' column
test.x<- test[,-1]
train.y <- train[,1] #label column
test.y<- test[,1]
train.x <- train.x/255
test.x <- test.x/255

Bagging

Fitting model

library(ipred)
credit.bag<- bagging(as.factor(default)~., data = credit.train, nbagg=100, coob=T, method="class")
credit.bag

## 
## Bagging classification trees with 100 bootstrap replications 
## 
## Call: bagging.data.frame(formula = as.factor(default) ~ ., data = credit.train, 
##     nbagg = 100, coob = T, method = "class")
## 
## Out-of-bag estimate of misclassification error:  0.1862

We can extract any single tree from bagging, and it is an rpart object.

length(credit.bag$mtrees)

## [1] 100

one_tree <- credit.bag$mtrees[[1]][[1]]  # rpart object

Prediction accuracy

predprob.test.bag<- predict(credit.bag, newdata = credit.test, type="prob")[,2]
library(ROCR)
pred <- prediction(predprob.test.bag, credit.test$default)
perf <- performance(pred, "tpr", "fpr")
plot(perf, colorize=TRUE)

#Get the AUC
unlist(slot(performance(pred, "auc"), "y.values"))

## [1] 0.7681437

The classification results are generated by specifying type="class".

predclass.test.bag<- predict(credit.bag, newdata = credit.test, type="class")
table(credit.test$default, predclass.test.bag, dnn = c("True", "Pred"))

##     Pred
## True    0    1
##    0 4418  265
##    1  827  490

## cost when FN:FP=5:1
cost51 <- function(true, pred){
  sum(true==1 & pred==0)*5+sum(true==0 & pred==1)*1
}
cost51(true=credit.test$default, pred=predclass.test.bag)

## [1] 4400

Or we can set our own threshold to make the classification. In this case, we need probability output from the model. This can address the asymetric cost issue, where the model by default uses majority vote which is equivalent to 0.5 cutoff probability.

predprob.test.bag<- predict(credit.bag, newdata = credit.test, type="prob")[,2]
predclass.test.bag2 <- (predprob.test.bag>mean(credit.train$default))*1
table(credit.test$default, predclass.test.bag2, dnn = c("True", "Pred"))

##     Pred
## True    0    1
##    0 3360 1323
##    1  428  889

cost51(true=credit.test$default, pred=predclass.test.bag2)

## [1] 3463

Comparing with a single tree. For fair comparison, we assume the loss is symmetric.

library(rpart)
credit.rpart <- rpart(formula = default ~ ., data = credit.train, 
                      method = "class", parms=list(loss=matrix(c(0,5,1,0), nrow=2)))
credit.test.pred.tree1<- predict(credit.rpart, credit.test, type="class")
table(credit.test$default, credit.test.pred.tree1, dnn=c("Truth","Predicted"))

##      Predicted
## Truth    0    1
##     0 2520 2163
##     1  283 1034

cost51(true=credit.test$default, pred=credit.test.pred.tree1)

## [1] 3578

Which one is better?

Note that parms=list(loss=matrix(c(0,10,1,0), nrow=2)) in bagging() does not work. As far as I know, bagging() cannot handle asymetric loss. There might be other packages can do so. However, we can write our own algorithm for bagging based on rpart, and it should be pretty straightforward.

Write our own bagging with Paralell computing

library(doParallel)  # for parallel backend to foreach
library(foreach)     # for parallel processing with for loops
detectCores()        # detect how many cores are on you PC

## [1] 24

cl <- makeCluster(8) # use 8 workers
registerDoParallel(cl) # register the parallel backend
ptc <- proc.time() # time stamp
mybag.predictions <- foreach(1:100, .packages = "rpart", .combine = cbind) %dopar% {
    # bootstrap copy of training data
    index <- sample(nrow(credit.train), replace = TRUE)
    credit_train_boot <- credit.train[index, ]  
    
    # fit tree to bootstrap copy
    bagged_tree <- rpart(formula = default ~ ., data = credit_train_boot, method = "class",
                         parms = list(loss=matrix(c(0,5,1,0), nrow = 2))) 
    
    predict(bagged_tree, newdata = credit.test, type="class")
}
stopCluster(cl)
proc.time()-ptc  # time stamp

##    user  system elapsed 
##    0.08    0.03   10.78

Let’s look at the results from paralell computing, and make it as what we want.

dim(mybag.predictions)

## [1] 6000  100

mybag.predictions[1:10,1:10]

##    result.1 result.2 result.3 result.4 result.5 result.6 result.7 result.8
## 1         1        1        2        1        1        2        1        1
## 2         2        2        2        2        2        2        2        2
## 3         1        1        1        1        1        1        1        1
## 4         1        1        1        1        1        1        1        1
## 5         2        2        1        2        2        2        1        1
## 6         2        2        2        2        2        2        2        2
## 7         1        1        1        1        1        1        1        1
## 8         2        2        2        2        2        2        2        2
## 9         2        2        2        2        2        2        2        2
## 10        2        2        2        2        2        2        2        1
##    result.9 result.10
## 1         1         1
## 2         2         2
## 3         1         1
## 4         1         1
## 5         2         1
## 6         2         2
## 7         1         2
## 8         2         2
## 9         2         2
## 10        2         2

mybag.predictions <- mybag.predictions-1
mybag.pred <- apply(mybag.predictions, 1, mean)
summary(mybag.pred)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.4000  0.4573  0.9900  1.0000

mybag.pred.class <- (mybag.pred>0.5)*1
table(credit.test$default, mybag.pred.class, dnn=c("Truth","Predicted"))

##      Predicted
## Truth    0    1
##     0 2849 1834
##     1  316 1001

cost51(true=credit.test$default, pred=mybag.pred.class)

## [1] 3414

Random Forests

Model fitting

library(randomForest)

## Warning: package 'randomForest' was built under R version 4.5.2

## randomForest 4.7-1.2

## Type rfNews() to see new features/changes/bug fixes.

## 
## Attaching package: 'randomForest'

## The following object is masked from 'package:dplyr':
## 
##     combine

## The following object is masked from 'package:ggplot2':
## 
##     margin

credit.rf<- randomForest(as.factor(default)~., data = credit.train, cutoff=c(5/6,1/6)) # classwt does not change the results
credit.rf

## 
## Call:
##  randomForest(formula = as.factor(default) ~ ., data = credit.train,      cutoff = c(5/6, 1/6)) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 4
## 
##         OOB estimate of  error rate: 37.48%
## Confusion matrix:
##       0    1 class.error
## 0 10993 7688   0.4115411
## 1  1307 4012   0.2457229

The argument cutoff specifies the cut off probability of each class (e.g. Y=0/1) for prediction. In our example, we specify (5/6, 1/6), meaning that the rule of majority vote for 0 needs to be at least 5/6 and for 1 needs to be 1/6, instead of half-half. This is equivalent to asymmetric cost of FN:FP=5:1. Therefore, the MR or OOB error rate may not be a good measure of prediction accuracy, and you have to calculate your own asymmetric cost

We can again easily plot the error rate vs. ntree.

plot(credit.rf, lwd=rep(2, 3))
legend("right", legend = c("OOB Error", "FPR", "FNR"), lwd=rep(2, 3), lty = c(1,2,3), col = c("black", "red", "green"))

You may try different specification of cutoff and see the change of confusion matrix. In addition, for prediction purpose, you can also apply your own cutoff on the probabilistic prediction as we learned in logistic regression.

credit.rf1<- randomForest(as.factor(default)~., data = credit.train) 
credit.rf1

## 
## Call:
##  randomForest(formula = as.factor(default) ~ ., data = credit.train) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 4
## 
##         OOB estimate of  error rate: 18.15%
## Confusion matrix:
##       0    1 class.error
## 0 17663 1018  0.05449387
## 1  3339 1980  0.62774958

plot(credit.rf1, lwd=rep(2, 3))
legend("right", legend = c("OOB Error", "FPR", "FNR"), lwd=rep(2, 3), lty = c(1,2,3), col = c("black", "red", "green"))

Variable importance

credit.rf$importance

##           MeanDecreaseGini
## LIMIT_BAL        419.97285
## SEX               76.88544
## EDUCATION        158.90498
## MARRIAGE          91.21211
## AGE              466.77659
## PAY_0            785.72758
## PAY_2            361.64651
## PAY_3            232.04314
## PAY_4            204.20056
## PAY_5            190.04321
## PAY_6            154.26134
## BILL_AMT1        481.95808
## BILL_AMT2        449.42511
## BILL_AMT3        427.43616
## BILL_AMT4        417.52508
## BILL_AMT5        413.95284
## BILL_AMT6        417.45805
## PAY_AMT1         417.65035
## PAY_AMT2         400.57410
## PAY_AMT3         377.77006
## PAY_AMT4         361.03641
## PAY_AMT5         357.32099
## PAY_AMT6         374.82836

We can visualize it with vip package.

library(vip)
vip(credit.rf, num_features = 15, geom = "point")

Prediction

credit.rf.pred<- predict(credit.rf, newdata=credit.test, type = "prob")[,2]
library(ROCR)
pred <- prediction(credit.rf.pred, credit.test$default)
perf <- performance(pred, "tpr", "fpr")
plot(perf, colorize=TRUE)

#Get the AUC
unlist(slot(performance(pred, "auc"), "y.values"))

## [1] 0.773028

Below is the confusion matrix based on class prediction

credit.rf.class.test<- predict(credit.rf, newdata=credit.test, type = "class")
table(credit.test$default, credit.rf.class.test, dnn = c("True", "Pred"))

##     Pred
## True    0    1
##    0 2779 1904
##    1  313 1004

cost51(true=credit.test$default, pred=credit.rf.class.test)

## [1] 3469

Below is the confusion matrix based on probabilistic prediction with user specified cutoff (overall default rate).

credit.rf.class.test<- (credit.rf.pred>mean(credit.train$default))*1
table(credit.test$default, credit.rf.class.test, dnn = c("True", "Pred"))

##     Pred
## True    0    1
##    0 3493 1190
##    1  454  863

cost51(true=credit.test$default, pred=credit.rf.class.test)

## [1] 3460

Boosting

Model fitting

library(gbm)

## Warning: package 'gbm' was built under R version 4.5.2

## Loaded gbm 2.2.2

## This version of gbm is no longer under development. Consider transitioning to gbm3, https://github.com/gbm-developers/gbm3

credit.boost<- gbm(default~., data = credit.train, distribution = "bernoulli", 
                   n.trees = 2000, cv.folds = 5, n.cores = 5)

The argument distribution is to specify the type of likelihood function, n.tree is to specify the maximum number of trees, cv.folds is optional. We use cross-validation to choose the best number of trees. Recall that Boosting can easily overfit the data. The last argument n.cores is to specify the number of cores to be used for parallel computing.

# relative influence
summary(credit.boost)

##                 var    rel.inf
## PAY_0         PAY_0 41.1042725
## PAY_2         PAY_2  7.6991207
## BILL_AMT1 BILL_AMT1  4.2492063
## BILL_AMT2 BILL_AMT2  3.8090360
## PAY_5         PAY_5  3.6937636
## LIMIT_BAL LIMIT_BAL  3.5849787
## PAY_3         PAY_3  3.2796699
## PAY_AMT2   PAY_AMT2  3.0550542
## AGE             AGE  3.0066809
## BILL_AMT3 BILL_AMT3  2.9676779
## PAY_4         PAY_4  2.9037943
## PAY_AMT3   PAY_AMT3  2.5976721
## BILL_AMT6 BILL_AMT6  2.4235600
## BILL_AMT5 BILL_AMT5  2.3060828
## BILL_AMT4 BILL_AMT4  2.2394675
## PAY_AMT1   PAY_AMT1  2.2084209
## PAY_AMT4   PAY_AMT4  2.0484894
## PAY_6         PAY_6  1.8122691
## PAY_AMT6   PAY_AMT6  1.6695578
## PAY_AMT5   PAY_AMT5  1.6219944
## MARRIAGE   MARRIAGE  0.7003176
## EDUCATION EDUCATION  0.6855226
## SEX             SEX  0.3333908

Below is the figure of cross-validation to choose the optimal number of trees. The black curve represents the training error and green the CV error.

best.iter <- (gbm.perf(credit.boost, method = "cv"))

best.iter

## [1] 461

Prediction

# predicted probability
pred.credit.boost<- predict(credit.boost, newdata = credit.test, n.trees = best.iter, type="response")
# AUC
pred <- prediction(pred.credit.boost, credit.test$default)
# perf <- performance(pred, "tpr", "fpr")
# plot(perf, colorize=TRUE)
unlist(slot(performance(pred, "auc"), "y.values"))

## [1] 0.7776554

The predicted value is the probability if we specify type="response", as in logistic regression.

Below we show the confusion matrix using the cutoff as the sample proportion.

pred.credit.boost.class<- (pred.credit.boost>mean(credit.train$default))*1
table(credit.test$default, pred.credit.boost.class, dnn = c("True", "Pred"))

##     Pred
## True    0    1
##    0 3807  876
##    1  518  799

cost51(true=credit.test$default, pred=pred.credit.boost.class)

## [1] 3466

Compare with logit model

#Fit logistic regression model
credit.glm<- glm(default~., data = credit.train, family=binomial)
pred.glm<- predict(credit.glm, credit.test, type="response")
# AUC
pred <- prediction(pred.glm, credit.test$default)
unlist(slot(performance(pred, "auc"), "y.values"))

## [1] 0.7131271

#Get binary prediction
credit.test.pred.glm<- (pred.glm>mean(credit.train$default))*1
#Confusion matrix
table(credit.test$default, credit.test.pred.glm, dnn=c("True","Pred"))

##     Pred
## True    0    1
##    0 3285 1398
##    1  488  829

cost51(true=credit.test$default, pred=credit.test.pred.glm)

## [1] 3838

XGboost

library(xgboost)
# Put data in DMatrix objects (XGBoost's internal format)
dtrain <- xgb.DMatrix(data = model.matrix(~., credit.train[,-24])[,-1], label = credit.train$default)
dtest <- xgb.DMatrix(data = model.matrix(~., credit.test[,-24])[,-1], label = credit.test$default)

params <- list(
  objective = "binary:logistic",  # logit loss for binary classification
  eta       = 0.05,                # learning rate
  max_depth = 4                  # tree depth
)

ptm <- proc.time()
fit.xgboost.class <- xgboost(
  data = dtrain,
  params = params,
  nrounds = 100,
  nthread = 6,
  verbose = 0
)
proc.time()-ptm

##    user  system elapsed 
##    1.97    0.39    0.55

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 <- 700   # upper bound; early stopping will usually stop earlier

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

# Best number of trees according to CV
best_nrounds <- xgboost.cv$best_iteration
best_nrounds

## [1] 122

# Plot test AUC over boosting iterations
plot(xgboost.cv$evaluation_log$iter,
     xgboost.cv$evaluation_log$test_auc_mean,
     type = "l",
     xlab = "Boosting iteration",
     ylab = "CV test AUC")
abline(v = best_nrounds, lty = 2)

Prediction

pred.credit.xgboost<- predict(fit.xgboost.class, newdata = model.matrix(~., credit.test[,-24])[,-1])
# AUC
pred <- prediction(pred.credit.xgboost, credit.test$default)
# perf <- performance(pred, "tpr", "fpr")
# plot(perf, colorize=TRUE)
unlist(slot(performance(pred, "auc"), "y.values"))

## [1] 0.7818113

Confusion matrix

Below we show the confusion matrix using the cutoff as the sample proportion.

pred.credit.xgboost.class<- (pred.credit.xgboost>mean(credit.train$default))*1
table(credit.test$default, pred.credit.xgboost.class, dnn = c("True", "Pred"))

##     Pred
## True    0    1
##    0 3743  940
##    1  493  824

cost51(true=credit.test$default, pred=pred.credit.xgboost.class)

## [1] 3405

XGboost for large dataset – MNIST data

We can’t see the advantage of XGboost for above example (small data). It even takes longer than gbm(). However, the fast computing speed can be seen for large scale data analysis.

ptm <- proc.time()
xgboost.digit.fit<- xgboost(data = train.x, label = train.y, eta = 0.2,  num_class=10,
                    nthread = 8, nrounds = 100, objective = "multi:softmax", verbose = 0)
proc.time()-ptm

##    user  system elapsed 
##  409.33    2.41   56.44

xgboost.digit.pred = predict(xgboost.digit.fit, test.x)
mean(xgboost.digit.pred==test.y)

## [1] 0.9674405

Over 95% accuracy rate!!!

Compare with `gbm()`.

It will take a very long time, could be hours depending on the computing power.

ptm <- proc.time()
gbm.digit.fit<- gbm(train.y~., data=data.frame(train.y, train.x), distribution = "multinomial", 
                    shrinkage = 0.1, n.trees = 500, cv.folds=5, n.cores = 8)
proc.time()-ptm

# performance
best.iter <- (gbm.perf(gbm.digit.fit, method = "cv"))

gbm.digit.pred <- predict(gbm.digit.fit, data.frame(test.y, test.x), n.trees = best.iter, type = "response")
gbm.digit.pred <- apply(gbm.digit.pred, 1, which.max)-1
mean(gbm.digit.pred==test.y)

Advanced Tree Models (for Classification) – Bagging, Random Forests, and Boosting

Shaobo Li

Bagging

Fitting model

Prediction accuracy

Write our own bagging with Paralell computing

Random Forests

Model fitting

Variable importance

Prediction

Boosting

Model fitting

Prediction

Compare with logit model

XGboost

Cross-validation and choosing the number of trees

Prediction

Confusion matrix

XGboost for large dataset – MNIST data

Compare with gbm().

Compare with `gbm()`.