Tree-Based Methods

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 ˆpm1. The x-axis should display ˆpm1, 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, ˆpm1 = 1 - ˆpm2. You could make this plot by hand, but it will be much easier in R.

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

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.

a.) Split the data

library(ISLR2)
library(tree)

data(Carseats)

set.seed(1)
train <- sample(1:nrow(Carseats), nrow(Carseats)/2)

carseats.train <- Carseats[train, ]
carseats.test <- Carseats[-train, ]

b.) Regression Tree

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

tree.carseats <- tree(Sales ~ ., data = carseats.train)
summary(tree.carseats)

Regression tree:
tree(formula = Sales ~ ., data = carseats.train)
Variables actually used in tree construction:
[1] "ShelveLoc"   "Price"       "Age"         "Advertising"
[5] "CompPrice"   "US"         
Number of terminal nodes:  18 
Residual mean deviance:  2.167 = 394.3 / 182 
Distribution of residuals:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-3.88200 -0.88200 -0.08712  0.00000  0.89590  4.09900 

plot(tree.carseats)
text(tree.carseats, pretty = 0, cex = .4)

pred <- predict(tree.carseats, newdata = carseats.test)
tree.test.mse <- mean((pred - carseats.test$Sales)^2)

cat("MSE:", tree.test.mse, "\n")

MSE: 4.922039 

The tree first splits on ShelveLoc and Price, indicating the two as the strongest predictors for Sales.

c.) Cross Validation

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

cv.carseats <- cv.tree(tree.carseats)
plot(cv.carseats$size, cv.carseats$dev,
     type = "b",
     xlab = "Tree Size",
     ylab = "CV Deviance")

best.size <- cv.carseats$size[which.min(cv.carseats$dev)]
cat("Best Size:", best.size, "\n") 

Best Size: 5 

prune.carseats <- prune.tree(tree.carseats, best = best.size)
plot(prune.carseats)
text(prune.carseats, pretty = 0, cex = .6)

yhat.prune <- predict(prune.carseats, newdata = carseats.test)
prune.mse <- mean((yhat.prune - carseats.test$Sales)^2)
cat("Prune MSE:", prune.mse, "\n")

Prune MSE: 5.186482 

The error rates for the pruned and unpruned trees are very similar, indicating that pruning the tree would not make a significant difference for this data set.

d.) Bagging Approach

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.

library(randomForest)

set.seed(1)
p <- ncol(carseats.train) - 1

bag.carseats <-  randomForest(Sales ~ ., 
                              data = carseats.train, 
                              mtry = p, importance = TRUE)
bag.carseats

Call:
 randomForest(formula = Sales ~ ., data = carseats.train, mtry = p,      importance = TRUE) 
               Type of random forest: regression
                     Number of trees: 500
No. of variables tried at each split: 10

          Mean of squared residuals: 2.889221
                    % Var explained: 63.26

yhat.bag <- predict(bag.carseats, newdata = carseats.test)
bag.test.MSE <- mean((yhat.bag - carseats.test$Sales)^2)
cat("Bagging MSE:", bag.test.MSE, "\n")

Bagging MSE: 2.605253 

importance(bag.carseats)

               %IncMSE IncNodePurity
CompPrice   24.8888481    170.182937
Income       4.7121131     91.264880
Advertising 12.7692401     97.164338
Population  -1.8074075     58.244596
Price       56.3326252    502.903407
ShelveLoc   48.8886689    380.032715
Age         17.7275460    157.846774
Education    0.5962186     44.598731
Urban        0.1728373      9.822082
US           4.2172102     18.073863

varImpPlot(bag.carseats)

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.

set.seed(1)
mtry.values <- c(2,3,4,5,6,7,p)
rf.test.mse <- rep(NA, length(mtry.values))

for (i in seq_along(mtry.values)) {
  rf.fit <- randomForest(Sales ~ .,
                         data = carseats.train,
                         mtry = mtry.values[i], 
                         importance = TRUE)
  yhat.rf <- predict(rf.fit, 
                     newdata = carseats.test)
  rf.test.mse[i] <- mean((yhat.rf - carseats.test$Sales)^2)
}

data.frame(mtry = mtry.values, test.MSE = rf.test.mse)

plot(mtry.values,
     rf.test.mse,
     type = "b",
     pch = 19,
     xlab = "mtry",
     ylab = "Test MSE",
     main = "Random Forest Test MSE v. # of split variables tested")

best.mtry <- mtry.values[which.min(rf.test.mse)]
rf.best <- randomForest(Sales ~., 
                        data = carseats.test,
                        mtry = best.mtry,
                        importance = TRUE)

importance(rf.best)

              %IncMSE IncNodePurity
CompPrice   17.368496    121.061683
Income       8.099290     77.099105
Advertising 21.295981    146.013186
Population   1.387677     60.365860
Price       53.825774    423.426286
ShelveLoc   65.551343    587.325928
Age         13.650662    106.605549
Education    2.561897     37.091145
Urban        1.161303      9.504955
US           1.318976      7.126472

varImpPlot(rf.best)

From the plot of MSE vs Number of split variables tests, one can gather that the Test error decreases as the number of variables in each split increases. It is noted that the most significant decrease occurred when m increased from 2 to 6, where after the MSE begins to plateau, with the best MSE (2.608) being at m = 10. From the rf.best plots, we conclude that ShelveLoc and Price are the most significant predictor variables, having both the highest increasing MSE and Node Purity. The removal of either variable would substantially decrease prediction accuracy.

f.) BART

Now analyze the data using BART, and report your results

library(BART)
x <- Carseats[, -which(names(Carseats) == "Sales")]

Error in .rs.exprMutatesPackageLibrary(part) : 
  argument "part" is missing, with no default

y <- Carseats$Sales

x <- model.matrix(Sales ~. -1, data = Carseats) |> as.data.frame()

xtrain <- x[train, ]
ytrain <- y[train]
xtest <- x[-train, ]
ytest <- y[-train]

set.seed(1)
bartfit <- gbart(xtrain,
                 ytrain, 
                 x.test = xtest)

*****Calling gbart: type=1
*****Data:
data:n,p,np: 200, 12, 200
y1,yn: 2.781850, 1.091850
x1,x[n*p]: 107.000000, 1.000000
xp1,xp[np*p]: 111.000000, 1.000000
*****Number of Trees: 200
*****Number of Cut Points: 63 ... 1
*****burn,nd,thin: 100,1000,1
*****Prior:beta,alpha,tau,nu,lambda,offset: 2,0.95,0.273474,3,0.23074,7.57815
*****sigma: 1.088371
*****w (weights): 1.000000 ... 1.000000
*****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,12,0
*****printevery: 100

MCMC
done 0 (out of 1100)
done 100 (out of 1100)
done 200 (out of 1100)
done 300 (out of 1100)
done 400 (out of 1100)
done 500 (out of 1100)
done 600 (out of 1100)
done 700 (out of 1100)
done 800 (out of 1100)
done 900 (out of 1100)
done 1000 (out of 1100)
time: 3s
trcnt,tecnt: 1000,1000

#test error
yhat.bart <- bartfit$yhat.test.mean
bart.test.mse <- mean((ytest - yhat.bart)^2)

cat("BART MSE:", bart.test.mse, "\n")

BART MSE: 1.432639 

#check # of times each variable appeared
ord <- order(bartfit$varcount.mean,
             decreasing = TRUE)
bartfit$varcount.mean[ord]

          Price   ShelveLocGood           USYes       CompPrice 
         27.531          21.300          20.896          20.731 
      Education ShelveLocMedium    ShelveLocBad          Income 
         20.059          19.842          18.918          18.831 
            Age      Population     Advertising        UrbanYes 
         18.465          17.913          17.501          17.403 

BART outperformed the other tree-based methods tested, since its approach uses a Bayesian framework in fitting many small trees while regularizing the model. This leads to a reduction in the models tendency to fit noise on the training data, therefore reducing the chance of overfitting. Hence, BART produced the lowest MSE when compared to the rest of the tree-based methods tried.

Question 9

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

data("OJ")

a.) Random training set

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

set.seed(1)
train.oj <- sample(1:nrow(OJ), 800)

oj.train <- OJ[train.oj, ]
oj.test <- OJ[-train.oj, ]

b.) 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?

tree.oj <- tree(Purchase ~ .,
                data = oj.train)
summary(tree.oj)

Classification tree:
tree(formula = Purchase ~ ., data = oj.train)
Variables actually used in tree construction:
[1] "LoyalCH"       "PriceDiff"     "SpecialCH"     "ListPriceDiff"
[5] "PctDiscMM"    
Number of terminal nodes:  9 
Residual mean deviance:  0.7432 = 587.8 / 791 
Misclassification error rate: 0.1588 = 127 / 800 

The variables used for the tree include: “LoyalCH”, “PriceDiff”, “SpecialCH”, “ListPriceDiff”, and “PctDiscMM”. The tree has a total of 9 terminal nodes. and an MSE of 15%.

c.) Terminal Nodes Interpretation

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.

tree.oj

node), split, n, deviance, yval, (yprob)
      * denotes terminal node

 1) root 800 1073.00 CH ( 0.60625 0.39375 )  
   2) LoyalCH < 0.5036 365  441.60 MM ( 0.29315 0.70685 )  
     4) LoyalCH < 0.280875 177  140.50 MM ( 0.13559 0.86441 )  
       8) LoyalCH < 0.0356415 59   10.14 MM ( 0.01695 0.98305 ) *
       9) LoyalCH > 0.0356415 118  116.40 MM ( 0.19492 0.80508 ) *
     5) LoyalCH > 0.280875 188  258.00 MM ( 0.44149 0.55851 )  
      10) PriceDiff < 0.05 79   84.79 MM ( 0.22785 0.77215 )  
        20) SpecialCH < 0.5 64   51.98 MM ( 0.14062 0.85938 ) *
        21) SpecialCH > 0.5 15   20.19 CH ( 0.60000 0.40000 ) *
      11) PriceDiff > 0.05 109  147.00 CH ( 0.59633 0.40367 ) *
   3) LoyalCH > 0.5036 435  337.90 CH ( 0.86897 0.13103 )  
     6) LoyalCH < 0.764572 174  201.00 CH ( 0.73563 0.26437 )  
      12) ListPriceDiff < 0.235 72   99.81 MM ( 0.50000 0.50000 )  
        24) PctDiscMM < 0.196196 55   73.14 CH ( 0.61818 0.38182 ) *
        25) PctDiscMM > 0.196196 17   12.32 MM ( 0.11765 0.88235 ) *
      13) ListPriceDiff > 0.235 102   65.43 CH ( 0.90196 0.09804 ) *
     7) LoyalCH > 0.764572 261   91.20 CH ( 0.95785 0.04215 ) *

Looking at terminal node 7, we note that there are 261 observations, where customers with LoyalCH great than .764 are classified as CH. Here the class probabilities for purchasing CH are 95.8% while there is a 4.2% probability of purchasing MM, with a low deviance of 91.20.

d.) Plot of Tree

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

plot(tree.oj)
text(tree.oj, pretty = 0, cex = .6)

The tree suggests that LoyalCH (customer brand loyalty for CH) has the most significance as a predictor variable. This indicates that customer loyalty is the primary factor in influencing brand choice, followed by PriceDiff, ListPriceDiff, SpecialCH, and PctDiscMM in the case of loyalty not being the singular determining factor.

e.) Confusion Matrix and Test Error

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?

tree.pred <- predict(tree.oj, 
                     newdata = oj.test,
                     type = "class")
conf.mat <- table(Predicted = tree.pred, 
                  Actual = oj.test$Purchase)
conf.mat

         Actual
Predicted  CH  MM
       CH 160  38
       MM   8  64

#test error
test.error <- 1 - sum(diag(conf.mat)) / sum(conf.mat)
cat("Test Error:", test.error*100, "%\n")

Test Error: 17.03704 %

f.) Cross Validation Tree

Apply the cv.tree() function to the training set in order to determine the optimal tree size.

set.seed(1)

cv.oj <- cv.tree(tree.oj, 
            FUN = prune.misclass)
cv.oj

$size
[1] 9 8 7 4 2 1

$dev
[1] 145 145 146 146 167 315

$k
[1]       -Inf   0.000000   3.000000   4.333333  10.500000 151.000000

$method
[1] "misclass"

attr(,"class")
[1] "prune"         "tree.sequence"

g.) Plot Tree Size v Error Rate

Produce a plot with tree size on th ex-axis and cross-validated classification error rate on the y-axis.

plot(cv.oj$size, 
     cv.oj$dev, 
     type = "b", 
     pch = 19,
     xlab = "Tree Size",
     ylab = "CV Error")

h.) Best Tree Size

Which tree size corresponds to the lowest cross-validated classification error rate?

best.size.oj <- cv.oj$size[which.min(cv.oj$dev)]
cat("Best Size:", best.size.oj, "\n")

Best Size: 9 

i.) Pruned Tree for Optimal Tree Size

Produce a pruned tree corresponding to the optimal tree size obtained using cross-validation. If cross-validation does not lead to selection of a pruned tree, then create a pruned tree with five terminal nodes.

size.to.use <- if(best.size.oj == max(cv.oj$size)) 5 else best.size.oj

prune.oj <- prune.misclass(tree.oj, best = size.to.use)

plot(prune.oj)
text(prune.oj, 
     pretty = 0,
     cex = .8)

j.) Pruned v Unpruned Tree Error Rates

Compare the training error rates between the pruned an unpruned trees. Which is higher?

pred.unpruned.train <- predict(tree.oj,
                               newdata = oj.train,
                               type = "class")
train.error.unpruned <- mean(pred.unpruned.train != oj.train$Purchase)

pred.pruned.train <- predict(prune.oj,
                             newdata = oj.train,
                             type = "class")
train.error.pruned <- mean(pred.pruned.train != oj.train$Purchase)

cat("Unpruned:",train.error.unpruned*100, "%\n", "Pruned:", train.error.pruned*100, "%\n")

Unpruned: 15.875 %
 Pruned: 16.25 %

k.) Comparining the Test Error Rates

Compare the test error rates between the pruned and unpruned trees. Which is higher?

Unpruned: 15.875 %
 Pruned: 16.25 %

When compared the test error results are close, with the unpruned test error rate being lower at 15.88% and the pruned test error rate being higher at 16.25%. This suggests that pruning the tree did not improve the prediction on the test set. However since the difference in the test error is so small, both models perform well, with pruning not providing any substatial difference to the test data.