---
title: "TMA4268 Statistical Learning V2018"
subtitle: "Module 6: Solution sketches"
author: "Thiago G. Martins, Department of Mathematical Sciences, NTNU"
date: "Spring 2018"
output: #3rd letter intentation hierarchy
#  beamer_presentation:
###    incremental: true # or >* for one at a time
#  slidy_presentation:
#    font_adjustment: +1  
   prettydoc::html_pretty:
    theme: architect
    highlight: github
#  pdf_document:
  #  toc: true
  #  toc_depth: 2

---

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

# Recommended exercise 1

Solution can be found [here](https://en.wikipedia.org/wiki/Bias%E2%80%93variance_tradeoff#Derivation).

# Recommended exercise 2

For the least square estimator, the solution can be found in the first session [here](https://en.wikipedia.org/wiki/Proofs_involving_ordinary_least_squares).

For the maximum likelihood estimator, the solution can be found [here](https://en.wikipedia.org/wiki/Proofs_involving_ordinary_least_squares#Maximum_likelihood_approach).

# Recommended exercise 3

```{r, message=FALSE, warning=FALSE}
library(ISLR) # Package with data for an Introduction to Statistical 
              # Learning with Applications in R
```

```{r, message=FALSE, warning=FALSE}
# Load Credit dataset
data(Credit)

# Check column names
names(Credit)

# Check dataset shape
dim(Credit)
```

```{r, message=FALSE, warning=FALSE}
head(Credit)
```

```{r, message=FALSE, warning=FALSE}
# Select variable to plot
pairwise_scatter_data <- Credit[,c("Balance", "Age", "Cards", "Education", 
                                   "Income", "Limit", "Rating")]
```

```{r, message=FALSE, warning=FALSE}
# Simplest possible pairwise scatter plot
pairs(pairwise_scatter_data)
```

```{r, message=FALSE, warning=FALSE, cache=TRUE}
# More interesting but slower pairwise plot from package GGally
library(GGally)
ggpairs(data=pairwise_scatter_data)
```

Check [here](https://tgmstat.wordpress.com/2013/11/13/plot-matrix-with-the-r-package-ggally/) for quick get started to ggpairs

# Recommended exercise 4

From a linear algebra perspective, when $n < p$ we have more unknowns than equations, leading to either no solution or multiple solutions to the regression problem. 

In other words, the design matrix $X$ do not have full rank, leading $X^TX$ not being invertible and therefore we are unable to directly compute 

$$ \hat{\beta} = (X^TX)^{-1} X^T y $$

# Recommended exercise 5

```{r, message=FALSE, warning=FALSE}
# Exclude 'ID' column
credit_data <- subset(Credit, select=-c(ID))

# Counting the dummy variables as well
credit_data_number_predictors <- 11

# Take a look at the data
head(credit_data)

# Summary statistics
summary(credit_data)

# Create train and test set indexes
set.seed(1)
train_perc <- 0.75
credit_data_train_index <- sample(1:nrow(credit_data), nrow(credit_data)*train_perc)
credit_data_test_index <- (-credit_data_train_index)

# Create train and test set
credit_data_training <- credit_data[credit_data_train_index, ]
credit_data_testing <- credit_data[credit_data_test_index, ]
```

```{r, message=FALSE, warning=FALSE}
library(leaps)

# Perform best subset selection using all the predictors and the training data
best_subset_method=regsubsets(Balance~.,credit_data_training,nvmax=credit_data_number_predictors)

# Save summary obj
best_subset_method_summary=summary(best_subset_method)
```

```{r, message=FALSE, warning=FALSE}
# Plot RSS, Adjusted R^2, C_p and BIC

par(mfrow=c(2,2))
plot(best_subset_method_summary$rss,xlab="Number of Variables",ylab="RSS",type="l")
plot(best_subset_method_summary$adjr2,xlab="Number of Variables",ylab="Adjusted RSq",type="l")
bsm_best_adjr2 = which.max(best_subset_method_summary$adjr2)

points(bsm_best_adjr2,best_subset_method_summary$adjr2[bsm_best_adjr2], col="red",cex=2,pch=20)
plot(best_subset_method_summary$cp,xlab="Number of Variables",ylab="Cp",type='l')
bsm_best_cp=which.min(best_subset_method_summary$cp)

points(bsm_best_cp,best_subset_method_summary$cp[bsm_best_cp],col="red",cex=2,pch=20)
bsm_best_bic=which.min(best_subset_method_summary$bic)

plot(best_subset_method_summary$bic,xlab="Number of Variables",ylab="BIC",type='l')
points(bsm_best_bic,best_subset_method_summary$bic[bsm_best_bic],col="red",cex=2,pch=20)
```

```{r, message=FALSE, warning=FALSE}
# Create a prediction function to make predictions
# for regsubsets with id predictors included
predict.regsubsets=function(object,newdata,id,...){
  form=as.formula(object$call[[2]])
  mat=model.matrix(form,newdata)
  coefi=coef(object,id=id)
  xvars=names(coefi)
  mat[,xvars]%*%coefi
}

# Create indexes to divide the data between folds
k=10
set.seed(1)
folds=sample(1:k,nrow(credit_data_training),replace=TRUE)
cv.errors=matrix(NA,k,credit_data_number_predictors, dimnames=list(NULL, paste(1:credit_data_number_predictors)))

# Perform CV
for(j in 1:k){
  best_subset_method=regsubsets(Balance~.,data=credit_data_training[folds!=j,],nvmax=credit_data_number_predictors)
  for(i in 1:credit_data_number_predictors){
    pred=predict(best_subset_method,credit_data_training[folds==j,],id=i)
    cv.errors[j,i]=mean( (credit_data_training$Balance[folds==j]-pred)^2)
    }
}

# Compute mean cv errors for each model size
mean.cv.errors=apply(cv.errors,2,mean)
mean.cv.errors

# Plot the mean cv errors
par(mfrow=c(1,1))
plot(mean.cv.errors,type='b')
```

```{r, message=FALSE, warning=FALSE}
# Fit the selected model using the whole training data
# and compute test error

# models selected
number_predictors_selected <- 4

# Create info for lm call
variables <- names(coef(best_subset_method,id=number_predictors_selected))
variables <- variables[!variables %in% "(Intercept)"]
bsm_formula <- as.formula(best_subset_method$call[[2]])
bsm_design_matrix <- model.matrix(bsm_formula,credit_data_training)[, variables]
bsm_data_train <- data.frame(Balance = credit_data_training$Balance, bsm_design_matrix)

# Fit a standard linear model using only the selected 
# predictors on the training data
model_best_subset_method <- lm(formula = bsm_formula, bsm_data_train)
summary(model_best_subset_method)

# Make predictions on the test set
bsm_design_matrix_test <- model.matrix(bsm_formula,credit_data_testing)[, variables]
bsm_predictions <- predict(object = model_best_subset_method, newdata = as.data.frame(bsm_design_matrix_test))

# Compute test squared errors
bsm_squared_errors <- (credit_data_testing$Balance-bsm_predictions)^2
squared_errors <- data.frame(bsm_squared_errors=bsm_squared_errors)


# test MSE
mean(bsm_squared_errors)
```

# Recommended exercise 6

Similar analysis as previous exercise, simply replace Best Subset Selection (`best_subset_method=regsubsets(Balance~.,credit_data,nvmax=credit_data_number_predictors)`) by Forward Stepwise Selection (`regfit.fwd=regsubsets(Balance~.,credit_data,nvmax=credit_data_number_predictors,method="forward")`), Backward Stepwise Selection (`regfit.fwd=regsubsets(Balance~.,credit_data,nvmax=credit_data_number_predictors,method="backward")`) and Hybrid Stepwise Selection 
(`regfit.fwd=regsubsets(Balance~.,credit_data,nvmax=credit_data_number_predictors,method="seqrep")`)

# Recommended exercise 7

```{r, message=FALSE, warning=FALSE}
library(glmnet) # Package Lasso and Elastic-Net Regularized 
                # Generalized Linear Models
```

```{r, message=FALSE, warning=FALSE}
x_train <- model.matrix(Balance~.,credit_data_training)[,-1]
y_train <- credit_data_training$Balance

x_test <- model.matrix(Balance~.,credit_data_testing)[,-1]
y_test <- credit_data_testing$Balance

ridge_mod <- glmnet(x_train,y_train,alpha=0)

set.seed(1)
cv.out=cv.glmnet(x_train, y_train,alpha=0)
plot(cv.out)

best_lambda_ridge <- cv.out$lambda.min
best_lambda_ridge

ridge_predictions = predict(ridge_mod,s=best_lambda_ridge,newx=x_test)
ridge_square_errors <- as.numeric((ridge_predictions-y_test)^2)
squared_errors <- data.frame(ridge_square_errors = ridge_square_errors, squared_errors)

```

# Recommended exercise 8

```{r, message=FALSE, warning=FALSE}
lasso_mod <- glmnet(x_train,y_train,alpha=1)

set.seed(1)
cv.out=cv.glmnet(x_train, y_train,alpha=1)
plot(cv.out)

best_lambda_lasso <- cv.out$lambda.min
best_lambda_lasso

lasso_predictions = predict(lasso_mod,s=best_lambda_lasso,newx=x_test)
lasso_square_errors <- as.numeric((lasso_predictions-y_test)^2)
squared_errors <- data.frame(lasso_square_errors = lasso_square_errors, squared_errors)
```

# Recommended exercise 9

```{r, message=FALSE, warning=FALSE}
standard_lm <- lm(formula = Balance~., data = credit_data_training)
standard_lm_predictions <- predict(object = standard_lm, newdata = credit_data_testing)
standard_lm_square_errors <- as.numeric((standard_lm_predictions-credit_data_testing$Balance)^2)
squared_errors <- data.frame(standard_lm_square_errors = standard_lm_square_errors, squared_errors)


# Comparison between methods as measures by computing MSE on a independent test set
library(reshape2)
library(ggplot2)
ggplot(melt(squared_errors)) + geom_boxplot(aes(variable, value))
```

# Recommended exercise 10

```{r, message=FALSE, warning=FALSE}
x <- model.matrix(Balance~.,credit_data)[,-1]

credit_pca <- prcomp(x, center = TRUE, scale. = TRUE) 

print(credit_pca)

plot(credit_pca, type = "l")

summary(credit_pca)
```

The first PC explain along $25\%$ of the variability in the data. Then the second PC explain an extra $15\%$ of the variability in the data. From the third PC until $8$th PC the extra variability explained per PC varies between $7.5\%$ to $10\%$, dropping to $3.6\%$ on the $9$th PCA. So I would likely use $8$ PCs for the Credit dataset.

# Recommended exercise 11

```{r, message=FALSE, warning=FALSE}
library(pls)

set.seed(1)

pcr_model <- pcr(Balance~., data=credit_data_training,scale=TRUE, validation="CV")
validationplot(pcr_model,val.type="MSEP")
```

```{r, message=FALSE, warning=FALSE}
pcr_predictions = predict(pcr_model,credit_data_testing,ncomp=10)
pcr_square_errors <- as.numeric((pcr_predictions-credit_data_testing$Balance)^2)
squared_errors <- data.frame(pcr_square_errors = pcr_square_errors, squared_errors)
mean(pcr_square_errors)
```

```{r, message=FALSE, warning=FALSE}
ggplot(melt(squared_errors)) + geom_boxplot(aes(variable, value))
```

# Recommended exercise 12

```{r, message=FALSE, warning=FALSE}
library(pls)

set.seed(1)

plsr_model <- plsr(Balance~., data=credit_data_training,scale=TRUE, validation="CV")
validationplot(plsr_model,val.type="MSEP")
```

```{r, message=FALSE, warning=FALSE}
plsr_predictions = predict(plsr_model,credit_data_testing,ncomp=3)
plsr_square_errors <- as.numeric((plsr_predictions-credit_data_testing$Balance)^2)
squared_errors <- data.frame(plsr_square_errors = plsr_square_errors, squared_errors)
mean(plsr_square_errors)
```

```{r, message=FALSE, warning=FALSE}
ggplot(melt(squared_errors)) + geom_boxplot(aes(variable, value))

colMeans(squared_errors) 
```
