We need the following packages:

```{r message=FALSE, warning=FALSE}
library(tidyverse)
library(leaps)
library(glmnet)
library(mlbench)
library(doMC)
library(caret)
library(class)
library(reshape2)
library(plotROC)
library(kernlab)
library(bestglm)
```


# Best subset selection 

## Exercise 1

We consider the **ozone** data set from the **mlbench** package

```{r}
data("Ozone")
summary(Ozone)
```

We first delete individuals with missing data and the three first variables which contains the date of the observations:

```{r}
donnees <- Ozone %>% na.omit()
donnees <- donnees[,-c(1:3)]
summary(donnees)
```

We then split the data into:

  * a training set of size 150;
  * a test set of size 53.
  
```{r}
set.seed(1234)
perm <- sample(nrow(donnees))
ind.train <- perm[1:150]
ind.test <- perm[-c(1:150)]
train <- donnees[ind.train,]
test <- donnees[ind.test,]
```

The problem is to explain the `Daily maximum one-hour-average ozone reading` (V4) by the other variables. We first consider the linear model

```{r}
linear.model <- lm(V4~.,data=train)
summary(linear.model)
```

It seems that some variables are not necessary in the model.

1. Explain the output of the following command

```{r}
mod.sel <- regsubsets(V4~.,data=train)
summary(mod.sel)
```


2. Select the model which optimize Mallows's $C_p$ and $BIC$ criteria.


3. We consider the quadratic risk for the three models:
$$E[(Y-\widehat m(X))^2].$$
This risk is estimated with the test set according to
$$\frac{1}{n_{test}}\sum_{i\in test}(Y_i-\widehat m(X_i))^2.$$
Compute the estimated risks for the three linear models.


## Exercise 2 (variable selection)

We  consider the model
$$Y=\beta_0+\beta_1X_1+\dots+\beta_pX_p+\varepsilon$$
where $X_1,\dots,X_p$ are independent variables with distribution $\mathcal N(0,1)$ and $\varepsilon$ has distribution $\mathcal N(0,0.5^2)$ ($X=(X_1,\dots,X_p)$ and $\varepsilon$ are independent).

We assume that $p=105$ and that

* $\beta_0=0$, $\beta_1=\dots=\beta_5=1$
* $\beta_6=\dots=\beta_{105}=0$.
So here, only variables $X_1,\dots,X_5$ are relevant to explain $Y$.

1. Generate $n=1000$ observations $(x_1,y_1),\dots,(x_n,y_n)$ according to this model. Put these observations in a **data.frame df** with dimension $1000\times 106$.

2. Fit a linear model (**lm** function) with *df* and print the estimator of $\beta_0,\dots,\beta_{105}$.


3. We propose to make a variable selection procedure with a backward selection approach using **BIC** criterion. You can use *nvmax=30* in **regsubsets**.


4. Generate a test dataset (**df.test**) $(x_{n+1},y_{n+1}),\dots,(x_{n+m},y_{n+m})$ with $m=500$ according to the same model. Compute the estimated mean square error
$$\frac{1}{m}\sum_{i\in test}(\widehat y_i-y_i)^2$$
for the two models (full and backward).


# Penalized regression

## Exercise 3 (Lasso and ridge for ozone data)

We consider the same problem as in exercise 1. We propose to fit ridge ans lasso estimates and to compare these models with those of exercise 1.

**glmnet** function requires the $X$ matrix and the $Y$ vector as inputs (we can not use a formula). To compute the $X$ matrix, we use the function **model.matrix**:
```{r}
train.X <- model.matrix(V4~.,data=train)
test.X <- model.matrix(V4~.,data=test)
```

1. Draw the coefficient paths for ridge and lasso.

2. Select the shrinkage parameter for lasso regression with **cv.glmnet**.

3. Do the same for ridge regression. Explain the problem and solve it.

4. Estimate the quadratic error for the selected ridge and lasso models. Compare this error with the error of the tree other models.

5. Conclude.



## Exercise 4

We consider the spam dataset from the **kernlab** package:

```{r}
data(spam)
```

We split the data into a training dataset of size 3000 and a test dataset of size 1601.

```{r}
set.seed(1234)
perm <- sample(nrow(spam))
train <- spam[perm[1:3000],]
test <- spam[perm[3001:4601],]
train.X <- model.matrix(type~.,data=train)
test.X <- model.matrix(type~.,data=test)
```

1. Fit a logistic model on the train dataset with all the variables.

2. We propose to make a variable selection procedure with a backward selection approach using **BIC** criterion (this approach will be described in details in other lectures). You just have to use the **step** function with the **direction="backward"** and **k=log(nrow(train))** options. We call it **mod.back** (it may takes few minutes).

3. Fit a logistic lasso model on the training data (select the shrinkage parameter with **cv.glmnet**).

4. Fit a logistic ridge model on the training data (select the shrinkage parameter with **cv.glmnet**).

5. Fit a nearest neigbor rule on the training data (use caret to select the number of nearest neigbor by 10 folds cross validation).

6. Make a comparison of the 5 methods with the error probability (estimated on the test dataset).

7. Make a comparison of the 5 methods with ROC curve and AUC (estimated on the test dataset).

