ISLR - Chapter 5

Resampling Methods - Labs

Cross Validation

---

We will utilize the Auto dataset of ISLR package to compare different linear methods of varying flexibility to model mpg as function of horespower using CV approaches

Leave One Out Cross Validation (LOOCV)

Loading the required packages

require(boot)

Loading required package: boot

Plot of mpg vs horsepower

Using glm() to model mpg as a function of horsepower

cv.glm(Auto, glm.fit)$delta

[1] 24.23151 24.23114

cv.glm is bit slow even for a linear model as it fits any method n times. However, just for linear model we have a direct simple formula for comuting the CV estimate of test error without having to fit the model n times.

This formula is given as -

LOOCV(n) = 1/n sum{[(y_i - y_hat_i)/(1-h_i)]^2}

We will write a user-defined function to compute test error estimate using this formula

(test_err <- loocv(glm.fit))

[1] 24.23151

In the above code blocks, we estimated the test error by fitting a linear model on mpg using linear horespower term.

Now we will use LOOCV to compare linear models of varying degrees.

{
  glm.fit <- glm(mpg ~ poly(horsepower, i))
  cv.error[i] <- loocv(glm.fit)
}

Error in eval(predvars, data, env) : object 'mpg' not found

Plot of LOOCV test error estimate against linear model of varying flexibilities

The above plot shows that linear model of degree 1 does a poor job. The test error estimates drops down for degree 2 model. Inclusion of higher degree terms do not make much difference.

Hence we can conclude that a linear model of degree 2 is suited for modeling the relationship between mpg and horsepower

---

K- fold Cross Validation (K = 10)

The concept remains the same as LOOCV, except that instead of leaving one observation out of the model fitting process, we leave out approx (n/k) observations. The model fitting process is repeated k times

We will reuse the above code block with a small change. We will set the argument “k” to 10 in the cv.glm()

{
  glm.fit <- glm(mpg ~ poly(horsepower, i), data = Auto)
  cv.error10[i] <- cv.glm(Auto, glm.fit, k = 10)$delta[1]
}

Error in cv.glm(Auto, glm.fit, k = 10) : unused argument (k = 10)

Comparing the plot of test error estimate from LOOCV and K-fold CV approach

As we observe there is not much difference in the outputs of the 2 methods.

Generally speaking, we use 10-fold (or 5-fold) CV approach as it is computationally less expensive than LOOCV and it tends to produce more stable measure.

---

---

Bootstrap

---

We will illustrate the use of Bootstrap to estimate the standard error of the coefficients estimated by the linear regression method

Here we use the bootstrap approach in order to assess the variability of the estimates for β0 and β1, the intercept and slope terms for the linear regression model that uses horsepower to predict mpg in the Auto data set.

  return(coef(lm(mpg ~ horsepower, data = data, subset = index)))

Error in as.data.frame.default(data, optional = TRUE) : 
  cannot coerce class ‘"function"’ to a data.frame

Calling boot() function to repeat the model fit process on 1000 bootstrapped samples

boot(data = Auto, boot.fn, R = 1000)

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Auto, statistic = boot.fn, R = 1000)


Bootstrap Statistics :
      original        bias    std. error
t1* 39.9358610  0.0040975524 0.848703651
t2* -0.1578447 -0.0001175595 0.007325433

The output indicates that the bootstrap estimate of the standard error of Beta0 = 0.848 and that of Beta1 = 0.0074.

We can compare these values with the one produced by lm function by executing summary function on its output

summary(lm.out)

Call:
lm(formula = mpg ~ horsepower, data = Auto)

Residuals:
     Min       1Q   Median       3Q      Max 
-13.5710  -3.2592  -0.3435   2.7630  16.9240 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 39.935861   0.717499   55.66   <2e-16 ***
horsepower  -0.157845   0.006446  -24.49   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.906 on 390 degrees of freedom
Multiple R-squared:  0.6059,    Adjusted R-squared:  0.6049 
F-statistic: 599.7 on 1 and 390 DF,  p-value: < 2.2e-16

We see that the 2 values are different.

Does this mean the our boostrap estimate of standard error of cofficients is incorrect? No, this is not the case. Infact the bootstrap estimate is more accurate as it does not depend on any assumptions.

To determine the SE of cofficients we need to estimate the unknown error(noise) variance. And estimation of this noise variance depends on the linear model being correct.

As we have seen that the relationship between mpg and horsepower is not linear. So the residuals from linear model will be inflated and so will be the estimated error variance.

To validate this, we will now fit a quardatic model and compare the estimates from 2 methods

boot(Auto, boot.fn2, R = 10000)

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Auto, statistic = boot.fn2, R = 10000)


Bootstrap Statistics :
      original        bias    std. error
t1*   23.44592 -0.0009480357   0.3932161
t2* -120.13774  0.1498780916   4.7499967
t3*   44.08953 -0.3022614152   4.7778677

summary(lm.out2)$coef

                       Estimate Std. Error   t value      Pr(>|t|)
(Intercept)            23.44592  0.2209163 106.13030 2.752212e-289
poly(horsepower, 2)1 -120.13774  4.3739206 -27.46683  4.169400e-93
poly(horsepower, 2)2   44.08953  4.3739206  10.08009  2.196340e-21