Assignment 7
Rex K.
7/23/2020
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 xaxis should display ˆpm1, ranging from 0 to 1, and the y-axis should display the value of the Gini index, classification error, and entropy.
x <- seq(0, 1, 0.01)
gini <- x * (1 - x) * 2
entropy <- -(x * log(x) + (1 - x) * log(1 - x))
class.err <- 1 - pmax(x, 1 - x)
matplot(x, cbind(gini, entropy, class.err), col = c("blue", "dark green", "black"))
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.
library(ISLR)
attach(Carseats)
set.seed(42)
str(Carseats)## 'data.frame': 400 obs. of 11 variables:
## $ Sales : num 9.5 11.22 10.06 7.4 4.15 ...
## $ CompPrice : num 138 111 113 117 141 124 115 136 132 132 ...
## $ Income : num 73 48 35 100 64 113 105 81 110 113 ...
## $ Advertising: num 11 16 10 4 3 13 0 15 0 0 ...
## $ Population : num 276 260 269 466 340 501 45 425 108 131 ...
## $ Price : num 120 83 80 97 128 72 108 120 124 124 ...
## $ ShelveLoc : Factor w/ 3 levels "Bad","Good","Medium": 1 2 3 3 1 1 3 2 3 3 ...
## $ Age : num 42 65 59 55 38 78 71 67 76 76 ...
## $ Education : num 17 10 12 14 13 16 15 10 10 17 ...
## $ Urban : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 2 2 1 1 ...
## $ US : Factor w/ 2 levels "No","Yes": 2 2 2 2 1 2 1 2 1 2 ...(a) Split the data set into a training set and a test set.
train = sample(dim(Carseats)[1], dim(Carseats)[1]/2)
Carseats.train = Carseats[train, ]
Carseats.test = Carseats[-train, ](b) Fit a regression tree to the training set. Plot the tree, and interpret the results. What test MSE do you obtain?
library(tree)
tree.carseats = tree(Sales ~ ., data = Carseats.train)
summary(tree.carseats)##
## Regression tree:
## tree(formula = Sales ~ ., data = Carseats.train)
## Variables actually used in tree construction:
## [1] "ShelveLoc" "Price" "Age" "Income" "Advertising"
## [6] "CompPrice" "Population" "Urban"
## Number of terminal nodes: 18
## Residual mean deviance: 2.266 = 412.4 / 182
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.46100 -1.01400 -0.09829 0.00000 1.06900 3.68600plot(tree.carseats)
text(tree.carseats, pretty = 0)
pred.carseats = predict(tree.carseats, Carseats.test)
mean((Carseats.test$Sales - pred.carseats)^2)## [1] 5.686401Test MSE for the first version is approximately \(5.6\).
(c) 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, FUN = prune.tree)
par(mfrow = c(1, 2))
plot(cv.carseats$size, cv.carseats$dev, type = "b")
plot(cv.carseats$k, cv.carseats$dev, type = "b")
It looks like the best size is about \(9\).
pruned.carseats = prune.tree(tree.carseats, best = 9)
par(mfrow = c(1, 1))
plot(pruned.carseats)
text(pruned.carseats, pretty = 0)
pred.pruned = predict(pruned.carseats, Carseats.test)
mean((Carseats.test$Sales - pred.pruned)^2)## [1] 5.271184Test MSE for this run is approximately \(5.2\)
(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.
library(randomForest)
bag.carseats = randomForest(Sales ~ ., data = Carseats.train, mtry = 10, ntree = 500,
importance = T)
bag.pred = predict(bag.carseats, Carseats.test)
mean((Carseats.test$Sales - bag.pred)^2)## [1] 2.360426Test MSE for this bagging run improved to approximately \(2.3\)
importance(bag.carseats)## %IncMSE IncNodePurity
## CompPrice 30.4417501 197.998343
## Income 9.2101020 104.302860
## Advertising 12.7498522 104.098442
## Population 2.4997582 62.703767
## Price 54.8337819 432.196870
## ShelveLoc 56.8325702 429.294512
## Age 10.4795981 133.921947
## Education -0.5607019 42.616731
## Urban 1.8401309 10.026923
## US 0.8045295 5.682508Price, ShelveLoc, and CompPrice are three most important predictors of Sale in this run.
(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.
rf.carseats = randomForest(Sales ~ ., data = Carseats.train, mtry = 15, ntree = 500,
importance = T)
rf.pred = predict(rf.carseats, Carseats.test)
mean((Carseats.test$Sales - rf.pred)^2)## [1] 2.400509importance(rf.carseats)## %IncMSE IncNodePurity
## CompPrice 31.1078667 203.322883
## Income 11.3993309 101.974302
## Advertising 14.1706506 108.355426
## Population 4.1880633 60.600259
## Price 51.2163567 433.626565
## ShelveLoc 58.7354903 421.414600
## Age 13.2550638 138.930890
## Education 1.0897721 41.170417
## Urban 2.0183686 9.737253
## US 0.5886619 4.957915Price, ShelveLoc, and CompPrice remain the three most important predictors of Sale but the MSE decrease slightly to approximately \(2.3\). Changing \(m\) varies test MSE between \(2.3\) and \(2.45\).
Version two of this comes from the Lab notes.
cs <- Carseats # creating another data set to use like the example show in the lab notes.
High=ifelse(cs$Sales<=8,"No","Yes")
cs1=data.frame(cs,High)(a) Split the data set into a training set and a test set.
train2 <- sample(1:nrow(cs1), 200)
Carseats.test <- cs1[-train2,]
High.test=High[-train2](b) Fit a regression tree to the training set. Plot the tree, and interpret the results. What test MSE do you obtain?
tree.carseats <- tree(High~.-Sales, cs1, subset=train2)
summary(tree.carseats)##
## Classification tree:
## tree(formula = High ~ . - Sales, data = cs1, subset = train2)
## Variables actually used in tree construction:
## [1] "ShelveLoc" "Price" "Income" "CompPrice" "Advertising"
## [6] "Age" "Population" "Education"
## Number of terminal nodes: 20
## Residual mean deviance: 0.4016 = 72.28 / 180
## Misclassification error rate: 0.095 = 19 / 200plot(tree.carseats)
text(tree.carseats,pretty=0)
tree.pred=predict(tree.carseats,Carseats.test,type="class")
table(tree.pred,High.test)## High.test
## tree.pred No Yes
## No 96 21
## Yes 29 54tp <- 91
tn <- 64
fp <- 20
fn <- 25
print("Accuracy")## [1] "Accuracy"acc <- (tp+tn)/(tp+tn+fp+fn)
print(acc)## [1] 0.775print("Precision")## [1] "Precision"prec <- tp/(tp+tn)
print(prec)## [1] 0.5870968print("Recall")## [1] "Recall"rec <- tp/(tp+fn)
print(rec)## [1] 0.7844828print("F-Score")## [1] "F-Score"fsc <- (2*prec*rec)/(prec+rec)
print(fsc)## [1] 0.6715867(c) Use cross-validation in order to determine the optimal level of tree complexity. Does pruning the tree improve the test MSE?
cs_cv <-cv.tree(tree.carseats, FUN = prune.misclass)
par(mfrow = c(1, 2))
plot(cs_cv$size, cs_cv$dev, type = "b")
plot(cs_cv$k, cs_cv$dev, type = "b")
It looks like the best size is about \(8\).
prune.carseats=prune.misclass(tree.carseats,best=8)
par(mfrow = c(1, 1))
plot(prune.carseats)
text(prune.carseats, pretty = 0)
tree.pred=predict(prune.carseats,Carseats.test,type="class")
table(tree.pred,High.test)## High.test
## tree.pred No Yes
## No 101 26
## Yes 24 49tp <- 101
tn <- 43
fp <- 24
fn <- 26
print("Accuracy")## [1] "Accuracy"acc <- (tp+tn)/(tp+tn+fp+fn)
print(acc)## [1] 0.742268print("Precision")## [1] "Precision"prec <- tp/(tp+tn)
print(prec)## [1] 0.7013889print("Recall")## [1] "Recall"rec <- tp/(tp+fn)
print(rec)## [1] 0.7952756print("F-Score")## [1] "F-Score"fsc <- (2*prec*rec)/(prec+rec)
print(fsc)## [1] 0.7453875This did marginally worse than initial results.
(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.carseats <- randomForest(High~.-Sales, cs1, subset=train2, mtry = 13, ntree = 500, importance = T)
bag.carseats##
## Call:
## randomForest(formula = High ~ . - Sales, data = cs1, mtry = 13, ntree = 500, importance = T, subset = train2)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 10
##
## OOB estimate of error rate: 20.5%
## Confusion matrix:
## No Yes class.error
## No 94 17 0.1531532
## Yes 24 65 0.2696629yhat.bag <- predict(bag.carseats,newdata=cs1[-train2,])
table(yhat.bag,High.test)## High.test
## yhat.bag No Yes
## No 107 22
## Yes 18 53tp <- 107
tn <- 53
fp <- 18
fn <- 22
print("Accuracy")## [1] "Accuracy"acc <- (tp+tn)/(tp+tn+fp+fn)
print(acc)## [1] 0.8print("Precision")## [1] "Precision"prec <- tp/(tp+tn)
print(prec)## [1] 0.66875print("Recall")## [1] "Recall"rec <- tp/(tp+fn)
print(rec)## [1] 0.8294574print("F-Score")## [1] "F-Score"fsc <- (2*prec*rec)/(prec+rec)
print(fsc)## [1] 0.7404844Looks like it improves and that we’re trading Precision and Recall here.
(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.
bag.carseats <- randomForest(High~.-Sales, cs1, subset=train2, mtry = 50, ntree = 700, importance = T)
bag.carseats##
## Call:
## randomForest(formula = High ~ . - Sales, data = cs1, mtry = 50, ntree = 700, importance = T, subset = train2)
## Type of random forest: classification
## Number of trees: 700
## No. of variables tried at each split: 10
##
## OOB estimate of error rate: 19.5%
## Confusion matrix:
## No Yes class.error
## No 94 17 0.1531532
## Yes 22 67 0.2471910yhat.bag <- predict(bag.carseats,newdata=cs1[-train2,])
table(yhat.bag,High.test)## High.test
## yhat.bag No Yes
## No 104 22
## Yes 21 53tp <- 104
tn <- 54
fp <- 21
fn <- 21
print("Accuracy")## [1] "Accuracy"acc <- (tp+tn)/(tp+tn+fp+fn)
print(acc)## [1] 0.79print("Precision")## [1] "Precision"prec <- tp/(tp+tn)
print(prec)## [1] 0.6582278print("Recall")## [1] "Recall"rec <- tp/(tp+fn)
print(rec)## [1] 0.832print("F-Score")## [1] "F-Score"fsc <- (2*prec*rec)/(prec+rec)
print(fsc)## [1] 0.7349823Again, it Looks like we’re trading Precision and Recall more than accuracy.
Problem 9:
This problem involves the OJ data set which is part of the ISLR package.
attach(OJ)
set.seed(42)(a) Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations.
train <- sample(dim(OJ)[1], 800)
oj_train <- OJ[train, ]
oj_test <- OJ[-train, ](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?
library(tree)
oj_tree <- tree(Purchase ~ ., data = oj_train)
summary(oj_tree)##
## Classification tree:
## tree(formula = Purchase ~ ., data = oj_train)
## Variables actually used in tree construction:
## [1] "LoyalCH" "SalePriceMM" "PriceDiff"
## Number of terminal nodes: 8
## Residual mean deviance: 0.7392 = 585.5 / 792
## Misclassification error rate: 0.1638 = 131 / 800This tree ran with three variables: LoyalCH, SalePriceMM, and PriceDiff and has \(8\) terminal nodes. The training misclassification error rate is \(0.1638\).
(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 1066.00 CH ( 0.61500 0.38500 )
## 2) LoyalCH < 0.48285 285 296.00 MM ( 0.21404 0.78596 )
## 4) LoyalCH < 0.064156 64 0.00 MM ( 0.00000 1.00000 ) *
## 5) LoyalCH > 0.064156 221 260.40 MM ( 0.27602 0.72398 )
## 10) SalePriceMM < 2.04 128 123.50 MM ( 0.18750 0.81250 ) *
## 11) SalePriceMM > 2.04 93 125.00 MM ( 0.39785 0.60215 ) *
## 3) LoyalCH > 0.48285 515 458.10 CH ( 0.83689 0.16311 )
## 6) LoyalCH < 0.753545 230 282.70 CH ( 0.69565 0.30435 )
## 12) PriceDiff < 0.265 149 203.00 CH ( 0.57718 0.42282 )
## 24) PriceDiff < -0.165 32 38.02 MM ( 0.28125 0.71875 ) *
## 25) PriceDiff > -0.165 117 150.30 CH ( 0.65812 0.34188 )
## 50) LoyalCH < 0.703993 105 139.60 CH ( 0.61905 0.38095 ) *
## 51) LoyalCH > 0.703993 12 0.00 CH ( 1.00000 0.00000 ) *
## 13) PriceDiff > 0.265 81 47.66 CH ( 0.91358 0.08642 ) *
## 7) LoyalCH > 0.753545 285 111.70 CH ( 0.95088 0.04912 ) *I will use terminal node labeled 11). The \(\tt{split}\) variable is SalePriceMM for this node. It has a split value of \(2.04\) and there are \(93\) points within the sub-tree. The \(\tt{deviance}\) for points contained within this node is \(125.00\). The presence of the * at the end denotes this as a terminal node. The \(\tt{prediction}\) for this node is Sales = \(\tt{MM}\). About \(40\)% of the point in this node are less than 2.04 and the other \(60\)% points are greater than 2.04.
(d) Create a plot of the tree, and interpret the results.
plot(oj_tree)
text(oj_tree, pretty = 0)
‘LoyalCH’ is the most important variable of the tree as observed in the first 3 nodes. If LoyalCH < 0.48 \(\tt{MM}\) and if LoyalCH > 0.48, the tree predicts \(\tt{CH}\). The intermediate values of LoyalCH also depends on the values of SalesPrice and PriceDiff in the deeper nodes.
(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?
oj_pred <- predict(oj_tree, oj_test, type = "class")
table(oj_test$Purchase, oj_pred)## oj_pred
## CH MM
## CH 125 36
## MM 15 94tp <- 125
tn <- 94
fp <- 15
fn <- 36
print("Accuracy")## [1] "Accuracy"acc <- (tp+tn)/(tp+tn+fp+fn)
print(acc)## [1] 0.8111111print("Precision")## [1] "Precision"prec <- tp/(tp+tn)
print(prec)## [1] 0.5707763print("Recall")## [1] "Recall"rec <- tp/(tp+fn)
print(rec)## [1] 0.7763975print("F-Score")## [1] "F-Score"fsc <- (2*prec*rec)/(prec+rec)
print(fsc)## [1] 0.6578947(f) Apply the cv.tree() function to the training set in order to determine the optimal tree size.
cv_oj <- cv.tree(oj_tree, FUN = prune.tree)(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 = "Deviance")
(h) Which tree size corresponds to the lowest cross-validated classification error rate?
A visual inspection shows a tree size of \(5\) to be the most optimal.
(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.
oj_pruned <- prune.tree(oj_tree, best = 5)(j) Compare the training error rates between the pruned and unpruned trees. Which is higher?
summary(oj_pruned)##
## Classification tree:
## snip.tree(tree = oj_tree, nodes = c(5L, 12L))
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff"
## Number of terminal nodes: 5
## Residual mean deviance: 0.7833 = 622.7 / 795
## Misclassification error rate: 0.1812 = 145 / 800This tree ran with two variables this time: LoyalCH, and PriceDiff and has \(5\) terminal nodes. The training misclassification error rate increased to \(0.1812\).
(k) Compare the test error rates between the pruned and unpruned trees. Which is higher?
pred.unpruned <- predict(oj_tree, oj_test, type = "class")
misclass.unpruned <- sum(oj_test$Purchase != pred.unpruned)
misclass.unpruned/length(pred.unpruned)## [1] 0.1888889pred.pruned <- predict(oj_pruned, oj_test, type = "class")
misclass.pruned <- sum(oj_test$Purchase != pred.pruned)
misclass.pruned/length(pred.pruned)## [1] 0.2185185Pruned \(0.189\) and unpruned \(0.22\).trees have different test error rates. This suggests the pruning took out some relevant infomration in regard to accuracy.