Assignment7

Problem 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.

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

• All three measures reach their maximum when p̂m1 = 0.5 (a perfectly impure node with equal class proportions) and their minimum at p̂m1 = 0 or 1 (a pure node). Entropy takes the largest values, followed by the Gini index, with classification error the smallest. The Gini index and entropy are smooth, differentiable curves while classification error is piecewise linear — which is one reason Gini and entropy are preferred for growing trees, as they are more sensitive to changes in node purity.

---

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

library(ISLR2)
library(tree)

set.seed(1)
n <- nrow(Carseats)
train_idx <- sample(1:n, n / 2)
carseats_train <- Carseats[train_idx, ]
carseats_test  <- Carseats[-train_idx, ]

(b) 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" "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 = 0.7)

tree_pred <- predict(tree_carseats, carseats_test)
mean((tree_pred - carseats_test$Sales)^2)

[1] 4.922039

• The tree splits first on ShelveLoc, it is the most important factor for sales. Price appears repeatedly throughout the tree, with lower prices leading to higher predicted sales. The test MSE is approximately 4.9, meaning predictions are off by about $2,200 on average.

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

library(rpart)
library(caret)
library(rattle)

#cost complexity = decrease the overall lack of fit by a factor .1 
tree.carseats <- rpart(Sales~.,data = carseats_train, method = "anova", control = rpart.control(minsplit = 15, cp=0.1)) 
summary(tree.carseats)

Call:
rpart(formula = Sales ~ ., data = carseats_train, method = "anova", 
    control = rpart.control(minsplit = 15, cp = 0.1))
  n= 200 

         CP nsplit rel error    xerror       xstd
1 0.2146655      0 1.0000000 1.0181022 0.09266075
2 0.1016064      1 0.7853345 0.8557999 0.07675451
3 0.1000000      2 0.6837281 0.8629270 0.07825594

Variable importance
ShelveLoc     Price CompPrice 
       63        31         6 

Node number 1: 200 observations,    complexity param=0.2146655
  mean=7.57815, MSE=7.863433 
  left son=2 (158 obs) right son=3 (42 obs)
  Primary splits:
      ShelveLoc   splits as  LRL,       improve=0.21466550, (0 missing)
      Price       < 89.5  to the right, improve=0.17140880, (0 missing)
      Age         < 61.5  to the right, improve=0.06863684, (0 missing)
      Advertising < 2.5   to the left,  improve=0.06010481, (0 missing)
      US          splits as  LR,        improve=0.04966938, (0 missing)
  Surrogate splits:
      Price < 163.5 to the left,  agree=0.795, adj=0.024, (0 split)

Node number 2: 158 observations,    complexity param=0.1016064
  mean=6.908291, MSE=6.105264 
  left son=4 (134 obs) right son=5 (24 obs)
  Primary splits:
      Price       < 94.5  to the right, improve=0.16565390, (0 missing)
      Age         < 50.5  to the right, improve=0.10951490, (0 missing)
      ShelveLoc   splits as  L-R,       improve=0.08044967, (0 missing)
      Advertising < 9.5   to the left,  improve=0.07610689, (0 missing)
      Income      < 61.5  to the left,  improve=0.06165666, (0 missing)
  Surrogate splits:
      CompPrice < 98.5  to the right, agree=0.88, adj=0.208, (0 split)

Node number 3: 42 observations
  mean=10.0981, MSE=6.439368 

Node number 4: 134 observations
  mean=6.482687, MSE=5.178796 

Node number 5: 24 observations
  mean=9.284583, MSE=4.619916 

tree.carseats

n= 200 

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

1) root 200 1572.6870  7.578150  
  2) ShelveLoc=Bad,Medium 158  964.6316  6.908291  
    4) Price>=94.5 134  693.9586  6.482687 *
    5) Price< 94.5 24  110.8780  9.284583 *
  3) ShelveLoc=Good 42  270.4534 10.098100 *

printcp(tree.carseats)

Regression tree:
rpart(formula = Sales ~ ., data = carseats_train, method = "anova", 
    control = rpart.control(minsplit = 15, cp = 0.1))

Variables actually used in tree construction:
[1] Price     ShelveLoc

Root node error: 1572.7/200 = 7.8634

n= 200 

       CP nsplit rel error  xerror     xstd
1 0.21467      0   1.00000 1.01810 0.092661
2 0.10161      1   0.78533 0.85580 0.076755
3 0.10000      2   0.68373 0.86293 0.078256

plotcp(tree.carseats)

fancyRpartPlot(tree.carseats)

#function to get 
#taking cptable and asking which has the min cross validation error
tree.carseats$cptable[which.min(tree.carseats$cptable[,"xerror"]),"CP"]

[1] 0.1016064

#prun
carseats.prun <- prune(tree.carseats,cp=tree.carseats$cptable[which.min(tree.carseats$cptable[,"xerror"]),"CP"])
fancyRpartPlot(carseats.prun)

# Original tree predictions
tree_pred <- predict(tree.carseats, newdata = carseats_test)
tree_mse <- mean((tree_pred - carseats_test$Sales)^2)

# Pruned tree predictions
pruned_pred <- predict(carseats.prun, newdata = carseats_test)
pruned_mse <- mean((pruned_pred - carseats_test$Sales)^2)

tree_mse

[1] 5.138384

pruned_mse

[1] 5.788948

• Cross-validation selects a tree close to the full size, so pruning does not improve the test MSE, it stays around the same value or slightly worse. Pruning mainly buys us a simpler, more interpretable tree rather than better predictions on this data set.

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

bag_mse

[1] 2.605253

• Bagging substantially improves the test MSE, cutting it down to roughly half of the single-tree error 2.6. The importance measures confirm that Price and ShelveLoc are by far the most important variables, followed by CompPrice and Age.

(e) 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)
rf_carseats <- randomForest(Sales ~ ., data = carseats_train, mtry = 3, importance = TRUE)
rf_pred <- predict(rf_carseats, carseats_test)
rf_mse  <- mean((rf_pred - carseats_test$Sales)^2)
rf_mse

[1] 2.960559

randomForest::importance(rf_carseats)

               %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

# Effect of m on test error
library(caret)
set.seed(1)

rf_grid <- expand.grid(mtry = 1:10)
rf_fit <- train(Sales ~ .,data = carseats_train,method = "rf",tuneGrid = rf_grid,
                trControl = trainControl(method = "cv", number = 10))

rf_fit

Random Forest 

200 samples
 10 predictor

No pre-processing
Resampling: Cross-Validated (10 fold) 
Summary of sample sizes: 180, 180, 180, 179, 180, 180, ... 
Resampling results across tuning parameters:

  mtry  RMSE      Rsquared   MAE     
   1    2.290542  0.5368420  1.892590
   2    2.011591  0.6187520  1.660789
   3    1.880137  0.6481762  1.540409
   4    1.812118  0.6567222  1.484342
   5    1.760494  0.6640174  1.442265
   6    1.729397  0.6702447  1.418660
   7    1.707664  0.6754862  1.402371
   8    1.697450  0.6719438  1.391994
   9    1.702145  0.6634968  1.396636
  10    1.692196  0.6684809  1.386273

RMSE was used to select the optimal model using the smallest value.
The final value used for the model was mtry = 10.

plot(rf_fit)

• The random forest again identifies Price and ShelveLoc as dominant. The plot of test MSE against m shows that very small values of m (like 1 or 2) hurt performance because the trees are too constrained, while performance improves and levels off as m grows. For this data set the best results come from relatively large m because there are only a couple of truly strong predictors, so decorrelating the trees by restricting m provides little benefit.

(f) Now analyze the data using BART, and report your results.

library(BART)

x_train <- carseats_train[, -1]
y_train <- carseats_train$Sales
x_test  <- carseats_test[, -1]

set.seed(1)
bart_fit  <- gbart(x_train, y_train, x.test = x_test)

*****Calling gbart: type=1
*****Data:
data:n,p,np: 200, 14, 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,14,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

bart_pred <- bart_fit$yhat.test.mean
bart_mse  <- mean((bart_pred - carseats_test$Sales)^2)
bart_mse

[1] 1.450842

• BART performs very well on this data, achieving a test MSE competitive with bagging and random forests. BART’s approach with built-in regularization handles this data set nicely with the default settings.

---

Problem 9

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

(a) Create a training set containing a random sample of 800 observations, and a test set 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) 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" "PctDiscMM"    
Number of terminal nodes:  9 
Residual mean deviance:  0.7432 = 587.8 / 791 
Misclassification error rate: 0.1588 = 127 / 800 

• The tree uses only a small handful of the available variables. The training error rate is approximately 15.9%, and the tree has around 9 terminal nodes. It’s notable that out of 17 predictors, loyalty and price capture nearly all of the signal.

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

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 ) *

• Consider one of the terminal nodes reached when LoyalCH is very low. For customers with extremely low Citrus Hill loyalty, the tree predicts Minute Maid with very high, nearly all training observations in that node purchased MM.

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

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

• The plot shows LoyalCH at the top splits of the tree, the first two levels split on loyalty alone. Customers with high CH loyalty almost always purchase CH, and those with low loyalty purchase MM. Price variables looks like it only matter for customers in the middle loyalty range.

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

tree_oj_pred <- predict(tree_oj, oj_test, type = "class")
table(tree_oj_pred, oj_test$Purchase)

            
tree_oj_pred  CH  MM
          CH 160  38
          MM   8  64

test_err <- mean(tree_oj_pred != oj_test$Purchase)
test_err

[1] 0.1703704

• The test error rate is approximately 17%, only slightly higher than the training error. The confusion matrix shows the tree does a somewhat better job classifying CH purchases than MM purchases.

(f) 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) 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 Errors")

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

best_oj_size <- cv_oj$size[which.min(cv_oj$dev)]
best_oj_size

[1] 9

• The lowest cross-validated error typically occurs around a tree size of 7-9 nodes. If several sizes tie, the smallest tree among them is preferred for simplicity.

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

prune_size <- ifelse(best_oj_size < summary(tree_oj)$size, best_oj_size, 5)
pruned_oj  <- prune.misclass(tree_oj, best = prune_size)
plot(pruned_oj)
text(pruned_oj, pretty = 0, cex = 0.8)

summary(pruned_oj)

Classification tree:
snip.tree(tree = tree_oj, nodes = c(4L, 10L))
Variables actually used in tree construction:
[1] "LoyalCH"       "PriceDiff"     "ListPriceDiff" "PctDiscMM"    
Number of terminal nodes:  7 
Residual mean deviance:  0.7748 = 614.4 / 793 
Misclassification error rate: 0.1625 = 130 / 800 

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

train_pred_full   <- predict(tree_oj, oj_train, type = "class")
train_pred_pruned <- predict(pruned_oj, oj_train, type = "class")

mean(train_pred_full != oj_train$Purchase)

[1] 0.15875

mean(train_pred_pruned != oj_train$Purchase)

[1] 0.1625

• The pruned tree has a slightly higher training error rate. This is expected, the unpruned tree has more terminal nodes and can fit the training data more closely. Removing nodes necessarily makes the fit to the training data the same or worse.

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

test_pred_pruned <- predict(pruned_oj, oj_test, type = "class")
mean(tree_oj_pred != oj_test$Purchase)

[1] 0.1703704

mean(test_pred_pruned != oj_test$Purchase)

[1] 0.162963

• The test error rates are very similar, with the pruned tree performing about the same as the unpruned tree. The takeaway is that pruning simplified the tree considerably without hurting test performance, a smaller tree with equivalent accuracy is the better model since it’s easier to interpret and less likely to overfit.