Question 3

p=seq(0,1,0.0001)
#Gini
G=2*p*(1-p)
#Classification Error
E=1-pmax(p,1-p)
#Entropy
D=-(p*log(p) + (1-p)*log(1-p))

plot(p,D, col="red",ylab="")
lines(p,E,col='green')
lines(p,G,col='blue')
legend(0.3,0.15,c("Entropy", "Missclassification","Gini"),lty=c(1,1,1),lwd=c(2.5,2.5,2.5),col=c('red','green','blue'))


Question 8

a

library(ISLR2)
set.seed(1)

train_idx <- sample(nrow(Carseats), nrow(Carseats) / 2)
train_set <- Carseats[train_idx, ]
test_set <- Carseats[-train_idx, ]

b

library(tree)

tree.fit <- tree(Sales ~ ., data = train_set)
plot(tree.fit)
text(tree.fit, pretty = 0)


tree.pred <- predict(tree.fit, test_set)
tree.mse <- mean((test_set$Sales - tree.pred)^2)
cat("Regression Tree Test MSE:", tree.mse, "\n")
Regression Tree Test MSE: 4.922039 

Test MSE is 4.92.

c

cv.carseats <- cv.tree(tree.fit)
plot(cv.carseats$size, cv.carseats$dev, type = "b", xlab = "Tree Size", ylab = "Deviance")


best_size <- cv.carseats$size[which.min(cv.carseats$dev)]
cat("Optimal Tree Size by CV:", best_size, "\n")
Optimal Tree Size by CV: 18 
prune.carseats <- prune.tree(tree.fit, best = best_size)
plot(prune.carseats)
text(prune.carseats, pretty = 0)


prune.pred <- predict(prune.carseats, test_set)
prune.mse <- mean((test_set$Sales - prune.pred)^2)
cat("Pruned Tree Test MSE:", prune.mse, "\n")
Pruned Tree Test MSE: 4.922039 

Cross-validation selects a tree size of 18. Since the optimal size is the full tree, pruning does not alter the tree and the test MSE remains 4.92.

d

library(randomForest)
set.seed(1)

bag.fit <- randomForest(Sales ~ ., data = train_set, mtry = 10, importance = TRUE)
bag.pred <- predict(bag.fit, test_set)
bag.mse <- mean((test_set$Sales - bag.pred)^2)
cat("Bagging Test MSE:", bag.mse, "\n")
Bagging Test MSE: 2.605253 
importance(bag.fit)
               %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

Test MSE with bagging is 2.61. Price and ShelveLoc are the two most important variables.

e

set.seed(1)

# Fit Random Forest (mtry = p/3 = 10/3 approx 3)
rf.fit <- randomForest(Sales ~ ., data = train_set, mtry = 3, importance = TRUE)
rf.pred <- predict(rf.fit, test_set)
rf.mse <- mean((test_set$Sales - rf.pred)^2)
cat("Random Forest Test MSE:", rf.mse, "\n")
Random Forest Test MSE: 2.960559 
importance(rf.fit)
               %IncMSE IncNodePurity
CompPrice   14.8840765     158.82956
Income       4.3293950     125.64850
Advertising  8.2215192     107.51700
Population  -0.9488134      97.06024
Price       34.9793386     385.93142
ShelveLoc   34.9248499     298.54210
Age         14.3055912     178.42061
Education    1.3117842      70.49202
Urban       -1.2680807      17.39986
US           6.1139696      33.98963

Test MSE with Random Forest (mtry = 3) is 2.96. Increasing the value of mtry towards 10 (bagging) decreases the test MSE from 2.96 to 2.61. Price and ShelveLoc remain the most important predictors.

f

library(gbm)
set.seed(1)

boost.fit <- gbm(Sales ~ ., data = train_set, distribution = "gaussian", n.trees = 1000, shrinkage = 0.01)
boost.pred <- predict(boost.fit, test_set, n.trees = 1000)
boost.mse <- mean((test_set$Sales - boost.pred)^2)
cat("Boosting Test MSE:", boost.mse, "\n")
Boosting Test MSE: 2.033312 

Test MSE with Boosting is 2.03.


Question 9

a

set.seed(1)

train_idx <- sample(nrow(OJ), 800)
train_set <- OJ[train_idx, ]
test_set <- OJ[-train_idx, ]

b

tree.fit <- tree(Purchase ~ ., data = train_set)
summary(tree.fit)

Classification tree:
tree(formula = Purchase ~ ., data = train_set)
Variables actually used in tree construction:
[1] "LoyalCH"       "PriceDiff"     "SpecialCH"     "ListPriceDiff" "PctDiscMM"    
Number of terminal nodes:  9 
Residual mean deviance:  0.7432 = 587.8 / 791 
Misclassification error rate: 0.1588 = 127 / 800 

The training error rate is 15.88% (127 misclassified out of 800) and the tree has 9 terminal nodes.

c

tree.fit
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 ) *

Looking at node 4 (split on LoyalCH < 0.280875): it has 177 observations, a deviance of 140.5 and 86.4% of the observations in this branch choose MM.

d

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

e

tree.pred <- predict(tree.fit, test_set, type = "class")
table(tree.pred, test_set$Purchase)
         
tree.pred  CH  MM
       CH 160  38
       MM   8  64
test_err <- mean(tree.pred != test_set$Purchase)
cat("Test Error Rate:", test_err, "\n")
Test Error Rate: 0.1703704 

f

set.seed(1)
cv.oj <- cv.tree(tree.fit, FUN = prune.misclass)
cv.oj
$size
[1] 9 8 7 4 2 1

$dev
[1] 145 145 146 146 167 315

$k
[1]       -Inf   0.000000   3.000000   4.333333  10.500000 151.000000

$method
[1] "misclass"

attr(,"class")
[1] "prune"         "tree.sequence"

g

plot(cv.oj$size, cv.oj$dev, type = "b", xlab = "Tree Size", ylab = "Deviance")

h

Tree size 8 corresponds to the lowest deviance.

i

prune.oj <- prune.tree(tree.fit, best = 6, method = "misclass")
plot(prune.oj)
text(prune.oj, pretty = 0)

j

unpruned_train_err <- summary(tree.fit)$misclass[1] / summary(tree.fit)$misclass[2]
pruned_train_err <- summary(prune.oj)$misclass[1] / summary(prune.oj)$misclass[2]

cat("Unpruned Training Error Rate:", unpruned_train_err, "\n")
Unpruned Training Error Rate: 0.15875 
cat("Pruned Training Error Rate:", pruned_train_err, "\n")
Pruned Training Error Rate: 0.1625 

k

prune.pred <- predict(prune.oj, test_set, type = "class")
prune_test_err <- mean(prune.pred != test_set$Purchase)

cat("Unpruned Test Error Rate:", test_err, "\n")
Unpruned Test Error Rate: 0.1703704 
cat("Pruned Test Error Rate:", prune_test_err, "\n")
Pruned Test Error Rate: 0.162963 
LS0tCnRpdGxlOiAiQ2hyaXMgU2VycmFubyAtIEFzc2lnbm1lbnQgNyIKb3V0cHV0OgogIGh0bWxfbm90ZWJvb2s6CiAgICB0b2M6IHRydWUKICAgIHRvY19mbG9hdDogdHJ1ZQogICAgZWNobzogdHJ1ZQotLS0KCiMjIFF1ZXN0aW9uIDMKCiMjIyAKCmBgYHtyfQpwPXNlcSgwLDEsMC4wMDAxKQojR2luaQpHPTIqcCooMS1wKQojQ2xhc3NpZmljYXRpb24gRXJyb3IKRT0xLXBtYXgocCwxLXApCiNFbnRyb3B5CkQ9LShwKmxvZyhwKSArICgxLXApKmxvZygxLXApKQoKcGxvdChwLEQsIGNvbD0icmVkIix5bGFiPSIiKQpsaW5lcyhwLEUsY29sPSdncmVlbicpCmxpbmVzKHAsRyxjb2w9J2JsdWUnKQpsZWdlbmQoMC4zLDAuMTUsYygiRW50cm9weSIsICJNaXNzY2xhc3NpZmljYXRpb24iLCJHaW5pIiksbHR5PWMoMSwxLDEpLGx3ZD1jKDIuNSwyLjUsMi41KSxjb2w9YygncmVkJywnZ3JlZW4nLCdibHVlJykpCmBgYAoKLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tCgojIyBRdWVzdGlvbiA4CgojIyMgYQoKYGBge3J9CmxpYnJhcnkoSVNMUjIpCnNldC5zZWVkKDEpCgp0cmFpbl9pZHggPC0gc2FtcGxlKG5yb3coQ2Fyc2VhdHMpLCBucm93KENhcnNlYXRzKSAvIDIpCnRyYWluX3NldCA8LSBDYXJzZWF0c1t0cmFpbl9pZHgsIF0KdGVzdF9zZXQgPC0gQ2Fyc2VhdHNbLXRyYWluX2lkeCwgXQpgYGAKCiMjIyBiCgpgYGB7cn0KbGlicmFyeSh0cmVlKQoKdHJlZS5maXQgPC0gdHJlZShTYWxlcyB+IC4sIGRhdGEgPSB0cmFpbl9zZXQpCnBsb3QodHJlZS5maXQpCnRleHQodHJlZS5maXQsIHByZXR0eSA9IDApCgp0cmVlLnByZWQgPC0gcHJlZGljdCh0cmVlLmZpdCwgdGVzdF9zZXQpCnRyZWUubXNlIDwtIG1lYW4oKHRlc3Rfc2V0JFNhbGVzIC0gdHJlZS5wcmVkKV4yKQpjYXQoIlJlZ3Jlc3Npb24gVHJlZSBUZXN0IE1TRToiLCB0cmVlLm1zZSwgIlxuIikKYGBgCgpUZXN0IE1TRSBpcyBgciByb3VuZCh0cmVlLm1zZSwgMilgLgoKIyMjIGMKCmBgYHtyfQpjdi5jYXJzZWF0cyA8LSBjdi50cmVlKHRyZWUuZml0KQpwbG90KGN2LmNhcnNlYXRzJHNpemUsIGN2LmNhcnNlYXRzJGRldiwgdHlwZSA9ICJiIiwgeGxhYiA9ICJUcmVlIFNpemUiLCB5bGFiID0gIkRldmlhbmNlIikKCmJlc3Rfc2l6ZSA8LSBjdi5jYXJzZWF0cyRzaXplW3doaWNoLm1pbihjdi5jYXJzZWF0cyRkZXYpXQpjYXQoIk9wdGltYWwgVHJlZSBTaXplIGJ5IENWOiIsIGJlc3Rfc2l6ZSwgIlxuIikKCnBydW5lLmNhcnNlYXRzIDwtIHBydW5lLnRyZWUodHJlZS5maXQsIGJlc3QgPSBiZXN0X3NpemUpCnBsb3QocHJ1bmUuY2Fyc2VhdHMpCnRleHQocHJ1bmUuY2Fyc2VhdHMsIHByZXR0eSA9IDApCgpwcnVuZS5wcmVkIDwtIHByZWRpY3QocHJ1bmUuY2Fyc2VhdHMsIHRlc3Rfc2V0KQpwcnVuZS5tc2UgPC0gbWVhbigodGVzdF9zZXQkU2FsZXMgLSBwcnVuZS5wcmVkKV4yKQpjYXQoIlBydW5lZCBUcmVlIFRlc3QgTVNFOiIsIHBydW5lLm1zZSwgIlxuIikKYGBgCgpDcm9zcy12YWxpZGF0aW9uIHNlbGVjdHMgYSB0cmVlIHNpemUgb2YgYHIgYmVzdF9zaXplYC4gU2luY2UgdGhlIG9wdGltYWwgc2l6ZSBpcyB0aGUgZnVsbCB0cmVlLCBwcnVuaW5nIGRvZXMgbm90IGFsdGVyIHRoZSB0cmVlIGFuZCB0aGUgdGVzdCBNU0UgcmVtYWlucyBgciByb3VuZChwcnVuZS5tc2UsIDIpYC4KCiMjIyBkCgpgYGB7cn0KbGlicmFyeShyYW5kb21Gb3Jlc3QpCnNldC5zZWVkKDEpCgpiYWcuZml0IDwtIHJhbmRvbUZvcmVzdChTYWxlcyB+IC4sIGRhdGEgPSB0cmFpbl9zZXQsIG10cnkgPSAxMCwgaW1wb3J0YW5jZSA9IFRSVUUpCmJhZy5wcmVkIDwtIHByZWRpY3QoYmFnLmZpdCwgdGVzdF9zZXQpCmJhZy5tc2UgPC0gbWVhbigodGVzdF9zZXQkU2FsZXMgLSBiYWcucHJlZCleMikKY2F0KCJCYWdnaW5nIFRlc3QgTVNFOiIsIGJhZy5tc2UsICJcbiIpCgppbXBvcnRhbmNlKGJhZy5maXQpCmBgYAoKVGVzdCBNU0Ugd2l0aCBiYWdnaW5nIGlzIGByIHJvdW5kKGJhZy5tc2UsIDIpYC4gUHJpY2UgYW5kIFNoZWx2ZUxvYyBhcmUgdGhlIHR3byBtb3N0IGltcG9ydGFudCB2YXJpYWJsZXMuCgojIyMgZQoKYGBge3J9CnNldC5zZWVkKDEpCgojIEZpdCBSYW5kb20gRm9yZXN0IChtdHJ5ID0gcC8zID0gMTAvMyBhcHByb3ggMykKcmYuZml0IDwtIHJhbmRvbUZvcmVzdChTYWxlcyB+IC4sIGRhdGEgPSB0cmFpbl9zZXQsIG10cnkgPSAzLCBpbXBvcnRhbmNlID0gVFJVRSkKcmYucHJlZCA8LSBwcmVkaWN0KHJmLmZpdCwgdGVzdF9zZXQpCnJmLm1zZSA8LSBtZWFuKCh0ZXN0X3NldCRTYWxlcyAtIHJmLnByZWQpXjIpCmNhdCgiUmFuZG9tIEZvcmVzdCBUZXN0IE1TRToiLCByZi5tc2UsICJcbiIpCgppbXBvcnRhbmNlKHJmLmZpdCkKYGBgCgpUZXN0IE1TRSB3aXRoIFJhbmRvbSBGb3Jlc3QgKG10cnkgPSAzKSBpcyBgciByb3VuZChyZi5tc2UsIDIpYC4gSW5jcmVhc2luZyB0aGUgdmFsdWUgb2YgbXRyeSB0b3dhcmRzIDEwIChiYWdnaW5nKSBkZWNyZWFzZXMgdGhlIHRlc3QgTVNFIGZyb20gYHIgcm91bmQocmYubXNlLCAyKWAgdG8gYHIgcm91bmQoYmFnLm1zZSwgMilgLiBQcmljZSBhbmQgU2hlbHZlTG9jIHJlbWFpbiB0aGUgbW9zdCBpbXBvcnRhbnQgcHJlZGljdG9ycy4KCiMjIyBmCgpgYGB7cn0KbGlicmFyeShnYm0pCnNldC5zZWVkKDEpCgpib29zdC5maXQgPC0gZ2JtKFNhbGVzIH4gLiwgZGF0YSA9IHRyYWluX3NldCwgZGlzdHJpYnV0aW9uID0gImdhdXNzaWFuIiwgbi50cmVlcyA9IDEwMDAsIHNocmlua2FnZSA9IDAuMDEpCmJvb3N0LnByZWQgPC0gcHJlZGljdChib29zdC5maXQsIHRlc3Rfc2V0LCBuLnRyZWVzID0gMTAwMCkKYm9vc3QubXNlIDwtIG1lYW4oKHRlc3Rfc2V0JFNhbGVzIC0gYm9vc3QucHJlZCleMikKY2F0KCJCb29zdGluZyBUZXN0IE1TRToiLCBib29zdC5tc2UsICJcbiIpCmBgYAoKVGVzdCBNU0Ugd2l0aCBCb29zdGluZyBpcyBgciByb3VuZChib29zdC5tc2UsIDIpYC4KCi0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQoKIyMgUXVlc3Rpb24gOQoKIyMjIGEKCmBgYHtyfQpzZXQuc2VlZCgxKQoKdHJhaW5faWR4IDwtIHNhbXBsZShucm93KE9KKSwgODAwKQp0cmFpbl9zZXQgPC0gT0pbdHJhaW5faWR4LCBdCnRlc3Rfc2V0IDwtIE9KWy10cmFpbl9pZHgsIF0KYGBgCgojIyMgYgoKYGBge3J9CnRyZWUuZml0IDwtIHRyZWUoUHVyY2hhc2UgfiAuLCBkYXRhID0gdHJhaW5fc2V0KQpzdW1tYXJ5KHRyZWUuZml0KQpgYGAKClRoZSB0cmFpbmluZyBlcnJvciByYXRlIGlzIDE1Ljg4JSAoMTI3IG1pc2NsYXNzaWZpZWQgb3V0IG9mIDgwMCkgYW5kIHRoZSB0cmVlIGhhcyA5IHRlcm1pbmFsIG5vZGVzLgoKIyMjIGMKCmBgYHtyfQp0cmVlLmZpdApgYGAKCkxvb2tpbmcgYXQgbm9kZSA0IChzcGxpdCBvbiBMb3lhbENIIFw8IDAuMjgwODc1KTogaXQgaGFzIDE3NyBvYnNlcnZhdGlvbnMsIGEgZGV2aWFuY2Ugb2YgMTQwLjUgYW5kIDg2LjQlIG9mIHRoZSBvYnNlcnZhdGlvbnMgaW4gdGhpcyBicmFuY2ggY2hvb3NlIE1NLgoKIyMjIGQKCmBgYHtyfQpwbG90KHRyZWUuZml0KQp0ZXh0KHRyZWUuZml0LCBwcmV0dHkgPSAwKQpgYGAKCiMjIyBlCgpgYGB7cn0KdHJlZS5wcmVkIDwtIHByZWRpY3QodHJlZS5maXQsIHRlc3Rfc2V0LCB0eXBlID0gImNsYXNzIikKdGFibGUodHJlZS5wcmVkLCB0ZXN0X3NldCRQdXJjaGFzZSkKCnRlc3RfZXJyIDwtIG1lYW4odHJlZS5wcmVkICE9IHRlc3Rfc2V0JFB1cmNoYXNlKQpjYXQoIlRlc3QgRXJyb3IgUmF0ZToiLCB0ZXN0X2VyciwgIlxuIikKYGBgCgojIyMgZgoKYGBge3J9CnNldC5zZWVkKDEpCmN2Lm9qIDwtIGN2LnRyZWUodHJlZS5maXQsIEZVTiA9IHBydW5lLm1pc2NsYXNzKQpjdi5vagpgYGAKCiMjIyBnCgpgYGB7cn0KcGxvdChjdi5vaiRzaXplLCBjdi5vaiRkZXYsIHR5cGUgPSAiYiIsIHhsYWIgPSAiVHJlZSBTaXplIiwgeWxhYiA9ICJEZXZpYW5jZSIpCmBgYAoKIyMjIGgKClRyZWUgc2l6ZSA4IGNvcnJlc3BvbmRzIHRvIHRoZSBsb3dlc3QgZGV2aWFuY2UuCgojIyMgaQoKYGBge3J9CnBydW5lLm9qIDwtIHBydW5lLnRyZWUodHJlZS5maXQsIGJlc3QgPSA2LCBtZXRob2QgPSAibWlzY2xhc3MiKQpwbG90KHBydW5lLm9qKQp0ZXh0KHBydW5lLm9qLCBwcmV0dHkgPSAwKQpgYGAKCiMjIyBqCgpgYGB7cn0KdW5wcnVuZWRfdHJhaW5fZXJyIDwtIHN1bW1hcnkodHJlZS5maXQpJG1pc2NsYXNzWzFdIC8gc3VtbWFyeSh0cmVlLmZpdCkkbWlzY2xhc3NbMl0KcHJ1bmVkX3RyYWluX2VyciA8LSBzdW1tYXJ5KHBydW5lLm9qKSRtaXNjbGFzc1sxXSAvIHN1bW1hcnkocHJ1bmUub2opJG1pc2NsYXNzWzJdCgpjYXQoIlVucHJ1bmVkIFRyYWluaW5nIEVycm9yIFJhdGU6IiwgdW5wcnVuZWRfdHJhaW5fZXJyLCAiXG4iKQpjYXQoIlBydW5lZCBUcmFpbmluZyBFcnJvciBSYXRlOiIsIHBydW5lZF90cmFpbl9lcnIsICJcbiIpCmBgYAoKIyMjIGsKCmBgYHtyfQpwcnVuZS5wcmVkIDwtIHByZWRpY3QocHJ1bmUub2osIHRlc3Rfc2V0LCB0eXBlID0gImNsYXNzIikKcHJ1bmVfdGVzdF9lcnIgPC0gbWVhbihwcnVuZS5wcmVkICE9IHRlc3Rfc2V0JFB1cmNoYXNlKQoKY2F0KCJVbnBydW5lZCBUZXN0IEVycm9yIFJhdGU6IiwgdGVzdF9lcnIsICJcbiIpCmNhdCgiUHJ1bmVkIFRlc3QgRXJyb3IgUmF0ZToiLCBwcnVuZV90ZXN0X2VyciwgIlxuIikKYGBgCg==