# Proposed solutions to exercise set 5 in STK2100/Spring 2024
#
# - These solutions are taken from the week 7 problems Spring 2020
# - See also  Ch 6.5.1 in James et al. (2013) for the code examples in problem 6.11

# -------------------------------
# 6.2 from ITSL
# -------------------------------

# -------------------------------
# a) Lasso

# iii. Less flexible and better predictions because of less variance, more bias

# -------------------------------
# b) Ridge Regression

# iii. Same as lasso


# -------------------------------
# c) Non-linear methods

# ii. More flexible, less bias, more variance


# -------------------------------
# 6.3 from ITSL
# -------------------------------

# -------------------------------
# a)

# (iv) Steadily decreases: As we increase s
# from 0, all β 's increase from 0 to their least square estimate values. Training error for 0 β
# s is the maximum and it steadily decreases to the Ordinary Least Square RSS

# -------------------------------
# b)

# (ii) Decrease initially, and then eventually start increasing in a U shape: When s=0
# , all β s are 0, the model is extremely simple and has a high test RSS. As we increase s, beta s
# assume non-zero values and model starts fitting well on test data and so test RSS decreases.
# Eventually, as beta s approach their full blown OLS values, they start overfitting to the training
# data, increasing test RSS.

# -------------------------------
# c)

# (iii) Steadily increase: When s=0
# , the model effectively predicts a constant and has almost no variance. As we increase s,
# the models includes more β s and their values start increasing. At this point, the values of β
# s become highly dependent on training data, thus increasing the variance.

# -------------------------------
# d)

# (iv) Steadily decrease: When s=0, the model effectively predicts a constant and hence the
# prediction is far from actual value. Thus bias is high. As s increases, more β s become
# non-zero and thus the model continues to fit training data better. And thus, bias decreases.

# -------------------------------
# e)
# (v) Remains constant: By definition, irreducible error is model independent and hence
# irrespective of the choice of s, remains constant.


# -------------------------------
# 6.4 from ITSL
# -------------------------------

# -------------------------------
# a)

# (iii) Steadily increase: As we increase λ from 0, all β 's decrease from their least square
# estimate values to 0. Training error for full-blown-OLS β s is the minimum and it steadily
# increases as β s are reduced to 0.

# -------------------------------
# b)

# (ii) Decrease initially, and then eventually start increasing in a U shape: When λ=0, all β s
# have their least square estimate values. In this case, the model tries to fit hard to training
# data and hence test RSS is high. As we increase λ, beta s start reducing to zero and some of the
# overfitting is reduced. Thus, test RSS initially decreases. Eventually, as beta s approach 0, the
# model becomes too simple and test RSS increases.

# -------------------------------
# c)

# (iv) Steadily decreases: When λ=0, the β s have their least square estimate values. The actual
# estimates heavily depend on the training data and hence variance is high. As we increase λ, β s start
# decreasing and model becomes simpler. In the limiting case of λ approaching infinity, all beta s reduce
# to zero and model predicts a constant and has no variance.

# -------------------------------
# d)

# (iii) Steadily increases: When λ=0, β s have their least-square estimate values and hence have the
# least bias. As λ increases, β s start reducing towards zero, the model fits less accurately to training
# data and hence bias increases. In the limiting case of λ approaching infinity, the model predicts a
# constant and hence bias is maximum.

# -------------------------------
# e)

# (v) Remains constant: By definition, irreducible error is model independent and hence
# irrespective of the choice of λ, remains constant.

# -------------------------------
# 6.11 from ITSL
# -------------------------------

# -------------------------------
# a)
rm(list=ls())

# Best Subset Selection

set.seed(1)
library(MASS)
library(leaps)
library(glmnet)
library(pls)


predict.regsubsets = function(object, newdata, id, ...) {
  form = as.formula(object$call[[2]])
  mat = model.matrix(form, newdata)
  coefi = coef(object, id = id)
  mat[, names(coefi)] %*% coefi
}

k = 10
p = ncol(Boston) - 1
folds = sample(rep(1:k, length = nrow(Boston)))
cv.errors = matrix(NA, k, p)
for (i in 1:k) {
  # returns best subset of each size 1, ..., p
  best.fit = leaps::regsubsets(crim ~ ., data = Boston[folds != i, ], nvmax = p)
  for (j in 1:p) {
    pred = predict(best.fit, Boston[folds == i, ], id = j)
    cv.errors[i, j] = mean((Boston$crim[folds == i] - pred)^2)
  }
}
rmse.cv = sqrt(apply(cv.errors, 2, mean))
plot(rmse.cv, pch = 19, type = "b")

which.min(rmse.cv)
rmse.cv[which.min(rmse.cv)] #6.543

# Lasso
x = model.matrix(crim ~ . - 1, data = Boston)
y = Boston$crim
cv.lasso = cv.glmnet(x, y, type.measure = "mse")
plot(cv.lasso)
coef(cv.lasso)

sqrt(cv.lasso$cvm[cv.lasso$lambda == cv.lasso$lambda.1se]) #7.921

# Ridge Regression
x = model.matrix(crim ~ . - 1, data = Boston)
y = Boston$crim
cv.ridge = cv.glmnet(x, y, type.measure = "mse", alpha = 0)
plot(cv.ridge)

coef(cv.ridge)

sqrt(cv.ridge$cvm[cv.ridge$lambda == cv.ridge$lambda.1se]) #7.669

# PCR
pcr.fit = pcr(crim ~ ., data = Boston, scale = TRUE, validation = "CV")
summary(pcr.fit)

# 13 component pcr fit has lowest CV/adjCV RMSEP. Lowest CV: 6.488

# -------------------------------
# b)

# See above answers for cross-validate mean squared errors of selected models.

# -------------------------------
# c)

# I would choose the 9 parameter best subset model because it had the best
# cross-validated RMSE, after PCR, but it was simpler model than the 13 component PCR model.