STAT-Assignment-7
Shalini Kurumathur
4/20/2022
R Markdown
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 x-axis 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)
#Classification error rate
class.error= 1-pmax(p,1-p)
#Gini index
gini= p*(1-p)*2
#Entropy
entropy = -(p*log(p) + (1-p)*log(1-p))
matplot(p, cbind(gini, entropy, class.error), col = c("green", "red", "purple"))
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)
Carseats <- read.csv("C:/Users/kksha/OneDrive/coursework/Statistics/assignment/Carseats.csv")
Carseats$Urban <-as.factor(Carseats$Urban)
Carseats$US <-as.factor(Carseats$US)
Carseats$ShelveLoc<-as.factor(Carseats$ShelveLoc)
Carseats$CompPrice<-as.numeric(Carseats$CompPrice)
Carseats$Income<-as.numeric(Carseats$Income)
Carseats$Advertising<-as.numeric(Carseats$Advertising)
Carseats$Population<-as.numeric(Carseats$Population)
Carseats$Price<-as.numeric(Carseats$Price)
Carseats$Age<-as.numeric(Carseats$Age)
Carseats$Education<-as.numeric(Carseats$Education)
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.
set.seed(1)
#Splitting into training and test dataset
train = sample(dim(Carseats)[1], dim(Carseats)[1]/2)
ctrain = Carseats[train, ]
ctest = 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(MASS)
library(tree)## Warning: package 'tree' was built under R version 4.1.3set.seed(1)
tree.car=tree(Sales~.,ctrain)
summary(tree.car)##
## Regression tree:
## tree(formula = Sales ~ ., data = ctrain)
## 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.09900Six variables have been used in the tree: “ShelveLoc”,“Price”,“Age”,“Advertising”,“CompPrice” and “US”. Residual mean devaiance is the sum of squared errors for the tree, which is 2.167.
Plotting the regression tree
plot(tree.car)
text(tree.car , pretty =0)
pred.car = predict(tree.car, ctest)
mean((ctest$Sales - pred.car)^2)## [1] 4.922039The mean squared error is 4.922039
(c) Use cross-validation in order to determine the optimal level of tree complexity. Does pruning the tree improve the test MSE?
#using cv.tree to check if pruning the tree will improve performance
cv.car=cv.tree(tree.car)
plot(cv.car$size ,cv.car$dev ,type='b')
prune.car=prune.tree(tree.car ,best=9)
plot(prune.car)
text(prune.car , pretty =0)
pred.pruned = predict(prune.car, ctest)
mean((ctest$Sales - pred.pruned)^2)## [1] 4.918134The mean squared error has now reduced to 4.918134.
(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)## Warning: package 'randomForest' was built under R version 4.1.3## randomForest 4.7-1## Type rfNews() to see new features/changes/bug fixes.set.seed(1)
bag.car= randomForest(Sales~.,data=ctrain,mtry=10,importance =TRUE)
bag.car##
## Call:
## randomForest(formula = Sales ~ ., data = ctrain, mtry = 10, importance = TRUE)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 10
##
## Mean of squared residuals: 2.889221
## % Var explained: 63.26bag.pred = predict(bag.car, ctest)
mean((ctest$Sales - bag.pred)^2)## [1] 2.605253Bagging reduced the MSE to 2.605253.
importance(bag.car)## %IncMSE IncNodePurity
## CompPrice 24.8888481 170.182937
## Income 4.7121131 91.264880
## Advertising 12.7692401 97.164338
## Population -1.8074075 58.244596
## Price 56.3326252 502.903407
## ShelveLoc 48.8886689 380.032715
## Age 17.7275460 157.846774
## Education 0.5962186 44.598731
## Urban 0.1728373 9.822082
## US 4.2172102 18.073863“Price”,“ShelveLoc”,“CompPrice” and “Age” are important variables.
(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.
#Random forest is same as bagging, but uses less number of predictors. Here we use mtry as 5
library( randomForest)
set.seed(1)
rf.car= randomForest(Sales~.,data=ctrain,mtry=5,importance =TRUE)
rf.car##
## Call:
## randomForest(formula = Sales ~ ., data = ctrain, mtry = 5, importance = TRUE)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 5
##
## Mean of squared residuals: 3.060785
## % Var explained: 61.08rf.pred = predict(rf.car, ctest)
mean((ctest$Sales - rf.pred)^2)## [1] 2.714168Random Forest worsens the test MSE compared to Bagging.
importance(rf.car)## %IncMSE IncNodePurity
## CompPrice 17.4126238 157.53631
## Income 2.9969399 110.40731
## Advertising 11.0485672 105.75049
## Population -1.5321044 80.73318
## Price 43.3572135 452.02367
## ShelveLoc 44.4474163 331.64508
## Age 14.5322339 176.64252
## Education 0.8237454 55.91141
## Urban -2.7805788 11.07321
## US 3.7773881 23.75322“Price”, “ShelveLoc” and “Age” are three important variables.
**9. This problem involves the OJ data set which is part of the ISLR package.
library(ISLR)
View(OJ)(a) Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations.
library(ISLR)
attach(OJ)
dim(OJ)## [1] 1070 18train = sample(dim(OJ)[1],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?
library(MASS)
library(tree)
set.seed(1)
tree.OJ=tree(Purchase~.,OJtrain)
summary(tree.OJ)##
## Classification tree:
## tree(formula = Purchase ~ ., data = OJtrain)
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff" "ListPriceDiff"
## Number of terminal nodes: 8
## Residual mean deviance: 0.7543 = 597.4 / 792
## Misclassification error rate: 0.1725 = 138 / 800The number of terminal nodes are 8. The training error rate is 0.1725 The tree has three variables: “LoyalCH”,“PriceDiff” and “ListPriceDiff”.
(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 1063.00 CH ( 0.61875 0.38125 )
## 2) LoyalCH < 0.5036 355 432.90 MM ( 0.29859 0.70141 )
## 4) LoyalCH < 0.275386 170 138.40 MM ( 0.14118 0.85882 )
## 8) LoyalCH < 0.0356415 59 10.14 MM ( 0.01695 0.98305 ) *
## 9) LoyalCH > 0.0356415 111 113.30 MM ( 0.20721 0.79279 ) *
## 5) LoyalCH > 0.275386 185 254.10 MM ( 0.44324 0.55676 )
## 10) PriceDiff < 0.05 77 83.74 MM ( 0.23377 0.76623 ) *
## 11) PriceDiff > 0.05 108 146.00 CH ( 0.59259 0.40741 ) *
## 3) LoyalCH > 0.5036 445 336.80 CH ( 0.87416 0.12584 )
## 6) LoyalCH < 0.764572 185 207.50 CH ( 0.75135 0.24865 )
## 12) PriceDiff < 0.265 111 148.20 CH ( 0.61261 0.38739 )
## 24) ListPriceDiff < 0.235 64 88.47 MM ( 0.46875 0.53125 ) *
## 25) ListPriceDiff > 0.235 47 45.91 CH ( 0.80851 0.19149 ) *
## 13) PriceDiff > 0.265 74 25.11 CH ( 0.95946 0.04054 ) *
## 7) LoyalCH > 0.764572 260 84.77 CH ( 0.96154 0.03846 ) *Let us take terminal node label (10). The splitting node label is ‘PriceDiff’. The splitting value of this node is ‘0.05’.There are 77 points in the subtree below this node. A * denotes that it is a terminal node. Prediction of this node is Sales - MM.About 23.377 % points in the node have ‘CH’ as the value for Sales, rest 76.623% have the value ‘MM’ for Sales.
(d) Create a plot of the tree, and interpret the results.
plot(tree.OJ)
text(tree.OJ, pretty =0)
‘LoyalCH’ is a most important variable in the tree. If ‘LoyalCH’ is below 0.0356415, it predicts Sales as ‘MM’. If ‘LoyalCH’ is above 0.764572, it predicts ‘CH’. For the intermediate values of ‘LoyalCH’, it prediction of Sales depends on ‘PriceDiff’.
(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?
set.seed(1)
pred.OJ = predict(tree.OJ, OJtest, type = "class")
pred.unpruned = pred.OJ
table(OJtest$Purchase,pred.OJ)## pred.OJ
## CH MM
## CH 138 20
## MM 29 83Test error rate is (20+29)/(20+29+138+83) = 0.1814
(f) Apply the cv.tree() function to the training set in order to determine the optimal tree size.
#using cv.tree to check if pruning the tree will improve performance
cv.OJ=cv.tree(tree.OJ)
cv.OJ## $size
## [1] 8 7 6 5 4 3 2 1
##
## $dev
## [1] 712.9841 714.8318 714.8318 720.2176 748.5357 771.2814 776.5165
## [8] 1066.6902
##
## $k
## [1] -Inf 13.81949 15.00203 24.33767 34.20403 40.37572 44.49681
## [8] 293.83140
##
## $method
## [1] "deviance"
##
## 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.
#using cv.tree to check if pruning the tree will improve performance
plot(cv.OJ$size ,cv.OJ$dev ,type='b')
(h) Which tree size corresponds to the lowest cross-validated classification error rate?
The size of 6 shows lowest cross-validation error.
(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(prune.OJ)##
## Classification tree:
## snip.tree(tree = tree.OJ, nodes = c(12L, 4L))
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff"
## Number of terminal nodes: 6
## Residual mean deviance: 0.7887 = 626.2 / 794
## Misclassification error rate: 0.1775 = 142 / 800The training error rate of the pruned data is 0.1775.The training error rate is reduced after pruning.
(k) Compare the test error rates between the pruned and unpruned trees. Which is higher?
pred.pruned1 = predict(prune.OJ, OJtest, type = "class")
misclass.pruned = sum(OJtest$Purchase != pred.pruned1)
misclass.pruned/length(pred.pruned1)## [1] 0.1814815pred.unpruned = predict(tree.OJ, OJtest, type = "class")
misclass.unpruned = sum(OJtest$Purchase != pred.unpruned)
misclass.unpruned/length(pred.unpruned)## [1] 0.1814815The test error rates for both pruned and unpruned trees is the same.