Question 3 : Consider the Gini index, classification error, and entropy in a simple classification setting with two classes.

p = seq(0, 1, 0.001)
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), ylab = "gini.index, class.error, cross.entropy", col = c("red", "blue", "green"))

Question 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

a) Split the data set into a training set and a test set.

library(ISLR2)

data("Carseats")
set.seed(1)
train = sample(1:nrow(Carseats), nrow(Carseats)/2)
strain = Carseats[train, ]
stest = 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(ISLR)
install.packages("tree", repos = "https://cloud.r-project.org")
## 
## The downloaded binary packages are in
##  /var/folders/43/50tq8sf173g3h31yf1vx8v5r0000gn/T//RtmpZdPX2u/downloaded_packages
library(tree)
tree.seats = tree(Sales ~ ., data = strain)
summary(tree.seats)
## 
## Regression tree:
## tree(formula = Sales ~ ., data = strain)
## Variables actually used in tree construction:
## [1] "ShelveLoc"   "Price"       "Age"         "Advertising" "CompPrice"  
## [6] "US"         
## Number of terminal nodes:  18 
## Residual mean deviance:  2.167 = 394.3 / 182 
## Distribution of residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -3.88200 -0.88200 -0.08712  0.00000  0.89590  4.09900
plot(tree.seats)
text(tree.seats, pretty = 0)

treeseat.pred = predict(tree.seats, newdata = stest)
mean((treeseat.pred - stest$Sales)^2)
## [1] 4.922039

c) Use cross-validation in order to determine the optimal level of tree complexity. Does pruning the tree improve the test MSE?

set.seed(1)
cv.seats = cv.tree(tree.seats)
plot(cv.seats$size, cv.seats$dev, type = "b")

prune.car = prune.tree(tree.seats, best = 10)
plot(prune.car)
text(prune.car,pretty=0)

treeseat.pred = predict(prune.car, newdata = stest)
mean((treeseat.pred - stest$Sales)^2)
## [1] 4.918134

Here, pruning did not improve the MSE of my data set as it went up from 4.93 to 5.08

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.

install.packages("randomForest", repos = "https://cloud.r-project.org")
## 
## The downloaded binary packages are in
##  /var/folders/43/50tq8sf173g3h31yf1vx8v5r0000gn/T//RtmpZdPX2u/downloaded_packages
library(randomForest)

set.seed(1)
bag.seats = randomForest(Sales~., data = strain, mtry = 10, ntree = 551, importance = TRUE)
bagseat.pred = predict(bag.seats, newdata = stest)
mean((bagseat.pred - stest$Sales)^2)
## [1] 2.631294
importance(bag.seats)
##                 %IncMSE IncNodePurity
## CompPrice   25.76010601    165.344436
## Income       4.30770157     93.248899
## Advertising 15.19995373    102.306849
## Population  -1.73412418     59.239991
## Price       60.98677002    506.990148
## ShelveLoc   48.51140673    370.824288
## Age         19.12912973    155.074163
## Education    0.04806214     44.856419
## Urban        0.25153531      9.365518
## US           4.71331539     15.679252

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(1)
rando.seats = randomForest(Sales~., data = strain, mtry = 10, importance = TRUE)
randseat.pred = predict(rando.seats, newdata = stest)
mean((randseat.pred - stest$Sales)^2)
## [1] 2.640028

This MSE of 2.640028 is lower than the initial MSE of 4.922039 we got.

importance(rando.seats)
##                %IncMSE IncNodePurity
## CompPrice   24.3747505    165.896498
## Income       3.6994834     93.164819
## Advertising 14.5556207    103.424421
## Population  -2.1869464     59.516894
## Price       57.5700828    505.663115
## ShelveLoc   46.2409020    371.326940
## Age         18.5373118    154.490686
## Education   -0.3462442     43.863808
## Urban        0.5652769      9.350136
## US           4.2140560     16.066279

Question 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(1)
train = sample(1:nrow(OJ), 800)
OJtrain = OJ[train, ]
OJtest = 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 = OJtrain)
summary(tree.OJ)
## 
## Classification tree:
## tree(formula = Purchase ~ ., data = OJtrain)
## 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

This tree has 9 terminal nodes and the training error rate is 0.1588

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

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 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.196196 55   73.14 CH ( 0.61818 0.38182 ) *
##         25) PctDiscMM > 0.196196 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 ) *

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

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

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

treeOJ.pred = predict(tree.OJ, newdata = OJtest, type = "class")
table(treeOJ.pred, OJtest$Purchase)
##            
## treeOJ.pred  CH  MM
##          CH 160  38
##          MM   8  64

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

OJcv = cv.tree(tree.OJ, FUN = prune.misclass)
OJcv
## $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) Produce a plot with tree size on the x-axis and cross-validated classification error rate on the y-axis.

plot(OJcv$size, OJcv$dev, type = "b", xlab = "Tree Size", ylab = "cv classification error rate")

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

–> The tree size of 7 corresponds to the lowest cross-validated classification error rate

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=7)

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

summary(tree.OJ)
## 
## Classification tree:
## tree(formula = Purchase ~ ., data = OJtrain)
## 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(10L, 4L))
## 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

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

treeOJ.pred = predict(tree.OJ, newdata = OJtest, type = "class")
table(treeOJ.pred, OJtest$Purchase)
##            
## treeOJ.pred  CH  MM
##          CH 160  38
##          MM   8  64
unprunedOJvalerr = (38 + 8) / 270
unprunedOJvalerr
## [1] 0.1703704
pruneOJ.pred = predict(prune.OJ, newdata = OJtest, type = "class")
table(pruneOJ.pred, OJtest$Purchase)
##             
## pruneOJ.pred  CH  MM
##           CH 160  36
##           MM   8  66
prunedOJvalerr = (36 + 8) / 270
prunedOJvalerr
## [1] 0.162963

Pruning slightly decreases the test error rate, but they are close.