1 Exercise 3 — Gini index, classification error, and entropy

For a two-class problem, if \(\hat p_{m1}\) is the proportion of class 1 in node \(m\), then \(\hat p_{m2} = 1 - \hat p_{m1}\), and:

  • Gini index: \(G = 2\,\hat p_{m1}(1-\hat p_{m1})\)
  • Classification error: \(E = 1 - \max(\hat p_{m1}, 1-\hat p_{m1})\)
  • Entropy: \(D = -\hat p_{m1}\log(\hat p_{m1}) - \hat p_{m2}\log(\hat p_{m2})\)
p      <- seq(0.001, 0.999, by = 0.001)
gini   <- 2 * p * (1 - p)
class.error <- 1 - pmax(p, 1 - p)
entropy <- -(p * log(p) + (1 - p) * log(1 - p))

matplot(p, cbind(gini, class.error, entropy), type = "l", lty = 1, lwd = 2,
        col = c("red", "darkgreen", "blue"),
        xlab = expression(hat(p)[m1]), ylab = "Value",
        main = "Gini Index, Classification Error, and Entropy")
legend("top", legend = c("Gini index", "Classification error", "Entropy"),
       col = c("red", "darkgreen", "blue"), lty = 1, lwd = 2, bty = "n")

All three measures are symmetric around \(\hat p_{m1} = 0.5\) (where node impurity is maximal) and all equal zero at \(\hat p_{m1} = 0\) or \(1\) (a pure node). Entropy is scaled the largest, Gini index sits in between, and classification error is the smallest and least sensitive curve — which is why Gini index and entropy, not classification error, are preferred for growing trees (they are more sensitive to node purity).


2 Exercise 8 — Regression trees on the Carseats data set

data(Carseats)

2.1 (a) Train/test split

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

2.2 (b) Fit a regression tree

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"  
## [6] "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)

yhat <- predict(tree.carseats, newdata = Carseats.test)
tree.test.mse <- mean((yhat - Carseats.test$Sales)^2)
tree.test.mse
## [1] 4.922039

The tree typically splits first on ShelveLoc (shelving location) and Price, the two strongest predictors of Sales, followed by variables such as Age, CompPrice, and Advertising. The unpruned tree test MSE is usually in the range of ~4–5.

2.3 (c) Cross-validation to select tree complexity

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

best.size <- cv.carseats$size[which.min(cv.carseats$dev)]
best.size
## [1] 18
prune.carseats <- prune.tree(tree.carseats, best = best.size)
plot(prune.carseats)
text(prune.carseats, pretty = 0, cex = 0.8)

yhat.prune <- predict(prune.carseats, newdata = Carseats.test)
prune.test.mse <- mean((yhat.prune - Carseats.test$Sales)^2)
prune.test.mse
## [1] 4.922039

Whether pruning helps depends on the CV curve: if cv.tree() shows deviance still decreasing at the largest size, then the full (unpruned) tree is selected as “best” and pruning to a smaller size will not reduce test MSE here; often it makes the test MSE about the same or slightly worse, since the original tree in (b) was already fairly small. In general, pruning helps most when the original tree was allowed to grow very large/overfit.

2.4 (d) Bagging

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)
bag.test.mse
## [1] 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)

Bagging substantially reduces the test MSE relative to a single tree (often to roughly half the single-tree MSE, e.g., around 2–2.5). importance() typically identifies ShelveLoc and Price as, by far, the two most important variables (highest %IncMSE and IncNodePurity).

2.5 (e) Random forests

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, TestMSE = rf.test.mse)
##   mtry  TestMSE
## 1    2 3.401523
## 2    3 3.019287
## 3    4 2.821882
## 4    5 2.714157
## 5    6 2.638459
## 6    7 2.619630
## 7   10 2.608388
plot(mtry.values, rf.test.mse, type = "b", pch = 19,
     xlab = "mtry (number of variables tried at each split)",
     ylab = "Test MSE",
     main = "Random Forest Test MSE vs. mtry")

best.mtry <- mtry.values[which.min(rf.test.mse)]
rf.best <- randomForest(Sales ~ ., data = Carseats.train,
                         mtry = best.mtry, importance = TRUE)
importance(rf.best)
##                %IncMSE IncNodePurity
## CompPrice   25.8833322    168.843159
## Income       4.3970260     91.116446
## Advertising 12.0202979    103.004860
## Population  -1.6008612     58.293306
## Price       56.8359105    506.561117
## ShelveLoc   49.8878164    385.186420
## Age         16.2440038    156.095735
## Education    0.5719519     44.781080
## Urban        0.4550507      9.266563
## US           4.5323141     15.926806
varImpPlot(rf.best)

As m (mtry) increases toward p (the bagging case), test MSE generally decreases, since Carseats has a small number of genuinely strong predictors (ShelveLoc, Price) — decorrelating trees with a very small mtry (e.g. 2) tends to hurt performance here because many trees are forced to split on weaker variables. ShelveLoc and Price remain the top two variables in importance regardless of mtry.

2.6 (f) BART

x <- Carseats[, -which(names(Carseats) == "Sales")]
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: 1s
## trcnt,tecnt: 1000,1000
yhat.bart <- bartfit$yhat.test.mean
bart.test.mse <- mean((ytest - yhat.bart)^2)
bart.test.mse
## [1] 1.428145

BART typically achieves the lowest test MSE of all the methods tried in this exercise (often noticeably lower than bagging or random forests), consistent with the general pattern that BART tends to be highly competitive on moderate-sized regression problems like Carseats.

Summary of test MSEs:

data.frame(
  Method = c("Single Tree", "Pruned Tree", "Bagging", "Random Forest (best mtry)", "BART"),
  TestMSE = c(tree.test.mse, prune.test.mse, bag.test.mse, min(rf.test.mse), bart.test.mse)
)
##                      Method  TestMSE
## 1               Single Tree 4.922039
## 2               Pruned Tree 4.922039
## 3                   Bagging 2.605253
## 4 Random Forest (best mtry) 2.608388
## 5                      BART 1.428145

3 Exercise 9 — Classification trees on the OJ data set

data(OJ)

3.1 (a) Train/test split

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

3.2 (b) Fit a classification tree

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 tree typically uses only a handful of variables (often LoyalCH at the top, plus variables like PriceDiff, SpecialCH, or ListPriceDiff), with somewhere around 7–9 terminal nodes, and a training error rate typically in the range of 15–17%. LoyalCH (customer brand loyalty toward Citrus Hill) is almost always the dominant/first splitting variable, since it is highly predictive of purchase behavior.

3.3 (c) Detailed text output for a terminal node

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

Picking one terminal node line from the output (e.g., the first LoyalCH split ending in an asterisk * denoting a terminal node): the line reports (i) the node number, (ii) the splitting criterion that defines the region (e.g., LoyalCH < 0.076...), (iii) the number of training observations in that node, (iv) the deviance for the node, (v) the predicted class label for observations in that region, and (vi) the fraction of observations in the node belonging to each class (e.g., ( 0.96 0.04 )). A terminal node with, say, 60 observations and probabilities (0.02, 0.98) for (CH, MM) indicates a nearly pure region almost entirely composed of MM purchases.

3.4 (d) Plot the tree

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

LoyalCH clearly dominates the top of the tree: very low customer loyalty routes to a strong MM prediction, very high loyalty routes to a strong CH prediction, and the tree only needs additional variables like PriceDiff to distinguish among the “middle” loyalty customers.

3.5 (e) Predict on the test data / confusion matrix

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.rate <- 1 - sum(diag(conf.mat)) / sum(conf.mat)
test.error.rate
## [1] 0.1703704

The test error rate is usually close to the training error rate (roughly 15–20%), which suggests the unpruned tree is not badly overfitting on this data set.

3.6 (f) Cross-validation to determine 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"

3.7 (g) Plot tree size vs. cross-validated error rate

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

3.8 (h) Optimal tree size

best.size.oj <- cv.oj$size[which.min(cv.oj$dev)]
best.size.oj
## [1] 9

The tree size that minimizes CV error is typically somewhere in the range of 5–9 terminal nodes (frequently the CV curve is fairly flat across several sizes, so several sizes may tie for lowest error — which.min returns the first/smallest such size).

3.9 (i) Produce the pruned tree

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 = 0.8)

3.10 (j) Compare training error: pruned vs. unpruned

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)

c(Unpruned = train.error.unpruned, Pruned = train.error.pruned)
## Unpruned   Pruned 
##  0.15875  0.16250

The unpruned tree’s training error rate is always less than or equal to the pruned tree’s, since the unpruned tree is more flexible/complex and fits the training data more closely by construction.

3.11 (k) Compare test error: pruned vs. unpruned

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

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

c(Unpruned = test.error.unpruned, Pruned = test.error.pruned)
##  Unpruned    Pruned 
## 0.1703704 0.1629630

On the test set, the pruned tree’s error rate is typically very close to, or slightly better than, the unpruned tree’s — the extra terminal nodes in the unpruned tree mostly capture noise rather than genuine structure, so pruning to the CV-selected size tends to generalize about as well (or marginally better) while producing a simpler, more interpretable tree.