library(tree)
library(randomForest)

# a)
carseats <- as.data.frame(unclass(read.csv("Downloads/Carseats.xls")), stringsAsFactors = T)
n <- nrow(carseats)

test_prop <- 0.2
n_test_obs <- round(test_prop * n)

test_indices <- sample(1:n, replace = F, size = n_test_obs)

carseats_train <- carseats[-test_indices, ]
carseats_test <- carseats[test_indices, ]

# b)
par(mfrow=c(1,1))
tree_fit <- tree(Sales~., data = carseats_train)

mean((predict(tree_fit, carseats_test[,-1]) - carseats_test[,1])^2)

plot(tree_fit)
text(tree_fit, main='Initial tree')

# c)
tree_cv <- cv.tree(tree_fit)

tree_cv
plot(tree_cv$size, tree_cv$dev, type='l', xlab="size", ylab="test error")

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

best_size
pruned_tree <- prune.tree(tree_fit, best=best_size)

plot(pruned_tree)
text(pruned_tree, main='Pruned tree')

mean((predict(pruned_tree, carseats_test[,-1]) - carseats_test[,1])^2)

# e) 
forest_model <- randomForest(Sales~., data=carseats_train, importance=T)
imp_table <- importance(forest_model)
imp_table

varImpPlot(forest_model)

mean((predict(forest_model, carseats_test[,-1]) - carseats_test[,1])^2)

# Explanation: Increase in MSE and node purity if we randomly shuffle a predictor in 
# OOB data

m_errors <- numeric(10)

# We try different number of variables considered at each split 
for(i in 1:10){
  new_forest_model <- randomForest(Sales~., data=carseats_train, importance=T, mtry=i)
  m_errors[i] <- mean((predict(new_forest_model, carseats_test[,-1]) - carseats_test[,1])^2)
}

plot(1:10, m_errors, type='l')