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.

Gini: 2p1(1-p1)

Classification error: 1-max(p1, 1-p1)

Entropy: -p1*log2(p1)-(1-p1)log2(1-p1)

p_m1 <- seq(0, 1, length.out = 100)
gini <- 2 * p_m1 * (1 - p_m1)
classification_error <- 1 - pmax(p_m1, 1 - p_m1)
entropy <- -p_m1 * ifelse(p_m1 == 0, 0, log2(p_m1)) - (1 - p_m1) * ifelse((1 - p_m1) == 0, 0, log2(1 - p_m1))

data <- data.frame(p_m1, gini, classification_error, entropy)

ggplot(data, aes(x = p_m1)) +
  geom_line(aes(y = gini, color = "Gini Index")) +
  geom_line(aes(y = classification_error, color = "Classification Error")) +
  geom_line(aes(y = entropy, color = "Entropy")) +
  labs(x = expression(hat(p)[m1]), y = "Metric Value", title = "Comparison of Gini Index, Classification Error, and Entropy") +
  scale_color_manual(name = "Metrics", values = c("Gini Index" = "blue", "Classification Error" = "orange", "Entropy" = "green")) +
  theme_minimal() +
  theme(legend.title = element_text(size = 12), legend.text = element_text(size = 10))

ggsave("comparison_plot.png", width = 10, height = 6, dpi = 300)

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.

  1. Split the data set into a training set and a test set.
data("Carseats")

head(Carseats)
##   Sales CompPrice Income Advertising Population Price ShelveLoc Age Education
## 1  9.50       138     73          11        276   120       Bad  42        17
## 2 11.22       111     48          16        260    83      Good  65        10
## 3 10.06       113     35          10        269    80    Medium  59        12
## 4  7.40       117    100           4        466    97    Medium  55        14
## 5  4.15       141     64           3        340   128       Bad  38        13
## 6 10.81       124    113          13        501    72       Bad  78        16
##   Urban  US
## 1   Yes Yes
## 2   Yes Yes
## 3   Yes Yes
## 4   Yes Yes
## 5   Yes  No
## 6    No Yes
set.seed(42)

trainIndex <- createDataPartition(Carseats$Sales, p = 0.8, list = FALSE)
trainSet <- Carseats[trainIndex, ]
testSet <- Carseats[-trainIndex, ]
  1. Fit a regression tree to the training set. Plot the tree, and interpret the results. What test MSE do you obtain?
reg_tree <- tree(Sales ~ ., data = trainSet)

plot(reg_tree)
text(reg_tree, pretty = 0)

tree_preds <- predict(reg_tree, newdata = testSet)
tree_mse <- mean((tree_preds - testSet$Sales)^2)

cat("Test MSE:", tree_mse, "\n")
## Test MSE: 4.419773
  1. Use cross-validation in order to determine the optimal level of tree complexity. Does pruning the tree improve the test MSE?
reg_tree <- tree(Sales ~ ., data = trainSet)


cv_reg_tree <- cv.tree(reg_tree, FUN = prune.tree)
plot(cv_reg_tree$size, cv_reg_tree$dev, type = "b", xlab = "Number of Terminal Nodes", ylab = "Cross-Validation Error")

opt_size <- cv_reg_tree$size[which.min(cv_reg_tree$dev)]
pruned_tree <- prune.tree(reg_tree, best = opt_size)
pruned_preds <- predict(pruned_tree, newdata = testSet)
pruned_mse <- mean((pruned_preds - testSet$Sales)^2)

cat("Test MSE of pruned tree:", pruned_mse, "\n")
## Test MSE of pruned tree: 4.419773

Pruned Tree’s test MSE from 9 terminal nodes is slightly lower compared to the full tree’s test MSE.

  1. Use the bagging approach in order to analyze this data. What test MSE do you obtain? Use the importance() function to de- termine which variables are most important.
set.seed(42)
bagging_model <- randomForest(Sales ~ ., data = trainSet, mtry = ncol(trainSet) - 1, importance = TRUE)
bagging_predictions <- predict(bagging_model, newdata = testSet)

bagging_mse <- mean((bagging_predictions - testSet$Sales)^2)
cat("Test MSE with Bagging:", bagging_mse, "\n")
## Test MSE with Bagging: 2.68682
variable_importance <- importance(bagging_model)
print(variable_importance)
##                %IncMSE IncNodePurity
## CompPrice   43.0233684     291.29959
## Income      11.7505104     145.39993
## Advertising 25.6479432     215.54631
## Population  -1.3936689      83.40399
## Price       72.0512155     730.35533
## ShelveLoc   79.5149049     783.79279
## Age         21.3819204     195.41739
## Education    2.5287423      62.13057
## Urban        0.4126764      11.83756
## US           1.2199693      10.89709

Test MSE (2.69) is a lot lower compared to pruned tree. Variables Price and ShelveLoc were the most important.

  1. Use random forests to analyze this data. What test MSE do you obtain? Use the importance() function to determine which vari- ables are most important. Describe the effect of m, the number of variables considered at each split, on the error rate obtained.
set.seed(42)
rf_model <- randomForest(Sales ~ ., data = trainSet, mtry = sqrt(ncol(trainSet) - 1), importance = TRUE)
rf_predictions <- predict(rf_model, newdata = testSet)
rf_mse <- mean((rf_predictions - testSet$Sales)^2)

cat("Test MSE with Random Forests:", rf_mse, "\n")
## Test MSE with Random Forests: 2.967599
variable_importance_rf <- importance(rf_model)
print(variable_importance_rf)
##               %IncMSE IncNodePurity
## CompPrice   19.259916     241.25181
## Income       5.858968     188.83700
## Advertising 18.073452     226.20344
## Population  -2.350612     148.94221
## Price       47.742745     583.02079
## ShelveLoc   55.249672     619.86429
## Age         14.475603     267.17355
## Education    1.158960     104.44468
## Urban       -2.093707      21.27270
## US           2.238425      33.67042
m_values <- 1:(ncol(trainSet) - 1)
error_rates <- numeric(length(m_values))

for (m in m_values) {
  temp_model <- randomForest(Sales ~ ., data = trainSet, mtry = m, ntree = 100)
  temp_preds <- predict(temp_model, newdata = testSet)
  error_rates[m] <- mean((temp_preds - testSet$Sales)^2)
}

plot(m_values, error_rates, type = 'b', col = 'blue', xlab = "Number of Variables (m)", ylab = "Test MSE", main = "Effect of m on Random Forest Error Rate")

Test MSE (2.97) was slightly higher than the bagging model. Variables Price and ShelveLoc were still considered important.

  1. Now analyze the data using BART, and report your results.
set.seed(42)
trainSet[] <- lapply(trainSet, function(x) if(is.character(x)) factor(x) else x)
testSet[] <- lapply(testSet, function(x) if(is.character(x)) factor(x) else x)

trainSet <- na.omit(trainSet)
testSet <- na.omit(testSet)

x.train <- model.matrix(~ . - Sales, data = trainSet)
y.train <- trainSet$Sales
x.test <- model.matrix(~ . - Sales, data = testSet)
y.test <- testSet$Sales

bart_fit <- dbarts::bart(x.train = x.train, y.train = y.train, x.test = x.test)
## 
## Running BART with numeric y
## 
## number of trees: 200
## number of chains: 1, number of threads 1
## tree thinning rate: 1
## Prior:
##  k prior fixed to 2.000000
##  degrees of freedom in sigma prior: 3.000000
##  quantile in sigma prior: 0.900000
##  scale in sigma prior: 0.000719
##  power and base for tree prior: 2.000000 0.950000
##  use quantiles for rule cut points: false
##  proposal probabilities: birth/death 0.50, swap 0.10, change 0.40; birth 0.50
## data:
##  number of training observations: 321
##  number of test observations: 79
##  number of explanatory variables: 12
##  init sigma: 0.988581, curr sigma: 0.988581
## 
## Cutoff rules c in x<=c vs x>c
## Number of cutoffs: (var: number of possible c):
## (1: 100) (2: 100) (3: 100) (4: 100) (5: 100) 
## (6: 100) (7: 100) (8: 100) (9: 100) (10: 100) 
## (11: 100) (12: 100) 
## Running mcmc loop:
## iteration: 100 (of 1000)
## iteration: 200 (of 1000)
## iteration: 300 (of 1000)
## iteration: 400 (of 1000)
## iteration: 500 (of 1000)
## iteration: 600 (of 1000)
## iteration: 700 (of 1000)
## iteration: 800 (of 1000)
## iteration: 900 (of 1000)
## iteration: 1000 (of 1000)
## total seconds in loop: 0.553183
## 
## Tree sizes, last iteration:
## [1] 2 2 3 2 2 2 2 2 2 4 2 2 2 2 3 2 2 3 
## 3 2 1 2 2 3 2 4 2 3 2 3 2 2 2 2 2 2 2 2 
## 2 2 3 3 3 3 3 2 2 2 2 2 3 2 2 2 2 2 2 3 
## 3 3 3 2 2 3 1 2 3 2 4 2 3 2 3 2 3 5 4 2 
## 3 3 2 2 3 3 2 3 2 4 2 3 2 1 4 3 2 2 2 4 
## 2 3 3 2 3 4 2 3 3 3 2 3 3 2 1 5 2 3 2 3 
## 2 4 2 3 3 2 5 2 4 2 3 2 2 3 2 1 3 2 3 3 
## 3 3 4 3 3 2 2 5 2 1 2 2 3 6 2 2 2 2 2 2 
## 1 2 2 3 3 2 2 2 2 2 3 3 3 2 2 2 2 1 3 2 
## 2 2 2 3 2 2 2 2 2 1 1 3 2 2 2 2 3 3 2 4 
## 1 3 
## 
## Variable Usage, last iteration (var:count):
## (1: 0) (2: 20) (3: 31) (4: 25) (5: 24) 
## (6: 40) (7: 25) (8: 26) (9: 20) (10: 27) 
## (11: 22) (12: 31) 
## DONE BART
predictions <- bart_fit$yhat.test.mean
mse <- mean((predictions - y.test)^2)

cat("Test MSE with BART using dbarts:", mse, "\n")
## Test MSE with BART using dbarts: 1.80434

Test MSE was 1.80 which was much lower than the bagging and the pruned tree models’ test MSEs.

QUESTION 9

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

  1. Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations.
data("OJ")
set.seed(42)

train_indices <- sample(1:nrow(OJ), 800)
train_set <- OJ[train_indices, ]
test_set <- OJ[-train_indices, ]
  1. 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?
set.seed(42)
tree_model <- tree(Purchase ~ ., data = train_set)

summary(tree_model)
## 
## Classification tree:
## tree(formula = Purchase ~ ., data = train_set)
## Variables actually used in tree construction:
## [1] "LoyalCH"     "SalePriceMM" "PriceDiff"  
## Number of terminal nodes:  8 
## Residual mean deviance:  0.7392 = 585.5 / 792 
## Misclassification error rate: 0.1638 = 131 / 800
train_pred <- predict(tree_model, train_set, type = "class")
train_error_rate <- mean(train_pred != train_set$Purchase)
cat("Training Error Rate:", train_error_rate, "\n")
## Training Error Rate: 0.16375
num_terminal_nodes <- length(tree_model$frame$var[tree_model$frame$var == "<leaf>"])
cat("Number of Terminal Nodes:", num_terminal_nodes, "\n")
## Number of Terminal Nodes: 8

Training error rate was 0.16375 with 8 terminal nodes. High residual mean deviance suggests poor model fit.

  1. 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_model
## node), split, n, deviance, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 800 1066.00 CH ( 0.61500 0.38500 )  
##    2) LoyalCH < 0.48285 285  296.00 MM ( 0.21404 0.78596 )  
##      4) LoyalCH < 0.064156 64    0.00 MM ( 0.00000 1.00000 ) *
##      5) LoyalCH > 0.064156 221  260.40 MM ( 0.27602 0.72398 )  
##       10) SalePriceMM < 2.04 128  123.50 MM ( 0.18750 0.81250 ) *
##       11) SalePriceMM > 2.04 93  125.00 MM ( 0.39785 0.60215 ) *
##    3) LoyalCH > 0.48285 515  458.10 CH ( 0.83689 0.16311 )  
##      6) LoyalCH < 0.753545 230  282.70 CH ( 0.69565 0.30435 )  
##       12) PriceDiff < 0.265 149  203.00 CH ( 0.57718 0.42282 )  
##         24) PriceDiff < -0.165 32   38.02 MM ( 0.28125 0.71875 ) *
##         25) PriceDiff > -0.165 117  150.30 CH ( 0.65812 0.34188 )  
##           50) LoyalCH < 0.703993 105  139.60 CH ( 0.61905 0.38095 ) *
##           51) LoyalCH > 0.703993 12    0.00 CH ( 1.00000 0.00000 ) *
##       13) PriceDiff > 0.265 81   47.66 CH ( 0.91358 0.08642 ) *
##      7) LoyalCH > 0.753545 285  111.70 CH ( 0.95088 0.04912 ) *

Node 10, 128 observations. Criterion condition is when price is below 2.04, customers likely to purchase MM. No further splits.

  1. Create a plot of the tree, and interpret the results.
plot(tree_model)
text(tree_model, pretty = 0)

5 variables used, LoyalCH (root node) is most important.

  1. 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?
test_pred <- predict(tree_model, newdata = test_set, type = "class")
confusion_matrix <- table(test_set$Purchase, test_pred)
print(confusion_matrix)
##     test_pred
##       CH  MM
##   CH 125  36
##   MM  15  94
test_error_rate <- 1 - sum(diag(confusion_matrix)) / sum(confusion_matrix)
cat("Test Error Rate:", test_error_rate, "\n")
## Test Error Rate: 0.1888889
  1. Apply the cv.tree() function to the training set in order to determine the optimal tree size.
set.seed(42)
cv_tree <- cv.tree(tree_model, FUN = prune.misclass)

optimal_size <- which.min(cv_tree$dev)
cat("Optimal size based on cross-validated classification error rate:", cv_tree$size[optimal_size], "\n")
## Optimal size based on cross-validated classification error rate: 5
pruned_tree <- prune.misclass(tree_model, best = optimal_size)

summary(pruned_tree)
## 
## Classification tree:
## snip.tree(tree = tree_model, nodes = 2:3)
## Variables actually used in tree construction:
## [1] "LoyalCH"
## Number of terminal nodes:  2 
## Residual mean deviance:  0.945 = 754.1 / 798 
## Misclassification error rate: 0.1812 = 145 / 800

Tree size 5 was the optimal size based on cv.tree().

  1. Produce a plot with tree size on the x-axis and cross-validated classification error rate on the y-axis.
plot(cv_tree$size, cv_tree$dev, type = "b",
     xlab = "Size of tree", ylab = "Cross-validated Classification Error Rate")

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

Tree size of 5 had the lowest classification error rate.

  1. 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.
set.seed(42)
pruned_tree <- prune.misclass(tree_model, best = 5)
  
plot(pruned_tree)
text(pruned_tree, pretty = 0)

summary(pruned_tree)
## 
## Classification tree:
## snip.tree(tree = tree_model, nodes = c(2L, 25L))
## Variables actually used in tree construction:
## [1] "LoyalCH"   "PriceDiff"
## Number of terminal nodes:  5 
## Residual mean deviance:  0.8096 = 643.6 / 795 
## Misclassification error rate: 0.1638 = 131 / 800
  1. Compare the training error rates between the pruned and unpruned trees. Which is higher?
summary(tree_model)
## 
## Classification tree:
## tree(formula = Purchase ~ ., data = train_set)
## Variables actually used in tree construction:
## [1] "LoyalCH"     "SalePriceMM" "PriceDiff"  
## Number of terminal nodes:  8 
## Residual mean deviance:  0.7392 = 585.5 / 792 
## Misclassification error rate: 0.1638 = 131 / 800
set.seed(42)
pruned_tree2 <- prune.misclass(tree_model, newdata = train_set, best = 5)
summary(pruned_tree2)
## 
## Classification tree:
## snip.tree(tree = tree_model, nodes = c(2L, 25L))
## Variables actually used in tree construction:
## [1] "LoyalCH"   "PriceDiff"
## Number of terminal nodes:  5 
## Residual mean deviance:  0.8096 = 643.6 / 795 
## Misclassification error rate: 0.1638 = 131 / 800

Both error rates on the training set were the same (0.1638).

  1. Compare the test error rates between the pruned and unpruned trees. Which is higher?
cat("Unpruned test Error Rate:", test_error_rate, "\n\n")
## Unpruned test Error Rate: 0.1888889
Pruned_Test_Set <- predict(pruned_tree, newdata = test_set, type = "class")
confusion_matrix2 <- table(test_set$Purchase, Pruned_Test_Set)
print(confusion_matrix2)
##     Pruned_Test_Set
##       CH  MM
##   CH 125  36
##   MM  15  94
test_error_rate2 <- 1 - sum(diag(confusion_matrix2)) / sum(confusion_matrix2)
cat("Pruned test Error Rate:", test_error_rate2)
## Pruned test Error Rate: 0.1888889

Both error rates on the test set were the same (0.1889), but pruned tree is more interpretable.