library(tidyverse)
library(knitr)
library(ggthemes)
library(ISLR)
library(ggforce) 
theme_set(theme_bw())

Updated Schedule

Ensemble Models: 8.2

Combine simple models in order to obtain a more powerful model

Today: Bagging

Taken from Professor Wells

Additional instructions: don’t scroll past the white space until we get there

Ensemble Models

Who Wants to Be a Millionaire?

  • Who Wants to Be a Millionaire is a television gameshow that debuted in the 1990s and in which contestants answer a series of increasingly difficult multiple choice questions in order to win the grand prize of $1,000,000.
  • The original show included 3 “lifeline” options contestants could use to answer questions:
    • 50:50: Two randomly selected incorrect answers are eliminated
    • Phone a Friend: The contestant calls a friend and is given 30 seconds to discuss
    • Ask the Audience: Audience members each vote on the answer they think is correct
  • Question 1 Which lifeline has the highest chance of producing the correct answer? Why?





































Ensembling

At the Grinnell Walmart, consider the grocery section (south side of the store, to the left when you walk in). Write your answer to the following question on a sheet of paper. Do not consult your neighbor: Note this may have changed with the remodel, but I’m using the old numbers

  • Question 2 part A How many grocery aisles are there?





































  • In fact, there are 28 grocery aisles.
  • A survey of 30 current STEM students in STA 209 were asked the same question.
    • The average value of their guess was 28.4, with standard deviation of 16.2
    • The rMSE for individual students was 15.9.
    • However, the rMSE for the class was only 0.4.
  • Question 2 part B What factors about the class would cause such a drastic reduction in rMSE?

Ensemble Methods

  • Suppose we have \(m\) different models to predict \(Y\) based on \(X_1, \dots, X_n\). Suppose \(\hat{Y}_i\) is the prediction made by the \(i\)th model.

  • A simple ensemble model makes a prediction \(\hat{Y}\) as the weighted average of the predictions from each model: \[ \hat{Y} = w_1 \hat{Y}_1 + \dots + w_m \hat{Y}_m \qquad \textrm{where }w_1 + \dots w_m = 1, \quad w_i \geq 0 \]

  • Advantages of ensemble models?

    • Significantly more flexible than a single model
    • More efficient than single model
    • More resilient against model-building bias

  • Disadvantages?

    • Making predictions is more computationally expensive
    • Favors models with low test time
    • Diminishing returns on the number models that can be incorporated in ensemble

Bagging

Bootstrapping Review

  • Many occasions where it would be useful to have multiple samples from a population:
    • Estimating standard errors of a statistic
    • Creating confidence intervals to estimate a parameter
    • Estimating variability in test RMSE
  • However, subsetting existing sample to create many smaller samples is problematic:
    • Diminishes ability to build accurate models, since less data is used
    • Variability in statistics depends on sample size
  • Instead, we create bootstrap samples by sampling with replacement from the original sample a number of times equal to the original sample size
    • Since we sample with replacement, some observations from original sample are included more than once while others are omitted, which introduces a source of variability
    • Bootstrap sample is same size as original, and so can be used to estimate variability in statistic

Bootstrap Sample Computations

Question 3 Suppose we have a sample with 4 unique observations \(\{x_1,x_2, x_3, x_4\}\) and create a bootstrap sample.

  • What is the probability that \(x_1\) is not the 1st element in the bootstrap sample? What is the probability that it is not the 2nd?
  • What is the probability that \(x_1\) is not in the bootstrap sample at all?
  • What is the probability that \(x_2\) is not in the bootstrap sample?

Question 4 Now suppose we have \(n\) unique observations \(\{x_1, \dots, x_n\}\) and we create a bootstrap sample.

  • What is the probability that \(x_1\) is not in the bootstrap sample?
  • What is the probability that an arbitrary observation \(x_i\) is not in the bootstrap sample?
  • What happens to the probability that an arbitrary observation is not in the sample, as \(n\) goes to infinity?

Bagging

Suppose we only have one training set, but still want to build an ensemble of regression tree models. How can we do it?

  • Bagging (Bootstrap aggregation) was one of the earliest ensemble techniques
  • To create a bagged model, create many bootstrap samples from the original training set, and fit a decision tree to each. Average the resulting predictions.
    • Why do it?
  • Recall that decision trees tend to have high variance. But averaging the results of independent (or weakly dependent) variables decreases variance
    • Think about the Central Limit Theorem:
    • For large \(n\), \(\bar X\) is approximately Normal, with mean \(\mu\) and standard deviation \(\frac{\sigma}{\sqrt{n}}\)
  • Unlike a single tree model, we do not prune (we instead control variance by averaging)

Test Error for Bagged Models

  • Previously, we showed that an individual observation has probability \(e^{-1} \approx 0.368\) of not appearing in a bootstrap sample.
  • For each bootstrap, approximately 1/3 of observations are not included (called out-of-bag observations)
  • The out-of-bag observations can be used as a natural validation set for the bootstrap model.
  • We get an overall estimate of test MSE for the bagged model by averaging the MSE of each bootstrap model on its out-of-bag observations

A Bag of Trees

We return to the pdxTrees data set, this time expanding both our data set size and number of predictors:

library(pdxTrees)
my_pdxTrees <- get_pdxTrees_parks(park = c("Berkeley Park", "Woodstock Park", "Westmoreland Park", "Mt Scott Park", "Powell Park", "Kennilworth Park", "Sellwood Park", "Crystal Springs Rhododendron Garden", "Laurelhurst Park"))%>%  mutate_if(is.character, as.factor) %>%  dplyr::select(DBH, Condition, Tree_Height, Crown_Width_NS, Crown_Width_EW, Crown_Base_Height, Functional_Type, Mature_Size, Carbon_Sequestration_lb) %>% drop_na()
names(my_pdxTrees)
## [1] "DBH"                     "Condition"              
## [3] "Tree_Height"             "Crown_Width_NS"         
## [5] "Crown_Width_EW"          "Crown_Base_Height"      
## [7] "Functional_Type"         "Mature_Size"            
## [9] "Carbon_Sequestration_lb"
dim(my_pdxTrees)
## [1] 3015    9
set.seed(1)
library(rsample)
my_pdxTrees_split <- initial_split(my_pdxTrees )
my_pdxTrees_train <- training(my_pdxTrees_split)
my_pdxTrees_test <- testing(my_pdxTrees_split)
  • Can we improve on our previous model predicting Carbon_Sequestration_lb, now using more data and more predictors?

Bagged pdXTrees

  • Let’s get a few bootstrap samples using rsample:
library(rsample)
set.seed(1115)
pdx_bootstrap <- bootstraps(my_pdxTrees_train, times = 4)
  • And now build trees on each:
library(rpart)
get_tree <- function(split){
  bootstrap_sample <- analysis(split)
  model <- rpart(Carbon_Sequestration_lb ~., data = bootstrap_sample)
}
pdx_bootstrap$model <- map(pdx_bootstrap$splits, get_tree)

A few trees

par(mfcol = c(2, 2), mar = c(1, 1, 1, 1))
library(rpart.plot)
rpart.plot(pdx_bootstrap$model[[1]], box.palette = "Greens")
rpart.plot(pdx_bootstrap$model[[2]], box.palette = "Greens")
rpart.plot(pdx_bootstrap$model[[3]], box.palette = "Greens")
rpart.plot(pdx_bootstrap$model[[4]], box.palette = "Greens")

Performance

  • Let’s get predictions for each bootstrap:
get_predictions <- function(model){
  predictions <- predict(model, my_pdxTrees_test)
  data.frame(obs = my_pdxTrees_test$Carbon_Sequestration_lb, preds = predictions)
}
pdx_bootstrap$predictions <- map(pdx_bootstrap$model, get_predictions)
  • And calculate rmse on each using yardstick
library(yardstick)
results <- map_dfr(pdx_bootstrap$predictions, rmse, obs, preds)
results
## # A tibble: 4 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard        14.0
## 2 rmse    standard        14.3
## 3 rmse    standard        14.3
## 4 rmse    standard        13.1
mean(results$.estimate)
## [1] 13.89024

Variation in Model Predictions

  • How do individual tree predictions compare?
some_preds <- tibble(
  tree1 =  pdx_bootstrap[[4]][[1]]$preds,
  tree2 = pdx_bootstrap[[4]][[2]]$preds,
  tree3 = pdx_bootstrap[[4]][[3]]$preds,
  tree4 = pdx_bootstrap[[4]][[4]]$preds
) %>% rowwise() %>% mutate(bagged = mean(c_across( )))
some_preds%>% head() 
## # A tibble: 6 × 5
## # Rowwise: 
##   tree1 tree2 tree3 tree4 bagged
##   <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1  49.0  30.4  52.2  32.5   41.0
## 2  34.7  38.0  43.3  32.8   37.2
## 3  56.8  84.0  67.6  72.9   70.3
## 4  30.6  46.6  38.8  37.8   38.4
## 5  56.8  65.7  67.6  72.9   65.7
## 6  56.8  84.0  89.1  72.9   75.7
  • How does the bagged model RMSE compare to each individual tree’s RMSE?
bagged_results <- data.frame(obs = my_pdxTrees_test$Carbon_Sequestration_lb, preds = some_preds$bagged)

rbind(results, rmse(bagged_results, truth = obs, estimate = preds )) %>% mutate(model = c("tree 1", "tree 2", "tree 3", "tree 4", "bagged")) %>% select(model, everything())
## # A tibble: 5 × 4
##   model  .metric .estimator .estimate
##   <chr>  <chr>   <chr>          <dbl>
## 1 tree 1 rmse    standard        14.0
## 2 tree 2 rmse    standard        14.3
## 3 tree 3 rmse    standard        14.3
## 4 tree 4 rmse    standard        13.1
## 5 bagged rmse    standard        12.3
  • Note that the RMSE for the bagged tree is NOT simply the average RMSE. It is significantly lower!

The more trees the merrier?

Question 4 If 4 trees improved performance over 1, what if we bagged 10 trees? 100?





































Solutions

WWTBAM

Audience

Bootstrap Sample Computations Solutions

Suppose we have a sample with 4 unique observations \(\{x_1,x_2, x_3, x_4\}\) and create a bootstrap sample.

  • What is the probability that \(x_1\) is not the 1st element in the bootstrap sample? What is the probability that it is not the 2nd?

  • What is the probability that \(x_1\) is not in the bootstrap sample at all?

  • What is the probability that \(x_2\) is not in the bootstrap sample?

Now suppose we have \(n\) unique observations \(\{x_1, \dots, x_n\}\) and we create a bootstrap sample.

  • What is the probability that \(x_1\) is not in the bootstrap sample?

  • What is the probability that an arbitrary observation \(x_i\) is not in the bootstrap sample?

What happens to the probability that an arbitrary observation is not in the sample, as \(n\) goes to infinity?

The more trees the merrier?

Question 4 If 4 trees improved performance over 1, what if we bagged 10 trees? 100?

set.seed(11)
library(randomForest)
rfmodels<-list()
tree_size <- c(1:5,10*1:5, 75,100, 150, 200)
for (i in tree_size){
  rfmodels[[i]]<-randomForest(Carbon_Sequestration_lb ~ ., data = my_pdxTrees_train, ntree = i , mtry = 8, na.action = na.roughfix )
}
results <- data.frame()
for (i in tree_size){
  preds <-predict(rfmodels[[i]], newdata = my_pdxTrees_test)
  obs <- my_pdxTrees_test$Carbon_Sequestration_lb
  results <- rbind(results, data.frame(n_trees = i, preds, obs))
}

bagged_RMSE <- results  %>% group_by(n_trees) %>% rmse(truth = obs, estimate = preds) 

bagged_RMSE %>% 
ggplot( aes( x = n_trees, y = .estimate))+geom_point()+theme_bw()+geom_line()+labs(x = "Number of Trees", y = "RMSE", title = "RMSE of Bagged Tree, Based on Number of Trees Included")

  • Greatest gains by adding a small number of additional trees

  • Moderately small gains thereafter