Assignment 7
Brett Stokes
4/22/2021
Problem 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 ˆpm1. The x-axis should display ˆpm1, 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, ˆpm1 = 1− ˆpm2. You could make this plot by hand, but it will be much easier to make in R.
p <- seq(0, 1, 0.01)
gini.index <- 2 * p * (1 - p)
class.error <- 1 - pmax(p, 1 - p)
cross.entropy <- - (p * log(p) + (1 - p) * log(1 - p))
matplot(p, cbind(gini.index, class.error, cross.entropy), col = c("#66C2A5" , "#FC8D62", "#E78AC3"))
Problem 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(ISLR)
set.seed(100)
data1=Carseats
index = sample(1:nrow(data1), nrow(data1)*0.8)
train = data1[index, ]
test = data1[-index, ](b) Fit a regression tree to the training set. Plot the tree, and interpret the results. What test MSE do you obtain?
library(tree)## Warning: package 'tree' was built under R version 4.0.4set.seed(100)
model <- tree(Sales ~ ., data = train)
summary(model)
plot(model)
text(model, pretty = 0)
preds=predict(model, test)
MSE=mean((preds-train$Sales)^2)
MSE##
## Regression tree:
## tree(formula = Sales ~ ., data = train)
## Variables actually used in tree construction:
## [1] "ShelveLoc" "Price" "CompPrice" "Age" "Advertising"
## Number of terminal nodes: 16
## Residual mean deviance: 2.564 = 779.4 / 304
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -4.56500 -1.09800 0.05679 0.00000 1.07600 4.07600
## [1] 12.19055Our tree uses the variables ShelveLoc, Price, CompPrice, Age and Advertising. ShelveLoc seems to be the most important variable based on minimal depth. Our test MSE is 12.19055.
(c) Use cross-validation in order to determine the optimal level of tree complexity. Does pruning the tree improve the test MSE?
set.seed(100)
cv <- cv.tree(model)
plot(cv$size, cv$dev, type = "b")
model.prune <- prune.tree(model, best = 6)
plot(model.prune)
text(model.prune, pretty = 0)
preds=predict(model.prune, test)
MSE=mean((preds-train$Sales)^2)
MSE## [1] 12.4532Based on the above graph generated from cross validation, our optimal tree size would include 6 variables. We retrain our tree model to include 6 variables. We notice that our test MSE is now 12.4532 which is higher than our original tree. Clearly pruning will not help us in this case.
(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.
set.seed(100)
model.bag <- randomForest(Sales ~ ., data = train, mtry = 10, ntree = 1000, importance = TRUE)
importance(model.bag)## %IncMSE IncNodePurity
## CompPrice 52.68420105 287.518558
## Income 13.70800224 137.793247
## Advertising 31.69155182 187.392967
## Population -0.01069956 76.549416
## Price 105.75950970 694.229625
## ShelveLoc 110.47066501 845.653586
## Age 34.32548379 248.312910
## Education 3.44852182 60.959720
## Urban -2.27891048 9.333176
## US 6.10311355 12.883038preds <- predict(model.bag, test)
MSE=mean((preds - test$Sales)^2)
MSE## [1] 2.567417We note that Price and ShelveLoc are the most important variables based on our bagged results. Also, our test MSE has decreased dramatically to 2.567417.
(e) 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(100)
mtry.values = seq(1,10,1)
hyper_grid = expand.grid(mtry = mtry.values)
test.MSE = c()
for (i in 1:nrow(hyper_grid)) {
model.rf = randomForest(Sales ~ ., data = train, mtry = hyper_grid$mtry[i], ntree = 1000, importance = TRUE)
preds=predict(model.rf,test)
test.MSE[i]=mean((preds-test$Sales)^2)
}
mse.df = cbind(m_number = mtry.values, testMSE = test.MSE)
mse.df## m_number testMSE
## [1,] 1 4.676336
## [2,] 2 3.406331
## [3,] 3 2.954867
## [4,] 4 2.743929
## [5,] 5 2.585772
## [6,] 6 2.555366
## [7,] 7 2.512130
## [8,] 8 2.534927
## [9,] 9 2.506986
## [10,] 10 2.511641As we can see from the above output, generally we are getting better test MSE values as we increase m. The best test MSE score occurs at m = 9 variables.
Problem 9
This problem involves the “OJ” data set which is part of the “ISLR” package.
(a) Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations.
set.seed(100)
data2 = OJ
index2 <- sample(1:nrow(data2), 800)
train2 <- data2[index2, ]
test2 <- data2[-index2, ](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?
set.seed(100)
model2 <- tree(Purchase ~ ., data = train2)
summary(model2)##
## Classification tree:
## tree(formula = Purchase ~ ., data = train2)
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff" "ListPriceDiff" "SalePriceMM"
## Number of terminal nodes: 8
## Residual mean deviance: 0.7235 = 573 / 792
## Misclassification error rate: 0.1588 = 127 / 800For the above tree: the training error rate is 0.1588 and the number of terminal nodes is 8.
(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.
model2## node), split, n, deviance, yval, (yprob)
## * denotes terminal node
##
## 1) root 800 1070.00 CH ( 0.61000 0.39000 )
## 2) LoyalCH < 0.5036 355 420.20 MM ( 0.27887 0.72113 )
## 4) LoyalCH < 0.280875 167 130.20 MM ( 0.13174 0.86826 )
## 8) LoyalCH < 0.0356415 59 10.14 MM ( 0.01695 0.98305 ) *
## 9) LoyalCH > 0.0356415 108 106.40 MM ( 0.19444 0.80556 ) *
## 5) LoyalCH > 0.280875 188 254.40 MM ( 0.40957 0.59043 )
## 10) PriceDiff < 0.05 81 74.58 MM ( 0.17284 0.82716 ) *
## 11) PriceDiff > 0.05 107 144.90 CH ( 0.58879 0.41121 ) *
## 3) LoyalCH > 0.5036 445 336.80 CH ( 0.87416 0.12584 )
## 6) LoyalCH < 0.764572 180 204.60 CH ( 0.74444 0.25556 )
## 12) ListPriceDiff < 0.18 49 66.27 MM ( 0.40816 0.59184 ) *
## 13) ListPriceDiff > 0.18 131 101.10 CH ( 0.87023 0.12977 )
## 26) SalePriceMM < 2.04 51 59.94 CH ( 0.72549 0.27451 ) *
## 27) SalePriceMM > 2.04 80 25.59 CH ( 0.96250 0.03750 ) *
## 7) LoyalCH > 0.764572 265 85.16 CH ( 0.96226 0.03774 ) *We choose to interpret terminal node (9). At this point the split criteria is LoyalCH greater than 0.0356415. This branch has 108 observations. Most of the observations in this branch take the value of MM.
(d) Create a plot of the tree, and interpret the results.
plot(model2)
text(model2, pretty = 0)
We can see that the most important predictor is LoyalCH from the above plot.
(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?
preds <- predict(model2, test2, type = "class")
caret::confusionMatrix(as.factor(preds),as.factor(test2$Purchase))## Confusion Matrix and Statistics
##
## Reference
## Prediction CH MM
## CH 142 36
## MM 23 69
##
## Accuracy : 0.7815
## 95% CI : (0.7274, 0.8293)
## No Information Rate : 0.6111
## P-Value [Acc > NIR] : 1.771e-09
##
## Kappa : 0.5297
##
## Mcnemar's Test P-Value : 0.1182
##
## Sensitivity : 0.8606
## Specificity : 0.6571
## Pos Pred Value : 0.7978
## Neg Pred Value : 0.7500
## Prevalence : 0.6111
## Detection Rate : 0.5259
## Detection Prevalence : 0.6593
## Balanced Accuracy : 0.7589
##
## 'Positive' Class : CH
## 1-0.7815## [1] 0.2185From the above table, we note that the test error rate is 0.2185.
(f) Apply the cv.tree() function to the training set in order to determine the optimal tree size.
model2.cv <- cv.tree(model2, FUN = prune.misclass)
model2.cv## $size
## [1] 8 6 4 2 1
##
## $dev
## [1] 141 141 148 161 312
##
## $k
## [1] -Inf 0.0 4.5 9.5 157.0
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"(g) Produce a plot with tree size on the x-axis and cross-validated classification error rate on the y-axis.
plot(model2.cv$size, model2.cv$dev, type = "b", xlab = "Tree size", ylab = "Deviance")
(h) Which tree size corresponds to the lowest cross-validated classification error rate?
According to the above graph, I would say that a 6 node tree would be optimal for this particular model.
(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.
model2.prune <- prune.misclass(model2, best = 6)
plot(model2.prune)
text(model2.prune, pretty = 0)
(j) Compare the training error rates between the pruned and unpruned trees. Which is higher?
summary(model2)
summary(model2.prune)##
## Classification tree:
## tree(formula = Purchase ~ ., data = train2)
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff" "ListPriceDiff" "SalePriceMM"
## Number of terminal nodes: 8
## Residual mean deviance: 0.7235 = 573 / 792
## Misclassification error rate: 0.1588 = 127 / 800
##
## Classification tree:
## snip.tree(tree = model2, nodes = c(4L, 13L))
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff" "ListPriceDiff"
## Number of terminal nodes: 6
## Residual mean deviance: 0.7585 = 602.2 / 794
## Misclassification error rate: 0.1588 = 127 / 800The training error is exactly the same in my case for both the pruned and unpruned tree.
(k) Compare the test error rates between the pruned and unpruned trees. Which is higher?
preds <- predict(model2, test2, type = "class")
caret::confusionMatrix(as.factor(preds),as.factor(test2$Purchase))## Confusion Matrix and Statistics
##
## Reference
## Prediction CH MM
## CH 142 36
## MM 23 69
##
## Accuracy : 0.7815
## 95% CI : (0.7274, 0.8293)
## No Information Rate : 0.6111
## P-Value [Acc > NIR] : 1.771e-09
##
## Kappa : 0.5297
##
## Mcnemar's Test P-Value : 0.1182
##
## Sensitivity : 0.8606
## Specificity : 0.6571
## Pos Pred Value : 0.7978
## Neg Pred Value : 0.7500
## Prevalence : 0.6111
## Detection Rate : 0.5259
## Detection Prevalence : 0.6593
## Balanced Accuracy : 0.7589
##
## 'Positive' Class : CH
## 1-0.7815## [1] 0.2185preds=predict(model2.prune,test2, type = "class")
caret::confusionMatrix(as.factor(preds),as.factor(test2$Purchase))## Confusion Matrix and Statistics
##
## Reference
## Prediction CH MM
## CH 142 36
## MM 23 69
##
## Accuracy : 0.7815
## 95% CI : (0.7274, 0.8293)
## No Information Rate : 0.6111
## P-Value [Acc > NIR] : 1.771e-09
##
## Kappa : 0.5297
##
## Mcnemar's Test P-Value : 0.1182
##
## Sensitivity : 0.8606
## Specificity : 0.6571
## Pos Pred Value : 0.7978
## Neg Pred Value : 0.7500
## Prevalence : 0.6111
## Detection Rate : 0.5259
## Detection Prevalence : 0.6593
## Balanced Accuracy : 0.7589
##
## 'Positive' Class : CH
## 1-0.7815## [1] 0.2185For my case the test error for the pruned and unpruned tree is exactly the same: 0.2185. Therefore, the pruned model has obtained the same amount of accuracy with fewer variables.