```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# *Tutorial Session : CV, subset selection and regularized regression*

# Chapter 5.3 Lab: Cross-validation

```{r}
library(ISLR) # We are going to use 'Auto' dateset from this package
```

## Sample a subset from a dataset

```{r}
set.seed(1) # For reproducibility.

# Sample a subset from the data Auto
num_samples <- nrow(Auto)
train <- sample(num_samples,196)  ##create indexes for my training data
```

## Train & Validation (using MSE)

#### Model 1

```{r}
lm.fit <- lm(mpg~horsepower,data=Auto,subset=train)
attach(Auto)

# (Cross Validation) Calculate the Mean Squared Error (MSE) on the validation set. 
mean((mpg-predict(lm.fit,Auto))[-train]^2)  #calculate MSE on the validation set
```

```{r}
# We can also create a function for repeated use
MSE <- function(y_pred, y){
  return(mean((y_pred-y)^2))
}
MSE(mpg[-train], predict(lm.fit,Auto[-train,]))
```

#### Model 2 & 3

```{r}
lm.fit2 <- lm(mpg~poly(horsepower,2),data=Auto,subset=train)
MSE(mpg[-train], predict(lm.fit2,Auto[-train,]))

lm.fit3 <- lm(mpg~poly(horsepower,3),data=Auto,subset=train)
MSE(mpg[-train], predict(lm.fit3,Auto[-train,]))
```

#### Repeat with another seed

```{r}
set.seed(2)
train <- sample(num_samples,196)

lm.fit <- lm(mpg~horsepower,data=Auto,subset=train)
mean((mpg-predict(lm.fit,Auto))[-train]^2)

lm.fit2 <- lm(mpg~poly(horsepower,2),data=Auto,subset=train)
mean((mpg-predict(lm.fit2,Auto))[-train]^2)

lm.fit3 <- lm(mpg~poly(horsepower,3),data=Auto,subset=train)
mean((mpg-predict(lm.fit3,Auto))[-train]^2)
```

## Leave-One-Out Cross-Validation

#### Intro of glm function

```{r}
# glm can also do regular linear regression
# we want to use it because it's handy for cross validation
glm.fit <- glm(mpg~horsepower,data=Auto)
coef(glm.fit)

lm.fit <- lm(mpg~horsepower,data=Auto)
coef(lm.fit) # result are the same
```

```{r}
library(boot)   # the package for cross-validation under glm models
```

#### Apply Leave-One-Out Cross-Validation

```{r}
glm.fit <- glm(mpg~horsepower,data=Auto)
cv.err <- cv.glm(Auto,glm.fit) #by default it does Leave-One-Out Cross-Validation

cv.err$delta  # the first value is the one we covered in the lecture. 
#the second one is a bias corrected version (in the book by Davison, A.C. and Hinkley, D.V.).
#in many cases, they are very close.
```

#### Using Leave-One-Out CV to compare models

```{r}
cv.error <- rep(0,5) #save CV error for each of the 5 models we are going to fit

for (i in 1:5){
 glm.fit <- glm(mpg~poly(horsepower,i),data=Auto) # polynomial of horsepower of i-th order
 cv.error[i] <- cv.glm(Auto,glm.fit)$delta[1]
 }
cv.error
plot(cv.error)
```

## k-Fold Cross-Validation

Leave-one-out is equivalent as k-fold CV where k = data size. But we can use other values! Picking k smaller can reduce computation cost.

```{r}
set.seed(17)
cv.error.10 <- rep(0,10)
for (i in 1:10){
 glm.fit <- glm(mpg~poly(horsepower,i),data=Auto)
 cv.error.10[i] <- cv.glm(Auto,glm.fit,K=10)$delta[1] #does 10-fold CV by setting K=10
 }
cv.error.10
plot(cv.error.10)
```

# Chapter 6.5 Lab 1: Subset Selection

```{r}
library(ISLR) # FYI: we get the "Hitters" dataset from ISLR (actually we already loaded this package earlier, no need to do again)
```

## Data Cleansing

#### Get a sense of the data

```{r}
dim(Hitters)
```

```{r}
summary(Hitters)
```

#### Get rid of NA values

```{r}
sum(is.na(Hitters$Salary))  #how many samples have missing salary
Hitters <- na.omit(Hitters)    #create a data set with completely observed samples
```

Check data again

```{r}
dim(Hitters)
sum(is.na(Hitters))
```

## Best subset selection (Global)

```{r}
library(leaps)   # package for performing subset selection
```

#### Perform best subset selection

```{r}
#Here by default it only produce results up to model with 8 parameters.
regfit.full <- regsubsets(Salary~.,Hitters) # best subset selection 
summary(regfit.full)
```

```{r}
#But we can tune 'nvmax' to change the number of output
regfit.full <- regsubsets(Salary~.,data=Hitters,nvmax=19) #all subset selection with max # of variables=19
reg.summary <- summary(regfit.full)
```

#### Measure how good the fit is

```{r}
names(reg.summary) # Look at what 'reg.summary' includes
```

```{r}
reg.summary$cp # cp: Mallows's Cp, which is a statistic used to measure the goodness of fit (the lower the better)
```

#### Plot the Result

```{r}
par(mfrow=c(2,2))

## First, Plot RSS
#Plot 1: RSS
plot(reg.summary$rss,xlab="Number of Variables",ylab="RSS",type="l") #Rss is always decreasing as we include more variables


## Then, we do another 3 plots that measure goodness of fit using different criteria
#Plot 2: Adjusted RSq
plot(reg.summary$adjr2,xlab="Number of Variables",ylab="Adjusted RSq",type="l")

which.max(reg.summary$adjr2) # Here, we find where the maximum 'Adjusted RSq' occurs
points(11,reg.summary$adjr2[11], col="red",cex=2,pch=20) # We label that (a indicator of the best model)

#Plot 3: Cp
plot(reg.summary$cp,xlab="Number of Variables",ylab="Cp",type='l')
which.min(reg.summary$cp) # Look for minimum, because, here, the lower the better
points(10,reg.summary$cp[10],col="red",cex=2,pch=20)

#Plot 4: BIC
plot(reg.summary$bic,xlab="Number of Variables",ylab="BIC",type='l')
which.min(reg.summary$bic)
points(6,reg.summary$bic[6],col="red",cex=2,pch=20)
```

## Forward and Backward Stepwise Selection

#### Forward Stepwise Selection

```{r}
regfit.fwd <- regsubsets(Salary~.,data=Hitters,nvmax=19,method="forward") 
summary(regfit.fwd)
```

#### Backward Stepwise Selection

```{r}
regfit.bwd <- regsubsets(Salary~.,data=Hitters,nvmax=19,method="backward")
summary(regfit.bwd)
```

## Comparison Between the Three Methods

Best subset selection

```{r}
coef(regfit.full,7)  # coefficient of best subset selection with 7 variables
```

Forward stepwise selection

```{r}
coef(regfit.fwd,7)
```

Backward stepwise selection

```{r}
coef(regfit.bwd,7)
```

## Choosing Among Models by CV

#### Create Train and test set

```{r}
set.seed(1)
train <- sample(c(TRUE,FALSE), nrow(Hitters),rep=TRUE) 
#you may also use train_index=sample(c(1:263),131) to create training data
test <- (!train)
```

#### Perform Best subset selection, save MSE for comparison

you may also use k-fold cross-validation to select best model

```{r}
regfit.best <- regsubsets(Salary~.,data=Hitters[train,],nvmax=19)
test.mat <- model.matrix(Salary~.,data=Hitters[test,])   # create an X matrix of test data

val.errors <- rep(NA,19)
for(i in 1:19){
  # First, we need to calculate the predicted values. Here we show how ot do it manually. But I guess we can also do it using 'predict.glm' function. 
   coefi <- coef(regfit.best,id=i) #pull coefficients from the fitted model (betas in regression)
   pred <- test.mat[,names(coefi)]%*%coefi #calculate prediction (i.e. under standard regression formula notations, this is X %*% beta)
   
   val.errors[i] <- mean((Hitters$Salary[test]-pred)^2) # calculate MSE
}
```

Display those values

```{r}
val.errors
```

Find where the minimum MSE occurs

```{r}
which.min(val.errors)
```

Look at that model

```{r}
coef(regfit.best, 7)
```

## Compute AIC, BIC, adjusted R\^2 in linear models

#### Using 'lm'

```{r}
lm1 <- lm(Fertility ~ . , data <- swiss)

AIC(lm1)
BIC(lm1)
summary(lm1)$adj.r.squared
```

#### Using 'glm'

In 'glm', there's no R\^2 calculated (because it's not defined in general for GLM)

```{r}
glm1 <- glm(Salary~.,data=Hitters)

AIC(glm1)
BIC(glm1)
```

# Chapter 6.6 Lab 2: perform ridge regression and the Lasso

We first create design matrix 'X' and the predicted variable 'y', for later use.

```{r}
x <- model.matrix(Salary~., Hitters)[,-1]  #remove the intercept
y <- Hitters$Salary
```

## Ridge Regression

```{r}
library(glmnet)  # the package for regularized linear regression
```

#### Perform ridge regression with differnt 'lambda' penalty

```{r}
grid <- 10^seq(10,-2,length=100)   #create a grid for \lambda
ridge.mod <- glmnet(x,y,alpha=0,lambda=grid) #alpha=0 is the ridge penalty, alpha=1 is the lasso penalty
```

Look at result shape (100 models, each with 20 parameters)

```{r}
dim(coef(ridge.mod))
```

#### Result comparison

Look at the 50th model

```{r}
ridge.mod$lambda[50] #the lambda value in the 50th model
```

```{r}
coef(ridge.mod)[,50]
```

```{r}
sqrt(sum(coef(ridge.mod)[-1,50]^2))  #calculuate L2 norm of beta
```

Look at the 60th model

```{r}
ridge.mod$lambda[60]
```

```{r}
coef(ridge.mod)[,60]
```

```{r}
sqrt(sum(coef(ridge.mod)[-1,60]^2)) #As we decrease lambda, L_2 norm of beta increase 
```

#### Summary / Simplified ridge regression process

In practice, we just need to do the following. First we fit model with different lambdas

```{r}
set.seed(1)
cv.out <- cv.glmnet(x,y,alpha=0) #10 fold cross validation, you can also add lambda=grid
bestlam <- cv.out$lambda.min #get best lambda
bestlam
```

Then we look at the best model

```{r}
ridge.mod <- glmnet(x,y,alpha=0,lambda=bestlam) #get the model under the best lambda value
coef(ridge.mod)[,1] #look at coefficients

pred.ridge <- predict(ridge.mod, s = bestlam, newx = x) # get predictions
```

#### The Lasso

Similar process: first we fit model with differnt lambdas

```{r}
set.seed(1)
cv.out <- cv.glmnet(x,y,alpha=1) #10 fold cross validation
bestlam <- cv.out$lambda.min
bestlam
```

Then we look at the best model

```{r}
lasso.mod <- glmnet(x,y,alpha=1,lambda=bestlam) #get the model under the best lambda value
lasso.coef <- coef(lasso.mod)[,1]
lasso.coef[lasso.coef!=0]
```

Look at predictions

```{r}
pred.lasso <- predict(lasso.mod, s = bestlam, newx = x)
pred.lasso
```

Look at the prediction for a specific person. For example: '-Wade Boggs'

```{r}
index_boggs <- which(names(pred.lasso[,1])=="-Wade Boggs")

c(pred.lasso[index_boggs], pred.ridge[index_boggs])
```
