Introduction

This assignment will cover tree-based method exercises from Chapter 8 of the ISLR2 textbook. The methods I work with are random forest, boosting, decision trees, and regression trees. I utilize plots, rpart plots, importance plots, confusion matrix’s, and a variety of other tools in this lab.


Question 3.

Consider the Gini index, classification error, and entropy in a simple classification setting with two classes. Create a single plot that displays each of these quantities as a function of pˆm1. The x-axis should display pˆm1, ranging from 0 to 1, and the y-axis should display the value of the Gini index, classification error, and entropy. Hint: In a setting with two classes, pˆm1 =1 − pˆm2. You could make this plot by hand, but it will be much easier to make in R.

# Code for question 3 provided by Proffessor Joseph Campbell
p=seq(0,1,0.0001)
#Gini
G=2*p*(1-p)
#Classification Error
E=1-pmax(p,1-p)
#Entropy
D=-(p*log(p) + (1-p)*log(1-p))

plot(p,D, col="red",ylab="")
lines(p,E,col='green')
lines(p,G,col='blue')
legend(0.3,0.15,c("Entropy", "Missclassification","Gini"),lty=c(1,1,1),lwd=c(2.5,2.5,2.5),col=c('red','green','blue'))

Question 8. Carseats Data

In the lab, a classification tree was applied to the Carseats data set after converting Sales into a qualitative response variable. Now we will seek to predict Sales using regression trees and related approaches, treating the response as a quantitative variable.

attach(Carseats)

# Visualize the response variable
z<-Carseats |> ggplot(
    aes(
        Sales)
    ) + 
    geom_histogram(
        col="gray", 
        fill ="navyblue",
        ) +
    labs(x = "Sales", 
         y = "Count", 
         title = "Distribution of Carseat Sales Units in Thousands (per location)")

ggplotly(z)

Initial Thoughts: The response variable Sales seems to be normally distributed with a high volume spike in Sales in the 5 - 9 (thousand) ranges


(a) Test/train split

Split the data set into a training set and a test set.

set.seed(1) 

# create the test/train split. Here I perform a 70/30 split
train <- createDataPartition(Carseats$Sales, p=0.7, list=FALSE)

train_data <- Carseats[train, ]
test_data <- Carseats[-train, ]

(b) Regression Trees

Fit a regression tree to the training set. Plot the tree, and interpret the results. What test MSE do you obtain?

# fit the regression tree
tree_fit <- tree(Sales ~., Carseats, subset = train)
summary(tree_fit)
## 
## Regression tree:
## tree(formula = Sales ~ ., data = Carseats, subset = train)
## Variables actually used in tree construction:
## [1] "ShelveLoc"   "Price"       "Income"      "Advertising" "CompPrice"  
## [6] "Age"        
## Number of terminal nodes:  18 
## Residual mean deviance:  2.22 = 583.8 / 263 
## Distribution of residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -3.64200 -0.96000  0.02174  0.00000  0.96000  3.40700
plot(tree_fit)
text(tree_fit, 
     pretty=0, 
     cex = .65, 
     digits =3)

Initial Thoughts: This tree tells us a lot about which variables do the best job of predicting Sales. For readability remember that Left = Yes (the condition is TRUE) and Right = No (the condition is FALSE). Starting at the top of the tree and consider the following:

  • ShelveLoc (the quality of the shelving location for the car seats at each site) is the most important predictor of Sales

  • When ShelveLoc is Good, Price is less than 105.5, Age is less than 49.5, income is less than 61, then the predicted value of Sales is 8.43 (in thousands)

  • When ShelveLoc is Bad/Medium and Price is less than 109.5 then the predicted value of sales is 12.20 (in thousands)

  • When SelveLoc is Bad/Medium, price is greater than 109.5, advertising is greater than 0.5, and income is less than 35 then predicted value of Sales is 7.55 (in thousands)

# Calculate predictions
yhat.reg <- predict(tree_fit, newdata = Carseats[-train, ])
Carseats.test <- Carseats[-train, "Sales"]

# plot the predictions
plot(yhat.reg, 
     Carseats.test,
     xlab = "Predicted Sales",
     ylab = "Actual Sales",
     main = "Single Regression Tree - Actual vs Predicted Sales")
abline(0,1, col = "red", lwd = 2)

# Calculate test MSE
mean((yhat.reg-Carseats.test)^2)
## [1] 4.469011
sqrt(mean((yhat.reg-Carseats.test)^2))
## [1] 2.114003

Initial Thoughts: The test MSE comes out to 4.47 if we take the square root we get 2.11. Meaning that the sales predictions of this model are usually within 2.11 units off from the actual Sales values.

(c) Does Pruning do anything?

Use cross-validation in order to determine the optimal level of tree complexity. Does pruning the tree improve the test MSE?

# use CV to see whether pruning the tree will improve it
cv.carseats <- cv.tree(tree_fit)

# using the RSS (dev) calculate MSE
mse <- cv.carseats$dev / length(train)

# plot MSE
plot(cv.carseats$size, 
     mse, 
     type = 'b',
     xlab = "Tree Node Size",
     ylab = "MSE",
     main = "The Affect of Tree Node Size on MSE")

Initial Thoughts: The optimal level of tree complexity would be 19 nodes. The MSE sits right at 4.526317 for a tree with 18 node levels, slightly larger than the overall test MSE of 4.21. So, in this case pruning the tree does NOT improve MSE.

(d) Bagging

Use the bagging approach in order to analyze this data. What test MSE do you obtain? Use the importance() function to determine which variables are most important.

set.seed(1)

#use bagging on Random Forest
bag.carseats <- randomForest(Sales ~ ., 
                             data = Carseats, 
                             subset = train, 
                             mtry = 10,
                             importance = TRUE)

# Make predictions on the test set
yhat.bag <- predict(bag.carseats, newdata = test_data)
plot(yhat.bag, 
     test_data$Sales,
     xlab = "Predicted Sales",
     ylab = "Actual Sales",
     main = "Bagging - Actual Sales vs Predicted Sales",
     abline(0,1, col = "red", lwd = 2))

# Obtain test MSE
mean((yhat.bag - test_data$Sales)^2)
## [1] 2.95016
sqrt(mean((yhat.bag - test_data$Sales)^2))
## [1] 1.717603
importance(bag.carseats)
##                %IncMSE IncNodePurity
## CompPrice   29.3462503    184.444145
## Income       2.1950934     84.971443
## Advertising 28.4461920    226.470679
## Population  -0.7596617     67.714326
## Price       68.7814941    651.996671
## ShelveLoc   71.7824230    628.164064
## Age         21.3597782    166.657371
## Education    1.4885394     55.640049
## Urban       -3.4002675      8.580206
## US           5.0762780     12.380845
varImpPlot(bag.carseats)

Initial Thoughts: The test MSE of the bagged model turned out to be 2.95. Smaller than the regression tree model. The predictions also seem close to the prediction line indicating the model is doing a fairly good job predicting Sales. The test predictions in this dataset would be off when predicting sales on average by 1.72 units of Sales.

Importance Plots:

  • IncMSA - This plot measures which variables are the most importance for making predictions. ShelveLoc would be the most important variale in predicting sales with a very high MSE

  • IncNodePurity - This plot measures which variables help create the best tree splits. Price would be the most important factor in creating a tree.

(e) Random Forest

Use random forests to analyze this data. What test MSE do you obtain? Use the importance() function to determine which variables are most important. Describe the effect of m, the number of variables considered at each split, on the error rate obtained.

fit_control <- trainControl(
  method = "cv",
  number = 5
  #,classProbs = TRUE,
  #summaryFunction = twoClassSummary
)

set.seed(42)

rf_fit <- train(
  Sales ~.,
  data = train_data,
  method = "rf",
  #metric "ROC",
  trControl = fit_control,
  tuneGrid = expand.grid(mtry = c(1, 2, 3)),
  ntree = 300
)

plot(varImp(rf_fit))

yhat.rf <- predict(rf_fit, newdata = test_data)
plot(yhat.rf, 
     test_data$Sales,
     xlab = "Predicted Sales",
     ylab = "Actual Sales",
     main = "Random Forest - Actual Sales vs Predicted Sales",
     abline(0,1, col = "red", lwd = 2))

# Obtain test MSE
mean((yhat.rf - test_data$Sales)^2)
## [1] 4.065536
sqrt(mean((yhat.rf - test_data$Sales)^2))
## [1] 2.016317

Initial Thoughts: The importance plot shows that Price and ShelveLocGood are the most important predictor variables in all the random forests. Notice how the plot split Shelveloc into Medium and Good because it was a categorical variable. Also, mtry = 3 was the most effective way to split the trees. Meaning that before making a split the trees consider 3 variables as apposed to 1 or 2. The RMSE of 1.76, R^2 of .66 and MAE of 1.41 prove this. Also, the test MSE of the Randomforest is the largest MSE so far of 4.07.

(f) BART

#Fit the BART model
x <- Carseats[, -1]
y <- Carseats[, "Sales"]
xtrain <- x[train, ]
ytrain <- y[train]
xtest <- x[-train, ]
ytest <- y[-train]
set.seed(1)
bartfit <- gbart(xtrain, ytrain, x.test = xtest)
# Calculate MSE
yhat.bart <- bartfit$yhat.test.mean
mean((ytest - yhat.bart)^2)
## [1] 1.73499
yhat.bart <- bartfit$yhat.test.mean

plot(yhat.bart, 
     test_data$Sales,
     xlab = "Predicted Sales",
     ylab = "Actual Sales",
     main = "BART - Actual Sales vs Predicted Sales",
     abline(0,1, col = "red", lwd = 2))

Initial Thoughts: The BART model had an MSE of 1.73 which has been the best MSE thus far. See the table below for comparisons.


Comparisons

single_reg<-mean((yhat.reg-Carseats.test)^2)

bag<-mean((yhat.bag - test_data$Sales)^2)

rf<-mean((yhat.rf - test_data$Sales)^2)

bart<-mean((ytest - yhat.bart)^2)

datatable(
  data.frame(
    Model = c("Single Regression Tree", 
              "Bagging", 
              "Random Forest",
              "BART"),
    MSE = c(single_reg, bag, rf, bart)
  ),
  caption = "Comparison of Tree-Based Model Performance",
  rownames = FALSE,
  options = list(dom = "t")
)

Question 9. OJ Data

This problem involves the OJ data set which is part of the ISLR2 package.

(a) Test/Train Split

Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations.

attach(OJ)

set.seed(1) 

# Here I perform the test/train split
oj_train <- createDataPartition(OJ$Purchase, p=0.7, list=FALSE)

oj_train_data <- OJ[oj_train, ]
oj_test_data <- OJ[-oj_train, ]

(b) Decision Tree

Fit a tree to the training data, with Purchase as the response and the other variables as predictors. Use the summary() function to produce summary statistics about the tree, and describe the results obtained. What is the training error rate? How many terminal nodes does the tree have?

# create the tree
full_tree <- rpart(Purchase~.,
                   data=oj_train_data, 
                   method = "class",
                   control=rpart.control())
summary(full_tree)
## Call:
## rpart(formula = Purchase ~ ., data = oj_train_data, method = "class", 
##     control = rpart.control())
##   n= 750 
## 
##           CP nsplit rel error    xerror       xstd
## 1 0.48287671      0 1.0000000 1.0000000 0.04573100
## 2 0.02568493      1 0.5171233 0.5479452 0.03842134
## 3 0.02397260      3 0.4657534 0.5650685 0.03885138
## 4 0.02054795      4 0.4417808 0.5513699 0.03850854
## 5 0.01369863      5 0.4212329 0.5342466 0.03806642
## 6 0.01027397      6 0.4075342 0.5102740 0.03742113
## 7 0.01000000      7 0.3972603 0.5239726 0.03779370
## 
## Variable importance
##        LoyalCH        StoreID      PriceDiff    SalePriceMM        PriceMM 
##             47              9              8              6              5 
##         DiscMM      PctDiscMM  ListPriceDiff        PriceCH WeekofPurchase 
##              5              5              4              3              3 
##    SalePriceCH          STORE         Store7 
##              2              2              1 
## 
## Node number 1: 750 observations,    complexity param=0.4828767
##   predicted class=CH  expected loss=0.3893333  P(node) =1
##     class counts:   458   292
##    probabilities: 0.611 0.389 
##   left son=2 (425 obs) right son=3 (325 obs)
##   Primary splits:
##       LoyalCH     < 0.5036    to the right, improve=123.09670, (0 missing)
##       StoreID     < 3.5       to the right, improve= 32.56409, (0 missing)
##       PriceDiff   < 0.31      to the right, improve= 19.87539, (0 missing)
##       SalePriceMM < 2.04      to the right, improve= 17.34553, (0 missing)
##       DiscCH      < 0.255     to the right, improve= 15.60231, (0 missing)
##   Surrogate splits:
##       StoreID        < 3.5       to the right, agree=0.648, adj=0.188, (0 split)
##       PriceMM        < 1.89      to the right, agree=0.596, adj=0.068, (0 split)
##       WeekofPurchase < 246.5     to the right, agree=0.595, adj=0.065, (0 split)
##       ListPriceDiff  < 0.035     to the right, agree=0.591, adj=0.055, (0 split)
##       PriceCH        < 1.72      to the right, agree=0.588, adj=0.049, (0 split)
## 
## Node number 2: 425 observations,    complexity param=0.0239726
##   predicted class=CH  expected loss=0.1388235  P(node) =0.5666667
##     class counts:   366    59
##    probabilities: 0.861 0.139 
##   left son=4 (404 obs) right son=5 (21 obs)
##   Primary splits:
##       PriceDiff   < -0.39     to the right, improve=12.310240, (0 missing)
##       LoyalCH     < 0.7645725 to the right, improve= 9.114571, (0 missing)
##       DiscMM      < 0.57      to the left,  improve= 8.216284, (0 missing)
##       PctDiscMM   < 0.264375  to the left,  improve= 8.216284, (0 missing)
##       SalePriceMM < 1.585     to the right, improve= 7.859655, (0 missing)
##   Surrogate splits:
##       DiscMM      < 0.72      to the left,  agree=0.976, adj=0.524, (0 split)
##       SalePriceMM < 1.435     to the right, agree=0.976, adj=0.524, (0 split)
##       PctDiscMM   < 0.3342595 to the left,  agree=0.976, adj=0.524, (0 split)
##       SalePriceCH < 2.075     to the left,  agree=0.958, adj=0.143, (0 split)
## 
## Node number 3: 325 observations,    complexity param=0.02568493
##   predicted class=MM  expected loss=0.2830769  P(node) =0.4333333
##     class counts:    92   233
##    probabilities: 0.283 0.717 
##   left son=6 (181 obs) right son=7 (144 obs)
##   Primary splits:
##       LoyalCH   < 0.282272  to the right, improve=14.082430, (0 missing)
##       PriceDiff < 0.31      to the right, improve= 8.237351, (0 missing)
##       SpecialCH < 0.5       to the right, improve= 5.851637, (0 missing)
##       DiscCH    < 0.255     to the right, improve= 5.607657, (0 missing)
##       PctDiscCH < 0.132882  to the right, improve= 5.607657, (0 missing)
##   Surrogate splits:
##       STORE       < 2.5       to the left,  agree=0.609, adj=0.118, (0 split)
##       PriceCH     < 1.875     to the left,  agree=0.594, adj=0.083, (0 split)
##       SalePriceCH < 1.875     to the left,  agree=0.594, adj=0.083, (0 split)
##       StoreID     < 3.5       to the right, agree=0.569, adj=0.028, (0 split)
##       PriceMM     < 2.205     to the left,  agree=0.566, adj=0.021, (0 split)
## 
## Node number 4: 404 observations
##   predicted class=CH  expected loss=0.1113861  P(node) =0.5386667
##     class counts:   359    45
##    probabilities: 0.889 0.111 
## 
## Node number 5: 21 observations,    complexity param=0.01027397
##   predicted class=MM  expected loss=0.3333333  P(node) =0.028
##     class counts:     7    14
##    probabilities: 0.333 0.667 
##   left son=10 (9 obs) right son=11 (12 obs)
##   Primary splits:
##       LoyalCH < 0.742157  to the right, improve=3.5000000, (0 missing)
##       STORE   < 1.5       to the right, improve=3.0476190, (0 missing)
##       StoreID < 1.5       to the right, improve=0.7619048, (0 missing)
##       Store7  splits as  LR, improve=0.7619048, (0 missing)
##       PriceCH < 1.975     to the right, improve=0.1696970, (0 missing)
##   Surrogate splits:
##       PriceCH     < 2.04      to the right, agree=0.714, adj=0.333, (0 split)
##       DiscMM      < 0.47      to the left,  agree=0.714, adj=0.333, (0 split)
##       SalePriceMM < 1.64      to the right, agree=0.714, adj=0.333, (0 split)
##       SalePriceCH < 2.04      to the right, agree=0.714, adj=0.333, (0 split)
##       PctDiscMM   < 0.2224545 to the left,  agree=0.714, adj=0.333, (0 split)
## 
## Node number 6: 181 observations,    complexity param=0.02568493
##   predicted class=MM  expected loss=0.4143646  P(node) =0.2413333
##     class counts:    75   106
##    probabilities: 0.414 0.586 
##   left son=12 (107 obs) right son=13 (74 obs)
##   Primary splits:
##       PriceDiff     < 0.05      to the right, improve=12.694000, (0 missing)
##       SalePriceMM   < 1.94      to the right, improve=10.544100, (0 missing)
##       DiscMM        < 0.22      to the left,  improve= 6.766872, (0 missing)
##       PctDiscMM     < 0.0729725 to the left,  improve= 5.775128, (0 missing)
##       ListPriceDiff < 0.235     to the right, improve= 4.614858, (0 missing)
##   Surrogate splits:
##       SalePriceMM   < 1.94      to the right, agree=0.950, adj=0.878, (0 split)
##       DiscMM        < 0.08      to the left,  agree=0.834, adj=0.595, (0 split)
##       PctDiscMM     < 0.038887  to the left,  agree=0.834, adj=0.595, (0 split)
##       ListPriceDiff < 0.135     to the right, agree=0.779, adj=0.459, (0 split)
##       PriceMM       < 2.04      to the right, agree=0.773, adj=0.446, (0 split)
## 
## Node number 7: 144 observations
##   predicted class=MM  expected loss=0.1180556  P(node) =0.192
##     class counts:    17   127
##    probabilities: 0.118 0.882 
## 
## Node number 10: 9 observations
##   predicted class=CH  expected loss=0.3333333  P(node) =0.012
##     class counts:     6     3
##    probabilities: 0.667 0.333 
## 
## Node number 11: 12 observations
##   predicted class=MM  expected loss=0.08333333  P(node) =0.016
##     class counts:     1    11
##    probabilities: 0.083 0.917 
## 
## Node number 12: 107 observations,    complexity param=0.02054795
##   predicted class=CH  expected loss=0.4299065  P(node) =0.1426667
##     class counts:    61    46
##    probabilities: 0.570 0.430 
##   left son=24 (53 obs) right son=25 (54 obs)
##   Primary splits:
##       STORE          < 1.5       to the left,  improve=3.442309, (0 missing)
##       PriceDiff      < 0.49      to the right, improve=1.997351, (0 missing)
##       PriceMM        < 2.205     to the left,  improve=1.879220, (0 missing)
##       SalePriceMM    < 2.205     to the left,  improve=1.879220, (0 missing)
##       WeekofPurchase < 249.5     to the right, improve=1.713074, (0 missing)
##   Surrogate splits:
##       StoreID     < 5.5       to the right, agree=0.794, adj=0.585, (0 split)
##       Store7      splits as  RL, agree=0.794, adj=0.585, (0 split)
##       SalePriceCH < 1.775     to the left,  agree=0.748, adj=0.491, (0 split)
##       PriceMM     < 2.155     to the left,  agree=0.720, adj=0.434, (0 split)
##       PriceCH     < 1.875     to the left,  agree=0.710, adj=0.415, (0 split)
## 
## Node number 13: 74 observations
##   predicted class=MM  expected loss=0.1891892  P(node) =0.09866667
##     class counts:    14    60
##    probabilities: 0.189 0.811 
## 
## Node number 24: 53 observations
##   predicted class=CH  expected loss=0.3018868  P(node) =0.07066667
##     class counts:    37    16
##    probabilities: 0.698 0.302 
## 
## Node number 25: 54 observations,    complexity param=0.01369863
##   predicted class=MM  expected loss=0.4444444  P(node) =0.072
##     class counts:    24    30
##    probabilities: 0.444 0.556 
##   left son=50 (18 obs) right son=51 (36 obs)
##   Primary splits:
##       LoyalCH        < 0.49      to the right, improve=1.500000, (0 missing)
##       SalePriceCH    < 1.775     to the right, improve=1.463019, (0 missing)
##       PriceDiff      < 0.31      to the left,  improve=1.121357, (0 missing)
##       WeekofPurchase < 269.5     to the right, improve=1.066667, (0 missing)
##       SalePriceMM    < 2.15      to the left,  improve=1.060652, (0 missing)
##   Surrogate splits:
##       StoreID        < 3.5       to the right, agree=0.704, adj=0.111, (0 split)
##       STORE          < 3.5       to the right, agree=0.704, adj=0.111, (0 split)
##       WeekofPurchase < 252.5     to the left,  agree=0.685, adj=0.056, (0 split)
##       PriceDiff      < 0.43      to the right, agree=0.685, adj=0.056, (0 split)
##       ListPriceDiff  < 0.38      to the right, agree=0.685, adj=0.056, (0 split)
## 
## Node number 50: 18 observations
##   predicted class=CH  expected loss=0.3888889  P(node) =0.024
##     class counts:    11     7
##    probabilities: 0.611 0.389 
## 
## Node number 51: 36 observations
##   predicted class=MM  expected loss=0.3611111  P(node) =0.048
##     class counts:    13    23
##    probabilities: 0.361 0.639

(c) Interpreting the results.

Type in the name of the tree object in order to get a detailed text output. Pick one of the terminal nodes, and interpret the information displayed.

Initial Thoughts: The tree contains 9 nodes and has a training error of 14.5%. Meaning that 14.6% of the observations in the training set are misclassified while the other 85.4% are accurate. The top tree section tells us that members with a loyalty score to great than of .48 generally purchase Citril Hill orange juice and members with a loyalty score less than that usually go with Minute Maid.

(d) Plot

Create a plot of the tree, and interpret the results.

# visualize the tree
rpart.plot(
  full_tree,
  type = 2,
  extra = 104,
  fallen.leaves = TRUE,
  box.palette = "Blues",
  shadow.col = "gray",
  nn = TRUE
)

Initial Thoughts: The following factors tend to influence members in predicting which brand they’re loyal to and purchasing behavior. Members with a loyalty score greater than .76, price difference between CH and MM greater than $0.02, and store ID greater than 4.

(e) Confusion Matrix and Final Thoughts

Predict the response on the test data, and produce a confusion matrix comparing the test labels to the predicted test labels. What is the test error rate?

# create the tree
pred <- predict(
  full_tree, 
  newdata = oj_test_data, 
  type = "class"
                )

# confusion matrix
confusion <- table(
  Actual = oj_test_data$Purchase,
  Predicted = pred
)

confusion
##       Predicted
## Actual  CH  MM
##     CH 168  27
##     MM  28  97
mean(pred != oj_test_data$Purchase)
## [1] 0.171875

Final Thoughts: The confusion matrix shows that the decision tree model predicted well on the test set. With a 13.3% error rate we can see that the vast majority of the predictions were accuately predicted.