Assignment-7
Sonali
2022-11-14
Problem 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 \(p^_m1\). 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)
gini.ind=2 * p * (1 - p)
classification.error=1 - pmax(p, 1 - p)
entropy=- (p * log(p) + (1 - p) * log(1 - p))
matplot(p, cbind(gini.ind, classification.error, entropy), pch=c(15,17,19) ,ylab = "gini.ind, classification.error, entropy",col = c("darkblue" , "yellow", "red"), type = 'b')
legend('bottom', inset=.01, legend = c('gini.ind', 'classification.error', 'entropy'), col = c("darkblue" , "yellow", "red"), pch=c(15,17,19))
Problem 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.
(a) Split the data set into a training set and a test set.
library(ISLR)
library(tree)
attach(Carseats)set.seed(1)
train = sample(dim(Carseats)[1], dim(Carseats)[1]/2)
Carseats.train = Carseats[train, ]
Carseats.test = Carseats[-train, ](b) Fit a regression tree to the training set. Plot the tree, and interpret the results. What test MSE do you obtain?
tree.carseats = tree(Sales ~ ., data = Carseats.train)
summary(tree.carseats)##
## Regression tree:
## tree(formula = Sales ~ ., data = Carseats.train)
## 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.09900plot(tree.carseats)
text(tree.carseats, pretty = 0)
pred.carseats = predict(tree.carseats, Carseats.test)
mean((Carseats.test$Sales - pred.carseats)^2)## [1] 4.922039Test MSE for the tree is 4.92
(c) Use cross-validation in order to determine the optimal level
of tree complexity. Does pruning the tree improve the test
MSE?
cv.carseats = cv.tree(tree.carseats, FUN = prune.tree)
par(mfrow = c(1, 2))
plot(cv.carseats$size, cv.carseats$dev, type = "b")
plot(cv.carseats$k, cv.carseats$dev, type = "b")
#Seems to be 17 is the best size
pruned.carseats = prune.tree(tree.carseats, best = 17)
par(mfrow = c(1, 1))
plot(pruned.carseats)
text(pruned.carseats, pretty = 0)
pred.pruned = predict(pruned.carseats, Carseats.test)
mean((Carseats.test$Sales - pred.pruned)^2)## [1] 4.827162OUr MSE has decreased to 4.82
(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)bag.carseats = randomForest(Sales ~ ., data = Carseats.train, mtry = 10, ntree = 500, importance = T)
bag.pred = predict(bag.carseats, Carseats.test)
mean((Carseats.test$Sales - bag.pred)^2)## [1] 2.657296importance(bag.carseats)## %IncMSE IncNodePurity
## CompPrice 23.07909904 171.185734
## Income 2.82081527 94.079825
## Advertising 11.43295625 99.098941
## Population -3.92119532 59.818905
## Price 54.24314632 505.887016
## ShelveLoc 46.26912996 361.962753
## Age 14.24992212 159.740422
## Education -0.07662320 46.738585
## Urban 0.08530119 8.453749
## US 4.34349223 15.157608By Using bagging approach the MSE has decreased to 2.65. Using
importance function we see ‘Price’, ‘Shelveloc’, ‘Age’and ’ComPrice’ are
the most important variables for sale.
(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.
rf.carseats = randomForest(Sales ~ ., data = Carseats.train, mtry = 5, ntree = 500, importance = T)
rf.pred = predict(rf.carseats, Carseats.test)
mean((Carseats.test$Sales - rf.pred)^2)## [1] 2.701665importance(rf.carseats)## %IncMSE IncNodePurity
## CompPrice 19.8160444 162.73603
## Income 2.8940268 106.96093
## Advertising 11.6799573 106.30923
## Population -1.6998805 79.04937
## Price 46.3454015 448.33554
## ShelveLoc 40.4412189 334.33610
## Age 12.5440659 169.06125
## Education 1.0762096 55.87510
## Urban 0.5703583 13.21963
## US 5.8799999 25.59797By using Random Forest the MSE has increased to 2.7 but the important
variables impacting the sales remains the same.
(f) Now analyze the data using BART, and report your
results.
library(BART)x=Carseats[,1:11]
y=Carseats[,"Sales"]
xtrain=x[train, ]
ytrain=y[train]
xtest=x[-train, ]
ytest=y[-train]
set.seed(5)
bart.fit=gbart(xtrain,ytrain,x.test=xtest)## *****Calling gbart: type=1
## *****Data:
## data:n,p,np: 200, 15, 200
## y1,yn: 2.781850, 1.091850
## x1,x[n*p]: 10.360000, 1.000000
## xp1,xp[np*p]: 11.220000, 1.000000
## *****Number of Trees: 200
## *****Number of Cut Points: 100 ... 1
## *****burn,nd,thin: 100,1000,1
## *****Prior:beta,alpha,tau,nu,lambda,offset: 2,0.95,0.273474,3,4.01382e-30,7.57815
## *****sigma: 0.000000
## *****w (weights): 1.000000 ... 1.000000
## *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,15,0
## *****printevery: 100
##
## MCMC
## done 0 (out of 1100)
## done 100 (out of 1100)
## done 200 (out of 1100)
## done 300 (out of 1100)
## done 400 (out of 1100)
## done 500 (out of 1100)
## done 600 (out of 1100)
## done 700 (out of 1100)
## done 800 (out of 1100)
## done 900 (out of 1100)
## done 1000 (out of 1100)
## time: 2s
## trcnt,tecnt: 1000,1000yhat.bart=bart.fit$yhat.test.mean
mean((ytest-yhat.bart)^2)## [1] 0.1100149By Using BART our test MSE is 0.11
detach(Carseats)Problem 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.
library(ISLR2)attach(OJ)
set.seed(5)
train=sample(dim(OJ)[1],800)
train.OJ=OJ[train,]
test.OJ=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?
oj.tree=tree(Purchase~.,data=train.OJ)
summary(oj.tree)##
## Classification tree:
## tree(formula = Purchase ~ ., data = train.OJ)
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff" "ListPriceDiff"
## Number of terminal nodes: 9
## Residual mean deviance: 0.7347 = 581.1 / 791
## Misclassification error rate: 0.1662 = 133 / 800The Summary statistics about this tree are: Misclassification error
rate is 0.1662, Number or terminal nodes is 9 and variables “LoyalCH”
PriceDiff” ListPriceDiff” are actually used in tree construction.
(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.
oj.tree## node), split, n, deviance, yval, (yprob)
## * denotes terminal node
##
## 1) root 800 1068.00 CH ( 0.61250 0.38750 )
## 2) LoyalCH < 0.5036 346 412.40 MM ( 0.28324 0.71676 )
## 4) LoyalCH < 0.280875 164 125.50 MM ( 0.12805 0.87195 )
## 8) LoyalCH < 0.0356415 56 10.03 MM ( 0.01786 0.98214 ) *
## 9) LoyalCH > 0.0356415 108 103.50 MM ( 0.18519 0.81481 ) *
## 5) LoyalCH > 0.280875 182 248.00 MM ( 0.42308 0.57692 )
## 10) PriceDiff < 0.05 71 67.60 MM ( 0.18310 0.81690 ) *
## 11) PriceDiff > 0.05 111 151.30 CH ( 0.57658 0.42342 ) *
## 3) LoyalCH > 0.5036 454 362.00 CH ( 0.86344 0.13656 )
## 6) PriceDiff < -0.39 31 40.32 MM ( 0.35484 0.64516 )
## 12) LoyalCH < 0.638841 10 0.00 MM ( 0.00000 1.00000 ) *
## 13) LoyalCH > 0.638841 21 29.06 CH ( 0.52381 0.47619 ) *
## 7) PriceDiff > -0.39 423 273.70 CH ( 0.90071 0.09929 )
## 14) LoyalCH < 0.705326 135 143.00 CH ( 0.77778 0.22222 )
## 28) ListPriceDiff < 0.255 67 89.49 CH ( 0.61194 0.38806 ) *
## 29) ListPriceDiff > 0.255 68 30.43 CH ( 0.94118 0.05882 ) *
## 15) LoyalCH > 0.705326 288 99.77 CH ( 0.95833 0.04167 ) *I am considering the node 29 where the root is split into nodes using
variable ListPriceDiff greater than .255 with 68 observation in it. The
prediction is made that Citrus Hill is being purchased with the deviance
of 30.43 and it is observed to be correct on 94.11% of the
observation.
(d) Create a plot of the tree, and interpret the
results.
plot(oj.tree)
text(oj.tree,pretty=0)
LoyalCH, which represents the brand loyalty of CH’s clients, is the
splitting variable at the top. As a result, it appears that LoyalCH is
the most crucial variable in this situation, followed by PriceDiff and
ListPriceDiff. Customers are more inclined to purchase Minute Maid if
they are less devoted to Citrus Hill. It takes a greater price gap to
persuade the more devoted Citrus Hill Customers to leave the
company.
(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?
oj.pred=predict(oj.tree,newdata=test.OJ,type="class")
table(test.OJ$Purchase,oj.pred)## oj.pred
## CH MM
## CH 148 15
## MM 32 75mean(oj.pred!=test.OJ$Purchase)## [1] 0.1740741The test error rate is 17.41%
(f) Apply the cv.tree() function to the training set in order to
determine the optimal tree size.
cv.oj=cv.tree(oj.tree,FUN=prune.misclass)
cv.oj## $size
## [1] 9 6 5 3 2 1
##
## $dev
## [1] 155 156 156 172 170 310
##
## $k
## [1] -Inf 0.0 1.0 8.5 9.0 150.0
##
## $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(cv.oj$size,cv.oj$dev,type="b",xlab="Tree Size",ylab="Deviance")
(h) Which tree size corresponds to the lowest cross-validated
classification error rate?
The lowest cross validated classification error rate is a tree with tree
size as 9 nodes.
(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.
oj.pruned=prune.misclass(oj.tree,best=9)
plot(oj.pruned)
text(oj.pruned,pretty=0)
(j) Compare the training error rates between the pruned and unpruned trees. Which is higher?
summary(oj.pruned)##
## Classification tree:
## tree(formula = Purchase ~ ., data = train.OJ)
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff" "ListPriceDiff"
## Number of terminal nodes: 9
## Residual mean deviance: 0.7347 = 581.1 / 791
## Misclassification error rate: 0.1662 = 133 / 800Both the trimmed and unpruned trees have the same test error rate.
The optimal size was found to be 9 terminal nodes using Cross validated
classification.
(k) Compare the test error rates between the pruned and unpruned
trees. Which is higher?
ojprune.pred=predict(oj.pruned, newdata=test.OJ, type="class")
mean(ojprune.pred != test.OJ$Purchase)## [1] 0.1740741As a result of the cross-validated classification finding that the unpruned tree was the best, we once again obtain the same test error rate.
detach(OJ)