library(tidyverse)
library(knitr)
library(caret)
times <- read_csv('../data/waiting_times.csv')
times %>% head(5) %>% kable()

wait_model <- lm(waiting ~ 0 + pub, data=times)

times %>% group_by(pub) %>% summarize(mean(waiting))

## # A tibble: 3 x 2
##   pub              `mean(waiting)`
##   <chr>                      <dbl>
## 1 The Berkeley                4.70
## 2 The Gallimaufry             1.40
## 3 The Milk Thistle            6.27

ggplot(times, aes(x=pub, y=waiting, color=pub)) + 
  geom_point()

library(rpart)
library(rpart.plot)
dtree <- rpart(waiting ~ pub, data=times, 
               control=list(minsplit=1, # override the default minimum split size
                            minbucket=1 # override the default minimum leaf size
                            )
               )
prp(dtree)

library(sf)

## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1

weca = read_sf("../data/weca.gpkg") %>% 
  st_drop_geometry() %>% 
  pivot_longer(price_dec_1995:price_dec_2018, 
               names_to=c(NA, 'quarter', 'year'), 
               values_to='price', names_sep='_')

dtree_weca <- rpart(price ~ year + quarter + la_name, 
                    data=weca %>% mutate(price = price/1000) # make the plot easier to read
                    )
prp(dtree_weca, faclen=20)

Overfitting and random forests

As you can tell, this procedure may make this model a bit too sensitive to the input data. The precision with which data can be sliced and diced may make for very predictions with very low prediction bias, sure.22 Remember: we mean bias in the sense that the predictions are fairly close to reality, not that the predictions are systematically over (or under) the correct value. But, the prediction variance can be immense. For example, in the weca tree above, just changing from 2000 to 2001 will nearly double the predicted house price!

So, random forests provide a way to collect together decision trees to reduce this prediction variance, while keeping the bias acceptably low. Basically, a random forest works by fitting many trees33 Hence “forest,” a good name for a collection of trees 😏 , and giving each tree a slightly different view of the data. They do this by randomizing the features that the decision trees “see.” For example, let’s think about a random forest trained to predict song genre using danceability, tempo, acousticness, instrumentalness, energy, and valence.

songs = read_csv("../data/midterm-songs.csv")

## Rows: 32833 Columns: 23

## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (10): track_id, track_name, track_artist, track_album_id, track_album_na...
## dbl (13): track_popularity, danceability, energy, key, loudness, mode, speec...

## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

dance_rf = train(playlist_genre ~ danceability + tempo + acousticness + instrumentalness + energy + valence, 
              data=songs, 
              method='ranger', # a faster random forest than method='rf'
              importance='impurity', # save measures of how important each feature is
              trControl= trainControl(method='repeatedcv', repeats=1, number=5))

## Growing trees.. Progress: 87%. Estimated remaining time: 4 seconds.
## Growing trees.. Progress: 60%. Estimated remaining time: 20 seconds.
## Growing trees.. Progress: 97%. Estimated remaining time: 1 seconds.
## Growing trees.. Progress: 66%. Estimated remaining time: 15 seconds.
## Growing trees.. Progress: 89%. Estimated remaining time: 3 seconds.
## Growing trees.. Progress: 64%. Estimated remaining time: 17 seconds.
## Growing trees.. Progress: 91%. Estimated remaining time: 3 seconds.
## Growing trees.. Progress: 53%. Estimated remaining time: 27 seconds.
## Growing trees.. Progress: 93%. Estimated remaining time: 4 seconds.

In this forest, each decision tree “sees” a different set of variables. Some trees may look at danceability, tempo, and energy (like we’ve seen before), but others may look at valence, instrumentalness, and acousticness. By ensuring that the trees see many different perspectives of the data, we can ensure that they are relatively independent. That is, they will not all use the same splits, nor even see the same basic picture of the dataset.44 li id="fn4">

This idea, that a collection of independent but “weak” prediction algorithms (like Decision Trees) can perform learning tasks well, is a central idea of ensemble learning, of which Random Forests are a central component.↩︎

This ensures that the trees won’t all use the same splits and (hopefully) will make it so that the decision trees don’t all make the same mistakes.55 li id="fn5">

In addition, another strategy that is used by default in Random Forest-trained models is called “Bootstrap Aggregation,” or “bagging,” for short. This is discussed in ISL 8.2.1, and (at its core) does a similar thing to the rows of the training data that the random forest does to the columns of a dataframe.↩︎

Importance

One very nice trait of random forests is that they allow you to get a measurement of the importance of a variable for the prediction. This is kind of like an effect size (from a regression framework). Usually, the “best” feature is rated 100. For any model fit in the caret framework, you can get feature importance using the varImp() function (short for variable Importance):

varImp(dance_rf)

## ranger variable importance
## 
##                  Overall
## tempo             100.00
## danceability       75.82
## energy             51.12
## valence            47.44
## acousticness       43.33
## instrumentalness    0.00

thus, you can see that the tempo variable is the most useful in predicting genre. The next most important feature is danceability. Energy, valence, and acousticness all hang around the same value. And, finally, we see that instrumentalness is not very useful at all in predicting the genre of a song’s playlist. A classic visualization of this data uses a lollipop plot:

varImp(dance_rf) %>% plot()

Your Turn

1. Using the earnings.csv data (which describes yearly earnings in thousands of dollars as a function of demographic variables), fit a single decision tree (using the rpart) to predict earnings using all of the variables in the data except age_band.77 li id="fn7">

   Note, you may need to drop NA values…↩︎

• Plot the structure of the decision tree as we did above.
• Reading the decision tree, what do you think is the most important variable in the data for predicting earnings? What do you think is the least important variable?
• Reading the tree, what is the predicted income of a university-educated 30 year old woman? How about a university-educated 30 year old man?
• Fit a linear model with the same specification as this decision tree. Which has a lower mean squared error?
• Challenge: plot the predicted earnings as age increases for men with a university degree, and for women with a university degree.88 li id="fn8">

  Remember: you can use predict(model, new_data) to make predictions about new data. And, you can set up a range of values from a to b, where each value is separated by step using seq(a,b,step)↩︎

2. Train a random forest with the same model specification as in Part 1 using whatever crossvalidation procedure you like.

• make a variable importance plot. Does this reflect the same structure as you saw in question 1?
• Does the random forest have lower mean squared error?
• Challenge: Why might the random forest with crossvalidation not be directly comparable to the tree you fit in question 1?99 li id="fn9">

  Think hard about what the MSE measures in part 1, versus what the MSE shown by caret measures.↩︎

• Challenge: compare this forest’s predictions to the tree’s predictions as a function of age for men or women with a university degree.