hw7
Bella Lin
2024-04-19
3. 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.
p <- seq(0, 1, 0.01)
gini <- p * (1 - p) * 2
entropy <- -(p * log(p) + (1 - p) * log(1 - p))
class.err <- 1 - pmax(p, 1 - p)
matplot(p, cbind(gini, entropy, class.err), col = c("pink", "skyblue", "mediumpurple"), type = "l",
xlab = "p̂ m1", ylab = "Value",
main = "Comparison of Gini Index, Entropy, and Classification Error")
legend("topright", legend = c("Gini Index", "Entropy", "Classification Error"), col = c("pink", "skyblue", "mediumpurple"), lty = 1)
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.
set.seed(123)
index_car <- sample(1:nrow(Carseats), 0.7 * nrow(Carseats))
train_car <- Carseats[index_car, ]
test_car <- Carseats[-index_car, ](b) Fit a regression tree to the training set. Plot the tree, and interpret the results. What test MSE do you obtain?
tree_model <- tree(Sales ~ ., data = train_car)
plot(tree_model)
text(tree_model, pretty = 0)
test_predictions <- predict(tree_model, newdata = test_car)test_mse <- mean((test_predictions - test_car$Sales)^2)
print(paste("Test MSE:", test_mse))## [1] "Test MSE: 3.60281802993858"The test MSE is 3.60.
(c) Use cross-validation in order to determine the optimal level of tree complexity. Does pruning the tree improve the test MSE?
set.seed(123)
cv.seats = cv.tree(tree_model)
plot(cv.seats$size, cv.seats$dev, type = "b")
optimal_tree_size <- cv.seats$size[which.min(cv.seats$dev)]
pruned_tree <- prune.tree(tree_model, best = optimal_tree_size)
plot(pruned_tree)
text(pruned_tree,pretty=0)
pruned_predictions = predict(pruned_tree, newdata = test_car)
pruned_test_mse <- mean((pruned_predictions - test_car$Sales)^2)
pruned_test_mse## [1] 4.353022According to the plot, we can see that the optimal level of tree complexity is 8. However, it did not improve the test MSE, MSE increase from 3.6 to 4.4.
(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.
bag_model <- randomForest(Sales ~ ., data = train_car, importance = TRUE)
bag_predictions <- predict(bag_model, newdata = test_car)
bag_test_mse <- mean((bag_predictions - test_car$Sales)^2)
bag_test_mse## [1] 2.727296var_importance <- importance(bag_model)
var_importance## %IncMSE IncNodePurity
## CompPrice 17.911554 239.39177
## Income 4.528663 164.30783
## Advertising 17.772111 181.38605
## Population -1.050638 125.99797
## Price 37.900069 478.57325
## ShelveLoc 47.483501 515.98209
## Age 17.369615 267.45670
## Education 3.283913 99.45872
## Urban 0.327023 20.45437
## US 4.037103 23.19051Using the bagging approach, we can see that the test MSE decrease from 3.6 to 2.7. Based on the variable importance measures: ShelveLoc, Price, and CompPrice are the most important predictors for predicting Sales in the random forest model, as they contribute the most to the reduction in mean squared error.
(e) Use random forests to analyze this data. What test MSE do you obtain? Use the importance() function to determine which vari- ables are most important. Describe the effect of m, the number of variables considered at each split, on the error rate obtained.
rf_model <- randomForest(Sales ~ ., data = train_car, mtry = 5, importance = TRUE)
rf_predictions <- predict(rf_model, newdata = test_car)
rf_test_mse <- mean((rf_predictions - test_car$Sales)^2)
rf_test_mse## [1] 2.357816rf_var_importance <- importance(rf_model)
rf_var_importance## %IncMSE IncNodePurity
## CompPrice 24.446299 245.02450
## Income 7.539530 138.70372
## Advertising 18.156128 175.24115
## Population -1.431855 92.96266
## Price 48.101500 531.42530
## ShelveLoc 58.605672 600.70245
## Age 19.370217 244.08801
## Education 3.042024 85.96575
## Urban 1.215012 13.42858
## US 2.431833 13.77914In this approach, the test MSE decrease from 2.7 to 2.35. Based on the variable importance measures for the random forest model: ShelveLoc, Price, and CompPrice are still the most important predictors for predicting Sales.
(f) Now analyze the data using BART, and report your results.
set.seed(123)
bart_model <- gbm(Sales ~ ., data = train_car, distribution = "gaussian", n.trees = 5000, interaction.depth = 4)
summary(bart_model)
## var rel.inf
## Price Price 26.5039858
## ShelveLoc ShelveLoc 23.0208745
## CompPrice CompPrice 14.1468503
## Age Age 11.0239146
## Advertising Advertising 8.7706646
## Income Income 7.8554956
## Population Population 5.9117306
## Education Education 2.1469674
## US US 0.3604688
## Urban Urban 0.2590478bart_predictions <- predict(bart_model, newdata = test_car, n.trees = bart_model$ntrees)## Using 5000 trees...bart_test_mse <- mean((bart_predictions - test_car$Sales)^2)
bart_test_mse## [1] 2.319632Lastly, using the BART approach, I acheieved the lowest test MSE of 2.32. According to the summary result, Price and ShelveLoc are the most important predictors in the BART model, with approximately 26.5% and 23.0% relative importance.
9. This problem involves the OJ data set which is part of the ISLR2 package.
(a) Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations.
set.seed(123)
index_OJ <- sample(1:nrow(OJ), 800)
train_OJ <- OJ[index_OJ, ]
test_OJ <- OJ[-index_OJ, ](b) 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_OJ)
summary(OJ_tree)##
## Classification tree:
## tree(formula = Purchase ~ ., data = train_OJ)
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff"
## Number of terminal nodes: 8
## Residual mean deviance: 0.7625 = 603.9 / 792
## Misclassification error rate: 0.165 = 132 / 800The summary indicates which variables were used in constructing the tree are LoyalCH and PriceDiff. The tree has 8 terminal nodes, the training error rate is 0.165, indicates high performance of this model.
(c) 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 1071.00 CH ( 0.60875 0.39125 )
## 2) LoyalCH < 0.5036 350 415.10 MM ( 0.28000 0.72000 )
## 4) LoyalCH < 0.276142 170 131.00 MM ( 0.12941 0.87059 )
## 8) LoyalCH < 0.0356415 56 10.03 MM ( 0.01786 0.98214 ) *
## 9) LoyalCH > 0.0356415 114 108.90 MM ( 0.18421 0.81579 ) *
## 5) LoyalCH > 0.276142 180 245.20 MM ( 0.42222 0.57778 )
## 10) PriceDiff < 0.05 74 74.61 MM ( 0.20270 0.79730 ) *
## 11) PriceDiff > 0.05 106 144.50 CH ( 0.57547 0.42453 ) *
## 3) LoyalCH > 0.5036 450 357.10 CH ( 0.86444 0.13556 )
## 6) PriceDiff < -0.39 27 32.82 MM ( 0.29630 0.70370 ) *
## 7) PriceDiff > -0.39 423 273.70 CH ( 0.90071 0.09929 )
## 14) LoyalCH < 0.705326 130 135.50 CH ( 0.78462 0.21538 )
## 28) PriceDiff < 0.145 43 58.47 CH ( 0.58140 0.41860 ) *
## 29) PriceDiff > 0.145 87 62.07 CH ( 0.88506 0.11494 ) *
## 15) LoyalCH > 0.705326 293 112.50 CH ( 0.95222 0.04778 ) *Looking at the 8th node, this terminal node identifies a group of customers with low loyalty to the Citrus Hill brand (LoyalCH < 0.0356415) and predicts that they are highly likely to purchase the Minute Maid (MM) brand (probability of approximately 98.2%).
(d) Create a plot of the tree, and interpret the results.
plot(OJ_tree)
text(OJ_tree, pretty = 0)
> According to the plot above, we can see there are 8 nodes. LoyalCH,
PriceDiff are the most important variables.
(e) 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?
test_predictions <- predict(OJ_tree, newdata = test_OJ, type = "class")
conf_matrix <- confusionMatrix(test_predictions, test_OJ$Purchase)
print(conf_matrix)## Confusion Matrix and Statistics
##
## Reference
## Prediction CH MM
## CH 150 34
## MM 16 70
##
## Accuracy : 0.8148
## 95% CI : (0.7633, 0.8593)
## No Information Rate : 0.6148
## P-Value [Acc > NIR] : 9.56e-13
##
## Kappa : 0.596
##
## Mcnemar's Test P-Value : 0.01621
##
## Sensitivity : 0.9036
## Specificity : 0.6731
## Pos Pred Value : 0.8152
## Neg Pred Value : 0.8140
## Prevalence : 0.6148
## Detection Rate : 0.5556
## Detection Prevalence : 0.6815
## Balanced Accuracy : 0.7883
##
## 'Positive' Class : CH
## test_error_rate <- 1 - conf_matrix$overall['Accuracy']
test_error_rate## Accuracy
## 0.1851852The overall accuracy of the model is approximately 81.48% with high true positive rate and low true negative rate. The test error rate, calculated as 1 - Accuracy, is 18.52%
(f) Apply the cv.tree() function to the training set in order to determine the optimal tree size.
cv_OJ <- cv.tree(OJ_tree)
cv_OJ## $size
## [1] 8 7 6 5 4 3 2 1
##
## $dev
## [1] 690.2346 692.8691 683.0423 717.4446 717.4446 743.3604 793.1266
## [8] 1073.0916
##
## $k
## [1] -Inf 12.03823 14.92474 25.76707 26.02613 38.91686 50.61655
## [8] 298.68751
##
## $method
## [1] "deviance"
##
## attr(,"class")
## [1] "prune" "tree.sequence"The tree sizes tested range from 1 to 8 terminal nodes. The corresponding cross-validated error rates range from approximately 707 to 1071.78. The tree size with the lowest cross-validated error rate is 6 nodes with error rate of 707.
(g) 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, type = "b", xlab = "Tree Size", ylab = "Cross-Validated Classification Error Rate")
(h) Which tree size corresponds to the lowest cross-validated classification error rate?
According to the plot created in (g), tree size of 6 corresponds to the lowest cross-validated classification error rate, which match the result from question (f).
(i) 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.
prune_OJ = prune.tree(OJ_tree, best = 6)(j) Compare the training error rates between the pruned and unpruned trees. Which is higher?
summary(prune_OJ)##
## Classification tree:
## snip.tree(tree = OJ_tree, nodes = c(4L, 14L))
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff"
## Number of terminal nodes: 6
## Residual mean deviance: 0.7945 = 630.9 / 794
## Misclassification error rate: 0.165 = 132 / 800summary(OJ_tree)##
## Classification tree:
## tree(formula = Purchase ~ ., data = train_OJ)
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff"
## Number of terminal nodes: 8
## Residual mean deviance: 0.7625 = 603.9 / 792
## Misclassification error rate: 0.165 = 132 / 800The training error rate for the pruned tree with 6 terminal nodes is 0.165, which is the same as the training error rate for the unpruned tree with 8 terminal nodes. There is no difference in the training error rates between the pruned and unpruned trees.
(k) Compare the test error rates between the pruned and unpruned trees. Which is higher?
pruned_pred <- predict(prune_OJ, newdata = test_OJ, type = "class")
pruned_conf <- confusionMatrix(pruned_pred, test_OJ$Purchase)
pruned_test_rate <- 1 - pruned_conf$overall['Accuracy']
pruned_test_rate## Accuracy
## 0.1851852Both the pruned and unpruned trees have test error rates of 0.185.