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:
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).
Carseats data setdata(Carseats)
set.seed(1)
train <- sample(1:nrow(Carseats), nrow(Carseats) / 2)
Carseats.train <- Carseats[train, ]
Carseats.test <- Carseats[-train, ]
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.
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.
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).
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.
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
OJ data setdata(OJ)
set.seed(1)
train.oj <- sample(1:nrow(OJ), 800)
OJ.train <- OJ[train.oj, ]
OJ.test <- OJ[-train.oj, ]
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.
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.
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.
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.
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"
plot(cv.oj$size, cv.oj$dev, type = "b", pch = 19,
xlab = "Tree Size", ylab = "CV Classification Error (deviance count)")
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).
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)
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.
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.