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 \(\hat{P}_{m1}\). The x-axis should display \(\hat{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, \(\hat{P}_{m1}\) = 1 − \(\hat{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.index = p * (1 - p) * 2
entropy = -(p * log(p) + (1 - p) * log(1 - p))
class.err = 1 - pmax(p, 1 - p)
matplot(p, cbind(gini.index, class.err, entropy), col = c("green", "purple", "blue"), pch=16, ylab = "gini.index, class.error, entropy")
legend('bottom', legend = c('gini.index', 'class.error', 'entropy'), col = c("green", "purple", "blue"), pch=c(16))

9. This problem involves the OJ data set which is part of the ISLR package.

attach(OJ)

(a) Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations.

set.seed(123)
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?

tree.oj = tree(Purchase~., data=OJ.train)
summary(tree.oj)
## 
## Classification tree:
## tree(formula = Purchase ~ ., data = OJ.train)
## 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 / 800

The tree has 8 terminal nodes and the two variables used are LoyalCh and PriceDiff. The training error rate is 0.165.

(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.

tree.oj
## 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 ) *

Terminal node 29 shows PriceDiff > 0.145 as the split criterion, there are 87 observations in that branch, the deviance equals 62.07, and the overall prediction in that branch is CH. 89% of the values in the branch are CH and 11% are MM.

(d) Create a plot of the tree, and interpret the results.

plot(tree.oj)
text(tree.oj, pretty=0)

If LoyalCH is less than 0.504, then MM is the more likely result. If LoyalCH is greater than 0.504, CH is the more likely result.

(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?

tree.pred=predict(tree.oj, newdata=OJ.test, type="class")
table(tree.pred,OJ.test$Purchase)
##          
## tree.pred  CH  MM
##        CH 150  34
##        MM  16  70

Error rate: 16+34/270 = 0.185

(f) Apply the cv.tree() function to the training set in order to determine the optimal tree size.

cv.oj = cv.tree(tree.oj, 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="CV Error")

(h) Which tree size corresponds to the lowest cross-validated classification error rate?

According to (g), 6 is the lowest tree size.

(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(tree.oj, best=6)
plot(prune.oj)
text(prune.oj, pretty=0)

(j) Compare the training error rates between the pruned and unpruned trees. Which is higher?

summary(tree.oj)
## 
## Classification tree:
## tree(formula = Purchase ~ ., data = OJ.train)
## 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 / 800
summary(prune.oj)
## 
## Classification tree:
## snip.tree(tree = tree.oj, 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 / 800

They are the same.

(k) Compare the test error rates between the pruned and unpruned trees. Which is higher?

prune.pred = predict(prune.oj, OJ.test, type = "class")
table(prune.pred, OJ.test$Purchase)
##           
## prune.pred  CH  MM
##         CH 150  34
##         MM  16  70

Error rate: 16+34/270 = 0.185. They are the same.