---
title: "Lab 5 Solutions"
output: pdf_document
---


### Validation set approach
#### 1)	Randomly pick half of the data as the training data. Remember to set a seed to make your result repeatable.

```{r}
library(ISLR)

set.seed(100)
train <- sample(392,196)
print(train)
```

#### 2)	Build a linear regression model based the training data.

```{r}
lm.fit.train <- lm(mpg~horsepower,data=Auto,subset=train)
print(lm.fit.train)
```

#### 3)	Estimate the test MSE based on the other half (as test data)

```{r}
error <- rep(0,3)
error[1] <- mean((Auto$mpg-predict(lm.fit.train, Auto))[-train]^2)
print(error[1])
```

#### 4)	Now try to build polynomial regression of degree 2 and 3 using `lm(y~poly(x,i))`, where `y` is the response variable, `x` is the predictor variable and `i` is the highest degree of `x`. Compute the test MSE for the two models.

```{r}
lm.fit.train2 <- lm(mpg~poly(horsepower,2),data=Auto,subset=train)
error[2] <- mean((Auto$mpg-predict(lm.fit.train2, Auto))[-train]^2)

lm.fit.train3 <- lm(mpg~poly(horsepower,3),data=Auto,subset=train)
error[3] <- mean((Auto$mpg-predict(lm.fit.train3, Auto))[-train]^2)


sapply(error,print)
```

or 

```{r}
for(i in 1:3){
  print(paste("CV error for degree", i, "is", error[i]))
}
```

#### 5)	What conclusion could we draw from the above comparison of degree 1 (linear) and degree 2 (quadratic) and degree 3 (cubic) regression models?

A model that predicts `mpg` using a quadratic function of `horsepower` performs better than a model that involves only a linear function of `horsepower`, and there is little evidence in favor of a model that uses a cubic function of `horsepower`.


#### 6)	Choose 10 different seeds. For each seed, calculate the test MSE for models of degree from 1 to 10. You may use a nested for-loop to do that. Plot the variability on the results. Can you obtain a similar plot as below.

The first snippet of code records all the result in a matrix. Each row of the matrix represents a seed. Each column of the matrix represents a degree of the polynomial. 

```{r}
set.seed(1)
train <- sample(392,196)
errors <- rep(0,10)
for(i in 1:10){
     lm.fit.train <- lm(mpg~poly(horsepower,i),data=Auto,subset=train)
      errors[i] <- mean((Auto$mpg-predict(lm.fit.train, Auto))[-train]^2)
}


plot(errors,
     col=1, pch=".", type="l",
     xlab="Degrees of Polynomial",ylab="Mean Squared Error",
     main="10 times random split",
     ylim = c(14,27), xlim= c(0,12))

errorMatrix <- matrix(nrow=10,ncol=10)
errorMatrix[1,] <- errors
for(i in 2 : 10){
  set.seed(i)
  train <- sample(392,196)
  for(j in 1:10){
    lm.fit.train <- lm(mpg~poly(horsepower,j),data=Auto,subset=train)
    errorMatrix[i,j] <- mean((Auto$mpg-predict(lm.fit.train,Auto))[-train]^2)
  }
  lines(errorMatrix[i,],col=i)
}

legend("bottomright",
       c("1","2","3","4","5","6","7","8","9","10"),
       lty=rep(1,10),col=1:10)

```

We can simplify the code above as follows:
```{r}

plot(1,
     col=1,pch=".",type="l",
     xlab="Degrees of Polynomial",ylab="Mean Squared Error",
     main="10 times random split",
     ylim = c(14,27), xlim= c(0,12))
#plot(1) creates a plot with one dot at (1,1). 
#Normally we use this to plot an "empty" plot. 
#We later use lines() to add lines to this plot.

errorMatrix <- matrix(nrow=10,ncol=10)
for(i in 1 : 10){ #i is the seed. 
  #You may change this to for example 
  #for(i in seq(100,by = 20, length.out=10))
  set.seed(i)
  train <- sample(392,196)
  for(j in 1:10){#j is the degree of the polynomial
    lm.fit.train <- lm(mpg~poly(horsepower,j),data=Auto,subset=train)
    errorMatrix[i,j] <- mean((Auto$mpg-predict(lm.fit.train,Auto))[-train]^2)
  }
  lines(errorMatrix[i,],col=i)
}

legend("bottomright",
       c("1","2","3","4","5","6","7","8","9","10"),
       lty=rep(1,10),col=1:10)

```

We can also use a single vector instead of a matrix to record and plot the test MSE. But the drawback is that the previous data will be overwritten. 
```{r}

plot(1,
     col=1,pch=".",type="l",
     xlab="Degrees of Polynomial",ylab="Mean Squared Error",
     main="10 times random split",
     ylim = c(14,27), xlim= c(0,12))
#plot(1) creates a plot with one dot at (1,1). 
#Normally we use this to plot an "empty" plot. 
#We later use lines() to add lines to this plot.

errors <- NULL #create an empty vector called errors
for(i in 1 : 10){ #i is the seed. 
  #You may change this to for example 
  #for(i in seq(100,by = 20, length.out=10))
  set.seed(i)
  train=sample(392,196)
  for(j in 1:10){#j is the degree of the polynomial
    lm.fit.train <- lm(mpg~poly(horsepower,j),data=Auto,subset=train)
    errors[j] <- mean((Auto$mpg-predict(lm.fit.train,Auto))[-train]^2)
  }
  lines(errors,col=i)
}

legend("bottomright",
       c("1","2","3","4","5","6","7","8","9","10"),
       lty=rep(1,10),col=1:10)

```

### (II) LOOCV

#### 7)	Experiment on the LOOCV for increasingly complex polynomial fits. More specifically, write a for-loop to increase the degree i, as in `lm(y~poly(x,i))`, from 1 to 10 and record the LOOCV estimate for the test error for each degree.

```{r}
library(boot)
cv.error <- rep(0,10)
for(i in 1:10){
     glm.fit <- glm(mpg~poly(horsepower,i),data=Auto)
     cv.error[i] <- cv.glm(Auto,glm.fit)$delta[1]
}
#the above for loop takes a while  
print(cv.error)
```

#### 8)	Plot the result from 7) where x-axis is the degree `i` and y-axis is the LOOCV estimate for the test error. 

```{r}
plot(cv.error,
     col="blue",pch=16,
     xlab="Degrees of Polynomial",ylab="Mean Squared Error",
     main="LOOCV",ylim=c(15,28))
lines(cv.error,col="blue")
```

or use the following to make the line broken:

```{r}
plot(cv.error,
     col="blue",pch=16,
     xlab="Degrees of Polynomial",ylab="Mean Squared Error",
     main="LOOCV",
     ylim=c(15,28),type="b") #type="b" means broken lines
```

### (III) K-fold CV

#### 9)	Set a seed. Write a for-loop to increase the degree `i`, as in `lm(y~poly(x,i))`, from 1 to 10 and record the 10-fold CV estimate for the test error for each degree.

#### 10)	Plot the result from 9) where x-axis is the degree `i` and y-axis is the 10-fold CV estimate for the test error.

#### 11)	Set 9 different seeds and repeat 9) and 10). Plot all the results into one plot.
```{r}
#plot(1) creates an empty plot as a base. 
plot(1,
     type="l",
     xlab="Degrees of Polynomial",ylab="Mean Squared Error",
     main="K-fold CV",
     ylim = c(14,27), xlim = c(0,12))

cv.errors <- rep(0,10)

for(i in 1:10){ #10 lines
        set.seed(i)
        for(j in 1:10){#10 degrees
          glm.fit <- glm(mpg~poly(horsepower,j),data=Auto)
          cv.errors[j] <- cv.glm(Auto,glm.fit,K=10)$delta[1]
    }
    lines(cv.errors,col=i)
}

legend(title = "degrees",
       "bottomright",
       c("1","2","3","4","5","6","7","8","9","10"),
       lty=rep(1,10),col=1:10)

```
