3.

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("red", "green", "blue"))

## 8.

(a)

library(MASS)
library(randomForest)
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
library(ISLR)
library(tree)
seatsDf <- Carseats
set.seed(3)
train <- sample(1:nrow(seatsDf), nrow(seatsDf) / 2)
seats.train <- seatsDf[train, ]
seats.test <- seatsDf[-train, ]

(b)

tree.seats <- tree(Sales ~ ., data = seats.train)
summary(tree.seats)
## 
## Regression tree:
## tree(formula = Sales ~ ., data = seats.train)
## Variables actually used in tree construction:
## [1] "ShelveLoc"   "Price"       "Age"         "CompPrice"   "Advertising"
## [6] "US"         
## Number of terminal nodes:  16 
## Residual mean deviance:  2.134 = 392.6 / 184 
## Distribution of residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -4.37400 -0.90790 -0.05181  0.00000  0.92840  3.82600
plot(tree.seats)
text(tree.seats, pretty = 0)

yhat <- predict(tree.seats, newdata = seats.test)
mean((yhat - seats.test$Sales)^2)
## [1] 4.784151

The test MSE from a single tree is 4.784.

(c)

cv.seats <- cv.tree(tree.seats)
plot(cv.seats$size, cv.seats$dev, type = "b")
tree.min <- which.min(cv.seats$dev)
points(tree.min, cv.seats$dev[tree.min], col = "red", cex = 2, pch = 20)

prune.seats <- prune.tree(tree.seats, best = 8)
plot(prune.seats)
text(prune.seats, pretty = 0)

yhat.prune <- predict(prune.seats, newdata = seats.test)
mean((yhat.prune - seats.test$Sales)^2)
## [1] 5.075903

In this case pruning the tree leads to a slight increase in error, with a new MSE of 5.076.

(d)

set.seed(6)
bag.seats <- randomForest(Sales ~ ., data = seats.train, mtry = 10, ntree = 1000, importance = TRUE)
yhat.bag <- predict(bag.seats, newdata = seats.test)
mean((yhat.bag - seats.test$Sales)^2)
## [1] 2.747505
importance(bag.seats)
##                %IncMSE IncNodePurity
## CompPrice   29.1872665    128.909712
## Income       3.4937567     63.072714
## Advertising 21.7891765    125.459028
## Population  -0.7957398     57.341794
## Price       69.9322087    404.579081
## ShelveLoc   80.4667295    507.639547
## Age         19.8312264    117.283535
## Education   -0.5767441     39.332251
## Urban       -1.2445711      8.559391
## US           8.0517710     11.619310

Using the bagging method yeilds the best test MSE so far of 2.748. Checking the importance of the variables show that ShelveLoc, Price, and CompPrice are the three most important variables.

(e)

set.seed(6)
rf.seats <- randomForest(Sales ~ ., data = seats.train, mtry = 3, ntree = 1000, importance = TRUE)
yhat.rf <- predict(rf.seats, newdata = seats.test)
mean((yhat.rf - seats.test$Sales)^2)
## [1] 3.362982
importance(rf.seats)
##                %IncMSE IncNodePurity
## CompPrice   14.6372150     125.13895
## Income       2.5131254     103.92417
## Advertising 19.6699694     164.97046
## Population  -0.1663633     102.79109
## Price       44.0660774     328.47456
## ShelveLoc   50.2666855     342.55792
## Age         10.8005302     134.97912
## Education    1.4084685      62.68963
## Urban       -0.8030655      14.38960
## US          10.3917946      34.80518

The random forest model perform better than the first two, but slightly worse than the bagged model, with a test MSE of 3.383. The most important variables in this model are ShelveLoc, Price, and Advertising. ## 9.

(a)

ojDf <- OJ
set.seed(1)
train <- sample(1:nrow(ojDf), 800)
OJ.train <- ojDf[train, ]
OJ.test <- ojDf[-train, ]

(b)

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"     "SpecialCH"     "ListPriceDiff"
## [5] "PctDiscMM"    
## Number of terminal nodes:  9 
## Residual mean deviance:  0.7432 = 587.8 / 791 
## Misclassification error rate: 0.1588 = 127 / 800

The training MSE is 0.7432, with 9 terminal nodes.

(c)

tree.oj
## node), split, n, deviance, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 800 1073.00 CH ( 0.60625 0.39375 )  
##    2) LoyalCH < 0.5036 365  441.60 MM ( 0.29315 0.70685 )  
##      4) LoyalCH < 0.280875 177  140.50 MM ( 0.13559 0.86441 )  
##        8) LoyalCH < 0.0356415 59   10.14 MM ( 0.01695 0.98305 ) *
##        9) LoyalCH > 0.0356415 118  116.40 MM ( 0.19492 0.80508 ) *
##      5) LoyalCH > 0.280875 188  258.00 MM ( 0.44149 0.55851 )  
##       10) PriceDiff < 0.05 79   84.79 MM ( 0.22785 0.77215 )  
##         20) SpecialCH < 0.5 64   51.98 MM ( 0.14062 0.85938 ) *
##         21) SpecialCH > 0.5 15   20.19 CH ( 0.60000 0.40000 ) *
##       11) PriceDiff > 0.05 109  147.00 CH ( 0.59633 0.40367 ) *
##    3) LoyalCH > 0.5036 435  337.90 CH ( 0.86897 0.13103 )  
##      6) LoyalCH < 0.764572 174  201.00 CH ( 0.73563 0.26437 )  
##       12) ListPriceDiff < 0.235 72   99.81 MM ( 0.50000 0.50000 )  
##         24) PctDiscMM < 0.196197 55   73.14 CH ( 0.61818 0.38182 ) *
##         25) PctDiscMM > 0.196197 17   12.32 MM ( 0.11765 0.88235 ) *
##       13) ListPriceDiff > 0.235 102   65.43 CH ( 0.90196 0.09804 ) *
##      7) LoyalCH > 0.764572 261   91.20 CH ( 0.95785 0.04215 ) *

Picking the last node, labeled 7, is populated by LoyalCh > 0.764572, contains 261 observations, has a deviance of 91.2, and is made up of about 96% CH and 4% MM.

(d)

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

This chart shows that LoyalCH is the most important split in determining Purchase.

(e)

tree.pred <- predict(tree.oj, OJ.test, type = "class")
table(tree.pred, OJ.test$Purchase)
##          
## tree.pred  CH  MM
##        CH 160  38
##        MM   8  64
(1 -(160 + 64)/270)
## [1] 0.1703704

The test error rate is approximately 17%. ### (f)

cv.oj <- cv.tree(tree.oj, FUN = prune.misclass)
cv.oj
## $size
## [1] 9 8 7 4 2 1
## 
## $dev
## [1] 150 150 149 158 172 315
## 
## $k
## [1]       -Inf   0.000000   3.000000   4.333333  10.500000 151.000000
## 
## $method
## [1] "misclass"
## 
## attr(,"class")
## [1] "prune"         "tree.sequence"

(g)

plot(cv.oj$size, cv.oj$dev, type = "b", xlab = "Tree size", ylab = "Deviance")

(h)

It appears that a 7-node tree has the lowest error rate.

(i)

prune.oj <- prune.misclass(tree.oj, best = 7)
plot(prune.oj)
text(prune.oj, pretty = 0)

(j)

summary(tree.oj)
## 
## Classification tree:
## tree(formula = Purchase ~ ., data = OJ.train)
## Variables actually used in tree construction:
## [1] "LoyalCH"       "PriceDiff"     "SpecialCH"     "ListPriceDiff"
## [5] "PctDiscMM"    
## Number of terminal nodes:  9 
## Residual mean deviance:  0.7432 = 587.8 / 791 
## Misclassification error rate: 0.1588 = 127 / 800
summary(prune.oj)
## 
## Classification tree:
## snip.tree(tree = tree.oj, nodes = c(4L, 10L))
## Variables actually used in tree construction:
## [1] "LoyalCH"       "PriceDiff"     "ListPriceDiff" "PctDiscMM"    
## Number of terminal nodes:  7 
## Residual mean deviance:  0.7748 = 614.4 / 793 
## Misclassification error rate: 0.1625 = 130 / 800

The misclassification rate is slightly higher for the pruned tree 16.25% to 15.88%.

(k)

prune.pred <- predict(prune.oj, OJ.test, type = "class")
table(prune.pred, OJ.test$Purchase)
##           
## prune.pred  CH  MM
##         CH 160  36
##         MM   8  66
(1 - (160 + 66)/270)
## [1] 0.162963

The test error rate differs from the train in that the misclassification rate is lower for the pruned tree, 16.3% to 17%.