Make Generation Reproducible

To make our random number generation reproducible, we set the seed of the generation algorithm.  This is done with the set.seed function, which takes a user-selected integer argument.  Let’s pick 1234 for this

argument and set the seed:

set .seed(1234)

Cross Validation

In this lab, we consider a class of methods that estimate the test error rate by holding out a subset of the training observations from the fitting process, and then applying the statistical learning method to those held out observations. We will go over three approaches of such to evaluate performance of statistical models/statistical learning methods.

1. Validation set approach: randomly divide the available set of observations into two parts, a training set and a validation set or hold-out set.  The model is fit on the training set, and the fitted model is used to predict the responses for the observations in the validation set.  The resulting validation set error rate-typically assessed using MSE in the case of a quantitative response-provides an estimate of the test error rate.

2. Leave-One-Out Cross Validation:

3. K-folds Cross Validation

We will use Auto Data set from the library ISLR. Load the data as follows

#  install .packages("ISLR")

#library(ISLR)

data (Auto)

head (Auto)

##

##  1

##  2

##  3

##  4

##  5

##  6

##

##  1

##  2

##  3

##  4

##  5

##  6

We will compare the following regression models:

M1 : mpg = β0 + β1 * horsepower + ϵ

M2 : mpg = β0 + β1 * horsepower + β2 * horsepower2 + β3 * horsepower3 + ϵ

Validation set approach

We will first divide the data set into two separate datasets one for training and one for validation as follows:

set .seed(1234)

train=sample(392 ,196)

Auto_train=Auto[train,]

Auto_val=Auto[-train,]

Note that function set.seed is used to obtain reproducible results while using a function returning non- deterministic output. Next, we fit the required models using the training set and calculate an estimate for test MSE using the validation set as follows

#M1

lm.fit=lm(mpg~horsepower,data=Auto_train)

lm.pred=predict(lm.fit ,Auto_val)

(val_err1=mean((Auto_val$mpg-lm .pred)ˆ2))

##  [1]  23 .93601

#M2

lm.fit2=lm(mpg~poly(horsepower,3),data=Auto_train)

lm.pred2=predict(lm.fit2 ,Auto_val)

(val_err2=mean((Auto_val$mpg-lm .pred2)ˆ2))

##  [1]  18 .20033

K-fold Cross Validation

In order to perform K-fold cross validation, we need to assign each data point to a fold at random. We will use the function createFolds in the package caret for this task as follows:

#  install .packages("caret")

#library(caret)

set .seed(1234)

nfolds=10  #Perform  10  fold  cross  validation

fold=createFolds (1:392 ,k=nfolds,list=FALSE)

fold  #check  if  assignment  looks  random

##      [1]    7    9    5    3    6    6    4    9    4  10    3    8    6  10    3  10    8    5    3    3    1    3    4    2    7 ##    [26]    6    7    4    1    9    4    2  10    4    8    5    8    1    9    4  10    3    7    2    8    2  10    7    8    8 ##    [51]    5    3    7    9    1    4    6    6    1  10    1    2    9    9    6    7    2    5    2  10    7    3    9    3    1 ##    [76]    5    5    5    5    1    1    6    1    7    6  10    8    4    2    2    9    8    6    8    2    4    7    9    8    7 ##  [101]    5    2    6  10    2    2    3    5    5    4    4    5  10    4    5    6    1    5    8    2    3    5    7    6    4 ##  [126]    6    9    8    6    3    7    3    8    2    9    2    2    6    4  10    1    3    7    4    2  10    2    9    5  10 ##  [151]    1    8    7    7    9  10    9    1    5  10    6    9    7    1    3    4    4    3    7    4    9    6    3    5    6 ##  [176]    9    1  10    1    7  10    8    2    1    8  10    8    9    8    3    1    9    6    7    8    1    6  10  10    9 ##  [201]    8    4    6    2    4    2    1  10    4    6    7    5    6    7    7    8    3    1    2    2    2    3    2    4    9 ##  [226]    5    1    6    8    2    7    8    6    7    9    3    4    5    1    2    5    7    6    4    5  10  10    5    5    8 ##  [251]    9  10    3    7  10    7    4    3    7    2  10    1    9    8    8    4    9  10    1    7    1    6    8    4    3 ##  [276]    3    6    1    5    8    3    1    1    8    9    6    5    5    9    3    2  10    4    9    9    9    6    4    3  10 ##  [301]    4    1    1    3    5    8    6    6    2    8    7    1    1    5  10    7    2    8    1    8  10    7    5    6    2 ##  [326]    7    7  10  10    4    7    9    9  10    5    9  10    2    7    6    3    6    3    4    8    8    5    5    1    9 ##  [351]    1    6    7    8    4    5    3    4    4    1    8    1    9    6  10    4    5    6    5    7    9    2  10    2    2 ##  [376]    4    3    5    1    4    2    2    3    3    8    6    9    2    3    8    9    3

#Number  of  data  points  in  each  fold  is  roughly  equal

table(fold)

##  fold

##    1    2    3    4    5    6    7    8    9  10

##  40  40  38  39  39  40  39  40  39  38

Now, we run a loop to train the model using all the data points, but those in one fold.  We calculate the estimate for test MSE using the rest. In the end, we take an average of test MSE estimates for each fold.

#initialize  error  to  0

kCV_err1=0

kCV_err2=0

for(i  in  1:nfolds)

{

m1=lm(mpg~horsepower,data=Auto[fold !=i,])

m2=lm(mpg~poly(horsepower,3),data=Auto[fold !=i,])

pred1=predict(m1,Auto[fold==i,])

pred2=predict(m2,Auto[fold==i,])

kCV_err1=kCV_err1  + mean ((Auto$mpg[fold==i]-pred1)ˆ2)

kCV_err2=kCV_err2  + mean ((Auto$mpg[fold==i]-pred2)ˆ2)

}

(kCV_err1=kCV_err1/nfolds)

##  [1]  24 .18192

(kCV_err2=kCV_err2/nfolds)

##  [1]  19 .54784

Leave-One-Out Cross Validation

For Leave-One-Out Cross Validation, we do not require to make folds as earlier, we just run a loop over all the data points leaving one point out of the training set at each step and calculating the estimate using the remaining point. In the end, we take an average of all the estimates.

#Initialize  error  to  0

LOO_err1=0

LOO_err2=0

for  (i  in  1:392)

{

m1=lm(mpg~horsepower,data=Auto[-i,])

m2=lm(mpg~poly(horsepower,3),data=Auto[-i,])

pred1=predict(m1,Auto[i,])

pred2=predict(m2,Auto[i,])

LOO_err1=LOO_err1  +  (Auto$mpg[i]-pred1) ˆ2

LOO_err2=LOO_err2  +  (Auto$mpg[i]-pred2) ˆ2

}

(LOO_err1=LOO_err1/392)

##                1

##  24 .23151

(LOO_err2=LOO_err2/392)

##                1

##  19 .33498

Using cv.glm for K-fold Cross Validation

For regression models, the function cv.glm in the package boot allows to calculate the K-fold cross validation error without manually writing a loop. It can be used as follows:

#library(boot)

glm.fit1=glm(mpg~horsepower,  data=Auto)

cv.err1=cv.glm(Auto,glm.fit1 ,K=10)

cv .err1$delta[ 1]

##  [1]  24 .31284