MReyes_Assignment7
Matthew Reyes
2026-07-17
Question 3
Consider the Gini index, classification error, and entropy in a simple classification setting with two classes. Create a single plot that dis- plays 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.
ANSWER: See below.
df <- tibble(
p = seq(0, 1, length.out = 500),
gini = 2 * p * (1 - p),
classification_error = 1 - pmax(p, 1 - p),
entropy = -(p * log2(p) + (1 - p) * log2(1 - p))
)
df$entropy[is.na(df$entropy)] <- 0
df_long <- df %>%
pivot_longer(cols = c(gini, classification_error, entropy),
names_to = "metric",
values_to = "value")
ggplot(df_long, aes(x = p, y = value, color = metric)) +
geom_line(linewidth = 1.2) +
labs(
title = "Gini Index, Classification Error, and Entropy",
x = expression(hat(p)[m1]),
y = "Value"
) +
theme_minimal() +
scale_color_manual(values = c(
"gini" = "blue",
"classification_error" = "red",
"entropy" = "green"
))
Question 8
- Split the data set into a training set and a test set.
ANSWER:
set.seed(123)
trainIndex <- createDataPartition(Carseats$Sales, p = 0.7, list = FALSE)
trainingData <- Carseats[trainIndex, ]
testingData <- Carseats[-trainIndex, ]
cat("Training data dimensions:", dim(trainingData), "\n")## Training data dimensions: 281 11cat("Testing data dimensions:", dim(testingData), "\n")## Testing data dimensions: 119 11trainingData$ShelveLoc <- factor(trainingData$ShelveLoc,
levels = levels(Carseats$ShelveLoc))
testingData$ShelveLoc <- factor(testingData$ShelveLoc,
levels = levels(Carseats$ShelveLoc))- Fit a regression tree to the training set. Plot the tree, and inter- pret the results. What test MSE do you obtain?
ANSWER: MSE is 5.717; RMSE is 2.391. It looks like the shelf location is perhaps the most important factor here, followed by price. If the shelf location is good, it likely has better sales regardless of price.
trControl <- trainControl(
method = "cv",
number = 10,
summaryFunction = defaultSummary # regression metrics: RMSE, Rsquared, MAE
)
set.seed(123)
model1 <- train(
Sales ~ .,
data = trainingData,
method = "rpart",
trControl = trControl,
metric = "RMSE" # FIXED: "MSE" is invalid
)## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo,
## : There were missing values in resampled performance measures.print(model1)## CART
##
## 281 samples
## 10 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 253, 253, 253, 253, 253, 253, ...
## Resampling results across tuning parameters:
##
## cp RMSE Rsquared MAE
## 0.04877915 2.28276 0.3396923 1.856057
## 0.09667155 2.36194 0.2949274 1.935341
## 0.26851707 2.56672 0.2177807 2.093814
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was cp = 0.04877915.# Plot the pruned tree (caret already pruned it)
rpart.plot(
model1$finalModel,
extra = 1,
fallen.leaves = TRUE,
type = 4,
tweak = 1.2
)
# Correct prediction using caret's built-in preprocessing
test_pred <- predict(model1, newdata = testingData)
# Test MSE
test_mse <- mean((test_pred - testingData$Sales)^2)
test_mse## [1] 5.717196# Test RMSE
test_rmse <- sqrt(test_mse)
test_rmse## [1] 2.391066- Use cross-validation in order to determine the optimal level of tree complexity. Does pruning the tree improve the test MSE?
ANSWER: When building the tree manually using cross-validation and selecting the optimal complexity, pruning accordingly didn’t help the MSE. In fact, it made it slightly worse.
full_tree <- rpart(
Sales ~ .,
data = trainingData,
method = "anova",
cp = 0 # grow a very large tree
)
printcp(full_tree)##
## Regression tree:
## rpart(formula = Sales ~ ., data = trainingData, method = "anova",
## cp = 0)
##
## Variables actually used in tree construction:
## [1] Advertising Age CompPrice Income Population Price
## [7] ShelveLoc
##
## Root node error: 2161.1/281 = 7.6908
##
## n= 281
##
## CP nsplit rel error xerror xstd
## 1 0.2685171 0 1.00000 1.00805 0.082881
## 2 0.0966716 1 0.73148 0.73838 0.058729
## 3 0.0487792 2 0.63481 0.70683 0.055831
## 4 0.0410861 3 0.58603 0.69813 0.055828
## 5 0.0334642 4 0.54495 0.66428 0.055514
## 6 0.0305383 6 0.47802 0.63987 0.052161
## 7 0.0303102 7 0.44748 0.61903 0.053080
## 8 0.0261170 9 0.38686 0.60958 0.051521
## 9 0.0249188 10 0.36074 0.57020 0.049867
## 10 0.0193015 11 0.33582 0.57333 0.049355
## 11 0.0132563 12 0.31652 0.57558 0.048428
## 12 0.0104637 13 0.30327 0.57284 0.052928
## 13 0.0076318 14 0.29280 0.56138 0.052557
## 14 0.0064169 15 0.28517 0.55628 0.052604
## 15 0.0063174 16 0.27875 0.55594 0.051926
## 16 0.0061120 17 0.27244 0.55933 0.052039
## 17 0.0038542 18 0.26632 0.55640 0.051910
## 18 0.0030167 19 0.26247 0.54852 0.052225
## 19 0.0020082 20 0.25945 0.55337 0.052457
## 20 0.0000000 21 0.25744 0.55334 0.052441best_cp <- full_tree$cptable[which.min(full_tree$cptable[,"xerror"]), "CP"]
best_cp## [1] 0.003016734pruned_tree <- prune(full_tree, cp = best_cp)
full_pred <- predict(full_tree, newdata = testingData)
full_mse <- mean((full_pred - testingData$Sales)^2)
pruned_pred <- predict(pruned_tree, newdata = testingData)
pruned_mse <- mean((pruned_pred - testingData$Sales)^2)
full_mse## [1] 4.403387pruned_mse## [1] 4.564439- Use the bagging approach in order to analyze this data. What test MSE do you obtain? Use the importance() function to de- termine which variables are most important.
ANSWER: MSE obtained after bagging is 2.5 (RMSE is 1.6).
set.seed(123)
p <- ncol(trainingData) - 1 # number of predictors
bag_model <- randomForest(
Sales ~ .,
data = trainingData,
mtry = p, # Bagging: use all predictors at each split
importance = TRUE
)
bag_model##
## Call:
## randomForest(formula = Sales ~ ., data = trainingData, mtry = p, importance = TRUE)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 10
##
## Mean of squared residuals: 2.605684
## % Var explained: 66.12bag_pred <- predict(bag_model, newdata = testingData)
bag_mse <- mean((bag_pred - testingData$Sales)^2)
bag_mse## [1] 2.585357bag_rmse <- sqrt(bag_mse)
bag_rmse## [1] 1.607905importance(bag_model)## %IncMSE IncNodePurity
## CompPrice 33.3349479 244.394085
## Income 5.8631159 102.875273
## Advertising 18.7483257 159.776528
## Population -0.3998659 76.500819
## Price 62.2476821 588.283743
## ShelveLoc 71.2961762 697.155907
## Age 15.6811225 158.021524
## Education 1.7440222 61.685166
## Urban -0.8846815 10.808252
## US 4.2855745 9.000114varImpPlot(bag_model)
- Use random forests to analyze this data. What test MSE do you obtain? Use the importance() function to determine which vari- ables are most important. Describe the effect of m, the number of variables considered at each split, on the error rate obtained.
ANSWER: Test MSE obtained is 2.94 (RMSE = 1.71). Fascinating that this appears to be slightly worse than the bagging approach. With bagging, m included all predictors, whereas here it
set.seed(123)
rf_model <- randomForest(
Sales ~ .,
data = trainingData,
importance = TRUE,
mtry = sqrt(p)
)
rf_model##
## Call:
## randomForest(formula = Sales ~ ., data = trainingData, importance = TRUE, mtry = sqrt(p))
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 3
##
## Mean of squared residuals: 2.869925
## % Var explained: 62.68rf_pred <- predict(rf_model, newdata = testingData)
rf_mse <- mean((rf_pred - testingData$Sales)^2)
rf_mse## [1] 2.82199rf_rmse <- sqrt(rf_mse)
rf_rmse## [1] 1.679878importance(rf_model)## %IncMSE IncNodePurity
## CompPrice 16.8460289 209.53729
## Income 3.6498653 160.61398
## Advertising 14.0779764 184.98404
## Population -0.6665519 134.46929
## Price 39.5357070 474.23443
## ShelveLoc 51.2651437 540.69756
## Age 10.9667323 201.35870
## Education -0.1728094 85.50893
## Urban -0.5513495 19.39221
## US 2.7594254 23.52524varImpPlot(rf_model)
- Now analyze the data using BART, and report your results.
ANSWER: The MSE is 1.7, which is significantly better than the previous MSEs.
set.seed(123)
train <- createDataPartition(Carseats$Sales, p = 0.7, list = FALSE)
x <- Carseats[, 2:10]
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: 281, 12, 119
## y1,yn: 3.727544, 2.217544
## x1,x[n*p]: 111.000000, 1.000000
## xp1,xp[np*p]: 138.000000, 1.000000
## *****Number of Trees: 200
## *****Number of Cut Points: 70 ... 1
## *****burn,nd,thin: 100,1000,1
## *****Prior:beta,alpha,tau,nu,lambda,offset: 2,0.95,0.276302,3,0.192449,7.49246
## *****sigma: 0.993970
## *****w (weights): 1.000000 ... 1.000000
## *****Dirichlet:sparse,theta,omega,a,b,rho,augment: 0,0,1,0.5,1,12,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: 1s
## trcnt,tecnt: 1000,1000yhat.bart <- bartfit$yhat.test.mean
mean((ytest - yhat.bart)^2)## [1] 1.705505Question 9
- Create a training set containing a random sample of 800 obser- vations, and a test set containing the remaining observations.
ANSWER: See below.
train_index <- createDataPartition(OJ$Purchase, p = 800/nrow(OJ), list = FALSE)
trainingData <- OJ[train_index, ]
testingData <- OJ[-train_index, ]
nrow(trainingData)## [1] 801nrow(testingData)## [1] 269- 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?
ANSWER: The tree model had an expected loss = 0.158 at the first node and 5 terminal nodes. The training error rate is about 15.4%.
tree_model1 <- rpart(
Purchase ~ .,
data = trainingData,
method = "class"
)
summary(tree_model1)## Call:
## rpart(formula = Purchase ~ ., data = trainingData, method = "class")
## n= 801
##
## CP nsplit rel error xerror xstd
## 1 0.53525641 0 1.0000000 1.0000000 0.04423447
## 2 0.02083333 1 0.4647436 0.5096154 0.03618222
## 3 0.01000000 4 0.3974359 0.4551282 0.03464304
##
## Variable importance
## LoyalCH PriceDiff SalePriceMM DiscMM PctDiscMM
## 63 7 7 5 5
## PriceMM WeekofPurchase StoreID ListPriceDiff PriceCH
## 4 4 2 2 1
## STORE Store7
## 1 1
##
## Node number 1: 801 observations, complexity param=0.5352564
## predicted class=CH expected loss=0.3895131 P(node) =1
## class counts: 489 312
## probabilities: 0.610 0.390
## left son=2 (504 obs) right son=3 (297 obs)
## Primary splits:
## LoyalCH < 0.482304 to the right, improve=144.79180, (0 missing)
## StoreID < 3.5 to the right, improve= 42.59574, (0 missing)
## PriceDiff < 0.015 to the right, improve= 26.12144, (0 missing)
## SalePriceMM < 1.84 to the right, improve= 23.67953, (0 missing)
## Store7 splits as RL, improve= 20.98707, (0 missing)
## Surrogate splits:
## StoreID < 3.5 to the right, agree=0.640, adj=0.030, (0 split)
## PriceMM < 1.89 to the right, agree=0.638, adj=0.024, (0 split)
## WeekofPurchase < 228.5 to the right, agree=0.635, adj=0.017, (0 split)
## DiscMM < 0.57 to the left, agree=0.634, adj=0.013, (0 split)
## SalePriceMM < 1.385 to the right, agree=0.634, adj=0.013, (0 split)
##
## Node number 2: 504 observations, complexity param=0.02083333
## predicted class=CH expected loss=0.1587302 P(node) =0.6292135
## class counts: 424 80
## probabilities: 0.841 0.159
## left son=4 (297 obs) right son=5 (207 obs)
## Primary splits:
## LoyalCH < 0.705699 to the right, improve=15.901960, (0 missing)
## PriceDiff < 0.015 to the right, improve=14.167510, (0 missing)
## SalePriceMM < 1.84 to the right, improve=13.289290, (0 missing)
## ListPriceDiff < 0.255 to the right, improve= 8.614616, (0 missing)
## DiscMM < 0.72 to the left, improve= 7.598318, (0 missing)
## Surrogate splits:
## WeekofPurchase < 237.5 to the right, agree=0.643, adj=0.130, (0 split)
## PriceCH < 1.755 to the right, agree=0.637, adj=0.116, (0 split)
## PriceMM < 2.04 to the right, agree=0.635, adj=0.111, (0 split)
## SalePriceMM < 1.64 to the right, agree=0.615, adj=0.063, (0 split)
## PctDiscMM < 0.1961965 to the left, agree=0.609, adj=0.048, (0 split)
##
## Node number 3: 297 observations
## predicted class=MM expected loss=0.2188552 P(node) =0.3707865
## class counts: 65 232
## probabilities: 0.219 0.781
##
## Node number 4: 297 observations
## predicted class=CH expected loss=0.05387205 P(node) =0.3707865
## class counts: 281 16
## probabilities: 0.946 0.054
##
## Node number 5: 207 observations, complexity param=0.02083333
## predicted class=CH expected loss=0.3091787 P(node) =0.258427
## class counts: 143 64
## probabilities: 0.691 0.309
## left son=10 (140 obs) right son=11 (67 obs)
## Primary splits:
## PriceDiff < 0.015 to the right, improve=16.414890, (0 missing)
## SalePriceMM < 1.84 to the right, improve=14.798450, (0 missing)
## ListPriceDiff < 0.255 to the right, improve=12.547930, (0 missing)
## PctDiscMM < 0.1155095 to the left, improve= 7.417401, (0 missing)
## DiscMM < 0.15 to the left, improve= 7.155280, (0 missing)
## Surrogate splits:
## SalePriceMM < 1.84 to the right, agree=0.952, adj=0.851, (0 split)
## PctDiscMM < 0.1155095 to the left, agree=0.894, adj=0.672, (0 split)
## DiscMM < 0.27 to the left, agree=0.874, adj=0.612, (0 split)
## WeekofPurchase < 234.5 to the right, agree=0.787, adj=0.343, (0 split)
## PriceMM < 2.04 to the right, agree=0.787, adj=0.343, (0 split)
##
## Node number 10: 140 observations
## predicted class=CH expected loss=0.1714286 P(node) =0.1747815
## class counts: 116 24
## probabilities: 0.829 0.171
##
## Node number 11: 67 observations, complexity param=0.02083333
## predicted class=MM expected loss=0.4029851 P(node) =0.08364544
## class counts: 27 40
## probabilities: 0.403 0.597
## left son=22 (24 obs) right son=23 (43 obs)
## Primary splits:
## ListPriceDiff < 0.235 to the right, improve=5.200046, (0 missing)
## PriceDiff < -0.165 to the right, improve=3.425619, (0 missing)
## DiscMM < 0.47 to the left, improve=2.395140, (0 missing)
## PctDiscMM < 0.227263 to the left, improve=2.395140, (0 missing)
## StoreID < 2.5 to the right, improve=1.464612, (0 missing)
## Surrogate splits:
## PriceDiff < -0.165 to the right, agree=0.746, adj=0.292, (0 split)
## StoreID < 5.5 to the right, agree=0.731, adj=0.250, (0 split)
## Store7 splits as RL, agree=0.731, adj=0.250, (0 split)
## STORE < 0.5 to the left, agree=0.731, adj=0.250, (0 split)
## PriceCH < 1.755 to the left, agree=0.701, adj=0.167, (0 split)
##
## Node number 22: 24 observations
## predicted class=CH expected loss=0.3333333 P(node) =0.02996255
## class counts: 16 8
## probabilities: 0.667 0.333
##
## Node number 23: 43 observations
## predicted class=MM expected loss=0.255814 P(node) =0.0536829
## class counts: 11 32
## probabilities: 0.256 0.744training_pred <- predict(tree_model1, trainingData, type = "class")
mean(training_pred != trainingData$Purchase)## [1] 0.1548065- 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.
ANSWER: Looking at the first split, we see that loyalty to CH is the most important predictor of which item a consumer purchases. Those with loyalty is greater than or equal to .48, they mostly buy CH. Otherwise, they go for MM.
tree_model1## n= 801
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 801 312 CH (0.61048689 0.38951311)
## 2) LoyalCH>=0.482304 504 80 CH (0.84126984 0.15873016)
## 4) LoyalCH>=0.705699 297 16 CH (0.94612795 0.05387205) *
## 5) LoyalCH< 0.705699 207 64 CH (0.69082126 0.30917874)
## 10) PriceDiff>=0.015 140 24 CH (0.82857143 0.17142857) *
## 11) PriceDiff< 0.015 67 27 MM (0.40298507 0.59701493)
## 22) ListPriceDiff>=0.235 24 8 CH (0.66666667 0.33333333) *
## 23) ListPriceDiff< 0.235 43 11 MM (0.25581395 0.74418605) *
## 3) LoyalCH< 0.482304 297 65 MM (0.21885522 0.78114478) *- Create a plot of the tree, and interpret the results.
ANSWER: The most significant predictor is loyalty to CH, as described above. Next to that, the second major predictor is the difference in price between MM and CH.
rpart.plot(tree_model1)
- 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?
ANSWER: The test error rate is about 21%.
test_pred <- predict(tree_model1, testingData, type = "class")
mean(test_pred != testingData$Purchase)## [1] 0.2081784confusionMatrix(test_pred, testingData$Purchase)## Confusion Matrix and Statistics
##
## Reference
## Prediction CH MM
## CH 129 21
## MM 35 84
##
## Accuracy : 0.7918
## 95% CI : (0.7383, 0.8387)
## No Information Rate : 0.6097
## P-Value [Acc > NIR] : 1.275e-10
##
## Kappa : 0.5728
##
## Mcnemar's Test P-Value : 0.08235
##
## Sensitivity : 0.7866
## Specificity : 0.8000
## Pos Pred Value : 0.8600
## Neg Pred Value : 0.7059
## Prevalence : 0.6097
## Detection Rate : 0.4796
## Detection Prevalence : 0.5576
## Balanced Accuracy : 0.7933
##
## 'Positive' Class : CH
## test_error <- mean(test_pred != testingData$Purchase)
test_error## [1] 0.2081784- Apply the cv.tree() function to the training set in order to determine the optimal tree size.
ANSWER: Optimal tree size is 8.
tree_model2 <- tree(Purchase ~ ., data = trainingData)
cv_oj <- cv.tree(tree_model2, FUN = prune.misclass)
best_size <- cv_oj$size[which.min(cv_oj$dev)]
best_size## [1] 8- Produce a plot with tree size on the x-axis and cross-validated classification error rate on the y-axis.
ANSWER:
plot(cv_oj$size,
cv_oj$dev,
type = "b",
xlab = "Tree Size (number of terminal nodes)",
ylab = "Cross‑Validated Classification Error",
main = "Cross‑Validation Results for Tree Size")
- Which tree size corresponds to the lowest cross-validated classi- fication error rate?
ANSWER: As shown above, it looks like 8 is the optimal size based on the cross-validated error rate. However, 5 also looks relatively good based on the depiction.
- 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.
ANSWER: See below.
pruned_tree <- prune.misclass(tree_model2, best = best_size)
plot(pruned_tree)
text(pruned_tree, pretty = 0)
- Compare the training error rates between the pruned and un- pruned trees. Which is higher?
ANSWER: The training error rates are identical… I have no idea why.
train_pred_unpruned <- predict(tree_model2, trainingData, type = "class")
train_error_unpruned <- mean(train_pred_unpruned != trainingData$Purchase)
train_error_unpruned## [1] 0.1573034train_pred_pruned <- predict(pruned_tree, trainingData, type = "class")
train_error_pruned <- mean(train_pred_pruned != trainingData$Purchase)
train_error_pruned## [1] 0.1573034- Compare the test error rates between the pruned and unpruned trees. Which is higher?
ANSWER: Again, the results are identical, and I have no idea why.
test_pred_unpruned <- predict(tree_model2, testingData, type = "class")
test_error_unpruned <- mean(test_pred_unpruned != testingData$Purchase)
test_error_unpruned## [1] 0.2007435test_pred_pruned <- predict(pruned_tree, testingData, type = "class")
test_error_pruned <- mean(test_pred_pruned != testingData$Purchase)
test_error_pruned## [1] 0.2007435