Assignment 7
Jordan Miller
2023-04-29
Question 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, ˆpm1 = 1 − ˆpm2. You could make this plot by hand, but it will be much easier to make in R.
pm1 <- seq(0, 1, 0.001)
pm2 <- 1 - pm1
error.rate <- 1 - pmax(pm1, pm2)
Gini <- pm1*pm2 + pm2*pm1
entropy <- -(pm1*log(pm1) + pm2*log(pm2))
data.frame(pm1, pm2, error.rate, Gini, entropy) %>%
pivot_longer(cols = c(error.rate, Gini, entropy), names_to = "error_method") %>%
ggplot(aes(x=pm1, y=value, color = error_method)) +
geom_line(linewidth = 1)## Warning: Removed 2 rows containing missing values (`geom_line()`).
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.
data("Carseats")8(a)
Split the data set into a training set and a test set.
# 60-40 Train-Test
train <- sample(400, 240)
test <- -train8(b)
Fit a regression tree to the training set. Plot the tree, and interpret the results. What test MSE do you obtain?
#train regression tree
CStree <- tree(Sales~., Carseats[train, ])
summary(CStree)##
## Regression tree:
## tree(formula = Sales ~ ., data = Carseats[train, ])
## Variables actually used in tree construction:
## [1] "ShelveLoc" "Price" "Income" "Advertising" "CompPrice"
## [6] "Age"
## Number of terminal nodes: 15
## Residual mean deviance: 2.259 = 508.2 / 225
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -4.4200 -1.0180 -0.0170 0.0000 0.9121 3.7990plot(CStree)
text(CStree, pretty = 0)
This regression tree produces 15 terminal nodes using 6 predictor variables.
yhat <- predict(CStree, newdata = Carseats[test, ])
CStree.test <- Carseats[test, "Sales"]
mean.CStree <- mean((yhat - CStree.test)^2)The test MSE of this regression tree is 5.0572.
8(c)
Use cross-validation in order to determine the optimal level of tree complexity. Does pruning the tree improve the test MSE?
cv.CStree <- cv.tree(CStree)
plot(cv.CStree$size , cv.CStree$dev, type = "b")
which.min(cv.CStree$size)## [1] 15Cross-validation does not prune the tree, keeping the number of terminal nodes at 15. However, improvements associated with increasing tree complexity beyond 8 terminal nodes appear modest.
We now investigate an 8-leaf model.
prune.CStree <- prune.tree(CStree, best = 8)
summary(prune.CStree)##
## Regression tree:
## snip.tree(tree = CStree, nodes = c(12L, 9L, 13L, 23L))
## Variables actually used in tree construction:
## [1] "ShelveLoc" "Price" "Income" "CompPrice"
## Number of terminal nodes: 8
## Residual mean deviance: 3.006 = 697.4 / 232
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -4.42000 -1.17000 -0.05784 0.00000 1.12000 4.08800plot(prune.CStree)
text(prune.CStree, pretty = 0)
This tree keeps only ShelveLoc, Price,
Income, and CompPrice as predictor
variables.
prune.yhat <- predict(prune.CStree, newdata=Carseats[test, ])
mean.prune <- mean((prune.yhat - CStree.test)^2)The test MSE of the 8-leaf tree is 5.44, which is not an improvement over the full, unpruned tree’s test error.
8(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.
bag.CStree <- randomForest(Sales~., data=Carseats, subset=train, mtry=10, importance=TRUE)
bag.CStree##
## Call:
## randomForest(formula = Sales ~ ., data = Carseats, mtry = 10, importance = TRUE, subset = train)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 10
##
## Mean of squared residuals: 2.679698
## % Var explained: 63.78yhat.bag <- predict(bag.CStree, newdata=Carseats[test, ])
mean.bag <- mean((yhat.bag - CStree.test)^2)
importance(bag.CStree)## %IncMSE IncNodePurity
## CompPrice 25.5820443 174.923376
## Income 5.6892468 89.385644
## Advertising 13.6945416 114.444139
## Population -0.1677107 72.058644
## Price 56.9047136 509.803832
## ShelveLoc 64.4515564 583.218798
## Age 15.6762128 132.136996
## Education -0.3552936 38.800815
## Urban -4.2025251 5.563751
## US 2.1276512 10.436122varImpPlot(bag.CStree)
The test MSE of the bagged tree model is 3.0006, which is a
significant improvement over the single regression tree. By far, the two
most important variables are ShelveLoc and
Price, which is consistent with the previous single tree
and pruned single tree models, although the apparent importance of the
remaining variables differ between bagged and single tree models.
8(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.CStree <- randomForest(Sales~., data=Carseats, subset=train, mtry=3, importance=TRUE)
RF.CStree##
## Call:
## randomForest(formula = Sales ~ ., data = Carseats, mtry = 3, importance = TRUE, subset = train)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 3
##
## Mean of squared residuals: 3.008398
## % Var explained: 59.34yhat.RF <- predict(RF.CStree, newdata=Carseats[test, ])
mean.RF <- mean((yhat.RF - CStree.test)^2)
importance(RF.CStree)## %IncMSE IncNodePurity
## CompPrice 15.1118874 158.91691
## Income 2.4510668 116.90314
## Advertising 15.0743360 168.71548
## Population 1.4381645 121.54422
## Price 33.7767316 389.78618
## ShelveLoc 40.5723337 398.40938
## Age 11.8226599 188.49007
## Education -0.1986517 77.44019
## Urban -0.1720398 14.60344
## US 6.0858179 26.96585varImpPlot(RF.CStree)
The test MSE of the random forest model is 3.6318, so this method
does not appear to be an improvement over bagging. However, the
importance of ShelveLoc and Price compared to
all other variables is corroborated by the random forest method.
RF.1 <- randomForest(Sales~., data=Carseats, subset=train, mtry=1, importance=TRUE)
yhat.RF1 <- predict(RF.1, newdata=Carseats[test, ])
mean.RF1 <- mean((yhat.RF1 - CStree.test)^2)
RF.5 <- randomForest(Sales~., data=Carseats, subset=train, mtry=5, importance=TRUE)
yhat.RF5 <- predict(RF.5, newdata=Carseats[test, ])
mean.RF5 <- mean((yhat.RF5 - CStree.test)^2)
RF.7 <- randomForest(Sales~., data=Carseats, subset=train, mtry=7, importance=TRUE)
yhat.RF7 <- predict(RF.7, newdata=Carseats[test, ])
mean.RF7 <- mean((yhat.RF7 - CStree.test)^2)
RF.9 <- randomForest(Sales~., data=Carseats, subset=train, mtry=9, importance=TRUE)
yhat.RF9 <- predict(RF.9, newdata=Carseats[test, ])
mean.RF9 <- mean((yhat.RF9 - CStree.test)^2)
m <- seq(1, 9, 2)
RF_error_rate <- c(mean.RF1, mean.RF, mean.RF5, mean.RF7, mean.RF9)
plot(m, RF_error_rate, type = "b")
As m increases, the error rate decreases.
8(f)
Now analyze the data using BART, and report your results.
x <- Carseats[, 2:11]
y <- Carseats[, "Sales"]
xtrain <- x[train, ]
ytrain <- y[train]
xtest <- x[-train, ]
ytest <- y[-train]
bartfit <- gbart(xtrain , ytrain , x.test = xtest)## *****Calling gbart: type=1
## *****Data:
## data:n,p,np: 240, 14, 160
## y1,yn: -1.026750, -2.736750
## x1,x[n*p]: 125.000000, 1.000000
## xp1,xp[np*p]: 138.000000, 1.000000
## *****Number of Trees: 200
## *****Number of Cut Points: 67 ... 1
## *****burn,nd,thin: 100,1000,1
## *****Prior:beta,alpha,tau,nu,lambda,offset: 2,0.95,0.281075,3,0.220944,7.54675
## *****sigma: 1.065015
## *****w (weights): 1.000000 ... 1.000000
## *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,14,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: 6s
## trcnt,tecnt: 1000,1000yhat.bart <- bartfit$yhat.test.mean
mean.bart <- mean((ytest - yhat.bart)^2)The BART test MSE is 1.4084, suggesting that BART is the decision tree method of choice for this problem.
Problem 9
This problem involves the OJ data set which is part
of the ISLR2 package.
data(OJ)9(a)
Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations.
OJtrain <- sample(1070, 800)
OJtest <- -OJtrain9(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?
OJtree <- tree(Purchase~., OJ[OJtrain, ])
summary(OJtree)##
## Classification tree:
## tree(formula = Purchase ~ ., data = OJ[OJtrain, ])
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff" "ListPriceDiff"
## Number of terminal nodes: 7
## Residual mean deviance: 0.739 = 586 / 793
## Misclassification error rate: 0.1638 = 131 / 800The misclassification error rate of this tree is 0.1638. The tree has 7 terminal nodes.
9(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.
OJtree## node), split, n, deviance, yval, (yprob)
## * denotes terminal node
##
## 1) root 800 1064.00 CH ( 0.61750 0.38250 )
## 2) LoyalCH < 0.5036 342 404.10 MM ( 0.27778 0.72222 )
## 4) LoyalCH < 0.276142 165 105.10 MM ( 0.09697 0.90303 ) *
## 5) LoyalCH > 0.276142 177 243.30 MM ( 0.44633 0.55367 )
## 10) PriceDiff < 0.05 76 80.79 MM ( 0.22368 0.77632 ) *
## 11) PriceDiff > 0.05 101 134.70 CH ( 0.61386 0.38614 ) *
## 3) LoyalCH > 0.5036 458 351.90 CH ( 0.87118 0.12882 )
## 6) LoyalCH < 0.705699 162 189.50 CH ( 0.72840 0.27160 )
## 12) ListPriceDiff < 0.255 80 110.10 CH ( 0.55000 0.45000 ) *
## 13) ListPriceDiff > 0.255 82 52.43 CH ( 0.90244 0.09756 ) *
## 7) LoyalCH > 0.705699 296 118.70 CH ( 0.94932 0.05068 )
## 14) PriceDiff < -0.39 12 16.30 CH ( 0.58333 0.41667 ) *
## 15) PriceDiff > -0.39 284 86.57 CH ( 0.96479 0.03521 ) *If the customer’s brand loyalty to Citrus Hill is less than 0.0276142, they are predicted to choose Minute Maid.
9(d)
Create a plot of the tree, and interpret the results.
plot(OJtree)
text(OJtree, pretty = 0)
9(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(OJtree, OJ[OJtest,], type = "class")
table(OJ.pred, OJ$Purchase[OJtest])##
## OJ.pred CH MM
## CH 143 44
## MM 16 67OJerror <- (16+44)/(143+16+44+67)The test error rate for this tree is 0.22222.
9(f)
Apply the cv.tree() function to the training set
in order to determine the optimal tree size.
cv.OJ <- cv.tree(OJtree, FUN = prune.misclass)
cv.OJ## $size
## [1] 7 4 2 1
##
## $dev
## [1] 139 139 163 306
##
## $k
## [1] -Inf 0.0 11.5 152.0
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"9(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")
9(h)
Which tree size corresponds to the lowest cross-validated classification error rate?
The 7-leaf tree (equivalent to unpruned) and the 4-leaf tree both have the lowest cross-validated error rate.
9(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.
Let us consider the 4-leaf tree:
prune.OJ <- prune.misclass(OJtree, best = 4)9(j)
Compare the training error rates between the pruned and unpruned trees. Which is higher?
summary(prune.OJ)##
## Classification tree:
## snip.tree(tree = OJtree, nodes = 3L)
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff"
## Number of terminal nodes: 4
## Residual mean deviance: 0.8448 = 672.5 / 796
## Misclassification error rate: 0.1638 = 131 / 800The training error rate of the 4-leaf tree is 0.1638, which is equivalent to the unpruned tree’s training error rate.
9(k)
Compare the test error rates between the pruned and unpruned trees. Which is higher?
prune.OJ.pred <- predict(prune.OJ, OJ[OJtest,], type = "class")
table(prune.OJ.pred, OJ$Purchase[OJtest])##
## prune.OJ.pred CH MM
## CH 143 44
## MM 16 67prune.OJerror <- (16+44)/(143+16+44+67)The test error rates of both the pruned and unpruned trees are 0.22222.