Penalized regressions

We need the following packages:

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

data("Ozone")
summary(Ozone)

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

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.

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

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

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:

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:

data(spam)

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

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).