The Goal

Last time, we talked about Best Subset Selection, which is an approach we use when we want to fit an LSLR model but we don’t know which of our X variables we want to include in that model. Best subset selection checks every possible subset of our X variables and finds the model that optimizes some metric, like the AIC or the \(R^2_{adj}\).

Today, we are are going to learn a different technique for building models when our Y is numeric. This technique is called a regression tree. Like BSS, regression trees are useful when we have a lot of potential X variables and we are not sure which we want to use. Let’s explore.

The Data

Instead of our typical data for this course, today we are going to work with a data set containing information on n = 333 penguins from three different species. To load the data, you need to use the following two lines of code:

library(palmerpenguins)
data(penguins)
penguins <- na.omit(penguins)

We have a client who is interested in building a model for Y = the body mass of a penguin in grams.

In addition to this response variable, we have information on 7 explanatory variables.

Building a Model

Our goal at the moment is to model body mass, which is a numeric response variable. This means we might tend to reach for a regression model. In our last class, we learned that Best Subset Selection (BSS) is a good tool to use when we want to build a regression model but we might not know which predictors (explanatory variables) that we want to use in our model.

However, BSS works best when (1) we do not have many strong relationship among our X variables and (2) our numeric X variables have a linear relationship with Y.

Let’s check for correlations among the X variables first. To do this, we need the following library:

suppressMessages(library(corrplot))

Once you have loaded the library, we need to create a correlation matrix, and then plot it.

# Create the correlation matrix
M <- cor(penguins[,-c(1:2,7)])
# Create the plot
corrplot(M, method = "number", type = "upper")

Let’s also check the shape of the relationship between bill depth and body mass.

ggplot(penguins, aes(x= bill_depth_mm, y = body_mass_g, color =species))+ 
  geom_point()+
  labs(x = "Bill Depth (in mm)", y = "Body Mass (in grams)", 
       title = "Figure 1", 
       caption = "A scatter plot of X = bill depth and Y = body mass ")

Do we notice in the picture that there are two defined groups in this graph? The Gentoo penguins seem to have a different group of points from the Adelie and Chinstrap penguins.

This grouping structure within the data, as well as the correlation among our potential X variables, indicates that these data might benefit from a type of model called a regression tree.

What is a Regresstion Tree?

A regression tree is a type of model we can build for a numeric Y variable. We can use numeric or categorical X variables to help model Y.

We start building a tree by looking at all the possible X variables are considering, just like we did in BSS. However, we don’t start building a linear model. Instead, we divide the data into two groups using one X variable. For instance, we might divide the data into male penguins and female penguins, or we might divide the group into penguins with bill length less than 20 mm and penguins with bill length greater than or equal to 20 mm.

How do we decide how to divide the data into these two groups? We choose the division (called a split in trees!) that gives us the smallest residual sum of squares (RSS or SSE). Just like in regular LSLR regression, our job is to minimize the residual sum of squares:

\[RSS = \sum_{i=1}^{n} (y_i -\hat{y}_i),\] where \(n\) is the number of rows in the data and \(\hat{y}_i\) is the value of \(y\) for row \(i\) that we would predict using our model.

Once we have found the first split that minimizes the RSS, we then try to add a second split. Can we make the RSS even smaller if we split our male penguins into those with bill depth greater than 5 mm and those with bill depth less than or equal to 5mm? And so on. We continue the process of dividing the data into groups into we no longer see a meaningful decrease in RSS (typically at least a 1% decrease).

Okay, that is the idea. Trees are really easiest to understand through an example, so let’s try one.

Building a Tree

Now that we have the idea of what we want to do, let’s build the tree! To build trees in R, there are two packages you need. The first allows you to actually build the tree. The second allows us to visualize the tree once we are done.

suppressMessages(library(rpart))
suppressMessages(library(rattle))

To build a tree in R, the function you need is rpart. Let’s stay I want to build a tree using island, species, bill length, and bill depth. The code I need for that is:

tree1 <- rpart(body_mass_g ~ island + species + bill_length_mm + 
                 bill_depth_mm, data = penguins)

To visualize the tree, use:

fancyRpartPlot(tree1, caption = "Tree 1")

You will notice that a a tree is built in layers. Each of the little boxes on the tree are called nodes. The very first node, labeled with a 1 in a square box, is called the root. This is where we start, with all n = 333 penguins in one group all together. Then, we split the penguins into two groups based on the value of one variable.

treeSplit1 <- rpart(body_mass_g ~ island + species + bill_length_mm +
                      bill_depth_mm, maxdepth = 1, data = penguins)
fancyRpartPlot(treeSplit1, caption = "Tree 1: First Split Only")

# Gentoo penguins
gentoo <-  penguins %>%
  filter(species == "Gentoo")

# Adelie and Chinstrap penguins
chinstrap_adelie <-  penguins %>%
  filter(species != "Gentoo")

# Male penguins
mean(gentoo$body_mass_g)

# Female penguins
mean(chinstrap_adelie$body_mass_g)

# Make predictions 
predictions_tree <- predict(treeSplit1, data = penguins)
# Obtain the RSS
sum( (penguins$body_mass_g - predictions_tree )^2)

# Obtain the TSS
sum( (penguins$body_mass_g - mean(penguins$body_mass_g))^2)

## [1] 215259666

treeSplit2 <- rpart(body_mass_g ~ island + species + bill_length_mm +
                      bill_depth_mm, maxdepth = 2, data = penguins)
fancyRpartPlot(treeSplit2, caption = "Tree 1: 3 Splits")

# Make predictions 
predictions_tree <- predict(treeSplit2, data = penguins)
# Obtain the RSS
sum( (penguins$body_mass_g - predictions_tree )^2)

Exploring the Whole Tree

Now that we have explored the basic structure of how a tree works, let’s see the entire tree.

tree1 <- rpart(body_mass_g ~ island + species + bill_length_mm +
                 bill_depth_mm, data = penguins)

fancyRpartPlot(tree1, caption = "Tree 1")

If your task is prediction, we have seen that we can use the tree to help us predict body mass. However, what if our goal is association? In other words, what if we are trying to understand the relationships in the data?

If your tree does not use all the possible X variables that you gave it, this means that adding more splits to the tree was no longer decreasing the residual sum of squares enough for that reduction to be meaningful. Once splitting no longer helps improve the residual sum of squares, we stop growing the tree.

Try it!