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

In a two-class setting \(\hat p_{m2} = 1 - \hat p_{m1}\), so writing \(p = \hat p_{m1}\):

  • Classification error: \(E = 1 - \max(p,\, 1-p)\)
  • Gini index: \(G = 2p(1-p)\)
  • Entropy: \(D = -\,p\log p - (1-p)\log(1-p)\)
p <- seq(0, 1, length = 500)

gini    <- 2 * p * (1 - p)
class.e <- 1 - pmax(p, 1 - p)
entropy <- -(p * log(p) + (1 - p) * log(1 - p))
entropy[is.nan(entropy)] <- 0          # limits at p = 0 and p = 1

plot(p, entropy, type = "l", col = "blue", lwd = 2, ylim = c(0, 0.75),
     xlab = expression(hat(p)[m1]), ylab = "Node impurity",
     main = "Gini index, classification error, and entropy")
lines(p, gini,    col = "red",       lwd = 2)
lines(p, class.e, col = "darkgreen", lwd = 2)
legend("topright", lwd = 2, col = c("red", "darkgreen", "blue"),
       legend = c("Gini index", "Classification error", "Entropy"))

All three measures are zero at \(p = 0\) and \(p = 1\) (a pure node) and reach their maximum at \(p = 0.5\) (a node split evenly between the two classes). The classification error is piecewise linear with a sharp kink at \(0.5\), whereas the Gini index and entropy are smooth and differentiable. Gini and entropy are also more sensitive to changes in node purity near the extremes — which is why they are preferred for growing trees (they reward splits that produce purer nodes), while classification error is typically used for pruning, when predictive accuracy is the goal. Entropy takes the largest values here (maximum \(\log 2 \approx 0.693\)) versus \(0.5\) for Gini and \(0.5\) for the error rate, but the three curves have very similar shape.

2 Exercise 8 — Regression trees on Carseats

2.1 (a) Split into a training and a test set

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

2.2 (b) Fit a regression tree

tree.carseats <- tree(Sales ~ ., data = Carseats, subset = train)
summary(tree.carseats)
## 
## Regression tree:
## tree(formula = Sales ~ ., data = Carseats, subset = 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)

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

Interpretation. The first split is on ShelveLoc (Good vs. Bad/Medium), confirming that shelving location is the single most important driver of Sales; subsequent splits use Price (lower price → higher sales), plus Age, CompPrice, Advertising and Income. The test MSE is about 4.92, i.e. an RMSE of roughly 2.22 thousand units.

2.3 (c) Cross-validation and pruning

set.seed(1)
cv.carseats <- cv.tree(tree.carseats)
plot(cv.carseats$size, cv.carseats$dev, type = "b",
     xlab = "Tree size (terminal nodes)", ylab = "CV deviance",
     main = "Cross-validation for pruning")
best.size <- cv.carseats$size[which.min(cv.carseats$dev)]
points(best.size, min(cv.carseats$dev), col = "red", pch = 19)

best.size
## [1] 18
prune.carseats <- prune.tree(tree.carseats, best = best.size)
plot(prune.carseats); text(prune.carseats, pretty = 0, cex = 0.7)

prune.pred <- predict(prune.carseats, Carseats.test)
prune.mse  <- mean((Carseats.test$Sales - prune.pred)^2)
c(unpruned_MSE = round(tree.mse, 3), pruned_MSE = round(prune.mse, 3))
## unpruned_MSE   pruned_MSE 
##        4.922        4.922

Cross-validation selects a tree with 18 terminal nodes, which is exactly the size of the unpruned tree. In other words, CV does not call for any pruning here: the CV deviance keeps decreasing out to the full tree, so the “pruned” tree is identical to the original one and the test MSE is unchanged (4.92). So on this split pruning does not improve the test MSE — the single tree is simply not flexible enough, and (as parts d–f show) the real gains come from ensembles rather than from tuning one tree’s size.

2.4 (d) Bagging

set.seed(1)
p <- ncol(Carseats) - 1        # 10 predictors
bag.carseats <- randomForest(Sales ~ ., data = Carseats, subset = train,
                             mtry = p, importance = TRUE)
bag.pred <- predict(bag.carseats, Carseats.test)
bag.mse  <- mean((Carseats.test$Sales - bag.pred)^2)
bag.mse
## [1] 2.605253
round(importance(bag.carseats), 2)
##             %IncMSE IncNodePurity
## CompPrice     24.89        170.18
## Income         4.71         91.26
## Advertising   12.77         97.16
## Population    -1.81         58.24
## Price         56.33        502.90
## ShelveLoc     48.89        380.03
## Age           17.73        157.85
## Education      0.60         44.60
## Urban          0.17          9.82
## US             4.22         18.07
varImpPlot(bag.carseats, main = "Bagging: variable importance")

Bagging cuts the test MSE to about 2.61, a large improvement over the single tree. The importance measures identify Price and ShelveLoc as by far the most important predictors, followed by CompPrice, Age and Advertising.

2.5 (e) Random forests and the effect of \(m\)

set.seed(1)
rf.carseats <- randomForest(Sales ~ ., data = Carseats, subset = train,
                            mtry = p / 3, importance = TRUE)
rf.pred <- predict(rf.carseats, Carseats.test)
rf.mse  <- mean((Carseats.test$Sales - rf.pred)^2)
rf.mse
## [1] 2.960559
round(importance(rf.carseats), 2)
##             %IncMSE IncNodePurity
## CompPrice     14.88        158.83
## Income         4.33        125.65
## Advertising    8.22        107.52
## Population    -0.95         97.06
## Price         34.98        385.93
## ShelveLoc     34.92        298.54
## Age           14.31        178.42
## Education      1.31         70.49
## Urban         -1.27         17.40
## US             6.11         33.99
set.seed(1)
m.vals <- 1:p
m.mse  <- sapply(m.vals, function(m) {
  fit <- randomForest(Sales ~ ., data = Carseats, subset = train, mtry = m)
  mean((Carseats.test$Sales - predict(fit, Carseats.test))^2)
})
plot(m.vals, m.mse, type = "b", xlab = "m (variables tried at each split)",
     ylab = "Test MSE", main = "Effect of m on random-forest test error")
points(which.min(m.mse), min(m.mse), col = "red", pch = 19)

data.frame(m = m.vals, test_MSE = round(m.mse, 3))
##     m test_MSE
## 1   1    4.800
## 2   2    3.469
## 3   3    3.036
## 4   4    2.805
## 5   5    2.687
## 6   6    2.599
## 7   7    2.671
## 8   8    2.580
## 9   9    2.620
## 10 10    2.601

Effect of \(m\). The test error is highest at \(m = 1\) (each split may only consider a single, possibly irrelevant, variable, so the trees are too weak and biased), falls rapidly as \(m\) increases, and then flattens out. The minimum here occurs at \(m = 8\), with test MSE 2.58. Because Carseats has only a few genuinely strong predictors (Price, ShelveLoc), larger \(m\) — up to bagging (\(m = p\)) — performs well; the usual variance-reduction benefit of a small \(m\) is muted since there are few correlated predictors to de-correlate.

2.6 (f) BART

x <- model.matrix(Sales ~ . - 1, data = Carseats)   # dummy-code the factors
y <- Carseats$Sales
set.seed(1)
bart.fit <- gbart(x[train, ], y[train], x.test = x[test, ])
bart.mse <- mean((y[test] - bart.fit$yhat.test.mean)^2)
bart.mse
## [1] 1.432639
# how often each variable was used in the ensemble's splitting rules
sort(round(colMeans(bart.fit$varcount / rowSums(bart.fit$varcount)), 4),
     decreasing = TRUE)[1:8]
##           Price   ShelveLocGood           USYes       CompPrice       Education 
##          0.1150          0.0890          0.0872          0.0866          0.0838 
## ShelveLocMedium    ShelveLocBad          Income 
##          0.0829          0.0791          0.0787
data.frame(
  method   = c("Single tree", "Pruned tree", "Bagging", "Random forest", "BART"),
  test_MSE = round(c(tree.mse, prune.mse, bag.mse, rf.mse, bart.mse), 3)
)
##          method test_MSE
## 1   Single tree    4.922
## 2   Pruned tree    4.922
## 3       Bagging    2.605
## 4 Random forest    2.961
## 5          BART    1.433

Summary. BART gives a test MSE of about 1.43 — the best of all five methods, roughly a 3.4-fold reduction relative to the single tree. BART’s most-used splitting variables are again Price and ShelveLoc, agreeing with the bagging and random-forest importance measures. Overall the ordering is BART < bagging < random forest \(\ll\) single/pruned tree: every ensemble clearly beats one tree, illustrating the standard trade-off in which we give up the interpretability of a single tree for substantially better predictive accuracy.

3 Exercise 9 — Classification trees on OJ

3.1 (a) Training set of 800 observations

set.seed(1)
train <- sample(nrow(OJ), 800)
OJ.train <- OJ[train, ]
OJ.test  <- OJ[-train, ]
dim(OJ.train); dim(OJ.test)
## [1] 800  18
## [1] 270  18

3.2 (b) Fit a tree and report the training error rate

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 summary reports a training misclassification error rate of 0.1588 (127 of 800 observations). Only a handful of the available variables are actually used, and the residual mean deviance is reported above.

3.3 (c) Plot the tree

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

n.leaves <- sum(tree.oj$frame$var == "<leaf>")
n.leaves
## [1] 9

The tree has 9 terminal nodes. LoyalCH (customer brand loyalty to Citrus Hill) is the dominant predictor — it forms the root split and several of the splits below it. Customers with high LoyalCH are predicted to buy CH; those with low loyalty are predicted to buy MM, with price-difference variables refining the middle range.

3.4 (d) Text summary of the tree

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

Interpreting a terminal node. Take the first terminal node listed (marked with *), node 8) LoyalCH < 0.0356415 59 10.14 MM ( 0.01695 0.98305 ) *. Reading the line: the split rule that defines the node is LoyalCH < 0.0356415; n = 59 training observations fall into it; its deviance is 10.14; the predicted class is MM; and the fitted class proportions are P(CH) = 0.017, P(MM) = 0.983. So of the 59 customers with almost no loyalty to Citrus Hill, 98.3% bought Minute Maid — a very pure node (hence the small deviance), and the tree confidently predicts MM for anyone landing there.

3.5 (e) Test-set predictions and confusion matrix

tree.test.pred <- predict(tree.oj, OJ.test, type = "class")
cm <- table(Predicted = tree.test.pred, Actual = OJ.test$Purchase)
cm
##          Actual
## Predicted  CH  MM
##        CH 160  38
##        MM   8  64
tree.test.err <- 1 - sum(diag(cm)) / sum(cm)
tree.test.err
## [1] 0.1703704

The test error rate is 0.1704 (about 17%).

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

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

plot(cv.oj$size, cv.oj$dev, type = "b",
     xlab = "Tree size (terminal nodes)",
     ylab = "CV misclassifications",
     main = "OJ: cross-validated error vs. tree size")
best.oj <- cv.oj$size[which.min(cv.oj$dev)]
points(best.oj, min(cv.oj$dev), col = "red", pch = 19)

3.8 (h) Tree size with the lowest CV error rate

best.oj
## [1] 9

The lowest cross-validated classification error corresponds to a tree with 9 terminal nodes. Note that the CV error is essentially tied across the larger sizes (sizes 9 and 8 both give 145 misclassifications, and size 7 is only one worse), so cross-validation does not clearly favour a smaller tree — it selects the full tree. Following the exercise’s instruction for exactly this situation, part (i) therefore prunes to five terminal nodes.

3.9 (i) The pruned tree

# if CV selects the full tree, use a 5-node tree as instructed
prune.size <- if (best.oj == n.leaves) 5 else best.oj
prune.oj <- prune.misclass(tree.oj, best = prune.size)
plot(prune.oj); text(prune.oj, pretty = 0, cex = 0.8)

prune.size
## [1] 5

3.10 (j) Training error rates: pruned vs. unpruned

train.err <- function(fit) {
  m <- summary(fit)$misclass
  m[1] / m[2]
}
c(unpruned = round(train.err(tree.oj), 4),
  pruned   = round(train.err(prune.oj), 4))
## unpruned   pruned 
##   0.1588   0.1625

The pruned tree has the higher (or equal) training error rate. This is expected: pruning removes splits, so the smaller tree can never fit the training data better than the larger tree from which it was derived — training error is monotonically non-increasing in tree size.

3.11 (k) Test error rates: pruned vs. unpruned

prune.test.pred <- predict(prune.oj, OJ.test, type = "class")
prune.test.err  <- mean(prune.test.pred != OJ.test$Purchase)
c(unpruned = round(tree.test.err, 4), pruned = round(prune.test.err, 4))
## unpruned   pruned 
##   0.1704   0.1630

The pruned tree has the lower test error rate (0.163 vs. 0.1704). This is the mirror image of part (j): pruning raises training error but can lower test error, because removing the deepest splits reduces overfitting. Here the 5-node tree is both simpler and at least as accurate on unseen data as the full 9-node tree, so it is clearly the better model.