---
title: "Cross-validation"
author: "Torben"
date: "August 20, 2018"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, comment = NA)
```

```{r libs, warning=FALSE, message=FALSE}
library(tidyverse)
theme_set(theme_bw())
```


# Generalisability of models

Let $y = f(x) + \varepsilon$ we our model, where $\varepsilon$ is a stochastic error with zero mean and variance $\sigma^2$.
Note, we don't assume anything about the distribution (e.g. not normality), only that the error is independent of $f(x)$.

Hence, we know that $$E[y] = E[f(x) + \varepsilon] = E[f(x)] + E[\varepsilon] = f(x)$$.

We do not necessarily know the *shape* of f(x) - but we wants to learn it. We may not even know which part of 
$x = (x_1,x_2,\dots,x_p)$ that affects $y$. We have just collected data on the phenomenon $y$ and hope the 
systematic components depends on the collected explanatory data, $x$.

## Example: Linear models

We know that linear models are models where $y = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p + \varepsilon$. Or in 
other words, $$y = f(x) + \varepsilon = X\beta + \varepsilon,$$ where $\varepsilon \sim N(0,\sigma^2I)$

### Example: `mtcars` and polynomial regression

```{r}
(mpg_hp_plot <- ggplot(mtcars, aes(y = mpg, x = hp)) + geom_point())
(mpg_hp_plot <- mpg_hp_plot + geom_smooth(method = "lm", formula = y ~ x))
(mpg_hp_plot <- mpg_hp_plot + geom_smooth(method = "lm", formula = y ~ x + I(x^2), colour = "red"))
(mpg_hp_plot <- mpg_hp_plot + geom_smooth(method = "lm", formula = y ~ x + I(x^2) + I(x^3), colour = "green"))
```

How many powers of `hp` is needed for a good fit?

What is a good fit? How is it measured?

The *root mean squared error* (RMSE) is one measure: 
$$RMSE = \sqrt{\frac{1}{n}\sum_{i=1}^n(y_i - \hat{y}_i)^2} = \sqrt{\frac{1}{n}\sum_{i=1}^nr_i^2} $$,
where $\hat{y}_i = \hat{f}(x_i)$ and $r_i$ is the $i$'th residual.

```{r}
rmse <- function(lm_obj){
  sqrt(mean(residuals(lm_obj)^2))
}
```

How does this perform?

```{r}
rmse(lm(mpg ~ hp, data = mtcars))
rmse(lm(mpg ~ poly(hp, 2), data = mtcars))
rmse(lm(mpg ~ poly(hp, 3), data = mtcars))
rmse(lm(mpg ~ poly(hp, 4), data = mtcars))
# and it goes on
```

How well does my model generalise? The above property only holds when *test* and *traning data* are the same.

# Test and traning data

We can split the data into *test* and *traning data* - however, a single split only gives us a single estimate of the
generalisation error. The solution is to do this $K$ times, by using $K$-fold cross validation

![10-fold cross validation](day-2-crossval.png)

## $K$-fold cross-validation

We divide the data into $K$ equal sized chunks - the first $K-1$ serves as training and the last as test. We permute the
role as test data $K$ times (each time the training is the non-test data).

```{r}
K <- 10
mtcars$cv_fold <- sample(K, size = nrow(mtcars), replace = TRUE)
head(mtcars)

cv_rmse <- function(power, K = 10){
  mtcars$pred <- NA
  for(k in 1:K){
    train_lm <- lm(mpg ~ poly(hp, degree = power), data = subset(mtcars, cv_fold != k))
    mtcars$pred[mtcars$cv_fold == k] <- predict(train_lm, newdata = subset(mtcars, cv_fold == k))
  }
  sqrt(mean((mtcars$mpg - mtcars$pred)^2))
}
```

So what happens when we run this?

```{r}
cv_rmse(power = 1)
cv_rmse(power = 2)
cv_rmse(power = 3)
cv_rmse(power = 4)
cv_rmse(power = 5)
```

We see that the cross-validated RMSE starts to increase after `power = 2` or `power = 3`, suggesting that 
we start to overfit to the traning data.

## CV is very powerful

No matter the type of model $f$ we are fitting to data, we can always do cross-validation. For some 
model types it is not possible to do hypothesis tests with the usual distributional assumptions as for `lm`,
hence we can utilise CV for the same type of questions..