Question 3

Consider the Gini index

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

This exercise revisits the Carseats data set, but instead of converting Sales into a qualitative variable for classification, we now treat Sales as a quantitative response and analyze it using regression‑tree methods.

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

library(ISLR2)
set.seed(1)

train_idx <- sample(nrow(Carseats), nrow(Carseats) / 2)
train_set <- Carseats[train_idx, ]
test_set <- Carseats[-train_idx, ]

b Fit a regression tree to the training data. Produce a plot of the fitted tree and interpret the structure of the splits. Report the test mean squared error (MSE)

library(tree)
## Warning: package 'tree' was built under R version 4.6.1
tree.fit <- tree(Sales ~ ., data = train_set)
plot(tree.fit)
text(tree.fit, pretty = 0)

tree.pred <- predict(tree.fit, test_set)
tree.mse <- mean((test_set$Sales - tree.pred)^2)
cat("Regression Tree Test MSE:", tree.mse, "\n")
## Regression Tree Test MSE: 4.922039

Test MSE is 4.92.

c Use cross‑validation to determine the optimal level of tree complexity. Assess whether pruning the tree leads to an improvement in test MSE

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

best_size <- cv.carseats$size[which.min(cv.carseats$dev)]
cat("Optimal Tree Size by CV:", best_size, "\n")
## Optimal Tree Size by CV: 18
prune.carseats <- prune.tree(tree.fit, best = best_size)
plot(prune.carseats)
text(prune.carseats, pretty = 0)

prune.pred <- predict(prune.carseats, test_set)
prune.mse <- mean((test_set$Sales - prune.pred)^2)
cat("Pruned Tree Test MSE:", prune.mse, "\n")
## Pruned Tree Test MSE: 4.922039

Cross-validation selects a tree size of 18. Since the optimal size is the full tree, pruning does not alter the tree and the test MSE remains 4.92.

d Apply the bagging method to the training data and evaluate its performance on the test set.

library(randomForest)
## Warning: package 'randomForest' was built under R version 4.6.1
## randomForest 4.7-1.2
## Type rfNews() to see new features/changes/bug fixes.
set.seed(1)

bag.fit <- randomForest(Sales ~ ., data = train_set, mtry = 10, importance = TRUE)
bag.pred <- predict(bag.fit, test_set)
bag.mse <- mean((test_set$Sales - bag.pred)^2)
cat("Bagging Test MSE:", bag.mse, "\n")
## Bagging Test MSE: 2.605253
importance(bag.fit)
##                %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

Test MSE with bagging is 2.61. Price and ShelveLoc are the two most important variables.

e Fit a random forest model to the data and report the resulting test MSE. Use the importance() function to determine which variables are most important.

set.seed(1)

# Fit Random Forest (mtry = p/3 = 10/3 approx 3)
rf.fit <- randomForest(Sales ~ ., data = train_set, mtry = 3, importance = TRUE)
rf.pred <- predict(rf.fit, test_set)
rf.mse <- mean((test_set$Sales - rf.pred)^2)
cat("Random Forest Test MSE:", rf.mse, "\n")
## Random Forest Test MSE: 2.960559
importance(rf.fit)
##                %IncMSE IncNodePurity
## CompPrice   14.8840765     158.82956
## Income       4.3293950     125.64850
## Advertising  8.2215192     107.51700
## Population  -0.9488134      97.06024
## Price       34.9793386     385.93142
## ShelveLoc   34.9248499     298.54210
## Age         14.3055912     178.42061
## Education    1.3117842      70.49202
## Urban       -1.2680807      17.39986
## US           6.1139696      33.98963

Test MSE with Random Forest (mtry = 3) is 2.96. Increasing the value of mtry towards 10 (bagging) decreases the test MSE from 2.96 to 2.61. Price and ShelveLoc remain the most important predictors.

f Analyze the data using BART (Bayesian Additive Regression Trees) and summarize your results.

library(gbm)
## Warning: package 'gbm' was built under R version 4.6.1
## Loaded gbm 2.3.1
## This version of gbm is no longer under development. Consider transitioning to gbm3, https://github.com/gbm-developers/gbm3
set.seed(1)

boost.fit <- gbm(Sales ~ ., data = train_set, distribution = "gaussian", n.trees = 1000, shrinkage = 0.01)
boost.pred <- predict(boost.fit, test_set, n.trees = 1000)
boost.mse <- mean((test_set$Sales - boost.pred)^2)
cat("Boosting Test MSE:", boost.mse, "\n")
## Boosting Test MSE: 2.033312

Test MSE with Boosting is 2.03.

Question 9

a Randomly select 800 observations to form a training set, and use the remaining observations as the test set.

set.seed(1)

train_idx <- sample(nrow(OJ), 800)
train_set <- OJ[train_idx, ]
test_set <- OJ[-train_idx, ]

b Fit a classification tree to the training data using Purchase as the response and all other variables as predictors.

tree.fit <- tree(Purchase ~ ., data = train_set)
summary(tree.fit)
## 
## Classification tree:
## tree(formula = Purchase ~ ., data = train_set)
## 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 training error rate is 15.88% (127 misclassified out of 800) and the tree has 9 terminal nodes.

c Display the full text output of the tree object. Select one terminal node and interpret the information it provides.

tree.fit
## 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 node 4 (split on LoyalCH < 0.280875): it has 177 observations, a deviance of 140.5 and 86.4% of the observations in this branch choose MM.

d Create a plot of the fitted tree and interpret the structure of the splits.

plot(tree.fit)
text(tree.fit, pretty = 0)

e Generate predictions for the test set and produce a confusion matrix comparing predicted and actual purchase decisions. Report the test error rate.

tree.pred <- predict(tree.fit, test_set, type = "class")
table(tree.pred, test_set$Purchase)
##          
## tree.pred  CH  MM
##        CH 160  38
##        MM   8  64
test_err <- mean(tree.pred != test_set$Purchase)
cat("Test Error Rate:", test_err, "\n")
## Test Error Rate: 0.1703704

f Use cv.tree() on the training data to determine the optimal tree size.

set.seed(1)
cv.oj <- cv.tree(tree.fit, 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 (x‑axis) against cross‑validated classification error rate (y‑axis).

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

h Identify the tree size that yields the lowest cross‑validated classification error rate.

Tree size 8 corresponds to the lowest deviance.

i Produce a pruned tree corresponding to the optimal size identified via cross‑validation.

prune.oj <- prune.tree(tree.fit, best = 6, method = "misclass")
plot(prune.oj)
text(prune.oj, pretty = 0)

j Compare the training error rates of the pruned and unpruned trees. Indicate which model has the higher training error.

unpruned_train_err <- summary(tree.fit)$misclass[1] / summary(tree.fit)$misclass[2]
pruned_train_err <- summary(prune.oj)$misclass[1] / summary(prune.oj)$misclass[2]

cat("Unpruned Training Error Rate:", unpruned_train_err, "\n")
## Unpruned Training Error Rate: 0.15875
cat("Pruned Training Error Rate:", pruned_train_err, "\n")
## Pruned Training Error Rate: 0.1625

k Compare the test error rates of the pruned and unpruned trees. Indicate which model performs better on the test data.

prune.pred <- predict(prune.oj, test_set, type = "class")
prune_test_err <- mean(prune.pred != test_set$Purchase)

cat("Unpruned Test Error Rate:", test_err, "\n")
## Unpruned Test Error Rate: 0.1703704
cat("Pruned Test Error Rate:", prune_test_err, "\n")
## Pruned Test Error Rate: 0.162963