library(leaps)

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

data.train <- prostate[prostate$train, ]
data.test <- prostate[!prostate$train, ]
n <- nrow(data.train)
p <- ncol(prostate) - 2 # number of covariates, remove response and train/test

# AIC and BIC
mod.bestSubset <- regsubsets(lcavol ~ lweight + age + lbph + svi + lcp + gleason + pgg45 + lpsa,
                             data = data.train)
summary(mod.bestSubset) # the difference is in the penalty, so given the number of variables,
# the best model is the same no matter the criterion used
summary(mod.bestSubset)$cp # cp is proportional to AIC
model.aic <- coef(mod.bestSubset, id = which.min(summary(mod.bestSubset)$cp))
summary(mod.bestSubset)$bic # BIC
round(exp(-0.5 * summary(mod.bestSubset)$bic) / sum(exp(-0.5 * summary(mod.bestSubset)$bic)), 2)
model.bic <- coef(mod.bestSubset, id = which.min(summary(mod.bestSubset)$bic))

# 5-fold cross-validation
totalNumberModels <- sum(sapply(1:p, function(x) choose(p, x))) # no null model
K <- 5
set.seed(20200908)
index <- sample(rep(1:K, ceiling(n/K)))[1:n]
cv.err.k <- matrix(NA, nrow = K, ncol = totalNumberModels)
for (k in 1:K)
{
  mod.bestSubset <- regsubsets(lcavol ~ lweight + age + lbph + svi + lcp + gleason + pgg45 + lpsa,
                               data = data.train[index != k, ], nbest = 2^p, really.big = TRUE)
  for (j in 1:totalNumberModels)
  {
    y.hat <- as.matrix(cbind(1, data.train[index == k, names(coef(mod.bestSubset, id = j))[-1] ])) %*%
               coef(mod.bestSubset, id = j)
    cv.err.k[k, j] <- mean((data.train[index == k, 'lcavol'] - y.hat)^2)
  }
  print(k)  
}
cv.err.5 <- apply(cv.err.k, 2, mean)

# 10-fold cross-validation
K <- 10
set.seed(20200908)
index <- sample(rep(1:K, ceiling(n/K))[1:n])
cv.err.k <- matrix(NA, nrow = K, ncol = totalNumberModels)
for (k in 1:K)
{
  mod.bestSubset <- regsubsets(lcavol ~ lweight + age + lbph + svi + lcp + gleason + pgg45 + lpsa,
                               data = data.train[index != k, ], nbest = 2^p, really.big = TRUE)
  for (j in 1:totalNumberModels)
  {
    y.hat <- as.matrix(cbind(1, data.train[index == k, names(coef(mod.bestSubset, id = j))[-1] ])) %*%
      coef(mod.bestSubset, id = j)
    cv.err.k[k, j] <- mean((data.train[index == k, 'lcavol'] - y.hat)^2)
  }
  print(k)  
}
cv.err.10 <- apply(cv.err.k, 2, mean)

# 0.632 bootstrap
B <- 20
set.seed(20200908)
bootstrap.err.b <- matrix(NA, nrow = B, ncol = totalNumberModels)
for (b in 1:B)
{
  index <- sample(1:n, replace = TRUE)
  mod.bestSubset <- regsubsets(lcavol ~ lweight + age + lbph + svi + lcp + gleason + pgg45 + lpsa,
                               data = data.train[index, ], nbest = 2^p, really.big = TRUE)
  for (j in 1:totalNumberModels)
  {
    y.hat <- as.matrix(cbind(1, data.train[-index, names(coef(mod.bestSubset, id = j))[-1] ])) %*%
      coef(mod.bestSubset, id = j)
    bootstrap.err.b[b, j] <- mean((data.train[-index, 'lcavol'] - y.hat)^2)
  }
  print(b)  
}
bootstrap.err <- apply(bootstrap.err.b, 2, mean)

train.err <- NULL
mod.bestSubset <- regsubsets(lcavol ~ lweight + age + lbph + svi + lcp + gleason + pgg45 + lpsa,
                             data = data.train, nbest = 2^p, really.big = TRUE)
for (j in 1:totalNumberModels)
{
  y.hat <- as.matrix(cbind(1, data.train[ , names(coef(mod.bestSubset, id = j))[-1] ])) %*%
    coef(mod.bestSubset, id = j)
  train.err[j] <- mean((data.train[ , 'lcavol'] - y.hat)^2)
}

bootstrap.err.632 <- 0.368 * train.err + 0.632 * bootstrap.err

order(cv.err.5)[1:10]
order(cv.err.10)[1:10]
order(bootstrap.err)[1:10]
order(bootstrap.err.632)[1:10]

# order selection models
model.cv5 <- coef(mod.bestSubset, id = which.min(cv.err.5))
model.cv10 <- coef(mod.bestSubset, id = which.min(cv.err.10))
model.bootstrap <- coef(mod.bestSubset, id = which.min(bootstrap.err))
model.bootstrap632 <- coef(mod.bestSubset, id = which.min(bootstrap.err.632))

# best models
models <- list(AIC = model.aic, BIC = model.bic, CV5 = model.cv5, CV10 = model.cv10,
               BOOTSTRAP = model.bootstrap, BOOTSTRAP_0.632 = model.bootstrap632)

test.error.exercise <- list(min(cv.err.5), min(cv.err.10), min(bootstrap.err),
                            min(bootstrap.err.632))
names(test.error.exercise) <- names(models)[-c(1:2)]
test.error.exercise

test.error <- lapply(models,
                     function(x, data.test)
                       mean((as.matrix(cbind(1, data.test[ , names(x)[-1] ])) %*% x - data.test[, 'lcavol'])^2),
                     data.test = data.test)
                       
test.error