STA-6543: Assignment 7
Austin Thrash - FFK221
2024-04-19
Table of Contents:
ISLR Chapter 8 Questions
Question 3: Gini Index
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.
proportion1 <- seq(0, 1, 0.001)
proportion2 <- 1 - proportion1
classification_error <- 1 - pmax(proportion1, proportion2)
gini <- proportion1 * (1-proportion1) + proportion2 * (1-proportion2)
entropy <- -proportion1 * log(proportion1) - proportion2 * log(proportion2)
data.frame(proportion1, proportion2, classification_error, gini, entropy) %>%
pivot_longer(cols = c(classification_error, gini, entropy), names_to = "metric") %>%
ggplot(aes(x = proportion1, y = value, col = factor(metric))) +
geom_line(linewidth = 1) +
scale_y_continuous(breaks = seq(0, 1, 0.1), minor_breaks = NULL) +
scale_color_hue(labels = c("Classification Error", "Entropy", "Gini")) +
labs(col = "Legend",
y = "Value",
x = "Proportion")## Warning: Removed 2 rows containing missing values (`geom_line()`).
Question 8: Regression trees
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.
8.(a) Splitting Data
Split the data set into a training set and a test set.
# Split the data
set.seed(12345)
train <- sample(1:nrow(Carseats), nrow(Carseats) / 2)8.(b) Regression Tree
Fit a regression tree to the training set. Plot the tree, and interpret the results. What test MSE do you obtain?
tree.Carseats <- tree(Sales ~ ., Carseats, subset = train)
summary(tree.Carseats)##
## Regression tree:
## tree(formula = Sales ~ ., data = Carseats, subset = train)
## Variables actually used in tree construction:
## [1] "ShelveLoc" "Price" "CompPrice" "Age" "US"
## [6] "Income" "Population" "Advertising"
## Number of terminal nodes: 19
## Residual mean deviance: 1.841 = 333.2 / 181
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -5.3080 -0.6933 0.1642 0.0000 0.8924 3.4020plot(tree.Carseats)
text(tree.Carseats, pretty = 0, cex = 0.5)
When looking at the tree, we can see where each split lies and at what value. The first split begins at ‘ShelveLoc’, the option of ‘Bad’ will lead you down the left side of the tree while the option ‘Medium’ will lead you down the right side of the tree. This process continues at each split of the tree, then finally arrives at a prediction for ‘Sales’ at the bottom of the tree.
yhat <- predict(tree.Carseats, newdata = Carseats[-train, ])
carseats.test <- Carseats[-train, "Sales"]
plot(yhat, carseats.test)
abline(0, 1)
mean((yhat - carseats.test)^2)## [1] 5.236977Based on the results of the model, it achieves a test MSE of 5.236977
8.(c) Cross-Validation
Use cross-validation in order to determine the optimal level of tree complexity. Does pruning the tree improve the test MSE?
cv.carseats <- cv.tree(tree.Carseats)
plot(cv.carseats$size, cv.carseats$dev, type = "b")
We can see that the optimal point seems to be 9, as it is the simplest and has a low deviance.
prune.carseats <- prune.tree(tree.Carseats, best = 9)
plot(prune.carseats)
text(prune.carseats, pretty = 0)
yhat <- predict(prune.carseats, newdata = Carseats[-train, ])
carseats.test <- Carseats[-train, "Sales"]
plot(yhat, carseats.test)
abline(0, 1)
mean((yhat - carseats.test)^2)## [1] 5.141645Pruning the decision tree down to 16 help decrease the test MSE from 5.236977 to 5.141645. Therefore, pruning the tree did help improve the test MSE.
8.(d) Bagging
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(12345)
bag.carseats <- randomForest(Sales ~ ., data = Carseats,
subset = train, mtry = 10, importance = TRUE)
bag.carseats##
## Call:
## randomForest(formula = Sales ~ ., data = Carseats, mtry = 10, importance = TRUE, subset = train)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 10
##
## Mean of squared residuals: 2.60257
## % Var explained: 63.73importance(bag.carseats)## %IncMSE IncNodePurity
## CompPrice 18.891751393 134.062837
## Income 2.451392651 70.848737
## Advertising 16.651146492 132.833223
## Population -1.441368637 59.141863
## Price 43.309642246 374.502459
## ShelveLoc 50.608908916 383.660180
## Age 19.903193742 192.799402
## Education 0.527296075 42.407971
## Urban -1.768327245 7.261618
## US -0.007330736 5.517387yhat <- predict(bag.carseats, newdata = Carseats[-train, ])
carseats.test <- Carseats[-train, "Sales"]
plot(yhat, carseats.test)
abline(0, 1)
mean((yhat - carseats.test)^2)## [1] 3.24628From the results we can see that ‘CompPrice’, ‘Advertising’, ‘Price’, ‘ShelveLoc’, and ‘Age’ seem to be the most important predictors when attempting to predict ‘Sales’. We can also notice that bagging has helped improve our test MSE and reduced it down to 3.24628.
8.(e) Random Forest
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.
set.seed(12345)
rf.carseats <- randomForest(Sales ~ ., data = Carseats,
subset = train, mtry = 9, importance = TRUE)
importance(rf.carseats)## %IncMSE IncNodePurity
## CompPrice 17.6902322 140.319225
## Income 2.3971098 73.749587
## Advertising 17.9166740 138.691129
## Population -0.9718496 56.918853
## Price 44.3251074 368.017147
## ShelveLoc 48.8718686 373.194651
## Age 18.8405766 187.013211
## Education -1.0019776 44.245674
## Urban -4.0700278 8.101418
## US 0.9250337 6.356382yhat.rf <- predict(rf.carseats, newdata = Carseats[-train, ])
mean((yhat.rf - carseats.test)^2)## [1] 3.235668We can see from the test MSE that Random Forest out performed the bagging method (3.24628 for bagging and 3.235668 for Random Forest), suggesting that Random Forest yields better results in this case. The importance of the predictors does not seem to vary that much from the bagging method. The effect of ‘m’ can play heavily into the accuracy of the model, a low ‘m’ can effect bias and cause the trees to be more correlated with one another, while a high ‘m’ will practically do the inverse effecting variance and lowering correlation between trees.
8.(f) Results
Now analyze the data using BART, and report your results.
set.seed(12345)
boost.carseats <- gbm(Sales ~ ., data = Carseats[train, ],
distribution = "gaussian", n.trees = 5000, interaction.depth = 4)
summary(boost.carseats)
## var rel.inf
## Price Price 24.0558395
## ShelveLoc ShelveLoc 20.0075801
## Age Age 14.2804000
## CompPrice CompPrice 12.1771605
## Advertising Advertising 9.0800697
## Income Income 8.9491531
## Population Population 7.2705051
## Education Education 3.0178871
## US US 0.8942359
## Urban Urban 0.2671691yhat.boost <- predict(boost.carseats, newdata = Carseats[-train, ], n.trees = 5000)
mean((yhat.boost - carseats.test)^2)## [1] 2.553012We can see that by using BART, the most influential predictors are the same ones found when bagging and Random Forest. The test MSE has decreased significantly when comparing BART to any of the other models tested.
Question 9
This problem involves the OJ data set which is part of the ISLR2 package.
9.(a) Splitting Data
Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations.
set.seed(12345)
train_index <- sample(1:nrow(OJ), 800)
train <- OJ[train_index, ]
test <- OJ[-train_index, ]9.(b) Fit a Tree
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?
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"
## Number of terminal nodes: 7
## Residual mean deviance: 0.7724 = 612.5 / 793
## Misclassification error rate: 0.17 = 136 / 800From the results of the classification tree, we can see that only 2 of the predictors were used when creating the splits, that being ‘LoyalCH’, ‘PriceDiff’. This model only has a missclassification error rate of 17% and only 7 terminal nodes.
9.(c) Interpret Tree
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.
oj.tree## node), split, n, deviance, yval, (yprob)
## * denotes terminal node
##
## 1) root 800 1059.000 CH ( 0.62500 0.37500 )
## 2) LoyalCH < 0.5036 349 425.200 MM ( 0.29799 0.70201 )
## 4) LoyalCH < 0.276142 170 138.400 MM ( 0.14118 0.85882 )
## 8) LoyalCH < 0.0356415 55 9.996 MM ( 0.01818 0.98182 ) *
## 9) LoyalCH > 0.0356415 115 115.100 MM ( 0.20000 0.80000 ) *
## 5) LoyalCH > 0.276142 179 246.100 MM ( 0.44693 0.55307 )
## 10) PriceDiff < 0.065 76 80.790 MM ( 0.22368 0.77632 ) *
## 11) PriceDiff > 0.065 103 137.600 CH ( 0.61165 0.38835 ) *
## 3) LoyalCH > 0.5036 451 334.500 CH ( 0.87805 0.12195 )
## 6) LoyalCH < 0.705699 138 164.300 CH ( 0.71739 0.28261 )
## 12) PriceDiff < 0.265 89 120.100 CH ( 0.59551 0.40449 ) *
## 13) PriceDiff > 0.265 49 22.570 CH ( 0.93878 0.06122 ) *
## 7) LoyalCH > 0.705699 313 126.300 CH ( 0.94888 0.05112 ) *Lets look at terminal node 9, which denoted by a “*”
If we start at the root node:
1) root 800 1059.000 CH ( 0.62500 0.37500 ) We can see that there are 800 observations, a deviance of 1059, with an overall prediction of CH(62.5% for CH and 37.5% for MM).
From that node we go down the following path:
LoyalCH = 0.5036
LoyalCH = 0.276142
LoyalCH = 0.0356415
Therefore, node 9 is a subset of purchases where ‘0.0356415 < LoyalCH < 0.276142’ and the overall prediction is MM. This node seems to be pretty good with a split of 80% MM vs 20% CH and 115 observations.
9.(d) Plot Tree
Create a plot of the tree, and interpret the results.
plot(oj.tree)
text(oj.tree, pretty = 0, cex = 0.5)
The plot matches the text output we did prior, both having a starting split at LoyalCH = 0.5036. We can see the node we selected in the bottom left corner of the plot, that node has the same prediction of ‘MM’
9.(e) Predictions
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?
predictions <- predict(oj.tree, test, type = "class")
table(predictions, test_labels = test$Purchase)## test_labels
## predictions CH MM
## CH 145 47
## MM 8 701 - mean(predictions == test$Purchase)## [1] 0.2037037The model achieves a test error rate of 20.37%
9.(f) cv.tree()
Apply the cv.tree() function to the training set in order to determine the optimal tree size.
cv.oj <- cv.tree(oj.tree, K = 10, FUN = prune.misclass)
cv.oj## $size
## [1] 7 4 2 1
##
## $dev
## [1] 160 160 169 300
##
## $k
## [1] -Inf 0.0 11.5 141.0
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"9.(g) Tree / Cross-Validation
Produce a plot with tree size on the x-axis and cross-validated classification error rate on the y-axis.
plot(cv.oj$size, cv.oj$dev / nrow(train), type = "b")
9.(h) Lowest Cross-Validated Classification Error Rate
Which tree size corresponds to the lowest cross-validated classification error rate?
From the plot, it seems as if a tree size of 4 corresponds to the lowest cross-validated classification error rate. So moving forward we will use a tree size of 4.
9.(i) Pruned Tree
Produce a pruned tree corresponding to the optimal tree size obtained using cross-validation. If cross-validation does not lead to selection of a pruned tree, then create a pruned tree with five terminal nodes.
pruned_oj <- prune.tree(oj.tree, best = 4)
pruned_oj## node), split, n, deviance, yval, (yprob)
## * denotes terminal node
##
## 1) root 800 1059.0 CH ( 0.62500 0.37500 )
## 2) LoyalCH < 0.5036 349 425.2 MM ( 0.29799 0.70201 )
## 4) LoyalCH < 0.276142 170 138.4 MM ( 0.14118 0.85882 ) *
## 5) LoyalCH > 0.276142 179 246.1 MM ( 0.44693 0.55307 ) *
## 3) LoyalCH > 0.5036 451 334.5 CH ( 0.87805 0.12195 )
## 6) LoyalCH < 0.705699 138 164.3 CH ( 0.71739 0.28261 ) *
## 7) LoyalCH > 0.705699 313 126.3 CH ( 0.94888 0.05112 ) *9.(j) Compare Training Error Rates
Compare the training error rates between the pruned and unpruned trees. Which is higher?
mean(predict(oj.tree, type = "class", newdata = train) != train$Purchase)## [1] 0.17mean(predict(pruned_oj, type = "class", newdata = train) != train$Purchase)## [1] 0.19875Our pruned tree has higher training error rates than the unpruned one. The unpruned tree seems to be out performing the pruned one.
9.(k) Compare Test Error Rates
Compare the test error rates between the pruned and unpruned trees. Which is higher?
mean(predict(oj.tree, type = "class", newdata = test) != test$Purchase)## [1] 0.2037037mean(predict(pruned_oj, type = "class", newdata = test) != test$Purchase)## [1] 0.2037037However, when obtaining the test error rates, both the pruned and unpruned trees seem to achieve the same value. Furthermore, changing the size of the trees has a negative impact on the accuracy of the pruned model.