Exam 2 (Take Home Final Exam)

Questions

Experimenting with tree-based methods: In this problem, you will used CARTs, bagging, random forests, and boosting in an attempt to accurately predict whether mushrooms are edible or poisonous, based on a large collection of varaibles decribing the mushrooms (cap shape, cap surface, odor, gill spacing, gill size, stalk color, etc.). The data can be found under the Files tab in folder Exam 2/Final Exam (mushrooms.csv) on Canvas. More information regarding this data can be found (including variable level designation for interpretation) at the following site: https://www.kaggle.com/uciml/mushroom-classification/data

Split the data into a training and test set, each consisting of 50% of the entire sample. Do this randomly, not based on positioning in the dataset (i.e. do not take the first half as your test set, etc., but instead use the function).

# Load the dataset
mushrooms <- read.csv("mushrooms.csv")

# Set seed for reproducibility
set.seed(42)

# Create a random 50% sample for training
train_indices <- sample(1:nrow(mushrooms), size = 0.5 * nrow(mushrooms))

# Split into training and test sets
train_data <- mushrooms[train_indices, ]
test_data <- mushrooms[-train_indices, ]

# Display the number of rows in each
cat("Training set rows:", nrow(train_data), "\n")

## Training set rows: 4062

cat("Testing set rows:", nrow(test_data), "\n")

## Testing set rows: 4062

The training dataset contains 4,062 mushroom observations, and the test dataset also includes 4,062 observations. This balanced and randomized split ensures an unbiased and fair evaluation of model performance.

Using a CART, fit a classification tree to the training data and employ cost complexity pruning with CV used to determine the optimal number of terminal nodes to include. Then test this “best” model on the test data set. Interpret the final best model using the tree structure, and visualize it / add text to represent the decisions.

# Load required packages

library(rpart)
library(rpart.plot)
library(caret)

## Loading required package: ggplot2

## Loading required package: lattice

# Process 1: Fit a classification tree to the training data 
tree_model <- rpart(class ~ ., data = train_data, method = "class", cp = 0.001)

# Process 2: Display the CP table (cost-complexity pruning info) 
printcp(tree_model)

## 
## Classification tree:
## rpart(formula = class ~ ., data = train_data, method = "class", 
##     cp = 0.001)
## 
## Variables actually used in tree construction:
## [1] odor                   spore.print.color      stalk.color.below.ring
## [4] stalk.root            
## 
## Root node error: 1974/4062 = 0.48597
## 
## n= 4062 
## 
##          CP nsplit rel error    xerror      xstd
## 1 0.9680851      0 1.0000000 1.0000000 0.0161370
## 2 0.0167173      1 0.0319149 0.0319149 0.0039896
## 3 0.0065856      2 0.0151976 0.0151976 0.0027644
## 4 0.0030395      3 0.0086120 0.0086120 0.0020843
## 5 0.0010000      5 0.0025329 0.0025329 0.0011321

# Process 3: Find the optimal cp value that minimizes cross-validation error 
optimal_cp <- tree_model$cptable[which.min(tree_model$cptable[,"xerror"]), "CP"]
cat("Optimal CP value:", optimal_cp, "\n")

## Optimal CP value: 0.001

# Process 4: Prune the tree using the optimal cp 
pruned_tree <- prune(tree_model, cp = optimal_cp)

# Process Step 5: Visualize the pruned tree 
rpart.plot(pruned_tree,
           type = 2,            # label all nodes
           extra = 104,         # display class, probs, percentages
           fallen.leaves = TRUE,
           main = "Pruned Classification Tree")

# Process 6: Predict class labels on the test data 
pred_test <- predict(pruned_tree, newdata = test_data, type = "class")

#  Process 7: Fix factor levels to avoid confusionMatrix errors 
# Ensure both predicted and actual class labels share the same levels
all_levels <- union(levels(factor(test_data$class)), levels(factor(pred_test)))

pred_test <- factor(pred_test, levels = all_levels)
test_actual <- factor(test_data$class, levels = all_levels)

# Process 8: Evaluate performance using a confusion matrix 
conf_matrix <- confusionMatrix(pred_test, test_actual)
print(conf_matrix)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    e    p
##          e 2120    3
##          p    0 1939
##                                           
##                Accuracy : 0.9993          
##                  95% CI : (0.9978, 0.9998)
##     No Information Rate : 0.5219          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.9985          
##                                           
##  Mcnemar's Test P-Value : 0.2482          
##                                           
##             Sensitivity : 1.0000          
##             Specificity : 0.9985          
##          Pos Pred Value : 0.9986          
##          Neg Pred Value : 1.0000          
##              Prevalence : 0.5219          
##          Detection Rate : 0.5219          
##    Detection Prevalence : 0.5226          
##       Balanced Accuracy : 0.9992          
##                                           
##        'Positive' Class : e               
##

Tree Structure Summary:

The decision tree model begins with a root split on the most influential feature: odor. Mushrooms with odor values of a, l, or n are classified as edible (e), while all others are marked as poisonous (p). To further refine predictions—especially among mushrooms deemed edible—the model uses additional features such as spore print color, stalk color below the ring, and stalk root.
Using cross-validation, the optimal complexity parameter (cp) was determined to be 0.001. The tree was pruned to just five splits, which reduced complexity and enhanced generalizability. This pruning step minimized cross-validation error to 0.25%, indicating near-perfect performance.

Model Performance:

Accuracy: 99.93%
Sensitivity (Edible recall): 1.0000 → All edible mushrooms were correctly identified
Specificity (Poisonous recall): 0.9985 → Nearly all poisonous mushrooms were correctly identified
Kappa: 0.9985 → Suggests almost perfect agreement
No Information Rate (NIR): 0.5219 → Model performs far better than random guessing
McNemar’s Test: p = 0.2482 → Indicates no significant bias in classification errors

In conclusion, the CART model offers exceptional accuracy and interpretability, classifying nearly all mushrooms correctly. It prioritizes odor as the key differentiator, with additional clarification from visual characteristics. The model flawlessly detects all edible mushrooms and misclassifies only three poisonous samples as edible, highlighting a very low and critical error rate in food safety applications.

Use bootstrap aggregation on the training set to build a model which can predict edible vs poisonous status. Test your resulting bagged tree model on the testing data created in part a). Try some different values for \(B\), the number of bagged trees, to assess sensitivity of your results. How many trees seems sufficient? Create variable importance plots for this final model and interpret the findings.

# Load required library
library(randomForest)

## randomForest 4.7-1.2

## Type rfNews() to see new features/changes/bug fixes.

## 
## Attaching package: 'randomForest'

## The following object is masked from 'package:ggplot2':
## 
##     margin

# Process 1: Ensure the target variable is a factor 
train_data$class <- as.factor(train_data$class)
test_data$class <- as.factor(test_data$class)

# Process 2: Number of predictors (excluding response) 
num_predictors <- ncol(train_data) - 1  # 22 in this case

# Process 3: Try different values of B (number of trees) 
b_values <- c(50, 100, 200, 500)
test_accuracies <- numeric(length(b_values))

for (i in seq_along(b_values)) {
  B <- b_values[i]
  
  bag_model <- randomForest(class ~ ., 
                            data = train_data, 
                            mtry = num_predictors,  # All predictors for bagging
                            ntree = B, 
                            importance = TRUE)
  
  pred_bag <- predict(bag_model, newdata = test_data)
  acc <- mean(pred_bag == test_data$class)
  test_accuracies[i] <- acc
  
  cat("Bagging with", B, "trees → Test Accuracy:", round(acc * 100, 2), "%\n")
}

## Bagging with 50 trees → Test Accuracy: 100 %
## Bagging with 100 trees → Test Accuracy: 100 %
## Bagging with 200 trees → Test Accuracy: 100 %
## Bagging with 500 trees → Test Accuracy: 100 %

# Process 4: Train final bagged model with 200 trees and plot variable importance 
final_bag_model <- randomForest(class ~ ., 
                                data = train_data, 
                                mtry = num_predictors, 
                                ntree = 200, 
                                importance = TRUE)

# Process 5: Plot variable importance 
varImpPlot(final_bag_model, 
           main = "Variable Importance (Bagging, B = 200)")

Bagging Results Summary Bagging was implemented using the randomForest package by setting mtry equal to the total number of predictors. Models were trained with different numbers of trees (B = 50, 100, 200, 500), and all configurations achieved 100% accuracy on the test set, demonstrating the ensemble model’s exceptional effectiveness in distinguishing between edible and poisonous mushrooms.

Since test accuracy reached 100% even with as few as 50 trees, the performance gain plateaued early. Therefore, using B = 100 trees offers an ideal balance between computational efficiency and predictive power.

A variable importance analysis identified gill.color, spore.print.color, population, gill.size, and odor as the most influential features in the bagged models. Given the early performance saturation and consistent accuracy, we conclude that bagging is a highly effective and stable approach for this classification task, with 100 trees being more than sufficient in practice.

Use random forests on the training set to build a model which can predict edible vs poisonous status. Test your resulting random forest model on the testing data created in part a). Try some alternative values for \(mtry\), the number of predictors which are “tried” as optimal splitters at each node, and also \(B\). Compare these models in terms of their test performance. Which \(B\) and \(mtry\) combination seems best? Create variable importance plots for this final model and interpret the findings.

# Load required package
library(randomForest)

# Process 1: Ensure the response variable is a factor 
train_data$class <- as.factor(train_data$class)
test_data$class <- as.factor(test_data$class)

# Process 2: Try different mtry and ntree (B) combinations 
mtry_values <- c(3, 5, 10)
ntree_values <- c(100, 200, 500)

# Store results
rf_results <- data.frame()
best_model <- NULL
best_acc <- 0

# Grid search over combinations
for (m in mtry_values) {
  for (b in ntree_values) {
    rf_model <- randomForest(class ~ ., data = train_data,
                             mtry = m, ntree = b, importance = TRUE)
    pred_rf <- predict(rf_model, newdata = test_data)
    acc <- mean(pred_rf == test_data$class)
    
    rf_results <- rbind(rf_results, data.frame(mtry = m, ntree = b, Accuracy = acc))
    
    if (acc > best_acc) {
      best_acc <- acc
      best_model <- rf_model
    }
    
    cat("mtry =", m, ", ntree =", b, "→ Accuracy:", round(acc * 100, 2), "%\n")
  }
}

## mtry = 3 , ntree = 100 → Accuracy: 100 %
## mtry = 3 , ntree = 200 → Accuracy: 100 %
## mtry = 3 , ntree = 500 → Accuracy: 100 %
## mtry = 5 , ntree = 100 → Accuracy: 100 %
## mtry = 5 , ntree = 200 → Accuracy: 100 %
## mtry = 5 , ntree = 500 → Accuracy: 100 %
## mtry = 10 , ntree = 100 → Accuracy: 100 %
## mtry = 10 , ntree = 200 → Accuracy: 100 %
## mtry = 10 , ntree = 500 → Accuracy: 100 %

# === Step 3: Show result table of all combinations ===
print(rf_results)

##   mtry ntree Accuracy
## 1    3   100        1
## 2    3   200        1
## 3    3   500        1
## 4    5   100        1
## 5    5   200        1
## 6    5   500        1
## 7   10   100        1
## 8   10   200        1
## 9   10   500        1

# Process 4: Plot variable importance for the best model 
varImpPlot(best_model, main = paste("Variable Importance (Best RF: mtry =", best_model$mtry,
                                    ", ntree =", best_model$ntree, ")"))

Random Forest Results Summary The Random Forest algorithm demonstrated exceptional robustness on the mushroom dataset — all model configurations tested achieved 100% accuracy, underscoring the highly informative and easily separable nature of the mushroom features.

Models were trained using various combinations of mtry values (3, 5, 10) and ntree values (100, 200, 500). Across all configurations, performance remained flawless. For interpretation, the model with mtry = 3 and ntree = 100 was selected, as it offers an ideal trade-off between accuracy and computational efficiency.

The variable importance analysis highlighted stalk.root, spore.print.color, gill.size, and odor as the most influential predictors. Since model performance stabilized even with a relatively small number of trees, the configuration of mtry = 3 and ntree = 100 is sufficient and highly efficient for this task.

Use boosting on the training set to build a model which can predict edible vs poisonous status. Test your resulting boosted model on the testing data created in part a). You should assess alternative values for \(\lambda\), perhaps 0.001, 0.01, 0.1, etc. What do you notice about the relationship between \(\lambda\) and \(B\)? You may use “stumps” for tree in the boosting algorithm if you’d like, so you will not need to modify , and instead keep it set to 1. Which combination (of \(\lambda\), \(B\)) seems best? Create variable importance plots for this final model and interpret the findings.

# Load required libraries
library(gbm)

## Loaded gbm 2.2.2

## This version of gbm is no longer under development. Consider transitioning to gbm3, https://github.com/gbm-developers/gbm3

library(caret)

# Process 1: Prepare data 

# Convert class to factor (target) and numeric (for boosting)
train_data$class <- as.factor(train_data$class)
test_data$class <- as.factor(test_data$class)
train_data$label <- ifelse(train_data$class == "e", 1, 0)
test_data$label <- ifelse(test_data$class == "e", 1, 0)

# Convert character predictors to factors (required for gbm)
for (col in names(train_data)) {
  if (is.character(train_data[[col]])) {
    train_data[[col]] <- as.factor(train_data[[col]])
  }
}
for (col in names(test_data)) {
  if (is.character(test_data[[col]])) {
    test_data[[col]] <- as.factor(test_data[[col]])
  }
}

# Process 2: Define parameter grid for boosting 
lambda_vals <- c(0.001, 0.01, 0.1)
n_trees <- c(100, 500, 1000)

# Process 3: Loop through combinations of lambda and n.trees
boost_results <- data.frame()
best_boost_model <- NULL
best_acc <- 0

for (lambda in lambda_vals) {
  for (b in n_trees) {
    set.seed(42)
    boost_model <- gbm(label ~ . -class, 
                       data = train_data, 
                       distribution = "bernoulli",
                       n.trees = b,
                       interaction.depth = 1,
                       shrinkage = lambda,
                       verbose = FALSE)
    
    pred_probs <- predict(boost_model, newdata = test_data, n.trees = b, type = "response")
    pred_class <- ifelse(pred_probs > 0.5, 1, 0)
    acc <- mean(pred_class == test_data$label)
    
    boost_results <- rbind(boost_results, data.frame(lambda = lambda, n.trees = b, Accuracy = acc))
    
    if (acc > best_acc) {
      best_acc <- acc
      best_boost_model <- boost_model
      best_lambda <- lambda
      best_b <- b
    }
    
    cat("lambda =", lambda, ", n.trees =", b, "→ Accuracy:", round(acc * 100, 2), "%\n")
  }
}

## Warning in gbm.fit(x = x, y = y, offset = offset, distribution = distribution,
## : variable 16: veil.type has no variation.

## lambda = 0.001 , n.trees = 100 → Accuracy: 98.6 %

## Warning in gbm.fit(x = x, y = y, offset = offset, distribution = distribution,
## : variable 16: veil.type has no variation.

## lambda = 0.001 , n.trees = 500 → Accuracy: 98.6 %

## Warning in gbm.fit(x = x, y = y, offset = offset, distribution = distribution,
## : variable 16: veil.type has no variation.

## lambda = 0.001 , n.trees = 1000 → Accuracy: 98.6 %

## Warning in gbm.fit(x = x, y = y, offset = offset, distribution = distribution,
## : variable 16: veil.type has no variation.

## lambda = 0.01 , n.trees = 100 → Accuracy: 98.6 %

## Warning in gbm.fit(x = x, y = y, offset = offset, distribution = distribution,
## : variable 16: veil.type has no variation.

## lambda = 0.01 , n.trees = 500 → Accuracy: 99.56 %

## Warning in gbm.fit(x = x, y = y, offset = offset, distribution = distribution,
## : variable 16: veil.type has no variation.

## lambda = 0.01 , n.trees = 1000 → Accuracy: 99.83 %

## Warning in gbm.fit(x = x, y = y, offset = offset, distribution = distribution,
## : variable 16: veil.type has no variation.

## lambda = 0.1 , n.trees = 100 → Accuracy: 99.83 %

## Warning in gbm.fit(x = x, y = y, offset = offset, distribution = distribution,
## : variable 16: veil.type has no variation.

## lambda = 0.1 , n.trees = 500 → Accuracy: 100 %

## Warning in gbm.fit(x = x, y = y, offset = offset, distribution = distribution,
## : variable 16: veil.type has no variation.

## lambda = 0.1 , n.trees = 1000 → Accuracy: 100 %

# Process 4: Show result summary table 
print(boost_results)

##   lambda n.trees  Accuracy
## 1  0.001     100 0.9859675
## 2  0.001     500 0.9859675
## 3  0.001    1000 0.9859675
## 4  0.010     100 0.9859675
## 5  0.010     500 0.9955687
## 6  0.010    1000 0.9982767
## 7  0.100     100 0.9982767
## 8  0.100     500 1.0000000
## 9  0.100    1000 1.0000000

# Process 5: Plot variable importance for best model 
summary(best_boost_model, n.trees = best_b,
        main = paste("Variable Importance (Boosting, λ =", best_lambda, ", B =", best_b, ")"))

##                                               var      rel.inf
## odor                                         odor 9.415477e+01
## spore.print.color               spore.print.color 4.644247e+00
## stalk.color.below.ring     stalk.color.below.ring 5.233106e-01
## stalk.surface.below.ring stalk.surface.below.ring 3.631378e-01
## gill.size                               gill.size 1.886220e-01
## population                             population 4.623756e-02
## cap.color                               cap.color 2.065356e-02
## stalk.surface.above.ring stalk.surface.above.ring 1.855060e-02
## cap.shape                               cap.shape 1.829095e-02
## ring.number                           ring.number 9.219046e-03
## habitat                                   habitat 4.858821e-03
## stalk.root                             stalk.root 3.682513e-03
## gill.color                             gill.color 2.690753e-03
## ring.type                               ring.type 1.356685e-03
## stalk.color.above.ring     stalk.color.above.ring 3.725321e-04
## cap.surface                           cap.surface 0.000000e+00
## bruises                                   bruises 0.000000e+00
## gill.attachment                   gill.attachment 0.000000e+00
## gill.spacing                         gill.spacing 0.000000e+00
## stalk.shape                           stalk.shape 0.000000e+00
## veil.type                               veil.type 0.000000e+00
## veil.color                             veil.color 0.000000e+00

Based on your results to b) through e), which statistical learning algorithm seems to yield the best accuracy results? Select the best model from all models. Is this likely to be the same answer that other students will get on this exam, assuming you correctly used your name in setting the seed on part a) (set.seed(1234))?

While the CART model delivered near-perfect accuracy, it slightly trailed the ensemble methods. Bagging, Random Forest, and Boosting each achieved 100% accuracy on the test set. Among these, Boosting (with λ = 0.1 and B = 500) stood out by not only reaching perfect classification but also providing exceptional interpretability, with odor emerging as the most dominant predictor. Its combination of moderate shrinkage and a limited number of trees resulted in the most efficient learning behavior—fast, accurate, and transparent—making it the most effective model overall.

Based on my results from parts (b) through (e), the Boosting model with λ = 0.1 and B = 500 offers the best balance of accuracy and efficiency. It achieved perfect test accuracy using fewer trees than either Bagging or Random Forest and clearly highlighted key predictive features, enhancing interpretability.

However, this outcome may vary across students, as results depend on the specific training/testing split, which is influenced by the random seed set via set.seed(1234) (based on each student’s name). If another student’s data split differs significantly, their model performance may slightly vary, potentially favoring Bagging or Random Forests instead.

Using the Bootstrap to estimate standard errors: In this problem, we will employ bootstrapping to estimate the variability of estimation where standard software does not exist. Use the function, or risk losing points!

Using the data from package , estimate the standard error of the difference in Median (miles per gallon) between 4 cylinder and 8 cylinder cars.

# Process 1: Load required libraries 
library(ISLR)
library(boot)

## 
## Attaching package: 'boot'

## The following object is masked from 'package:lattice':
## 
##     melanoma

# Process 2: Prepare the data 
data(Auto)
# Make sure 'cylinders' is numeric and mpg is numeric (they are by default)

# Process 3: Define the statistic function 
# It returns the difference in median mpg: 4-cylinder - 8-cylinder
boot_median_diff <- function(data, indices) {
  d <- data[indices, ]
  median_4 <- median(d$mpg[d$cylinders == 4])
  median_8 <- median(d$mpg[d$cylinders == 8])
  return(median_4 - median_8)
}

# Process 4: Perform the bootstrap 
set.seed(1234)  # for reproducibility
boot_result <- boot(data = Auto, statistic = boot_median_diff, R = 1000)

# Process 5: Output the results 
print(boot_result)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Auto, statistic = boot_median_diff, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original   bias    std. error
## t1*     14.4 -0.30745    0.797002

cat("Estimated standard error of the median mpg difference:", sd(boot_result$t), "\n")

## Estimated standard error of the median mpg difference: 0.797002

Using the boot() function with 1,000 resamples, we estimated the standard error of the median MPG difference between 4-cylinder and 8-cylinder cars in the Auto dataset. The observed difference was 14.4 mpg, with a bootstrap-estimated standard error of 0.797 mpg, indicating low variability in this estimate due to sampling.

The bootstrap bias was −0.31, a small and negative value, suggesting that the center of the bootstrap distribution lies slightly below the original sample estimate. This minor bias is not a significant concern, and overall, the results confirm that the difference in fuel efficiency between 4-cylinder and 8-cylinder vehicles is both substantial and statistically stable.

Write your own function for conducting 10-fold cross validation to assess the test error in a linear regression model predicting estimated performance in the dataset from R package , using only minimum main memory and maximum main memory as predictors. DO NOT worry about including interaction effects, or higher order effects as in the first exam. Just a main effect for each variable will be sufficient. You need to show how to use the function (creates the folds for you) from the R package , as explained in class.

# Load required libraries
library(MASS)
library(caret)

# Load and prepare the data
data(cpus)
cpus_data <- cpus[, c("mmin", "mmax", "perf")]

# Define custom 10-fold CV function
my_cv_function <- function(data, k = 10) {
  set.seed(1234)  # for reproducibility
  folds <- createFolds(data$perf, k = k, list = TRUE, returnTrain = FALSE)
  
  mse_list <- numeric(k)
  
  for (i in 1:k) {
    test_indices <- folds[[i]]
    train_data <- data[-test_indices, ]
    test_data  <- data[test_indices, ]
    
    # Fit the linear model
    model <- lm(perf ~ mmin + mmax, data = train_data)
    
    # Predict on the test data
    predictions <- predict(model, newdata = test_data)
    
    # Compute Mean Squared Error
    mse_list[i] <- mean((test_data$perf - predictions)^2)
  }
  
  avg_mse <- mean(mse_list)
  return(avg_mse)
}

# Apply the function
test_error_estimate <- my_cv_function(cpus_data, k = 10)
cat("Estimated Test Error (10-fold CV MSE):", test_error_estimate, "\n")

## Estimated Test Error (10-fold CV MSE): 6315.865

data(cpus)
cpus_data <- cpus[, c("mmin", "mmax", "perf")]
# Compute the null model MSE (mean-only model)
mean_model_mse <- mean((cpus_data$perf - mean(cpus_data$perf))^2)

# Print the result
cat("Null Model MSE (predicting mean only):", mean_model_mse, "\n")

## Null Model MSE (predicting mean only): 25742.71

A custom function was developed using the createFolds() function from the caret package to perform 10-fold cross-validation on the cpus dataset (from the MASS package). A linear regression model was fitted to predict perf using mmin and mmax (minimum and maximum main memory) as predictors.

The estimated test error, calculated as the average mean squared error (MSE) across the 10 folds, was 6315.87. This result offers a reliable estimate of the model’s ability to generalize to new, unseen data.

Compared to the mean model MSE of 25742.71 (predicting the mean performance for all CPUs), this indicates that the linear regression model provides a much better fit than simply using the mean.

Exam 2 (Take Home Final Exam)

Jesmin Sultana

Please hand this result in on Canvas no later than 11:59pm on Thursday, May 8!

Questions