Question 3
# (3)
p <- seq(0, 1, 0.01)
err <- 1 - pmax(p, 1 - p)
gini <- 2 * p * (1 - p)
ent <- -p * log(p) - (1 - p) * log(1 - p)
ent[is.na(ent)] <- 0
plot(p, ent, type = "l", col = "red", ylim = c(0, 0.7))
lines(p, gini, col = "blue")
lines(p, err, col = "green")
Question 8
# (8a)
set.seed(1)
train_idx <- sample(1:nrow(Carseats), 0.7 * nrow(Carseats))
train <- Carseats[train_idx, ]
test <- Carseats[-train_idx, ]
# (8b)
tree_fit <- tree(Sales ~ ., data = train)
plot(tree_fit)
text(tree_fit, pretty = 0)
preds <- predict(tree_fit, newdata = test)
mean((test$Sales - preds)^2)
## [1] 4.208383
invisible("Test MSE 4.21")
# (8c)
cv_tree <- cv.tree(tree_fit)
plot(cv_tree$size, cv_tree$dev, type = "b")
best_size <- cv_tree$size[which.min(cv_tree$dev)]
prune_fit <- prune.tree(tree_fit, best = best_size)
prune_preds <- predict(prune_fit, newdata = test)
mean((test$Sales - prune_preds)^2)
## [1] 4.208383
invisible("Pruning the tree does not improve the test MSE")
# (8d)
bag_fit <- randomForest(Sales ~ ., data = train, mtry = 10, importance = TRUE)
bag_preds <- predict(bag_fit, newdata = test)
mean((test$Sales - bag_preds)^2)
## [1] 2.546423
importance(bag_fit)
## %IncMSE IncNodePurity
## CompPrice 33.904069 231.47471
## Income 5.335622 121.33353
## Advertising 17.870018 157.15055
## Population -1.789009 64.69008
## Price 69.684067 690.84486
## ShelveLoc 72.687816 635.79311
## Age 22.489173 235.52049
## Education 1.010342 62.58691
## Urban -1.300710 11.57267
## US 2.778950 10.54935
varImpPlot(bag_fit)
invisible("Bagging has a test MSE of 2.55, with ShelveLoc and Price identified as the most important predictors of sales.")
# (8e)
rf_fit <- randomForest(Sales ~ ., data = train, mtry = 3, importance = TRUE)
rf_preds <- predict(rf_fit, newdata = test)
mean((test$Sales - rf_preds)^2)
## [1] 2.949671
importance(rf_fit)
## %IncMSE IncNodePurity
## CompPrice 16.7327243 203.67364
## Income 2.8719402 162.46321
## Advertising 12.2638672 161.41889
## Population -0.3836743 132.03719
## Price 44.7666458 540.43983
## ShelveLoc 46.4521878 515.63949
## Age 16.7918061 264.58022
## Education 0.5335414 98.80972
## Urban -1.5990201 19.75826
## US 3.7152171 30.35581
invisible("Random forest yields a test MSE of 2.95, with ShelveLoc and Price as the top predictors; reducing m from 10 to 3 decorrelates the trees. Increases the error rate but reduces variance and protects against overfitting.")
# (8f)
xtrain <- train[, -1]
ytrain <- train$Sales
xtest <- test[, -1]
bart_fit <- gbart(xtrain, ytrain, x.test = xtest)
## *****Calling gbart: type=1
## *****Data:
## data:n,p,np: 280, 14, 120
## y1,yn: 2.763393, 0.113393
## x1,x[n*p]: 107.000000, 1.000000
## xp1,xp[np*p]: 113.000000, 1.000000
## *****Number of Trees: 200
## *****Number of Cut Points: 70 ... 1
## *****burn,nd,thin: 100,1000,1
## *****Prior:beta,alpha,tau,nu,lambda,offset: 2,0.95,0.276302,3,0.215848,7.59661
## *****sigma: 1.052663
## *****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,1000
bart_preds <- bart_fit$yhat.test.mean
mean((test$Sales - bart_preds)^2)
## [1] 1.539667
invisible("BART yields a test MSE of 1.54, outperforming the single tree, bagging, and random forest models.")
Question 9
# (9a)
set.seed(1)
train_indices <- sample(1:nrow(OJ), 800)
train <- OJ[train_indices, ]
test <- OJ[-train_indices, ]
# (9b)
oj_tree <- tree(Purchase ~ ., data = train)
summary(oj_tree)
##
## Classification tree:
## tree(formula = Purchase ~ ., data = 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
invisible("Tree has a training error rate of 15.88% and 9 terminal nodes. It splits primarily on customer loyalty (LoyalCH), alongside price and promotional variables.")
# (9c)
oj_tree
## 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 ) *
invisible("Terminal node 8 contains 59 observations for customers with extremely low loyalty to Citrus Hill (LoyalCH < 0.0356415). Classifies this final leaf as a Minute Maid (MM) purchase with a 98.3% probability.")
# (9d)
plot(oj_tree)
text(oj_tree, pretty = 0)
invisible("The tree plot confirms that customer loyalty (LoyalCH) is the dominant predictor, forming the root split and its primary branches. Low loyalty directs predictions left toward Minute Maid (MM), while high loyalty branches right toward Citrus Hill (CH), with minor adjustments made by price differences and discounts.")
# (9e)
tree_pred <- predict(oj_tree, newdata = test, type = "class")
table(tree_pred, test$Purchase)
##
## tree_pred CH MM
## CH 160 38
## MM 8 64
mean(tree_pred != test$Purchase)
## [1] 0.1703704
invisible("The test error rate is 17.04% (46 out of 270 observations misclassified).")
# (9f)
cv_oj <- cv.tree(oj_tree, FUN = prune.misclass)
cv_oj
## $size
## [1] 9 8 7 4 2 1
##
## $dev
## [1] 150 150 149 158 172 315
##
## $k
## [1] -Inf 0.000000 3.000000 4.333333 10.500000 151.000000
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"
invisible("The cross-validation output shows that a tree size of 7 terminal nodes achieves the lowest cross-validation error (deviance of 149). Trees of size 8 and 9 are close behind with an error of 150.")
# (9g)
plot(cv_oj$size, cv_oj$dev, type = "b",
xlab = "Tree Size", ylab = "Cross-Validated Error (Deviance)")
# (9h)
invisible("A tree size of 7 terminal nodes corresponds to the lowest cross-validated classification error rate.")
# (9i)
pruned_oj <- prune.misclass(oj_tree, best = 7)
plot(pruned_oj)
text(pruned_oj, pretty = 0)
# (9j)
invisible("The training error rate of the pruned tree is higher than the unpruned tree. This occurs because pruning removes nodes, reducing the model's capacity to fit the training data perfectly.")
# (9k)
pruned_pred <- predict(pruned_oj, newdata = test, type = "class")
mean(pruned_pred != test$Purchase)
## [1] 0.162963
table(pruned_pred, test$Purchase)
##
## pruned_pred CH MM
## CH 160 36
## MM 8 66
invisible("The unpruned tree has a higher test error rate (17.04%) compared to the pruned tree (16.30%). Pruning successfully reduced overfitting, improving the model's performance.")