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. Hint: In a setting with two classes, pˆm1 = 1 − 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 <- 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"))
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.
Split the data set into a training set and a test set.
library(ISLR)
attach(Carseats)
set.seed(1)
subset<-sample(nrow(Carseats),nrow(Carseats)*0.7)
car.train<-Carseats[subset,]
car.test<-Carseats[-subset,]
Fit a regression tree to the training set. Plot the tree, and interpret the results. What test MSE do you obtain?
library(tree)
car.model.train<-tree(Sales~.,car.train)
summary(car.model.train)
##
## Regression tree:
## tree(formula = Sales ~ ., data = car.train)
## Variables actually used in tree construction:
## [1] "ShelveLoc" "Price" "Age" "Income" "CompPrice"
## [6] "Advertising"
## Number of terminal nodes: 18
## Residual mean deviance: 2.409 = 631.1 / 262
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -4.77800 -0.96100 -0.08865 0.00000 1.01800 4.14100
plot(car.model.train)
text(car.model.train,pretty=0)
tree.prediction<-predict(car.model.train,newdata=car.test)
tree.mse<-mean((car.test$Sales-tree.prediction)^2)
tree.mse
## [1] 4.208383
The test mse is at 4.2.
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.car<-cv.tree(car.model.train)
plot(cv.car$size,cv.car$dev,xlab = "Size of Tree",ylab = "Deviance",type = "b")
prune.car<-prune.tree(car.model.train,best=6)
plot(prune.car)
text(prune.car,pretty=0)
prune.predict<-predict(prune.car,car.test)
mean((prune.predict-car.test$Sales)^2)
## [1] 5.118217
With the pruned tree we get mse of about 5.11
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(caret)
## Loading required package: lattice
## Loading required package: ggplot2
bag.car<-train(Sales~.,car.train,method='rf',importance=TRUE)
#importance(bag.car)
varImp(bag.car)
## rf variable importance
##
## Overall
## Price 100.000
## ShelveLocGood 95.884
## CompPrice 45.236
## Age 35.085
## ShelveLocMedium 29.344
## Advertising 27.625
## Income 11.135
## USYes 7.139
## Education 5.290
## Population 1.012
## UrbanYes 0.000
bag.car.predict<-predict(bag.car,car.test)
mean((bag.car.predict-car.test$Sales)^2)
## [1] 2.652016
With randomForest, mse is 2.64 which is less mse than the pruned tree method.
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.car<-train(Sales~.,car.train,method='rf',importance=TRUE)
#importance(rf.car)
varImp(rf.car)
## rf variable importance
##
## Overall
## ShelveLocGood 100.000
## Price 99.035
## CompPrice 42.773
## Age 34.101
## ShelveLocMedium 31.497
## Advertising 28.167
## Income 9.834
## USYes 5.848
## Education 5.658
## Population 1.649
## UrbanYes 0.000
rf.car.predict<-predict(rf.car,car.test)
mean((rf.car.predict-car.test$Sales)^2)
## [1] 2.608392
At 2.61 mse, randomForest method gives use the lowest number. This is the best method for this data.
Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations.
library(ISLR)
attach(OJ)
train<-sample(1:nrow(OJ),800)
oj.train<-OJ[train,]
oj.test<-OJ[-train,]
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?
oj.tree<-tree(Purchase~.,oj.train)
fulltree.summary<-summary(oj.tree)
fulltree.summary
##
## Classification tree:
## tree(formula = Purchase ~ ., data = oj.train)
## Variables actually used in tree construction:
## [1] "LoyalCH" "STORE" "PriceDiff" "ListPriceDiff"
## Number of terminal nodes: 9
## Residual mean deviance: 0.7194 = 569 / 791
## Misclassification error rate: 0.1575 = 126 / 800
fulltree.misclass<-fulltree.summary$misclass[1]*100/fulltree.summary$misclass[2]
fulltree.misclass
## [1] 15.75
Full tree classification error rate is at 16.5%, and the terminal nodes came out to be 8.
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 1077.000 CH ( 0.600000 0.400000 )
## 2) LoyalCH < 0.5036 360 415.500 MM ( 0.263889 0.736111 )
## 4) LoyalCH < 0.280875 177 120.700 MM ( 0.107345 0.892655 )
## 8) STORE < 1.5 54 59.610 MM ( 0.240741 0.759259 ) *
## 9) STORE > 1.5 123 47.950 MM ( 0.048780 0.951220 ) *
## 5) LoyalCH > 0.280875 183 248.400 MM ( 0.415301 0.584699 )
## 10) PriceDiff < 0.065 77 83.740 MM ( 0.233766 0.766234 ) *
## 11) PriceDiff > 0.065 106 146.000 CH ( 0.547170 0.452830 ) *
## 3) LoyalCH > 0.5036 440 331.600 CH ( 0.875000 0.125000 )
## 6) PriceDiff < 0.31 303 284.000 CH ( 0.821782 0.178218 )
## 12) LoyalCH < 0.753545 134 166.600 CH ( 0.686567 0.313433 )
## 24) PriceDiff < -0.165 27 32.820 MM ( 0.296296 0.703704 )
## 48) ListPriceDiff < 0.115 11 14.420 CH ( 0.636364 0.363636 ) *
## 49) ListPriceDiff > 0.115 16 7.481 MM ( 0.062500 0.937500 ) *
## 25) PriceDiff > -0.165 107 111.400 CH ( 0.785047 0.214953 ) *
## 13) LoyalCH > 0.753545 169 86.610 CH ( 0.928994 0.071006 ) *
## 7) PriceDiff > 0.31 137 11.830 CH ( 0.992701 0.007299 ) *
Create a plot of the tree, and interpret the results.
plot(oj.tree)
text(oj.tree,pretty=0)
LoyalCH appears to be the most influential as it takes the top 3 nodes.
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.predict.test<-predict(oj.tree,newdata = oj.test,type = "class")
table(oj.test$Purchase,oj.predict.test,dnn = c("Actual","Predicted"))
## Predicted
## Actual CH MM
## CH 159 14
## MM 32 65
oj.testerror<-(21+37)/nrow(oj.test)
round(oj.testerror,3)
## [1] 0.215
Test error rate for the full grown tree is 21.5%
Apply the cv.tree() function to the training set in order to determine the optimal tree size
oj.cv.train<-cv.tree(oj.tree,FUN = prune.misclass)
oj.cv.train
## $size
## [1] 9 8 7 4 2 1
##
## $dev
## [1] 155 154 156 156 157 320
##
## $k
## [1] -Inf 0.000000 3.000000 3.666667 5.000000 170.000000
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"
Produce a plot with tree size on the x-axis and cross-validated classification error rate on the y-axis
plot(oj.cv.train$size,oj.cv.train$dev,xlab="Size of the Tree",ylab="CV Deviance",type = "b")
points(4,min(oj.cv.train$dev),col="red")
Which tree size corresponds to the lowest cross-validated classification error rate?
The deviance is minimum on a cross validated dateset for a tree size of 4 as seen above.
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.prune<-prune.misclass(oj.tree,best=4)
plot(oj.prune)
text(oj.prune,pretty=0)
Compare the training error rates between the pruned and unpruned trees. Which is higher?
prune.summary<-summary(oj.prune)
prune.summary
##
## Classification tree:
## snip.tree(tree = oj.tree, nodes = 4:3)
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff"
## Number of terminal nodes: 4
## Residual mean deviance: 0.8568 = 682 / 796
## Misclassification error rate: 0.175 = 140 / 800
prune.misclas<-prune.summary$misclass[1]*100/prune.summary$misclass[2]
Pruning did not seem to help as it was previously 16.5% for a full tree, and 17.5% as seen above.
Compare the test error rates between the pruned and unpruned trees. Which is higher?
prune.predict.test<-predict(oj.prune,newdata = oj.test,type="class")
table(oj.test$Purchase,prune.predict.test,dnn = c("Actual","Predicted"))
## Predicted
## Actual CH MM
## CH 161 12
## MM 39 58
oj.testerror.prune<-(40+15)/nrow(oj.test)
round(oj.testerror.prune,3)
## [1] 0.204
Full tree is at 21.5% misclassifications error rate, and pruned tree at 20.4% seen above. We see that pruning is less, meaning that over fitting is occurring for a full tree since it gave a lower training misclassification error and a higher test misclassificaiton error.