Assignment_7
Cassie Pruski
2026-07-17
3 Gini index, classification error, and entropy
p <- seq(0, 1, length.out = 200)
gini <- 2 * p * (1 - p)
entropy <- -(p * log(p) + (1 - p) * log(1 - p))
entropy[is.nan(entropy)] <- 0
class.error <- 1 - pmax(p, 1 - p)
plot(p, gini, type = "l", col = "red", lwd = 2, ylim = c(0, 0.75),
xlab = expression(hat(p)[m1]), ylab = "Value",
main = "Gini, Entropy, and Classification Error")
lines(p, entropy, col = "blue", lwd = 2)
lines(p, class.error, col = "darkgreen", lwd = 2)
legend("topright", legend = c("Gini index", "Entropy", "Classification error"),
col = c("red", "blue", "darkgreen"), lwd = 2)
8a
## Warning: package 'ISLR2' was built under R version 4.5.3## Warning: package 'tree' was built under R version 4.5.3## Warning: package 'randomForest' was built under R version 4.5.3## randomForest 4.7-1.2## Type rfNews() to see new features/changes/bug fixes.train <- sample(1:nrow(Carseats), nrow(Carseats) / 2)
Carseats.train <- Carseats[train, ]
Carseats.test <- Carseats[-train, ]8b
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.09900plot(tree.carseats)
text(tree.carseats, pretty = 0, cex = 0.7)
pred <- predict(tree.carseats, newdata = Carseats.test)
tree.mse <- mean((Carseats.test$Sales - pred)^2)
tree.mse## [1] 4.9220398c
set.seed(1)
cv.carseats <- cv.tree(tree.carseats)
plot(cv.carseats$size, cv.carseats$dev, type = "b",
xlab = "Tree size", ylab = "CV deviance")
best.size <- cv.carseats$size[which.min(cv.carseats$dev)]
best.size## [1] 18prune.carseats <- prune.tree(tree.carseats, best = best.size)
plot(prune.carseats)
text(prune.carseats, pretty = 0, cex = 0.7)
pred.prune <- predict(prune.carseats, newdata = Carseats.test)
prune.mse <- mean((Carseats.test$Sales - pred.prune)^2)
prune.mse## [1] 4.9220398d
set.seed(1)
bag.carseats <- randomForest(Sales ~ ., data = Carseats.train,
mtry = ncol(Carseats.train) - 1,
importance = TRUE)
pred.bag <- predict(bag.carseats, newdata = Carseats.test)
bag.mse <- mean((Carseats.test$Sales - pred.bag)^2)
bag.mse## [1] 2.605253importance(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.073863varImpPlot(bag.carseats)
8e
set.seed(1)
rf.carseats <- randomForest(Sales ~ ., data = Carseats.train,
mtry = (ncol(Carseats.train) - 1) / 3,
importance = TRUE)
pred.rf <- predict(rf.carseats, newdata = Carseats.test)
rf.mse <- mean((Carseats.test$Sales - pred.rf)^2)
rf.mse## [1] 2.960559importance(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.98963varImpPlot(rf.carseats)
mtry.vals <- 1:10
mse.vals <- rep(NA, length(mtry.vals))
for (m in mtry.vals) {
set.seed(1)
fit <- randomForest(Sales ~ ., data = Carseats.train, mtry = m)
pred <- predict(fit, newdata = Carseats.test)
mse.vals[m] <- mean((Carseats.test$Sales - pred)^2)
}
plot(mtry.vals, mse.vals, type = "b", xlab = "mtry", ylab = "Test MSE")
The random forest with m = p/3 ≈ 3 achieved a test MSE of 2.961, which
is better than the single pruned tree (4.922) but slightly worse than
full bagging with m = p = 10 (2.605). importance() shows Price and
ShelveLoc are by far the most important predictors, followed by
CompPrice and Age, with Population and Urban not contributing much. Test
MSE generally decreases as m increases. Restricting the number of
candidate variables at each split hurts more than it helps, because
Price and ShelveLoc are very dominant, so excluding them from many
splits costs more in individual tree quality than is gained from
decorrelating the trees.
8f
library(BART)## Warning: package 'BART' was built under R version 4.5.3## Loading required package: nlme## Loading required package: survivalx <- Carseats[, -1]
y <- Carseats[, "Sales"]
xtrain <- x[train, ]
ytrain <- y[train]
xtest <- x[-train, ]
ytest <- y[-train]
set.seed(1)
bart.fit <- gbart(xtrain, ytrain, x.test = xtest)## *****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,1000pred.bart <- bart.fit$yhat.test.mean
bart.mse <- mean((ytest - pred.bart)^2)
bart.mse## [1] 1.450842Fitting BART with gbart() yielded a test MSE of 1.451, which is the best result. Rather than growing many independent trees on bootstrap samples, BART fits an ensemble of intentionally weak, shallow trees, where each tree explains only the residual left by the others. It also discourages any single tree from becoming too influential, so this produces a more efficient, better-generalizing model here.
9a Split into training and test sets
data(OJ)
set.seed(1)
train <- sample(1:nrow(OJ), 800)
OJ.train <- OJ[train, ]
OJ.test <- OJ[-train, ]9b Fit a regression tree, plot it, compute test MSE
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 / 8009c Cross-validation to determine optimal tree complexity via pruning
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 ) *9d
plot(tree.oj)
text(tree.oj, pretty = 0, cex = 0.7)
##9e
tree.pred <- predict(tree.oj, newdata = OJ.test, type = "class")
table(Predicted = tree.pred, Actual = OJ.test$Purchase)## Actual
## Predicted CH MM
## CH 160 38
## MM 8 64test.error <- mean(tree.pred != OJ.test$Purchase)
test.error## [1] 0.17037049f & g
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",
xlab = "Tree size", ylab = "CV classification errors")
9h
Cross-validation showed sizes 9 and 8 tied for the lowest classification error of 145 misclassifications, with all smaller sizes performing worse. Size 8 was chosen for pruning since it achieves the same CV performance as the full tree with one fewer terminal node.
best.size <- 8
best.size## [1] 89i
prune.oj <- prune.misclass(tree.oj, best = best.size)
plot(prune.oj)
text(prune.oj, pretty = 0, cex = 0.7)
9j
The pruned 8-node tree’s training error of 0.15875 is identical to the unpruned 9-node tree’s training error of 0.15875.
pruned.train.pred <- predict(prune.oj, newdata = OJ.train, type = "class")
pruned.train.error <- mean(pruned.train.pred != OJ.train$Purchase)
pruned.train.error## [1] 0.15875unpruned.train.pred <- predict(tree.oj, newdata = OJ.train, type = "class")
unpruned.train.error <- mean(unpruned.train.pred != OJ.train$Purchase)
unpruned.train.error## [1] 0.158759k
Likewise, the pruned tree’s test error of 0.1703704 equals the unpruned tree’s test error of 0.1703704.
pruned.test.pred <- predict(prune.oj, newdata = OJ.test, type = "class")
pruned.test.error <- mean(pruned.test.pred != OJ.test$Purchase)
pruned.test.error## [1] 0.1703704unpruned.test.pred <- predict(tree.oj, newdata = OJ.test, type = "class")
unpruned.test.error <- mean(unpruned.test.pred != OJ.test$Purchase)
unpruned.test.error## [1] 0.1703704