You will need to load the core library for the course textbook. We will also be using various random number generation functions in today's seminar, so it would also be worthwhile setting the seed as shown in the code block below.

```{r}
library(ISLR)
set.seed(12345)
```


## Exercise 7.1 - Leave One Out Cross Validation

In this exercise, we will be using the `Weekly` dataset, which can be loaded with the following code:

```{r}
data("Weekly")
```

This data includes weekly percentage returns for the S&P 500 stock index between 1990 and 2010. To find out more about the variables included in this data, you can look at the help file, which can be accessed using `?Weekly`. We will be interested in predicting the `Direction` variable, which indicates whether the S&P 500 went up or down in a given week.

(a) Start by using the `summary()` command to explore the data. In how many weeks did the market move up?

(b) Fit a logistic regression model that predicts `Direction` using `Lag1` and `Lag2`. Recall that to estimate a logistic regression, you will need to use the `glm()` function, with `family = binomial` as an argument to that function. Summarise the results. Are either of the two lag variables significant? 

(c) Use the regression that you just estimated to calculate the predicted probability of the outcome for every observation in the data. To do so, you can simply use the `predict()` function on the estimated model object, with `type = "response"`. This will produce a vector of probabilities which will be larger when an observation is more likely to have moved up rather than down, given its covariate values. For how many observations is the probability of the S&P going up greater than .5?

(d) Create a new vector which is equal to `"Up"` if the model predicts upward movement in the stock market and `"Down"` if the model predicts downward movement in the stock market (i.e. use the vector of probabilities you created in part (c) and code which observations are greater than .5). Use this new variable to construct a confusion matrix with the original outcome variable, `Direction`. How many false positives are there? How many false negatives are there? What is the training set error rate?

    N.b. Remember that, in the context of classification problems, such as the one considered here, the error rate is $\frac{1}{n}\sum_{i=1}^N Err_i$. That is, for each observation, we predict the outcome and compare it to the observed outcome for that observation. If the predicted outcome is different from the actual outcome, then $Err_i = 1$, otherwise it is equal to 0. The error-rate is then just the average of those errors across all observations.

(e) As we discussed in the lecture, the training set error rate may underestimate the error rate for the test set. We will now use leave-one-out cross validation to estimate the test-set error rate. In the textbook sections for this topic (particularly sections 5.3.2 and 5.3.3 in James et al.), you will see that the `cv.glm()` function can be used in order to compute the LOOCV test error estimate. Alternatively, one can compute those quantities using just the `glm()` and `predict()` functions, and a `for` loop. Although it can be more laborious to do this using a loop, it also helps to clarify what is going on in this process!

    Start by re-fitting the logistic regression model that you fit in answer to part (b) (predicting `Direction` using `Lag1` and `Lag2`), but here use *all observations except the first observation to fit the model*. Remember, you can subset a `data.frame` by using the square brackets. So, here, you will want something like `my_data[-1,]`. How do the coefficients you estimate in this model differ from those in the model in part (b)?

(f) Use the model from (e) to predict the outcome of only the first observation. You can do this by predicting that the first observation will go up if the fitted value for that observation is greater than 0.5. Remember that to calculate predicted probabilities (rather than predicted log-odds) you will need to add `type = "response"` as an argument to the `predict()` function. Compare your prediction to the true value of `Direction` for the first observation. Was this observation correctly classified?

(g) To estimate the LOOCV error, we need to now replicate the process you followed above for *every* observation in the dataset. Fortunately, rather than writing hundreds of lines of code, we can accomplish this using a `for` loop. Write a `for` loop that loops over from `i=1` to `i=n`, where n is the number of observations in the data set (i.e. you could use `nrow(Weekly)`). At each iteration, your loop should performs each of the following steps:

    i. Fit a logistic regression model using all but the i-th observation to predict `Direction` using `Lag1` and `Lag2` (i.e. as in your answer to (e) above).
    ii. Calculate the probability of the market moving up for the i-th observation (this requires using the same code as you used for your answer to part (f) above.). 
    iii. Use the probability for the i-th observation in order to predict whether or not the market moves up (again, as you did in part (f)). 
    iv. Determine whether or not an error was made in predicting the direction for the i-th observation. If an error was made, then indicate this as a 1, and otherwise indicate it as a 0. Make sure that you store these results!

    The code provided below can be used to help with this question:

```{r, echo = TRUE, eval = FALSE}

# Set up a vector in which you can store the errors for each observation
errors <- rep(0, nrow(Weekly))

# Loop over each observation in the data
for (i in 1:nrow(Weekly)) {
  
  # Fit the logistic regression model
   ## YOUR CODE GOES HERE
   
   # Calculate the predicted probability, and check whether this is greater than .5
     ## YOUR CODE GOES HERE
   
   # Calculate the true outcome for the i-th observation
    ## YOUR CODE GOES HERE
   
   # Assign the errors vector to be equal to 1 if a prediction error was made (the vector is all 0s otherwise)
    ## YOUR CODE GOES HERE
}

# Count up the number of errors made

```


(h) Take the average of the `n` numbers obtained from the loop above in order to calculate the LOOCV estimate for the test error. Comment on the results. How does this compare to the estimate of the training set error calculated in part (d)?

## Exercise 7.2 - K-Fold Cross Validation

Running the leave-one-out procedure here does not take too long because the `Weekly` dataset has only `r nrow(Weekly)` observations, and this is a simple logistic regression model. However, often LOOCV is too computationally burdensome to be useful. In those cases, K-fold cross-validation is normally used instead. Though it is possible to code cross-validation manually, using only the functions we used above, in this case it will suffice to use some functions designed for the purpose.

(a) First, load the boot library (you may need to `install.packages("boot")` first):

```{r}

library(boot)

```

(b) Now fit the same logistic regression model as in part (b) in the question above - i.e. on the whole data

(c) Now use the `cv.glm()` function from the boot package to estimate the cross-validation error. Look at the help file for this function and you will see that you need to provide 4 arguments to make the function work: 

1. the data for estimating the model
2. the estimated model object from part (b) above
3. a "cost" function, which we have copied for you below (this just tells the function that you want to do cross-validation for a logit model)
4. a value for `K`, the number of folds to use. 

    Use this function to implement 10-fold cross validation and report the cross-validation error (this information is stored in the first element of `delta` component of the estimated model object). How does this estimate compare to the LOOCV reported in the question above?

```{r}
cost <- function(r, pi) mean(abs(r-pi)> 0.5)

```

(d) Experiment with some different values for `K` in the `cv.glm` function. How does this affect the cross-validation error?

## Exercise 7.3 - Lasso and ridge regression

In this exercise, we will predict the number of applications (`Apps`) received by US colleges using the other variables in the `College` data set. Again, if you would like to read more about the variables in the data, just type `?College`.

(b) Fit a linear model with `Apps` as the outcome, and all other variables in the data as predictors. Report the estimated coefficients from the model.

(c) Using the same outcome and the same predictors, fit a ridge regression model on the `College` data.  To get you started, the code below creates a matrix of predictors using the `model.matrix()` function, which makes working with `glmnet()` a little more straightforward. 

```{r, eval = FALSE}

library(glmnet) # Load the relevent package for estimating ridge and lasso regressions

train.mat <-  model.matrix(Apps ~ . , data = College)[,-1]

```

To estimate the ridge regression, you will need to use the `glmnet` function from the `library(glmnet)` package. This function takes four main arguments:

1. `x` - a matrix of predictor variables
2. `y` - a vector of outcomes
3. `alpha` - a parameter determining whether you want to run a ridge regression (`alpha = 0`) or a lasso regression (`alpha = 1`)
4. `lambda` - the penalty parameter applied to the ridge regression

    For the purposes of this exercise, use `train.mat` for `x`, and the outcome from the `College` data for `y`. Set `alpha = 1` and `lambda = 1000` (a very large penalty parameter).  Report the coefficients (`coef()`) estimated by this model. How do they compare to the coefficients that you estimated using OLS in part (b) above?

(d) Experiment by lowering the value of the `lambda` parameter. How does this affect the estimated coefficients? What happens when `lambda = 0`?

(e) To estimate a lasso regression, you can use exactly the same code as for the ridge regression, only changing the `alpha` parameter to be equal to 1 rather than 0. Try this now, and compare the OLS coefficients to the lasso coefficients when `lambda = 10`.

## Exercise 7.3 - Out-of-sample error for lasso and ridge regression (hard question)

*For this question, you might find it helpful to refer to pages 274-275 of the James et. al. textbook.*

(a) In this question we will use the validation set approach in order to estimate the test error rate using lasso and ridge regression. First, split the data set into a training set and a validation set. You may find the `sample()` function helpful for selecting observations at random to include in the two sets.

(b) Use least squares to estimate the model on the training set, and then report the test error obtained for the validation set. To calculate the test error, you will need to use the formula $\frac{1}{n}\sum_{i=1}^n (y_i - \hat{y_i})^2$, where $y_i$ is the observed value of the outcome (`Apps`) for observation $i$ and $\hat{y_i}$ is the fitted value for observation $i$ from the regression model.


(c) Using the same outcome and the same predictors, fit a ridge regression model on the training data that you just created. For the purposes of this question, set `lambda = 0.1`. Remember, you will have to convert the training data predictors using the `model.matrix()` function, as in part (c) of the previous question, in order for the `glmnet()` function to work. Use the `predict()` function to calculate fitted values for the test set by setting the `newx` argument to be equal to your test data. Report the test set MSE from this regression. How does it compare to the test set error you calculated in part (b)?


(d) Repeat the process you followed in part (c), but this time using a lasso regression (remember, a lasso is estimated by setting `alpha = 1` in `glmnet()`). Again, set `lambda = 0.1`. Report the test error obtained. How does this compare to the test set errors for the OLS and ridge models? 


(e) Our results so far have relied on the arbitrary choice of $\lambda = 0.1$. Let's try to do a bit better than that. Use 10-fold cross-validation to select the $\lambda$ parameter for both the ridge and lasso models. To do so, you will need to use the `cv.glmnet()` function, which takes the following main arguments:

|Argument|Value|
|:--------|:----------|
|`x`|A matrix of predictor variables|
|`y`|A response variable|
|`family`|A character string representing the family of the outcome variable. Use `"gaussian"` here for a numeric outcome.|
|`nfolds`|Number of folds to use in cross validation|
|`alpha`|Parameter expressing the type of penalised model you would like to fit. `alpha = 1` is the lasso penalty, and `alpha = 0` the ridge penalty.|

Implement the cross-validation using this function.

(f) Calculate fitted values for the test set you created in part a above using the lasso and ridge regressions that you fit in answer to parts c and d above, but this time using the best performing $\lambda$ value from the cross-validation you just conducted in part e. `cv.glmnet()` stores the optimal choice for $\lambda$ in the output of the fitted object. Extract it using `estimated_cv_object$lambda.min` and use as input to the `s` argument of the `predict()` function. (You can also see pages 276-277 of the James et. al. book if you are struggling here.) Report the test error obtained for both models.



(g) Does this cross-validation exercise help us to predict the outcome relative to the implementations above? 

