#install.packages("randomForest")
library(randomForest)

rm(list = ls())
set.seed(20221025)
# load the spam data
spam <- read.table('https://web.stanford.edu/~hastie/ElemStatLearn/datasets/spam.data', header = FALSE)
head(spam)

y <- spam$V58
X <- as.matrix(spam[, -58])
summary(X)

# split in training and test set
traintest <- read.table('https://web.stanford.edu/~hastie/ElemStatLearn/datasets/spam.traintest',
                        header = FALSE)
mean(traintest[, 1])
X.train <- X[traintest[, 1] == 0, ]
X.test <- X[traintest[, 1] == 1, ]
y.train <- as.factor(y[traintest[, 1] == 0])
y.test <- as.factor(y[traintest[, 1] == 1])


#example
rf1 <- randomForest(x = X.train, y = y.train, xtest = X.test, ytest = y.test, mtry = sqrt(ncol(X.train)))
rf1

oob_error <- NULL
test_error <- NULL
for(m in 1:57)
{
  rf <- randomForest(x = X.train, y = y.train, xtest = X.test, ytest = y.test, mtry = m, ntree = 200)
  oob_error[m] <- (rf$confusion[1, 2] + rf$confusion[2, 1]) / sum(rf$confusion)
  test_error[m] <- (rf$test$confusion[1, 2] + rf$test$confusion[2, 1]) / sum(rf$confusion)
  print(m)
}
  
plot(1:m, oob_error, type = 'l', ylim = c(0.02, 0.1), xlab="m", ylab="prediction error")
lines(1:m, test_error, col = 2)
legend("topright", legend=c("OOB error", "Test error"),
       col=c(1, 2), lty=c(1,1), cex=0.8)

# OOB errors are higher than test errors but have a similar trend

which(oob_error == min(oob_error))
which(test_error == min(test_error))

test_error[5]
test_error[8]
test_error[22]