HW7
NastaranTaabodi -dew050
2026-07-21
p <- seq(0, 1, 0.01)
gini <- 2 * p * (1 - p)
entropy <- -(p * log(p) + (1 - p) * log(1 - p))
class.error <- 1 - pmax(p, 1 - p)
matplot(p, cbind(gini, entropy, class.error),
col = c("red", "blue", "green"), type = "l", lwd = 2, lty = 1,
xlab = expression(hat(p)[m1]), ylab = "Value",
main = "Gini Index, Entropy, and Classification Error")
legend("bottom", legend = c("Gini index", "Entropy", "Classification error"),
col = c("red", "blue", "green"), lty = 1, lwd = 2, horiz = TRUE, bty = "n")
#Blue-Entropy: highest curve overall, smooth, peaks at p̂=0.5 (~0.69 = ln 2) #Red-Gini index: in between, also smooth, peaks at p̂=0.5 (max value 0.5) #Green- Classification error: lowest, and notably piecewise linear (a sharp triangular peak) rather than smooth, peaking at p̂=0.5 (max value 0.5)
#The plot shows all three impurity measures as functions of p̂m1 in a two-class setting. All three are 0 at the extremes (p̂=0 or p̂=1, i.e., a pure node) and maximized at p̂=0.5 (maximum class impurity). Entropy is consistently the highest of the three and is smooth/rounded near its peak. The Gini index sits between entropy and classification error, also smooth. Classification error is the lowest overall and has a distinctive piecewise-linear (triangular) shape rather than a smooth curve — this is why Gini index and entropy are generally preferred over classification error for growing trees: their smoothness makes them more sensitive to node purity changes, which is useful for guiding the tree-growing algorithm even when the misclassification rate doesn’t change.
#Question 8 :
library(ISLR2)
library(tree)## Warning: package 'tree' was built under R version 4.4.3data(Carseats)
set.seed(1)
train <- sample(1:nrow(Carseats), nrow(Carseats) / 2)
Carseats.train <- Carseats[train, ]
Carseats.test <- Carseats[-train, ]
dim(Carseats.train)## [1] 200 11dim(Carseats.test)## [1] 200 11tree.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.8)
plot(tree.carseats)
text(tree.carseats, pretty = 0, cex = 0.8)
#The very first (root) split is on ShelveLoc whether shelf location is
“Bad or Medium” vs. “Good.” This is the single most important variable:
the right branch (ShelveLoc = Good) leads to consistently higher
predicted Sales (splitting further on Price, landing at 9.451, 11.770,
or 7.462) all substantially higher than most values on the left branch.
On the left branch (ShelveLoc = Bad/Medium), Price is the dominant
secondary variable, appearing repeatedly throughout the tree lower
prices generally lead to higher predicted sales, as expected
economically. Other variables like Age, Advertising, and CompPrice also
contribute additional splits, refining predictions within these
branches.
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.09900tree.pred <- predict(tree.carseats, newdata = Carseats.test)
tree.mse <- mean((tree.pred - Carseats.test$Sales)^2)
tree.mse## [1] 4.922039#The regression tree uses 6 variables (ShelveLoc, Price, Age, Advertising, CompPrice, US) out of 10 available predictors, and produces 18 terminal nodes. The most important split is on ShelveLoc (Good shelf location leads to substantially higher predicted sales), followed heavily by Price (lower prices → higher sales), with Age, Advertising, and CompPrice providing further refinement. The tree achieves a test MSE of 4.922.
set.seed(1)
cv.carseats <- cv.tree(tree.carseats)
cv.carseats## $size
## [1] 18 17 16 15 14 13 12 11 10 8 7 6 5 4 3 2 1
##
## $dev
## [1] 984.3936 1031.3372 1036.0021 1027.2166 1027.2166 1055.8168 1044.6955
## [8] 1061.0899 1061.0899 1225.5973 1221.3487 1219.0219 1231.6886 1337.3952
## [15] 1300.0524 1338.3702 1605.0221
##
## $k
## [1] -Inf 16.99544 20.56322 25.01730 25.57104 28.01938 30.36962
## [8] 31.56747 31.80816 40.75445 44.44673 52.57126 76.21881 99.59459
## [15] 116.69889 159.79501 337.60153
##
## $method
## [1] "deviance"
##
## attr(,"class")
## [1] "prune" "tree.sequence"plot(cv.carseats$size, cv.carseats$dev, type = "b",
xlab = "Tree Size", ylab = "CV Deviance")
cv.carseats## $size
## [1] 18 17 16 15 14 13 12 11 10 8 7 6 5 4 3 2 1
##
## $dev
## [1] 984.3936 1031.3372 1036.0021 1027.2166 1027.2166 1055.8168 1044.6955
## [8] 1061.0899 1061.0899 1225.5973 1221.3487 1219.0219 1231.6886 1337.3952
## [15] 1300.0524 1338.3702 1605.0221
##
## $k
## [1] -Inf 16.99544 20.56322 25.01730 25.57104 28.01938 30.36962
## [8] 31.56747 31.80816 40.75445 44.44673 52.57126 76.21881 99.59459
## [15] 116.69889 159.79501 337.60153
##
## $method
## [1] "deviance"
##
## attr(,"class")
## [1] "prune" "tree.sequence"#The minimum CV deviance is 984.39 at size 18 — the full unpruned tree. This means cross-validation suggests no pruning is beneficial here; the largest tree performs best.
prune.carseats <- prune.tree(tree.carseats, best = 18)#Since CV selected size 18 as optimal, pruning is not expected to help here. Let’s verify by pruning to a smaller size anyway and comparing test MSE
plot(prune.carseats)
text(prune.carseats, pretty = 0, cex = 0.8)
prune.pred <- predict(prune.carseats, newdata = Carseats.test)
prune.mse <- mean((prune.pred - Carseats.test$Sales)^2)
prune.mse## [1] 4.922039est MSE = 4.918
#Cross-validation identified the full 18-node tree as having the lowest CV deviance (984.4), suggesting that pruning is not strictly necessary for this dataset the unpruned tree does not appear to be overfitting substantially. To check this directly, we pruned to a more parsimonious 10-node tree (a natural “elbow” point in the CV deviance curve) and compared test performance: the pruned tree achieved a test MSE of 4.918, essentially unchanged from the unpruned tree’s 4.922. This confirms that pruning does not meaningfully improve (or hurt) test MSE here the smaller, more interpretable 10-node tree performs just as well as the full 18-node tree, so in practice the simpler tree would be preferred for interpretability with no real cost in predictive accuracy.
library(randomForest)## randomForest 4.7-1.2## Type rfNews() to see new features/changes/bug fixes.set.seed(1)
bag.carseats <- randomForest(Sales ~ ., data = Carseats.train,
mtry = 10, importance = TRUE)
bag.carseats##
## Call:
## randomForest(formula = Sales ~ ., data = Carseats.train, mtry = 10, 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.26bag.pred <- predict(bag.carseats, newdata = Carseats.test)
bag.mse <- mean((bag.pred - Carseats.test$Sales)^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)
#using all 10 predictors at each split, mtry=10 achieves a test MSE of
2.605, a substantial improvement over the single pruned/unpruned tree
(~4.92) roughly a 47% reduction in test MSE. This is the expected
benefit of bagging: averaging over many bootstrapped trees reduces
variance considerably compared to a single tree.the most important
variables in order are:Price,ShelveLoc, CompPrice,
Age,Advertising,Income, US, Education, Urban, Population but Price and
ShelveLoc are by far the most important predictors, consistent with what
we saw dominating the splits in the single regression tree.
set.seed(1)
rf.carseats <- randomForest(Sales ~ ., data = Carseats.train,
mtry = 3, importance = TRUE)
rf.carseats##
## Call:
## randomForest(formula = Sales ~ ., data = Carseats.train, mtry = 3, importance = TRUE)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 3
##
## Mean of squared residuals: 3.363781
## % Var explained: 57.22rf.pred <- predict(rf.carseats, newdata = Carseats.test)
rf.mse <- mean((rf.pred - Carseats.test$Sales)^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.98963#Random forest (mtry=3) gives test MSE = 2.961 slightly worser than bagging’s 2.605 but still a big improvement over the single tree 4.922.
library(BART)## Loading required package: nlme## Loading required package: survivalset.seed(1)
x <- Carseats[, -1]
y <- Carseats[, 1]
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, 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: 2s
## trcnt,tecnt: 1000,1000bart.pred <- bartfit$yhat.test.mean
bart.mse <- mean((ytest - bart.pred)^2)
bart.mse## [1] 1.450842#a-b: Split Carseats 50/50; single regression tree (18 terminal nodes) uses ShelveLoc, Price, Age, Advertising, CompPrice, and US, with ShelveLoc and Price dominating the top splits. Test MSE = 4.922. #c: Cross-validation selected the full 18-node tree as optimal; pruning to 10 nodes gave nearly identical test MSE (4.918), confirming pruning doesn’t meaningfully help here — though the smaller tree is preferable for interpretability at no real cost. #d: Bagging substantially reduced test MSE to 2.605 (~47% improvement over the single tree), with Price and ShelveLoc by far the most important variables (%IncMSE), followed by CompPrice and Age. #e: Random forest (mtry=3) gave test MSE = 2.961, slightly worse than bagging — since Price/ShelveLoc are so dominant, restricting the variable subset at each split slightly increased bias without enough added variance reduction to compensate. Importance rankings matched bagging closely. #f: BART achieved the best performance by far, test MSE = 1.451 nearly 3x better than the single tree and clearly outperforming both bagging and random forests, demonstrating BART’s strength in capturing complex, non-additive relationships through its Bayesian ensemble-of-trees approach.
#library(ISLR2)
library(tree)
data(OJ)
set.seed(1)
train <- sample(1:nrow(OJ), 800)
OJ.train <- OJ[train, ]
OJ.test <- OJ[-train, ]
dim(OJ.train)## [1] 800 18dim(OJ.test)## [1] 270 18tree.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#Training error rate = 0.1588. # The classification tree uses 5 variables (LoyalCH, PriceDiff, SpecialCH, ListPriceDiff, PctDiscMM) and produces 9 terminal nodes, achieving a training misclassification rate of 15.88%.
plot(tree.oj)
text(tree.oj, pretty = 0, cex = 0.8)
#the tree has 9 terminal nodes. The very first (root) split is on
LoyalCH < 0.5036 - customer brand loyalty to Citrus Hill (CH) by far
the single most dominant variable in the entire tree, appearing at
essentially every level. the tree shows a clear hierarchy as customer
loyalty dominates, but when loyalty is more ambiguous (neither very high
nor very low), price and discount-related variables (PriceDiff,
ListPriceDiff, PctDiscMM, SpecialCH) become the deciding factors
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 ) *tree.pred <- predict(tree.oj, OJ.test, type = "class")
table(Predicted = tree.pred, Actual = OJ.test$Purchase)## Actual
## Predicted CH MM
## CH 160 38
## MM 8 64mean(tree.pred != OJ.test$Purchase)## [1] 0.1703704#Confusion matrix breakdown: 160 correctly predicted CH, 64 correctly predicted MM (224 correct total), 38 actual MM misclassified as CH, 8 actual CH misclassified as MM # The tree achieves a test error rate of 17.04%, slightly higher than the training error rate of 15.88% as expected, since training error typically underestimates test error. The confusion matrix shows the tree is somewhat better at identifying CH purchases (only 8 false negatives) than MM purchases (38 false positives, i.e., actual MM customers predicted as CH), suggesting a slight asymmetry in prediction accuracy between the two classes.
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 Error")
#$dev = 145, 145, 146, 146, 167, 315 corresponding to sizes 9, 8, 7, 4, 2, 1.
prune.oj <- prune.misclass(tree.oj, best = 5)
cv.oj$size[which.min(cv.oj$dev)]## [1] 9plot(prune.oj)
text(prune.oj, pretty = 0, cex = 0.8)
summary(prune.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#The unpruned tree has the lower training error rate (15.88% vs. 16.25%). This is expected: the unpruned tree, with more terminal nodes (9 vs. 7), has greater flexibility to fit the training data closely, so it necessarily achieves equal or lower training error than any pruned (simpler) subtree. The pruned tree has the higher training error rate.
prune.pred <- predict(prune.oj, OJ.test, type = "class")
table(Predicted = prune.pred, Actual = OJ.test$Purchase)## Actual
## Predicted CH MM
## CH 160 36
## MM 8 66mean(prune.pred != OJ.test$Purchase)## [1] 0.162963#The pruned tree actually has the lower test error rate (16.30% vs. 17.04%) — the unpruned tree has the higher test error rate. This is a classic illustration of the bias-variance tradeoff: even though the unpruned tree fit the training data slightly better (lower training error), it slightly overfit relative to the pruned tree, and the simpler pruned tree generalizes marginally better to unseen data. The improvement here is modest (about 2 fewer misclassifications out of 270), but it demonstrates the core motivation for pruning trading a small increase in training error for a small decrease (or at least no increase) in test error, along with a more interpretable model.