This assignment will cover tree-based method exercises from Chapter 8 of the ISLR2 textbook. The methods I work with are random forest, boosting, decision trees, and regression trees. I utilize plots, rpart plots, importance plots, confusion matrix’s, and a variety of other tools in this lab.
Consider the Gini index, classification error, and entropy in a simple classification setting with two classes. Create a single plot that displays each of these quantities as a function of pˆm1. The x-axis should display pˆm1, ranging from 0 to 1, and the y-axis should display the value of the Gini index, classification error, and entropy. Hint: In a setting with two classes, pˆm1 =1 − pˆm2. You could make this plot by hand, but it will be much easier to make in R.
# Code for question 3 provided by Proffessor Joseph Campbell
p=seq(0,1,0.0001)
#Gini
G=2*p*(1-p)
#Classification Error
E=1-pmax(p,1-p)
#Entropy
D=-(p*log(p) + (1-p)*log(1-p))
plot(p,D, col="red",ylab="")
lines(p,E,col='green')
lines(p,G,col='blue')
legend(0.3,0.15,c("Entropy", "Missclassification","Gini"),lty=c(1,1,1),lwd=c(2.5,2.5,2.5),col=c('red','green','blue'))
Carseats DataIn 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.
attach(Carseats)
# Visualize the response variable
z<-Carseats |> ggplot(
aes(
Sales)
) +
geom_histogram(
col="gray",
fill ="navyblue",
) +
labs(x = "Sales",
y = "Count",
title = "Distribution of Carseat Sales Units in Thousands (per location)")
ggplotly(z)
Initial Thoughts: The response variable
Sales seems to be normally distributed with a high volume
spike in Sales in the 5 - 9 (thousand) ranges
Split the data set into a training set and a test set.
set.seed(1)
# create the test/train split. Here I perform a 70/30 split
train <- createDataPartition(Carseats$Sales, p=0.7, list=FALSE)
train_data <- Carseats[train, ]
test_data <- Carseats[-train, ]
Fit a regression tree to the training set. Plot the tree, and interpret the results. What test MSE do you obtain?
# fit the regression tree
tree_fit <- tree(Sales ~., Carseats, subset = train)
summary(tree_fit)
##
## Regression tree:
## tree(formula = Sales ~ ., data = Carseats, subset = train)
## Variables actually used in tree construction:
## [1] "ShelveLoc" "Price" "Income" "Advertising" "CompPrice"
## [6] "Age"
## Number of terminal nodes: 18
## Residual mean deviance: 2.22 = 583.8 / 263
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.64200 -0.96000 0.02174 0.00000 0.96000 3.40700
plot(tree_fit)
text(tree_fit,
pretty=0,
cex = .65,
digits =3)
Initial Thoughts: This tree tells us a lot about which variables do the best job of predicting Sales. For readability remember that Left = Yes (the condition is TRUE) and Right = No (the condition is FALSE). Starting at the top of the tree and consider the following:
ShelveLoc (the quality of the shelving location for the car seats at each site) is the most important predictor of Sales
When ShelveLoc is Good, Price is less than 105.5, Age is less than 49.5, income is less than 61, then the predicted value of Sales is 8.43 (in thousands)
When ShelveLoc is Bad/Medium and Price is less than 109.5 then the predicted value of sales is 12.20 (in thousands)
When SelveLoc is Bad/Medium, price is greater than 109.5, advertising is greater than 0.5, and income is less than 35 then predicted value of Sales is 7.55 (in thousands)
# Calculate predictions
yhat.reg <- predict(tree_fit, newdata = Carseats[-train, ])
Carseats.test <- Carseats[-train, "Sales"]
# plot the predictions
plot(yhat.reg,
Carseats.test,
xlab = "Predicted Sales",
ylab = "Actual Sales",
main = "Single Regression Tree - Actual vs Predicted Sales")
abline(0,1, col = "red", lwd = 2)
# Calculate test MSE
mean((yhat.reg-Carseats.test)^2)
## [1] 4.469011
sqrt(mean((yhat.reg-Carseats.test)^2))
## [1] 2.114003
Initial Thoughts: The test MSE comes out to 4.47 if we take the square root we get 2.11. Meaning that the sales predictions of this model are usually within 2.11 units off from the actual Sales values.
Use cross-validation in order to determine the optimal level of tree complexity. Does pruning the tree improve the test MSE?
# use CV to see whether pruning the tree will improve it
cv.carseats <- cv.tree(tree_fit)
# using the RSS (dev) calculate MSE
mse <- cv.carseats$dev / length(train)
# plot MSE
plot(cv.carseats$size,
mse,
type = 'b',
xlab = "Tree Node Size",
ylab = "MSE",
main = "The Affect of Tree Node Size on MSE")
Initial Thoughts: The optimal level of tree complexity would be 19 nodes. The MSE sits right at 4.526317 for a tree with 18 node levels, slightly larger than the overall test MSE of 4.21. So, in this case pruning the tree does NOT improve MSE.
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.
set.seed(1)
#use bagging on Random Forest
bag.carseats <- randomForest(Sales ~ .,
data = Carseats,
subset = train,
mtry = 10,
importance = TRUE)
# Make predictions on the test set
yhat.bag <- predict(bag.carseats, newdata = test_data)
plot(yhat.bag,
test_data$Sales,
xlab = "Predicted Sales",
ylab = "Actual Sales",
main = "Bagging - Actual Sales vs Predicted Sales",
abline(0,1, col = "red", lwd = 2))
# Obtain test MSE
mean((yhat.bag - test_data$Sales)^2)
## [1] 2.95016
sqrt(mean((yhat.bag - test_data$Sales)^2))
## [1] 1.717603
importance(bag.carseats)
## %IncMSE IncNodePurity
## CompPrice 29.3462503 184.444145
## Income 2.1950934 84.971443
## Advertising 28.4461920 226.470679
## Population -0.7596617 67.714326
## Price 68.7814941 651.996671
## ShelveLoc 71.7824230 628.164064
## Age 21.3597782 166.657371
## Education 1.4885394 55.640049
## Urban -3.4002675 8.580206
## US 5.0762780 12.380845
varImpPlot(bag.carseats)
Initial Thoughts: The test MSE of the bagged model turned out to be 2.95. Smaller than the regression tree model. The predictions also seem close to the prediction line indicating the model is doing a fairly good job predicting Sales. The test predictions in this dataset would be off when predicting sales on average by 1.72 units of Sales.
Importance Plots:
IncMSA - This plot measures which variables are the most importance for making predictions. ShelveLoc would be the most important variale in predicting sales with a very high MSE
IncNodePurity - This plot measures which variables help create the best tree splits. Price would be the most important factor in creating a tree.
Use random forests to analyze this data. What test MSE 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.
fit_control <- trainControl(
method = "cv",
number = 5
#,classProbs = TRUE,
#summaryFunction = twoClassSummary
)
set.seed(42)
rf_fit <- train(
Sales ~.,
data = train_data,
method = "rf",
#metric "ROC",
trControl = fit_control,
tuneGrid = expand.grid(mtry = c(1, 2, 3)),
ntree = 300
)
plot(varImp(rf_fit))
yhat.rf <- predict(rf_fit, newdata = test_data)
plot(yhat.rf,
test_data$Sales,
xlab = "Predicted Sales",
ylab = "Actual Sales",
main = "Random Forest - Actual Sales vs Predicted Sales",
abline(0,1, col = "red", lwd = 2))
# Obtain test MSE
mean((yhat.rf - test_data$Sales)^2)
## [1] 4.065536
sqrt(mean((yhat.rf - test_data$Sales)^2))
## [1] 2.016317
Initial Thoughts: The importance plot shows that
Price and ShelveLocGood are the most important predictor variables in
all the random forests. Notice how the plot split Shelveloc into Medium
and Good because it was a categorical variable. Also,
mtry = 3 was the most effective way to split the trees.
Meaning that before making a split the trees consider 3 variables as
apposed to 1 or 2. The RMSE of 1.76, R^2 of .66 and MAE of 1.41 prove
this. Also, the test MSE of the Randomforest is the largest MSE so far
of 4.07.
#Fit the BART model
x <- Carseats[, -1]
y <- Carseats[, "Sales"]
xtrain <- x[train, ]
ytrain <- y[train]
xtest <- x[-train, ]
ytest <- y[-train]
set.seed(1)
bartfit <- gbart(xtrain, ytrain, x.test = xtest)
# Calculate MSE
yhat.bart <- bartfit$yhat.test.mean
mean((ytest - yhat.bart)^2)
## [1] 1.73499
yhat.bart <- bartfit$yhat.test.mean
plot(yhat.bart,
test_data$Sales,
xlab = "Predicted Sales",
ylab = "Actual Sales",
main = "BART - Actual Sales vs Predicted Sales",
abline(0,1, col = "red", lwd = 2))
Initial Thoughts: The BART model had an MSE of 1.73 which has been the best MSE thus far. See the table below for comparisons.
single_reg<-mean((yhat.reg-Carseats.test)^2)
bag<-mean((yhat.bag - test_data$Sales)^2)
rf<-mean((yhat.rf - test_data$Sales)^2)
bart<-mean((ytest - yhat.bart)^2)
datatable(
data.frame(
Model = c("Single Regression Tree",
"Bagging",
"Random Forest",
"BART"),
MSE = c(single_reg, bag, rf, bart)
),
caption = "Comparison of Tree-Based Model Performance",
rownames = FALSE,
options = list(dom = "t")
)
OJ DataThis problem involves the OJ data set which is part of the ISLR2 package.
Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations.
attach(OJ)
set.seed(1)
# Here I perform the test/train split
oj_train <- createDataPartition(OJ$Purchase, p=0.7, list=FALSE)
oj_train_data <- OJ[oj_train, ]
oj_test_data <- OJ[-oj_train, ]
Fit a tree to the training data, with Purchase as the response and the other variables as predictors. Use the summary() function to produce summary statistics about the tree, and describe the results obtained. What is the training error rate? How many terminal nodes does the tree have?
# create the tree
full_tree <- rpart(Purchase~.,
data=oj_train_data,
method = "class",
control=rpart.control())
summary(full_tree)
## Call:
## rpart(formula = Purchase ~ ., data = oj_train_data, method = "class",
## control = rpart.control())
## n= 750
##
## CP nsplit rel error xerror xstd
## 1 0.48287671 0 1.0000000 1.0000000 0.04573100
## 2 0.02568493 1 0.5171233 0.5479452 0.03842134
## 3 0.02397260 3 0.4657534 0.5650685 0.03885138
## 4 0.02054795 4 0.4417808 0.5513699 0.03850854
## 5 0.01369863 5 0.4212329 0.5342466 0.03806642
## 6 0.01027397 6 0.4075342 0.5102740 0.03742113
## 7 0.01000000 7 0.3972603 0.5239726 0.03779370
##
## Variable importance
## LoyalCH StoreID PriceDiff SalePriceMM PriceMM
## 47 9 8 6 5
## DiscMM PctDiscMM ListPriceDiff PriceCH WeekofPurchase
## 5 5 4 3 3
## SalePriceCH STORE Store7
## 2 2 1
##
## Node number 1: 750 observations, complexity param=0.4828767
## predicted class=CH expected loss=0.3893333 P(node) =1
## class counts: 458 292
## probabilities: 0.611 0.389
## left son=2 (425 obs) right son=3 (325 obs)
## Primary splits:
## LoyalCH < 0.5036 to the right, improve=123.09670, (0 missing)
## StoreID < 3.5 to the right, improve= 32.56409, (0 missing)
## PriceDiff < 0.31 to the right, improve= 19.87539, (0 missing)
## SalePriceMM < 2.04 to the right, improve= 17.34553, (0 missing)
## DiscCH < 0.255 to the right, improve= 15.60231, (0 missing)
## Surrogate splits:
## StoreID < 3.5 to the right, agree=0.648, adj=0.188, (0 split)
## PriceMM < 1.89 to the right, agree=0.596, adj=0.068, (0 split)
## WeekofPurchase < 246.5 to the right, agree=0.595, adj=0.065, (0 split)
## ListPriceDiff < 0.035 to the right, agree=0.591, adj=0.055, (0 split)
## PriceCH < 1.72 to the right, agree=0.588, adj=0.049, (0 split)
##
## Node number 2: 425 observations, complexity param=0.0239726
## predicted class=CH expected loss=0.1388235 P(node) =0.5666667
## class counts: 366 59
## probabilities: 0.861 0.139
## left son=4 (404 obs) right son=5 (21 obs)
## Primary splits:
## PriceDiff < -0.39 to the right, improve=12.310240, (0 missing)
## LoyalCH < 0.7645725 to the right, improve= 9.114571, (0 missing)
## DiscMM < 0.57 to the left, improve= 8.216284, (0 missing)
## PctDiscMM < 0.264375 to the left, improve= 8.216284, (0 missing)
## SalePriceMM < 1.585 to the right, improve= 7.859655, (0 missing)
## Surrogate splits:
## DiscMM < 0.72 to the left, agree=0.976, adj=0.524, (0 split)
## SalePriceMM < 1.435 to the right, agree=0.976, adj=0.524, (0 split)
## PctDiscMM < 0.3342595 to the left, agree=0.976, adj=0.524, (0 split)
## SalePriceCH < 2.075 to the left, agree=0.958, adj=0.143, (0 split)
##
## Node number 3: 325 observations, complexity param=0.02568493
## predicted class=MM expected loss=0.2830769 P(node) =0.4333333
## class counts: 92 233
## probabilities: 0.283 0.717
## left son=6 (181 obs) right son=7 (144 obs)
## Primary splits:
## LoyalCH < 0.282272 to the right, improve=14.082430, (0 missing)
## PriceDiff < 0.31 to the right, improve= 8.237351, (0 missing)
## SpecialCH < 0.5 to the right, improve= 5.851637, (0 missing)
## DiscCH < 0.255 to the right, improve= 5.607657, (0 missing)
## PctDiscCH < 0.132882 to the right, improve= 5.607657, (0 missing)
## Surrogate splits:
## STORE < 2.5 to the left, agree=0.609, adj=0.118, (0 split)
## PriceCH < 1.875 to the left, agree=0.594, adj=0.083, (0 split)
## SalePriceCH < 1.875 to the left, agree=0.594, adj=0.083, (0 split)
## StoreID < 3.5 to the right, agree=0.569, adj=0.028, (0 split)
## PriceMM < 2.205 to the left, agree=0.566, adj=0.021, (0 split)
##
## Node number 4: 404 observations
## predicted class=CH expected loss=0.1113861 P(node) =0.5386667
## class counts: 359 45
## probabilities: 0.889 0.111
##
## Node number 5: 21 observations, complexity param=0.01027397
## predicted class=MM expected loss=0.3333333 P(node) =0.028
## class counts: 7 14
## probabilities: 0.333 0.667
## left son=10 (9 obs) right son=11 (12 obs)
## Primary splits:
## LoyalCH < 0.742157 to the right, improve=3.5000000, (0 missing)
## STORE < 1.5 to the right, improve=3.0476190, (0 missing)
## StoreID < 1.5 to the right, improve=0.7619048, (0 missing)
## Store7 splits as LR, improve=0.7619048, (0 missing)
## PriceCH < 1.975 to the right, improve=0.1696970, (0 missing)
## Surrogate splits:
## PriceCH < 2.04 to the right, agree=0.714, adj=0.333, (0 split)
## DiscMM < 0.47 to the left, agree=0.714, adj=0.333, (0 split)
## SalePriceMM < 1.64 to the right, agree=0.714, adj=0.333, (0 split)
## SalePriceCH < 2.04 to the right, agree=0.714, adj=0.333, (0 split)
## PctDiscMM < 0.2224545 to the left, agree=0.714, adj=0.333, (0 split)
##
## Node number 6: 181 observations, complexity param=0.02568493
## predicted class=MM expected loss=0.4143646 P(node) =0.2413333
## class counts: 75 106
## probabilities: 0.414 0.586
## left son=12 (107 obs) right son=13 (74 obs)
## Primary splits:
## PriceDiff < 0.05 to the right, improve=12.694000, (0 missing)
## SalePriceMM < 1.94 to the right, improve=10.544100, (0 missing)
## DiscMM < 0.22 to the left, improve= 6.766872, (0 missing)
## PctDiscMM < 0.0729725 to the left, improve= 5.775128, (0 missing)
## ListPriceDiff < 0.235 to the right, improve= 4.614858, (0 missing)
## Surrogate splits:
## SalePriceMM < 1.94 to the right, agree=0.950, adj=0.878, (0 split)
## DiscMM < 0.08 to the left, agree=0.834, adj=0.595, (0 split)
## PctDiscMM < 0.038887 to the left, agree=0.834, adj=0.595, (0 split)
## ListPriceDiff < 0.135 to the right, agree=0.779, adj=0.459, (0 split)
## PriceMM < 2.04 to the right, agree=0.773, adj=0.446, (0 split)
##
## Node number 7: 144 observations
## predicted class=MM expected loss=0.1180556 P(node) =0.192
## class counts: 17 127
## probabilities: 0.118 0.882
##
## Node number 10: 9 observations
## predicted class=CH expected loss=0.3333333 P(node) =0.012
## class counts: 6 3
## probabilities: 0.667 0.333
##
## Node number 11: 12 observations
## predicted class=MM expected loss=0.08333333 P(node) =0.016
## class counts: 1 11
## probabilities: 0.083 0.917
##
## Node number 12: 107 observations, complexity param=0.02054795
## predicted class=CH expected loss=0.4299065 P(node) =0.1426667
## class counts: 61 46
## probabilities: 0.570 0.430
## left son=24 (53 obs) right son=25 (54 obs)
## Primary splits:
## STORE < 1.5 to the left, improve=3.442309, (0 missing)
## PriceDiff < 0.49 to the right, improve=1.997351, (0 missing)
## PriceMM < 2.205 to the left, improve=1.879220, (0 missing)
## SalePriceMM < 2.205 to the left, improve=1.879220, (0 missing)
## WeekofPurchase < 249.5 to the right, improve=1.713074, (0 missing)
## Surrogate splits:
## StoreID < 5.5 to the right, agree=0.794, adj=0.585, (0 split)
## Store7 splits as RL, agree=0.794, adj=0.585, (0 split)
## SalePriceCH < 1.775 to the left, agree=0.748, adj=0.491, (0 split)
## PriceMM < 2.155 to the left, agree=0.720, adj=0.434, (0 split)
## PriceCH < 1.875 to the left, agree=0.710, adj=0.415, (0 split)
##
## Node number 13: 74 observations
## predicted class=MM expected loss=0.1891892 P(node) =0.09866667
## class counts: 14 60
## probabilities: 0.189 0.811
##
## Node number 24: 53 observations
## predicted class=CH expected loss=0.3018868 P(node) =0.07066667
## class counts: 37 16
## probabilities: 0.698 0.302
##
## Node number 25: 54 observations, complexity param=0.01369863
## predicted class=MM expected loss=0.4444444 P(node) =0.072
## class counts: 24 30
## probabilities: 0.444 0.556
## left son=50 (18 obs) right son=51 (36 obs)
## Primary splits:
## LoyalCH < 0.49 to the right, improve=1.500000, (0 missing)
## SalePriceCH < 1.775 to the right, improve=1.463019, (0 missing)
## PriceDiff < 0.31 to the left, improve=1.121357, (0 missing)
## WeekofPurchase < 269.5 to the right, improve=1.066667, (0 missing)
## SalePriceMM < 2.15 to the left, improve=1.060652, (0 missing)
## Surrogate splits:
## StoreID < 3.5 to the right, agree=0.704, adj=0.111, (0 split)
## STORE < 3.5 to the right, agree=0.704, adj=0.111, (0 split)
## WeekofPurchase < 252.5 to the left, agree=0.685, adj=0.056, (0 split)
## PriceDiff < 0.43 to the right, agree=0.685, adj=0.056, (0 split)
## ListPriceDiff < 0.38 to the right, agree=0.685, adj=0.056, (0 split)
##
## Node number 50: 18 observations
## predicted class=CH expected loss=0.3888889 P(node) =0.024
## class counts: 11 7
## probabilities: 0.611 0.389
##
## Node number 51: 36 observations
## predicted class=MM expected loss=0.3611111 P(node) =0.048
## class counts: 13 23
## probabilities: 0.361 0.639
Type in the name of the tree object in order to get a detailed text output. Pick one of the terminal nodes, and interpret the information displayed.
Initial Thoughts: The tree contains 9 nodes and has a training error of 14.5%. Meaning that 14.6% of the observations in the training set are misclassified while the other 85.4% are accurate. The top tree section tells us that members with a loyalty score to great than of .48 generally purchase Citril Hill orange juice and members with a loyalty score less than that usually go with Minute Maid.
Create a plot of the tree, and interpret the results.
# visualize the tree
rpart.plot(
full_tree,
type = 2,
extra = 104,
fallen.leaves = TRUE,
box.palette = "Blues",
shadow.col = "gray",
nn = TRUE
)
Initial Thoughts: The following factors tend to influence members in predicting which brand they’re loyal to and purchasing behavior. Members with a loyalty score greater than .76, price difference between CH and MM greater than $0.02, and store ID greater than 4.
Predict the response on the test data, and produce a confusion matrix comparing the test labels to the predicted test labels. What is the test error rate?
# create the tree
pred <- predict(
full_tree,
newdata = oj_test_data,
type = "class"
)
# confusion matrix
confusion <- table(
Actual = oj_test_data$Purchase,
Predicted = pred
)
confusion
## Predicted
## Actual CH MM
## CH 168 27
## MM 28 97
mean(pred != oj_test_data$Purchase)
## [1] 0.171875
Final Thoughts: The confusion matrix shows that the decision tree model predicted well on the test set. With a 13.3% error rate we can see that the vast majority of the predictions were accuately predicted.