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.

set.seed(1)

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.

library(ISLR2)
## Warning: package 'ISLR2' was built under R version 4.6.1
library(tree)
## Warning: package 'tree' was built under R version 4.6.1
library(randomForest)
## Warning: package 'randomForest' was built under R version 4.6.1
## randomForest 4.7-1.2
## Type rfNews() to see new features/changes/bug fixes.
library(dbarts)
## Warning: package 'dbarts' was built under R version 4.6.1
data(Carseats)
  1. Split the data set into a training set and a test set.
set.seed(1)
n <- nrow(Carseats)
train_index <- sample(1:n, n / 2)
train <- Carseats[train_index,]
test <- Carseats[-train_index,]
  1. Fit a regression tree to the training set. Plot the tree, and interpret the results. What test MSE do you obtain? - 4.922039 TEST MSE result.
tree_fit <- tree(Sales ~ ., data = train)
summary(tree_fit)
## 
## Regression tree:
## tree(formula = Sales ~ ., data = 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_fit)
text(tree_fit, pretty = 0)

tree_pred <- predict(tree_fit, test)
tree_test_error <- mean((tree_pred - test$Sales)^2)
tree_test_error
## [1] 4.922039
  1. Use cross-validation in order to determine the optimal level of tree complexity. Does pruning the tree improve the test MSE? - Pruning still resulted the same test of ‘4.922039’ and did not improve. Cross-validation was used with an expectation of out of sample prediction error.
set.seed(1)
cv_tree <- cv.tree(tree_fit)
plot(cv_tree$size, cv_tree$dev, type = "b")

best_size <- cv_tree$size[which.min(cv_tree$dev)]
best_size
## [1] 18
pruned_tree <- prune.tree(tree_fit, best = best_size)
plot(pruned_tree)
text(pruned_tree, pretty = 0)

pruned_pred <- predict(pruned_tree, test)
pruned_test_error <- mean((pruned_pred - test$Sales)^2)
pruned_test_error
## [1] 4.922039
  1. 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. -Test MSE of ‘2.605253’.
set.seed(1)
bag_fit <- randomForest(Sales ~ ., data = train, mtry = ncol(train) - 1, importance = TRUE)
bag_pred <- predict(bag_fit, test)
bag_test_error <- mean((bag_pred - test$Sales)^2)
bag_test_error
## [1] 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
varImpPlot(bag_fit)

  1. 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. The error rate obtained was ‘2.960559’. The most important variables were Price and ShelveLoc due to having the highest %IncMSE. As m increased from a small number toward p, the test MSE decreased.
set.seed(1)
rf_fit <- randomForest(Sales ~ ., data = train, importance = TRUE)
rf_pred <- predict(rf_fit, test)
rf_test_error <- mean((rf_pred - test$Sales)^2)
rf_test_error
## [1] 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
varImpPlot(rf_fit)

mtry_vals <- 1:(ncol(train) - 1)
rf_errors <- sapply(mtry_vals, function(m) {
  set.seed(1)
  fit <- randomForest(Sales ~ ., data = train, mtry = m)
  pred <- predict(fit, test)
  mean((pred - test$Sales)^2)
})

names(rf_errors) <- mtry_vals
rf_errors
##        1        2        3        4        5        6        7        8 
## 4.799716 3.521097 3.014166 2.835478 2.710950 2.679895 2.642999 2.634485 
##        9       10 
## 2.598166 2.592505
plot(mtry_vals, rf_errors, type = "b", xlab = "mtry", ylab = "test MSE")

  1. Now analyze the data using BART, and report your results.
x_train <- train[, -which(names(train) == "Sales")]
y_train <- train$Sales
x_test <- test[, -which(names(test) == "Sales")]
y_test <- test$Sales

bart_fit <- 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.431428
## 
## Tree sizes, last iteration:
## [1] 2 3 2 4 2 2 2 2 3 3 2 2 1 4 4 2 4 2 
## 2 2 2 2 3 2 2 3 2 3 2 3 3 4 3 2 3 2 2 3 
## 3 2 1 2 3 3 2 2 2 2 2 2 4 2 2 1 2 2 3 1 
## 4 2 2 2 2 2 2 2 2 3 2 2 2 2 3 2 2 3 2 2 
## 3 2 2 3 2 2 2 2 1 1 3 2 2 2 3 2 2 3 2 2 
## 2 3 2 2 2 3 4 2 2 2 2 2 2 2 3 2 2 3 2 3 
## 2 2 2 2 3 2 2 2 2 4 3 3 2 3 2 2 3 3 3 3 
## 3 3 3 2 2 1 4 3 2 3 4 2 3 2 2 2 4 3 2 2 
## 3 2 3 2 2 3 2 2 3 2 3 2 2 4 2 4 3 2 1 2 
## 2 3 2 2 4 3 2 2 3 2 3 4 2 2 2 2 3 2 2 2 
## 2 2 
## 
## Variable Usage, last iteration (var:count):
## (1: 35) (2: 17) (3: 23) (4: 26) (5: 29) 
## (6: 24) (7: 29) (8: 24) (9: 14) (10: 21) 
## (11: 19) (12: 17) 
## DONE BART
bart_pred <- bart_fit$yhat.test.mean
bart_test_error <- mean((bart_pred - y_test)^2)
bart_test_error
## [1] 1.435738

Question 9

We will now consider the Boston housing data set, from the ISLR2 library.

library(ISLR2)
library(boot)

data(Boston)
  1. Based on this data set, provide an estimate for the population mean of medv. Call this estimate ˆμ. - estimate = 22.53281
mu_hat <- mean(Boston$medv)
mu_hat
## [1] 22.53281
  1. Provide an estimate of the standard error of ˆμ. Interpret this result. Hint: We can compute the standard error of the sample mean by dividing the sample standard deviation by the square root of the number of observations. - SE = .4088611, tells us the sample mean would expect to vary from sample to sample.
se_formula <- sd(Boston$medv) / sqrt(nrow(Boston))
se_formula
## [1] 0.4088611
  1. Now estimate the standard error of ˆμ using the bootstrap. How does this compare to your answer from (b)? -.4106622, similar to B.
set.seed(1)
boot_mean_fn <- function(data, index) {
  mean(data[index])
}
boot_mean <- boot(Boston$medv, boot_mean_fn, R = 1000)
boot_mean
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Boston$medv, statistic = boot_mean_fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original      bias    std. error
## t1* 22.53281 0.007650791   0.4106622
  1. Based on your bootstrap estimate from (c), provide a 95 % confidence interval for the mean of medv. Compare it to the results obtained using t.test(Boston$medv). Hint: You can approximate a 95 % confidence interval using the formula [ˆμ − 2SE(ˆμ), ˆμ + 2SE(ˆμ)].
se_boot_mean <- sd(boot_mean$t)
ci_boot <- c(mu_hat - 2 * se_boot_mean, mu_hat + 2 * se_boot_mean)
ci_boot
## [1] 21.71148 23.35413
t.test(Boston$medv)
## 
##  One Sample t-test
## 
## data:  Boston$medv
## t = 55.111, df = 505, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  21.72953 23.33608
## sample estimates:
## mean of x 
##  22.53281
  1. Based on this data set, provide an estimate, ˆμmed, for the median value of medv in the population.
mu_med_hat <- median(Boston$medv)
mu_med_hat
## [1] 21.2
  1. We now would like to estimate the standard error of ˆμmed. Unfortunately, there is no simple formula for computing the standard error of the median. Instead, estimate the standard error of the median using the bootstrap. Comment on your findings.
set.seed(1)
boot_median_fn <- function(data, index) {
  median(data[index])
}

boot_median <- boot(Boston$medv, boot_median_fn, R = 1000)
boot_median
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Boston$medv, statistic = boot_median_fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*     21.2 0.02295   0.3778075
  1. Based on this data set, provide an estimate for the tenth percentile of medv in Boston census tracts. Call this quantity ˆμ0.1. (You can use the quantile() function.)
mu_0.1_hat <- quantile(Boston$medv, 0.1)
mu_0.1_hat
##   10% 
## 12.75
  1. Use the bootstrap to estimate the standard error of ˆμ0.1. Comment on your findings. - SE of the 10th percentile is .4767526 and is the largest of the three standard errors (mean,median,10th percentile).
set.seed(1)
boot_percentile_fn <- function(data, index) {
  quantile(data[index], 0.1)
}

boot_percentile <- boot(Boston$medv, boot_percentile_fn, R = 1000)
boot_percentile
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Boston$medv, statistic = boot_percentile_fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    12.75  0.0339   0.4767526