Assignment 7

ISLR Chapter 08 (page 361): 3, 8, 9

Q3. Gini index, classification error, and entropy

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 to make 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'))

Q8. Predicting quantitative Sales in the Carseats dataset using regression trees, bagging, random forests, and BART.

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 set into a training set and a test set.

library(ISLR)
library(caret)

data(Carseat)

After loading necessary libraries and our carseat dataset, now we can split our data into train and test.

set.seed(210)

train_index <- createDataPartition(Carseats$Sales, p=0.5, list = FALSE)

train_set<- Carseats[train_index,]
test_set <- Carseats[-train_index, ]

dim(train_set)
## [1] 201  11
dim(test_set)
## [1] 199  11

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

library(tree)
## Warning: package 'tree' was built under R version 4.5.3
tree_model <- tree(Sales ~ ., data = train_set)
summary(tree_model)
## 
## Regression tree:
## tree(formula = Sales ~ ., data = train_set)
## Variables actually used in tree construction:
## [1] "ShelveLoc"   "Price"       "CompPrice"   "Advertising" "Age"        
## [6] "Income"      "Education"  
## Number of terminal nodes:  17 
## Residual mean deviance:  2.049 = 377 / 184 
## Distribution of residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -4.56500 -0.91250 -0.01929  0.00000  1.06100  3.78800
plot(tree_model)
text(tree_model, pretty = 0, cex = 0.7)

> Our tree model reveals that shelf location is the single most critical factor driving sales; locations with “Good” shelving placement are split off immediately to the right and yield much higher baseline sales than “Bad” or “Medium” locations on the left. Within both branches, price acts as the secondary dominant driver, showing that lower prices consistently boost sales volumes, while secondary factors like advertising budgets, competitor pricing, customer age, income levels, and education refine the deeper branches of the model.

pred_tree <- predict(tree_model, newdata = test_set)
mse_tree  <- mean((pred_tree - test_set$Sales)^2)
cat("Test MSE:", mse_tree, "\n")
## Test MSE: 5.187061
  1. Use cross-validation in order to determine the optimal level of tree complexity. Does pruning the tree improve the test MSE?
cv_tree <- cv.tree(tree_model)

plot(cv_tree$size, cv_tree$dev, type = "b", xlab = "Tree Size", ylab = "Deviance")

best_size <- cv_tree$size[which.min(cv_tree$dev)]
cat("Optimal Tree Size determined by CV:", best_size, "\n")
## Optimal Tree Size determined by CV: 11

Our optimal Tree Size is determined y cross-validation method showing 17, next we will adjust our tree model with prunning.

pruned_tree <- prune.tree(tree_model, best = best_size)

plot(pruned_tree)
text(pruned_tree, pretty = 0, cex = 0.7)

Now we can predict using our prunned model and observe the new Test MSE.

pred_pruned <- predict(pruned_tree, newdata = test_set)
mse_pruned  <- mean((pred_pruned - test_set$Sales)^2)
cat("Pruned Tree Test MSE:", mse_pruned, "\n")
## Pruned Tree Test MSE: 5.220978

Prunning the previous tree model did not improve the Test MSE yet it showed identical results.

  1. 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(210)
bag_model <- randomForest(Sales ~ ., data = train_set, mtry = 10, importance = TRUE)
bag_model
## 
## Call:
##  randomForest(formula = Sales ~ ., data = train_set, mtry = 10,      importance = TRUE) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 10
## 
##           Mean of squared residuals: 2.722852
##                     % Var explained: 63.5
pred_bag <- predict(bag_model, newdata = test_set)
mse_bag  <- mean((pred_bag - test_set$Sales)^2)
cat("Bagging Test MSE:", mse_bag, "\n")
## Bagging Test MSE: 2.768769

Bagging method improved our model Test MSE by almost cutting in half.

importance(bag_model)
##               %IncMSE IncNodePurity
## CompPrice   26.447822    179.793092
## Income       6.082824     87.005853
## Advertising 19.814043    145.077596
## Population  -1.500549     49.700263
## Price       55.362102    441.393933
## ShelveLoc   45.731828    339.744246
## Age         17.387162    150.602827
## Education    3.099620     44.343231
## Urban       -2.078941      7.867465
## US           1.003939      9.021848
varImpPlot(bag_model)

> Based on our model output, Price and ShelveLoc are the two most critical predictors of car seat sales by a wide margin.

  1. 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(210)
rf_model <- randomForest(Sales ~ ., data = train_set, mtry = 3, importance = TRUE)
rf_model
## 
## Call:
##  randomForest(formula = Sales ~ ., data = train_set, mtry = 3,      importance = TRUE) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 3
## 
##           Mean of squared residuals: 3.193399
##                     % Var explained: 57.19

Overall, this model explains 57.19% of the variation in car seat sales and maintains an average squared error of 3.19 on the test data.

pred_rf <- predict(rf_model, newdata = test_set)
mse_rf  <- mean((pred_rf - test_set$Sales)^2)
cat("Random Forest Test MSE:", mse_rf, "\n")
## Random Forest Test MSE: 3.308276

Random forest model performed better than previous models, however test MSE is higher than bagging method.

importance(rf_model)
##                 %IncMSE IncNodePurity
## CompPrice    9.34514938     140.13391
## Income       6.09792629     122.70316
## Advertising 14.98808150     140.68426
## Population  -1.99147676      99.52757
## Price       37.19383666     350.85843
## ShelveLoc   32.18620982     268.28621
## Age         13.21162501     176.66196
## Education    0.08429424      70.67185
## Urban       -2.18030644      13.40461
## US           1.43363795      18.83175
varImpPlot(rf_model)

> Price and ShelveLoc remain the top variables to drive this model, but their dominance is slightly reduced as other variables are included.

  1. Now analyze the data using BART, and report your results.
library(BART)
## Warning: package 'BART' was built under R version 4.5.3
## Loading required package: nlme
## Loading required package: survival
## 
## Attaching package: 'survival'
## The following object is masked from 'package:caret':
## 
##     cluster
X_train <- model.matrix(Sales ~ . - 1, data = train_set)
y_train <- train_set$Sales

X_test  <- model.matrix(Sales ~ . - 1, data = test_set)
set.seed(210)
bart_model <- gbart(X_train, y_train, x.test = X_test)
## *****Calling gbart: type=1
## *****Data:
## data:n,p,np: 201, 12, 199
## y1,yn: 2.013284, 2.223284
## x1,x[n*p]: 138.000000, 1.000000
## xp1,xp[np*p]: 111.000000, 1.000000
## *****Number of Trees: 200
## *****Number of Cut Points: 62 ... 1
## *****burn,nd,thin: 100,1000,1
## *****Prior:beta,alpha,tau,nu,lambda,offset: 2,0.95,0.2512,3,0.214239,7.48672
## *****sigma: 1.048731
## *****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: 4s
## trcnt,tecnt: 1000,1000
pred_bart <- bart_model$yhat.test.mean
mse_bart  <- mean((pred_bart - test_set$Sales)^2)
cat("BART Test MSE:", mse_bart, "\n")
## BART Test MSE: 1.624166

Among all the models tested, the BART method performed the best, achieving the lowest test error with a Test MSE of 1.62. In comparison, the single decision tree performed the worst MSE (5.19), while bagging method showed Test MSE (2.77) and random forests (3.31).

Q9. Building and Tuning Classification Trees to Predict Brand Purchases in the OJ Dataset

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

library(ISLR2)
attach(OJ)
str(OJ)
## 'data.frame':    1070 obs. of  18 variables:
##  $ Purchase      : Factor w/ 2 levels "CH","MM": 1 1 1 2 1 1 1 1 1 1 ...
##  $ WeekofPurchase: num  237 239 245 227 228 230 232 234 235 238 ...
##  $ StoreID       : num  1 1 1 1 7 7 7 7 7 7 ...
##  $ PriceCH       : num  1.75 1.75 1.86 1.69 1.69 1.69 1.69 1.75 1.75 1.75 ...
##  $ PriceMM       : num  1.99 1.99 2.09 1.69 1.69 1.99 1.99 1.99 1.99 1.99 ...
##  $ DiscCH        : num  0 0 0.17 0 0 0 0 0 0 0 ...
##  $ DiscMM        : num  0 0.3 0 0 0 0 0.4 0.4 0.4 0.4 ...
##  $ SpecialCH     : num  0 0 0 0 0 0 1 1 0 0 ...
##  $ SpecialMM     : num  0 1 0 0 0 1 1 0 0 0 ...
##  $ LoyalCH       : num  0.5 0.6 0.68 0.4 0.957 ...
##  $ SalePriceMM   : num  1.99 1.69 2.09 1.69 1.69 1.99 1.59 1.59 1.59 1.59 ...
##  $ SalePriceCH   : num  1.75 1.75 1.69 1.69 1.69 1.69 1.69 1.75 1.75 1.75 ...
##  $ PriceDiff     : num  0.24 -0.06 0.4 0 0 0.3 -0.1 -0.16 -0.16 -0.16 ...
##  $ Store7        : Factor w/ 2 levels "No","Yes": 1 1 1 1 2 2 2 2 2 2 ...
##  $ PctDiscMM     : num  0 0.151 0 0 0 ...
##  $ PctDiscCH     : num  0 0 0.0914 0 0 ...
##  $ ListPriceDiff : num  0.24 0.24 0.23 0 0 0.3 0.3 0.24 0.24 0.24 ...
##  $ STORE         : num  1 1 1 1 0 0 0 0 0 0 ...

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

set.seed(210)

train_index_oj <- sample(1:nrow(OJ), 800)

train_set_oj<- OJ[train_index_oj,]
test_set_oj <- OJ[-train_index_oj, ]

dim(train_set_oj)
## [1] 800  18
dim(test_set_oj)
## [1] 270  18

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

oj_tree <- tree(Purchase ~ ., data = train_set_oj)
summary(oj_tree)
## 
## Classification tree:
## tree(formula = Purchase ~ ., data = train_set_oj)
## Variables actually used in tree construction:
## [1] "LoyalCH"   "PriceDiff"
## Number of terminal nodes:  8 
## Residual mean deviance:  0.7314 = 579.3 / 792 
## Misclassification error rate: 0.1638 = 131 / 800

Our classification tree model creates exactly 8 terminal nodes using only two predictive variables: customer brand loyalty to Citrus Hill LoyalCH and the price difference between the brands PriceDiff. This structural simplicity results in a training misclassification error rate of 16.38%, meaning the model correctly identifies the purchase choice for 669 out of the 800 training observations while missing 131.

(c) 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.

oj_tree
## node), split, n, deviance, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 800 1073.00 CH ( 0.60625 0.39375 )  
##    2) LoyalCH < 0.5036 353  413.10 MM ( 0.27195 0.72805 )  
##      4) LoyalCH < 0.035047 53    0.00 MM ( 0.00000 1.00000 ) *
##      5) LoyalCH > 0.035047 300  376.10 MM ( 0.32000 0.68000 )  
##       10) PriceDiff < 0.05 125  116.30 MM ( 0.17600 0.82400 ) *
##       11) PriceDiff > 0.05 175  238.40 MM ( 0.42286 0.57714 )  
##         22) LoyalCH < 0.260429 64   64.60 MM ( 0.20312 0.79688 ) *
##         23) LoyalCH > 0.260429 111  152.80 CH ( 0.54955 0.45045 ) *
##    3) LoyalCH > 0.5036 447  345.00 CH ( 0.87025 0.12975 )  
##      6) LoyalCH < 0.764572 180  208.80 CH ( 0.73333 0.26667 )  
##       12) PriceDiff < 0.265 105  142.80 CH ( 0.58095 0.41905 )  
##         24) PriceDiff < -0.165 28   33.50 MM ( 0.28571 0.71429 ) *
##         25) PriceDiff > -0.165 77   95.55 CH ( 0.68831 0.31169 ) *
##       13) PriceDiff > 0.265 75   31.23 CH ( 0.94667 0.05333 ) *
##      7) LoyalCH > 0.764572 267   85.31 CH ( 0.96255 0.03745 ) *

Terminal Node Interpretation (Node 4): This node represents 53 customers with extremely low brand loyalty to Citrus Hill (score below 0.035). Because this group is completely uniform, the error rate here is zero. The model predicts a Minute Maid (MM) purchase with 100% certainty, as nobody in this group bought Citrus Hill.

(d) Create a plot of the tree, and interpret the results.

plot(oj_tree)
text(oj_tree, pretty = 0, cex = 0.7)

> Our classification tree shows that customer loyalty LoyalCH is the most important factor in deciding which orange juice to buy, forming the first split at the top. Customers with high brand loyalty track to the right and are generally predicted to buy Citrus Hill (CH), while those with lower loyalty track to the left and mostly choose Minute Maid (MM). Further down the tree, smaller splits use price differences PriceDiff and tighter loyalty ranges to fine-tune the final customer predictions.

(e) 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?

pred_oj <- predict(oj_tree, newdata = test_set_oj, type = "class")

conf_matrix <- table(Predicted = pred_oj, Actual = test_set_oj$Purchase)
print(conf_matrix)
##          Actual
## Predicted  CH  MM
##        CH 151  30
##        MM  17  72
test_error <- mean(pred_oj != test_set_oj$Purchase)
cat("Test Error Rate:", test_error, "\n")
## Test Error Rate: 0.1740741

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

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

(g) Produce a plot with tree size on the x-axis and cross-validated classification error rate on the y-axis.

plot(cv_oj$size, cv_oj$dev, type = "b", xlab = "Tree Size", ylab = "CV Classification Error Rate")

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

best_oj_size <- cv_oj$size[which.min(cv_oj$dev)]
cat("Optimal tree size with lowest CV error rate:", best_oj_size, "\n")
## Optimal tree size with lowest CV error rate: 8

(i) 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.

pruned_oj <- prune.misclass(oj_tree, best = 8)

plot(pruned_oj)
text(pruned_oj, pretty = 0, cex = 0.7)

(j) Compare the training error rates between the pruned and unpruned trees. Which is higher?

summary(oj_tree)$misclass
## [1] 131 800
summary(pruned_oj)$misclass
## [1] 131 800

The training error rates for the unpruned and pruned trees are identical at 16.38% (131/800). Neither is higher because the cross-validation process determined that the original tree size of 8 terminal nodes was already the optimal size, meaning no branches were removed during pruning.

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

pred_pruned_oj <- predict(pruned_oj, newdata = test_set_oj, type = "class")
pruned_test_error <- mean(pred_pruned_oj != test_set_oj$Purchase)

cat("Unpruned Test Error Rate:", test_error, "\n")
## Unpruned Test Error Rate: 0.1740741
cat("Pruned Test Error Rate:", pruned_test_error, "\n")
## Pruned Test Error Rate: 0.1740741

Again, our test error rates are identical at 17.41%. Because the pruned tree is identical to the unpruned tree, their predictions on the test dataset are exactly the same, resulting in the same misclassification rate.