The Goal

Today, we are focused on classification trees. We are going to deepen our understanding of how trees work, and we are going to fit and visualize tree models in R.

Formatting your lab

When you open your Markdown file, your first chunk likely looks like:

knitr::opts_chunk$set(echo = TRUE)

Change it to:

knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE, fig.asp =.5)

The Data Set

We are going to stick with our cell phone review data from the last lab. Recall that our client provided us a data set called train with 1000 phones that are labelled with their brand. They ask us to use those phones to build a model to (1) predict which brand a phone is and (2) describe any key traits in reviews that might suggest a phone is from a certain brand.

They then also provide us a data set called test with 500 phones to use to assess the predictive abilities of our model.

Let’s start off by loading our libraries. Most of these will look familiar, but a few are new. The new ones are the libraries that we need for classification trees. REMEMBER: If you get an error that says your computer does not have or cannot find a package, go to the top of your screen (the right of your camera for my Apple folks) and choose Tools. From there, choose Install Packages. From there, type in the name of the packages you need and install!

# Load the libraries
library(dplyr)
library(stringr)
library(ggplot2)
library(tidyr)
library(tidytext)
library(nnet)
library(tm)

# New - needed for trees
library(rpart)
library(rattle)
library(rpart.plot)

Once you have loaded the libraries, you can use the following code to load the data.

# Load the data
train <- read.csv("https://www.dropbox.com/scl/fi/y9srlcqu36qh4x2gxx1mf/CellPhonetrain.csv?rlkey=qc1bgktk8kdg8xo0uwhx9gbiv&st=bp2sdljz&dl=1")
test  <- read.csv("https://www.dropbox.com/scl/fi/i875r45iwefli11b8bhbq/CellPhonetest.csv?rlkey=ux8hcgfrb7dagm1lirmhszlxz&st=4gb6yu3q&dl=1")
train$brand <-as.factor(train$brand)
test$brand <-as.factor(test$brand)

Just as we did in Lab 9, our goal is to use the 118 features from LIWC features to (1) predict which phone brand each review is about and (2) to identify and discuss key traits in the text that identify reviews from each brand.

Today, we are going to do this using classification trees.

Classification Tree: One Split

Each step we use to build a tree is a question. These questions are thing like “If the word count is less than 50, is a review more likely to have been written about Apple, Samsung, or Motorola?” This means that trees make no assumptions about the shape of the relationship between X and Y. Instead, trees assume that we can use the features to group the data into authors using a series of splitting rules.

To see this, let’s build a tree with only one split.

tree_onesplit <- rpart( brand ~ . , data = train, method = "class", maxdepth = 1)

This looks a lot like the code we have been using all semester with:

• brand: our \(Y\) variable
• .: meaning we use all other variables as features
• method= "class": tells R to fit a classification tree.
• maxdepth = 1: Tells R right now we only want one splitting rule (so 2 leaves)

To visualize the tree, we can use

rpart.plot(tree_onesplit, sub = "" )

Here, Perception measures the percent of words in a review that are related to perception. These are words about how a phone feels, sounds, looks, etc. In other words, these are words about how a person physically interacts with the device.

Right now, this tree is very clear. Only one feature is being used, and if we needed to chat with clients, we could easily explain how to read the tree, both for predictions and for explaining the relationships highlighted in the tree.

However…there were 118 different features to choose from. Why did we choose perception words?? And why 11%??

rpart_splits <- function(fit, digits = getOption("digits")) {
  splits <- fit$splits
  if (!is.null(splits)) {
    ff <- fit$frame
    is.leaf <- ff$var == "<leaf>"
    n <- nrow(splits)
    nn <- ff$ncompete + ff$nsurrogate + !is.leaf
    ix <- cumsum(c(1L, nn))
    ix_prim <- unlist(mapply(ix, ix + c(ff$ncompete, 0), FUN = seq, SIMPLIFY = F))
    type <- rep.int("surrogate", n)
    type[ix_prim[ix_prim <= n]] <- "primary"
    type[ix[ix <= n]] <- "main"
    left <- character(nrow(splits))
    side <- splits[, 2L]
    for (i in seq_along(left)) {
      left[i] <- if (side[i] == -1L)
        paste("<", format(signif(splits[i, 4L], digits)))
      else if (side[i] == 1L)
        paste(">=", format(signif(splits[i, 4L], digits)))
      else {
        catside <- fit$csplit[splits[i, 4L], 1:side[i]]
        paste(c("L", "-", "R")[catside], collapse = "", sep = "")
      }
    }
    cbind(data.frame(var = rownames(splits),
                     type = type,
                     node = rep(as.integer(row.names(ff)), times = nn),
                     ix = rep(seq_len(nrow(ff)), nn),
                     left = left),
          as.data.frame(splits, row.names = F))
  }
}

gini <- function(tree,data,Y){
  
  variable <-rpart_splits(tree)[1,"var"]
  value <-rpart_splits(tree)[1,"index"]
  #print(length(rpart_splits(tree)[1,"left"]))
  if(length(rpart_splits(tree)[1,"left"]) >0 ){
    sign <- unlist(strsplit(rpart_splits(tree)[1,"left"], " "))[1]
    
    if(sign == ">=" ){
      Leaf <- ifelse( data[,variable] >= value, "Leaf 1", "Leaf 2")
    }
    
    if(sign == ">" ){
      Leaf <- ifelse( data[,variable] > value, "Leaf 1", "Leaf 2")
    }
    
    if(sign == "<=" ){
      Leaf <- ifelse( data[,variable] <= value, "Leaf 1", "Leaf 2")
    }
    
    if(sign == "<" ){
      Leaf <- ifelse( data[,variable] < value, "Leaf 1", "Leaf 2")
    }
    
    weights <- data.frame(table(Leaf)/length(Leaf))
    
    Leaf1 <- subset(data, Leaf== "Leaf 1")
    Leaf2 <- subset(data, Leaf== "Leaf 2")
    
    props1 <- data.frame(table(Leaf1[,Y])/nrow(Leaf1))
    props2 <- data.frame(table(Leaf2[,Y])/nrow(Leaf2))
    
    Impurity1 <- 1 -(sum(props1$Freq^2))
    Impurity2 <- 1 -(sum(props2$Freq^2))
    
    out <- weights$Freq[1]*Impurity1 + weights$Freq[2]*Impurity2
  }
  
  if(length(rpart_splits(tree)[1,"left"]) ==0){
    out <- NULL
  }
  
  out
  
  
}

gini(tree_onesplit, train,1)

Classification Tree 2: More Splits

When we grew our tree with one split, we used:

tree_onesplit <- rpart( brand ~ . , data = train, method = "class", maxdepth = 1)

The maxdepth = 1 argument is the part that tells R to stop after one split. If we remove it, the tree grows more!

tree_full <- rpart( brand ~ . , data = train, method = "class")

You will notice that this tree did not use all the features. Why did it stop growing???

Trees stop growing when they hit what is called a stopping rule. These are a series of rules that we choose that helps us keep trees from getting too large. Super larger trees (meaning trees with many leaves) are hard to interpret and generally do not predict well on test data…which is usually the whole point of prediction.

I don’t see any stopping rules in our code, though. Where did we specify these rules?

It turns out that R has done this for us.

tree_full$control

## $minsplit
## [1] 20
## 
## $minbucket
## [1] 7
## 
## $cp
## [1] 0.01
## 
## $maxcompete
## [1] 4
## 
## $maxsurrogate
## [1] 5
## 
## $usesurrogate
## [1] 2
## 
## $surrogatestyle
## [1] 0
## 
## $maxdepth
## [1] 30
## 
## $xval
## [1] 10

See all of this output?? All of these are default stopping rules in R. What do they all mean?

Here are the important stopping rules:

• minsplit= the minimum number of rows that must exist in a leaf in order for a split to be attempted.
• minbucket = the minimum number of rows allowed in any leaf.
• cp = the percent improvement in the Gini Index we require in order to accept a split. Right now, this is .01 = 1%. This means that if we cannot improve the Gini Index by at least 1%, we do not split a leaf (meaning we stop growing).
• maxdepth = the maximum depth of any leaf, with the root counted as depth 0. Basically, this is how we control how many splitting rules you go through to get from the root to a leaf.

Having default stopping rules is great, but we do need to think about whether or not they make sense for our data. Suppose I had 3 million rows of data. Do we really want to allow leaves as small as 7 rows??? The tree would be huge!!

If we want to change one of these stopping rules when growing a tree, we can adapt our tree growing code. For example, if we want to change the maxdepth stopping rule to 5, we use:

tree_new <- rpart( brand ~ . , data = train, maxdepth = 5)

For the rest of the lab, use tree_50 for all questions.

Note: If you are finding your tree hard to read, consider:

prp(tree_50)

As we can see, classification trees have a few nice properties! Even though we started with 118 features, only a few of them are actually used in our model. This makes it fairly easy to visualize and interpret the relationships in the data.

Predictions

Now that we have built and interpreted the tree, let’s use it for making predictions. We can (1) make predictions on the training data to see how well our model fits the data, and (2) we make predictions on test data. Let’s start with the training data.

To make predictions using a tree, we use the code:

yhat_test <- predict(tree_50, type = "class", newdata = test)

To make our confusion matrix, we use the same code as we did for multinomial:

holder <- table("Prediction"=yhat_test,"Actual" = test$brand)

# Label the rows 
row.names(holder) <- c("Predicted: Apple", "Predicted: Motorola", "Predicted: Samsung")

# Label the columns
colnames(holder) <- c("True: Apple", "True: Motorola", "True: Samsung")

# Format the confusion matrix
knitr::kable(holder)

References