Assignment #7

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.

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

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(ISLR2)
library(tree)

set.seed(1)
train <- sample(1:nrow(Carseats), nrow(Carseats) / 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.09900 

plot(tree_carseats)
text(tree_carseats, pretty = 0)

tree_pred <- predict(tree_carseats, Carseats_test)
tree_test_error <- mean((tree_pred - Carseats_test$Sales)^2)
tree_test_error

[1] 4.922039

The MSE is 4.922.

(c)

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_carseats <- cv.tree(tree_carseats)
cv_carseats

$size
 [1] 18 17 16 15 14 13 12 11 10  8  7  6  5  4  3  2  1

$dev
 [1]  984.3936 1031.3372 1036.0021 1027.2166 1027.2166 1055.8168
 [7] 1044.6955 1061.0899 1061.0899 1225.5973 1221.3487 1219.0219
[13] 1231.6886 1337.3952 1300.0524 1338.3702 1605.0221

$k
 [1]      -Inf  16.99544  20.56322  25.01730  25.57104  28.01938
 [7]  30.36962  31.56747  31.80816  40.75445  44.44673  52.57126
[13]  76.21881  99.59459 116.69889 159.79501 337.60153

$method
[1] "deviance"

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

plot(cv_carseats$size, cv_carseats$dev, type = "b")

best_size <- cv_carseats$size[which.min(cv_carseats$dev)]
best_size

[1] 18

pruned_tree <- prune.tree(tree_carseats, best = best_size)
plot(pruned_tree)
text(pruned_tree, pretty = 0)

pruned_pred <- predict(pruned_tree, Carseats_test)
pruned_test_error <- mean((pruned_pred - Carseats_test$Sales)^2)
pruned_test_error

[1] 4.922039

There is no improvement to MSE from pruning.

(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)

set.seed(1)
bag_carseats <- randomForest(Sales ~ ., data = Carseats_train, mtry = ncol(Carseats_train) - 1, importance = TRUE)
bag_pred <- predict(bag_carseats, Carseats_test)
bag_test_error <- mean((bag_pred - Carseats_test$Sales)^2)
bag_test_error

[1] 2.605253

The MSE obtained is 2.605 which is a lot lower than the previous MSE of 4.922.

importance(bag_carseats)

               %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

The importance measures show that Price and ShelveLoc are the most important, followed by CompPrice and Age.

varImpPlot(bag_carseats)

(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.

set.seed(1)
rf_carseats <- randomForest(Sales ~ ., data = Carseats_train, importance = TRUE)
rf_pred <- predict(rf_carseats, Carseats_test)
rf_test_error <- mean((rf_pred - Carseats_test$Sales)^2)
rf_test_error

[1] 2.960559

importance(rf_carseats)

               %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

varImpPlot(rf_carseats)

The MSE is 2.961, which is not slightly higher compared to the previous MSE of 2.605. Price and ShelveLoc are again the most important variables, followed by Age and CompPrice.

As m or mtry increases from small values towards p, the MSE decreases since Price, ShelveLoc, CompPrice, and Age drives Sales.

(f)

Now analyze the data using BART, and report your results.

library(dbarts)


x_train <- Carseats_train[, -which(names(Carseats_train) == "Sales")]
y_train <- Carseats_train$Sales
x_test <- Carseats_test[, -which(names(Carseats_test) == "Sales")]
y_test <- Carseats_test$Sales

set.seed(1)
bart_carseats <- bart(x_train, y_train, x.test = x_test)

Running BART with numeric y

number of trees: 200
number of chains: 1, default number of threads 1
tree thinning rate: 1
Prior:
    k prior fixed to 2.000000
    degrees of freedom in sigma prior: 3.000000
    quantile in sigma prior: 0.900000
    scale in sigma prior: 0.000964
    power and base for tree prior: 2.000000 0.950000
    use quantiles for rule cut points: false
    proposal probabilities: birth/death 0.50, swap 0.10, change 0.40; birth 0.50
data:
    number of training observations: 200
    number of test observations: 200
    number of explanatory variables: 12
    init sigma: 1.088371, curr sigma: 1.088371

Cutoff rules c in x<=c vs x>c
Number of cutoffs: (var: number of possible c):
(1: 100) (2: 100) (3: 100) (4: 100) (5: 100) 
(6: 100) (7: 100) (8: 100) (9: 100) (10: 100) 
(11: 100) (12: 100) 
Running mcmc loop:
iteration: 100 (of 1000)
iteration: 200 (of 1000)
iteration: 300 (of 1000)
iteration: 400 (of 1000)
iteration: 500 (of 1000)
iteration: 600 (of 1000)
iteration: 700 (of 1000)
iteration: 800 (of 1000)
iteration: 900 (of 1000)
iteration: 1000 (of 1000)
total seconds in loop: 0.517027

Tree sizes, last iteration:
[1] 2 2 2 3 2 2 1 2 2 2 2 3 2 2 2 1 2 2 
3 3 3 2 2 3 5 2 2 2 2 2 2 2 1 3 2 4 5 1 
2 3 4 2 3 2 2 2 3 2 2 2 2 3 3 4 2 2 4 3 
3 2 2 2 2 3 2 2 2 3 2 2 2 2 2 2 3 2 2 2 
3 3 3 1 2 2 2 2 2 2 2 2 2 3 2 4 2 2 3 3 
2 1 3 3 1 2 2 3 3 2 2 2 2 2 3 2 2 1 2 1 
2 2 3 3 2 2 2 2 2 2 2 3 2 2 2 4 2 3 2 2 
2 2 2 2 2 2 2 5 2 2 2 3 2 2 3 2 3 2 3 2 
2 3 3 2 3 4 3 5 2 3 2 2 2 1 2 3 2 2 2 3 
1 2 2 3 3 3 2 5 3 2 2 2 2 4 2 2 4 2 2 2 
2 4 

Variable Usage, last iteration (var:count):
(1: 21) (2: 28) (3: 32) (4: 18) (5: 35) 
(6: 24) (7: 18) (8: 15) (9: 24) (10: 21) 
(11: 20) (12: 15) 
DONE BART

bart_pred <- bart_carseats$yhat.test.mean
bart_test_error <- mean((bart_pred - y_test)^2)
bart_test_error

[1] 1.470076

BART gives a MSE of 1.470 which is the best so far as it outperforms the the single tree (4.922), bagging (2.605), and random forest (2.961).

Question 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.

summary(OJ)

 Purchase WeekofPurchase     StoreID        PriceCH         PriceMM     
 CH:653   Min.   :227.0   Min.   :1.00   Min.   :1.690   Min.   :1.690  
 MM:417   1st Qu.:240.0   1st Qu.:2.00   1st Qu.:1.790   1st Qu.:1.990  
          Median :257.0   Median :3.00   Median :1.860   Median :2.090  
          Mean   :254.4   Mean   :3.96   Mean   :1.867   Mean   :2.085  
          3rd Qu.:268.0   3rd Qu.:7.00   3rd Qu.:1.990   3rd Qu.:2.180  
          Max.   :278.0   Max.   :7.00   Max.   :2.090   Max.   :2.290  
     DiscCH            DiscMM         SpecialCH        SpecialMM     
 Min.   :0.00000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
 1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
 Median :0.00000   Median :0.0000   Median :0.0000   Median :0.0000  
 Mean   :0.05186   Mean   :0.1234   Mean   :0.1477   Mean   :0.1617  
 3rd Qu.:0.00000   3rd Qu.:0.2300   3rd Qu.:0.0000   3rd Qu.:0.0000  
 Max.   :0.50000   Max.   :0.8000   Max.   :1.0000   Max.   :1.0000  
    LoyalCH          SalePriceMM     SalePriceCH      PriceDiff      
 Min.   :0.000011   Min.   :1.190   Min.   :1.390   Min.   :-0.6700  
 1st Qu.:0.325257   1st Qu.:1.690   1st Qu.:1.750   1st Qu.: 0.0000  
 Median :0.600000   Median :2.090   Median :1.860   Median : 0.2300  
 Mean   :0.565782   Mean   :1.962   Mean   :1.816   Mean   : 0.1465  
 3rd Qu.:0.850873   3rd Qu.:2.130   3rd Qu.:1.890   3rd Qu.: 0.3200  
 Max.   :0.999947   Max.   :2.290   Max.   :2.090   Max.   : 0.6400  
 Store7      PctDiscMM        PctDiscCH       ListPriceDiff  
 No :714   Min.   :0.0000   Min.   :0.00000   Min.   :0.000  
 Yes:356   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.140  
           Median :0.0000   Median :0.00000   Median :0.240  
           Mean   :0.0593   Mean   :0.02731   Mean   :0.218  
           3rd Qu.:0.1127   3rd Qu.:0.00000   3rd Qu.:0.300  
           Max.   :0.4020   Max.   :0.25269   Max.   :0.440  
     STORE      
 Min.   :0.000  
 1st Qu.:0.000  
 Median :2.000  
 Mean   :1.631  
 3rd Qu.:3.000  
 Max.   :4.000  

set.seed(1)
train <- sample(1:nrow(OJ), 800)
OJ_train <- OJ[train, ]
OJ_test <- 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?

tree_oj <- tree(Purchase ~ ., data = OJ_train)
summary(tree_oj)

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

The classification tree shows 9 terminal nodes and a misclassification error rate of 0.1588, where 127 of the 800 are misclassified. It uses five predictors of LoyalCH, PriceDiff, SpecialCH, ListPriceDiff, and PctDiscMM, and LoyalCH. The residual mean deviance of 0.7432 is the training fit.

(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 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 ) *

Terminal node 8 of LoyalCH < 0.0356 has 59 observations with a deviance of 10.14 shows customers with no brand loyalty towards Citrus Hill. The tree predicts Minute Maid where customers with lower loyalty are likely to buy Minute Maid instead.

(d)

Create a plot of the tree, and interpret the results.

plot(tree_oj)
text(tree_oj, pretty = 0)

The split on LoyalCH < 0.5036 shows that it is a leading predictor while splitting customers into Citrus Hill and Minute Maid.

(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?

tree_pred <- predict(tree_oj, OJ_test, type = "class")
table(OJ_test$Purchase, tree_pred)

    tree_pred
      CH  MM
  CH 160   8
  MM  38  64

mean(tree_pred != OJ_test$Purchase)

[1] 0.1703704

Out of the 270 test observations, the tree correctly classifies 160 true CH purchases and 64 true MM purchases and misclassifies 8 actual CH buyers as MM and 38 actual MM buyers as CH. This results in an error rate of 0.1704.

(f)

Apply the cv.tree() function to the training set in order to determine the optimal tree size.

set.seed(1)
cv_oj <- cv.tree(tree_oj, 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)

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 = "Classification Error Rate")

(h)

Which tree size corresponds to the lowest cross-validated classification error rate?

best_size <- cv_oj$size[which.min(cv_oj$dev)]
best_size

[1] 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.

prune_oj <- prune.misclass(tree_oj, best = 5)
summary(prune_oj)

Classification tree:
snip.tree(tree = tree_oj, nodes = c(4L, 10L))
Variables actually used in tree construction:
[1] "LoyalCH"       "PriceDiff"     "ListPriceDiff" "PctDiscMM"    
Number of terminal nodes:  7 
Residual mean deviance:  0.7748 = 614.4 / 793 
Misclassification error rate: 0.1625 = 130 / 800 

plot(prune_oj)
text(prune_oj, pretty = 0)

Since finding the best size using cross-validated classification resulted in 9, it did not lead to the selection of a pruned tree. Therefore, a pruned tree with five terminal nodes is performed. However, as seen while performing the cv.tree() function, a tree size of 5 is not shown in the list, so the terminal node is defaulted to 7.

(j)

Compare the training error rates between the pruned and unpruned trees. Which is higher?

Comparing the two trees, the unpruned tree with 9 terminal nodes has a misclassification error of 0.1588. However, the pruned tree with 7 terminal nodes has a misclassification error of 0.1625, which is higher.

(k)

Compare the test error rates between the pruned and unpruned trees. Which is higher?

prune_pred <- predict(prune_oj, OJ_test, type = "class")
table(OJ_test$Purchase, prune_pred)

    prune_pred
      CH  MM
  CH 160   8
  MM  36  66

mean(prune_pred != OJ_test$Purchase)

[1] 0.162963

The unpruned tree (tree_pred) has a higher error rate of 0.1704 than the pruned tree of 0.1630. This is the opposite result we saw for the training error where the pruned tree was higher.