Assignment #8.
Tree-based Methods
Richmond Owusu Duah
Question 8
In the lab, a classification tree was applied to the Carseats data set after converting Sales
into a qualitative response variable. Now we will seek to predict Sales using regression
trees and related approaches, treating the response as a quantitative variable.
a. Split the data set into a training set and a test set.
library(ISLR2)
## Warning: package 'ISLR2' was built under R version 4.4.2
library(tree)
## Warning: package 'tree' was built under R version 4.4.2
library(randomForest)
## Warning: package 'randomForest' was built under R version 4.4.2
## randomForest 4.7-1.2
## Type rfNews() to see new features/changes/bug fixes.
library(BART)
## Warning: package 'BART' was built under R version 4.4.2
## Loading required package: nlme
## Loading required package: survival
head(Carseats)
## Sales CompPrice Income Advertising Population Price ShelveLoc Age Educat
ion
## 1 9.50 138 73 11 276 120 Bad 42
17
## 2 11.22 111 48 16 260 83 Good 65
10
## 3 10.06 113 35 10 269 80 Medium 59
12
## 4 7.40 117 100 4 466 97 Medium 55
14
## 5 4.15 141 64 3 340 128 Bad 38
13
## 6 10.81 124 113 13 501 72 Bad 78
16
## Urban US
## 1 Yes Yes
## 2 Yes Yes
## 3 Yes Yes
## 4 Yes Yes
## 5 Yes No
## 6 No Yes
data(Carseats)
[Link](42)
train <- sample(c(TRUE, FALSE), nrow(Carseats), replace = TRUE)
b. Fit a regression tree to the training set. Plot the tree, and interpret the
results. What test MSE do you obtain?
tr <- tree(Sales ~ ., data = Carseats[train, ])
summary(tr)
##
## Regression tree:
## tree(formula = Sales ~ ., data = Carseats[train, ])
## Variables actually used in tree construction:
## [1] "ShelveLoc" "Price" "Income" "Advertising" "CompPrice"
## [6] "Age"
## Number of terminal nodes: 16
## Residual mean deviance: 2.356 = 424.1 / 180
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -4.54900 -0.82980 0.03075 0.00000 0.89250 4.83100
plot(tr)
text(tr, pretty = 0, digits = 2, cex = 0.8)
carseats_mse <- function(model) {
p <- predict(model, newdata = Carseats[!train, ])
mean((p - Carseats[!train, "Sales"])^2)
}
carseats_mse(tr)
## [1] 4.559764
Comment: We may conclude that the Test MSE is about 4.56
c. Use cross-validation in order to determine the optimal level of tree
complexity. Does pruning the tree improve the test MSE?
res <- [Link](tr)
plot(res$size, res$dev, type = "b", xlab = "Tree size", ylab = "Deviance")
min <- [Link](res$dev)
abline(v = res$size[min], lty = 2, col = "red")
ptr <- [Link](tr, best = res$size[min])
plot(ptr)
text(ptr, pretty = 0, digits = 2, cex = 0.8)
carseats_mse(ptr)
## [1] 5.023125
Comment: Yes, pruning the tree increases the Test MSE to 5.02
d. Use the bagging approach in order to analyze this data. What test MSE
do you obtain? Use the importance() function to determine which
variables are most important.
bagged <- randomForest(Sales ~ ., data = Carseats[train, ], mtry = 10,
ntree = 200, importance = TRUE)
carseats_mse(bagged)
## [1] 2.762861
Comment: Bagging decreases the Test MSE to 2.8
importance(bagged)
## %IncMSE IncNodePurity
## CompPrice 11.2608998 104.474222
## Income 5.0953983 73.275066
## Advertising 12.9011190 125.886762
## Population 3.4071044 60.095200
## Price 34.6904380 450.952728
## ShelveLoc 33.7059874 374.808575
## Age 7.9101141 143.652934
## Education -2.1154997 32.712444
## Urban 0.9604097 7.029648
## US 3.1336559 6.287048
Comment: “Price” and “ShelveLoc” are the two most important variables.
e. Use random forests to analyze this data. What test error rate do you
obtain? Use the importance() function to determine which variables are
most important. Describe the effect of m, the number of variables
considered at each split, on the error rate obtained.
rf <- randomForest(Sales ~ ., data = Carseats[train, ], mtry = 3,
ntree = 500, importance = TRUE)
carseats_mse(rf)
## [1] 3.439357
Comment: Test MSE of 3.4
importance(rf)
## %IncMSE IncNodePurity
## CompPrice 8.5717587 122.75189
## Income 2.8955756 116.33951
## Advertising 13.0681194 128.13563
## Population 2.0475415 104.03803
## Price 34.7934136 342.84663
## ShelveLoc 39.0704834 292.56638
## Age 7.7941744 135.69061
## Education 0.8770806 64.67614
## Urban -0.3301478 13.83594
## US 6.2716539 22.07306
Comment: In this case also, “Price” and “ShelveLoc” are the two most important
variables.
f. Now analyze the data using BART, and report your results.
custom_predict_bart <- function(model, ...) model$[Link]
x_train <- Carseats[train, -which(names(Carseats) == "Sales")]
y_train <- Carseats[train, "Sales"]
x_test <- Carseats[!train, -which(names(Carseats) == "Sales")]
bartfit <- gbart(x_train, y_train, [Link] = x_test)
## *****Calling gbart: type=1
## *****Data:
## data:n,p,np: 196, 14, 204
## y1,yn: 2.070867, 2.280867
## x1,x[n*p]: 138.000000, 1.000000
## xp1,xp[np*p]: 141.000000, 1.000000
## *****Number of Trees: 200
## *****Number of Cut Points: 58 ... 1
## *****burn,nd,thin: 100,1000,1
## *****Prior:beta,alpha,tau,nu,lambda,offset: 2,0.95,0.287616,3,0.21118,7.42
913
## *****sigma: 1.041218
## *****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: 5s
## trcnt,tecnt: 1000,1000
carseats_mse <- function(model) {
mean((custom_predict_bart(model) - Carseats[!train, "Sales"])^2)
}
carseats_mse(bartfit)
## [1] 1.631285
Comment: The test MSE for BART is 1.6, which is an improvement over random forest and
bagging.