```{r echo = FALSE, results = 'hide', message = FALSE}
library(tidyverse)
library(mosaic)

theme_set(theme_bw())                 # change theme for ggplot2
trellis.par.set(theme=theme.mosaic()) # change default color scheme for lattice

knitr::opts_chunk$set(
    message = FALSE,
    warning = FALSE,
    tidy    = FALSE)
```

```{r}
set.seed(5678) # Initialize the random-number generator for reproducibility (pick any number)
### Load some helper functions
source("http://colindawson.net/stat213/code/helper_functions.R")
```

## Generating Synthetic Data From a Known Population Process

First we will generate synthetic data, again from a known model, so we have something to judge the performance of our metrics against.

The 'ground truth' model is the following **cubic polynomial**:

$$Y_i = 1 + 0.5 X_i - 0.5 X_i^2 + 0.25 X_i^3 + \varepsilon$$

where there is a single predictor variable X, and the residuals satisfy the usual regression conditions:

$$\varepsilon_i \sim N(0, 4)$$

The following script will generate the data as FakePolyData and plot the true function through the noised data

```{r}
source("http://colindawson.net/stat213/code/gen_poly_data.R")
```


## Fitting Polynomial Models to the Generated Data

Now we'll fit a series of 10 polynomial models to this data, ranging from degree 1 to degree 9.

The code below uses a function called `fit_polynomial_model()` that I wrote for this demo (the source code is in the `helper_functions.R` script that we loaded with `source()` above)

```{r}
m1 <- fit_polynomial_model(degree = 1, linear_formula = y ~ x, data = FakePolyData)
m2 <- fit_polynomial_model(degree = 2, linear_formula = y ~ x, data = FakePolyData)
m3 <- fit_polynomial_model(degree = 3, linear_formula = y ~ x, data = FakePolyData)
m4 <- fit_polynomial_model(degree = 4, linear_formula = y ~ x, data = FakePolyData)
m5 <- fit_polynomial_model(degree = 5, linear_formula = y ~ x, data = FakePolyData)
m6 <- fit_polynomial_model(degree = 6, linear_formula = y ~ x, data = FakePolyData)
m7 <- fit_polynomial_model(degree = 7, linear_formula = y ~ x, data = FakePolyData)
m8 <- fit_polynomial_model(degree = 8, linear_formula = y ~ x, data = FakePolyData)
m9 <- fit_polynomial_model(degree = 9, linear_formula = y ~ x, data = FakePolyData)
```

Let's look at the coefficients from one of the models:

```{r}
m3 %>% coef() %>% round(2)
```

The 'ground truth' model is the following **cubic polynomial**:

$$Y_i = 1 + 0.5 X_i - 0.5 X_i^2 + 0.25 X_i^3 + \varepsilon_i$$

The fitted cubic model is:

$$Y_i = 0.96 + 2.70 X_i - 0.59 X_i^2 + 0.13 X_i^3 + \varepsilon_i$$

which looks like this, plotted against the data:

```{r}
f_hat_3 <- makeFun(m3)
gf_point(y ~ x, data = FakePolyData) +
  stat_function(fun = f_hat_3)
```

**Aside:** We could also have done this in one line using R's handy `lapply()` function:

```{r}
modelList <- lapply(
  1:9, 
  FUN            = fit_polynomial_model, 
  linear_formula = y ~ x,
  data           = FakePolyData)
```

(the code above calls `fit_polynomial_model()`, setting the first argument first to 1, then to 2, then to 3, etc., keeping the `linear_formula=` and `data=` arguments constant, and returns the results of all nine calls of the function in a list)

Let's plot the nine fitted model curves on the data. The function `plot_data_with_polynomial_model()` is also a custom one I wrote for this demo, and is defined in `helper_functions.R`. Again, for conciseness I'm calling the function ten times all at once by using `lapply()`, but we could call it one by one.

```{r, results = 'hide'}
lapply(
  modelList, 
  FUN  = plot_data_with_polynomial_model, 
  data = FakePolyData,
  xlim = c(-5,5),
  ylim = c(-60,40))
```

## Evaluating the Models Using In-Sample Fit Metrics

Next let's calculate some fit metrics, such as $R^2$ and adjusted $R^2$, for these ten models. The `rsquared_adj()` function is defined in `helper_functions.R`.

```{r}
## Compute metrics for each model and assemble the results into a data table
fit_metrics <- tibble(
  ModelID                 = factor(1:9),
  R_squared               = sapply(modelList, FUN = rsquared),
  adjusted_R_squared      = sapply(modelList, FUN = rsquared_adj))
```

Let's see how regular and adjusted $R^2$ compare:

```{r}
### Plot the R^2 results
fit_metrics %>%
  pivot_longer(
    cols      = c(R_squared, adjusted_R_squared),
    names_to  = "FitMeasure", 
    values_to = "Value") %>%  
  gf_point(
    Value ~ ModelID,
    color = ~FitMeasure,
    xlab  = "Polynomial Order",
    ylab  = "Proportion of Variability Accounted For (higher is better)") +
  scale_color_discrete(
    name = "Fit Measure",
    labels = list(
      adjusted_R_squared = expression(paste("Adjusted ", R^2)),
      R_squared          = expression(R^2)))
```

Unadjusted $R^2$ of course gets higher and higher as we add complexity. Adjusted $R^2$ maxes out round polynomial degree 5 or 6 and then falls off. 

Recall that the true generating function was a _third_ order polynomial though, so this is still overfitting to some extent.

## Evaluating The Models Using Out of Sample Validation

We're going to build up to cross-validation, but we'll go one step at a time so you can see what we're doing at each step.

**The code below is more step-by-step than you'd really need** -- I'll show some 'all-in-one' code at the end, but it's useful to go through it step-by-step to understand what's going on.

### Simple Validation

Before we do **cross**-validation, let's try just doing **validation** with a single training set and a single validation set.

First, we split the data randomly, taking 10% case as the validation set (since this dataset only has 12 datapoints, the validation set will only have 1).

```{r}
N <- nrow(FakePolyData) # sample size

ValidationSet <- FakePolyData %>% 
  slice_sample(prop = 1/10)

TrainSet <- setdiff(FakePolyData, ValidationSet)
```

Now that we've split the data, let's fit some models on the training set.

```{r}
### Now fit some polynomials on the training set only
modelsTrain <- lapply(
  1:9, 
  FUN            = fit_polynomial_model,
  linear_formula = y ~ x,
  data           = TrainSet)
names(modelsTrain) <- 1:9
```

and plot them with the _full_ dataset:

```{r, results = 'hide'}
lapply(
  modelsTrain, 
  FUN  = plot_data_with_polynomial_model,
  data = FakePolyData,
  xlim = c(-5, 5),
  ylim = c(-60, 40))
```

Can you guess which point was held out?

Let's plot the data with just the validation point:

```{r, results = 'hide'}
lapply(
  modelsTrain, 
  FUN  = plot_data_with_polynomial_model,
  data = ValidationSet,
  xlim = c(-5, 5),
  ylim = c(-60,40))
```

Now we'll compute the mean squared prediction error for the held out point for each of the models. The `sapply()` function below does the same thing as `lapply()`, it just "simplifies" the list of results to a format that fits in a column of a data table. We're appending these generalization error measures to the table of other fit metrics we calculated before.

```{r}
fit_metrics <- fit_metrics %>%
  mutate(
    generalization_mse  = sapply(modelsTrain, FUN = mean_squared_error, dataset = ValidationSet),
    generalization_rmse = sapply(modelsTrain, FUN = root_mse,           dataset = ValidationSet))
  
```

Here's the Mean Squared (Generalization) Error for the 9 models:

```{r, results = 'hide'}
fit_metrics %>%
  gf_point(
    generalization_mse ~ ModelID, 
    xlab  = "Polynomial Order",
    ylab  = "Validation Mean Squared Error (lower is better)")
```

As we saw from the plots, the 9th order model is doing terribly.

The scale of the other values is tight, so let's plot the square root of the mean squared error instead.

```{r}
fit_metrics %>%
  gf_point(
    generalization_rmse ~ ModelID, 
    xlab  = "Polynomial Order",
    ylab  = "Validation Root MSE (lower is better)")
```

So for this particular validation fold, the 7th order polynomial generalizes the best, though the 4th order gets very close.

Compare this to the mean squared prediction error on the training set itself:

```{r}
fit_metrics <- fit_metrics %>%
  mutate(
    train_rmse = sapply(modelsTrain, root_mse, dataset = TrainSet))
```

```{r, results = 'hide'}
fit_metrics %>%
  select(ModelID, generalization_rmse, train_rmse) %>%
  pivot_longer(
    cols      = c(generalization_rmse, train_rmse),
    names_to  = "Metric",
    values_to = "Value") %>%
  gf_point(
    Value ~ ModelID, 
    color = ~Metric,
    xlab  = "Polynomial Order",
    ylab  = "Error Metric (lower is better)")
```


## Cross-Validation

The above is based on just one fixed validation point. What about cross-validation, where we let each point have a turn as the validation set?

The `do_k_fold_cv()` function will split the dataset you supply into `K` folds (for the value of `K` you supply as an argument), and do what we did above with each fold taking its turn as the validation set. It will calculate the error metric you supply for each of the validation folds, and return the average result.

Here are the results for model 5:

```{r}
MSE_CV_m5 <- do_k_fold_cv(
    model_formula  = formula(m5), 
    dataset        = FakePolyData, 
    K              = 12, 
    error_function = mean_squared_error)
RMSE_CV_m5 = sqrt(MSE_CV_m5)
MSE_CV_m5
RMSE_CV_m5
```

We can compute this for all of our different models and plot the results:

```{r}
formula_list <- sapply(modelList, formula)  ## Extracts the model formula from each of our models
fit_metrics <- fit_metrics %>% 
  mutate(
    MSE_CV = sapply(  ## Calling do_k_fold_cv() with each of the nine model formulas
      formula_list,
      FUN            = do_k_fold_cv,
      dataset        = FakePolyData,
      K              = N,
      error_function = mean_squared_error),
    RMSE_CV = sqrt(MSE_CV))
```

Now let's plot the resulting error measures.

```{r, results = 'hide'}
fit_metrics %>%
  gf_point(
    RMSE_CV ~ ModelID, 
    data  = results, 
    xlab  = "Polynomial Order",
    ylab  = "Leave-One-Out Root Mean Squared Error (lower is better)")
```

It's kind of hard to see the differences among the smaller values due to how massively larger the errors are for the most complex models, so let's just chop those off...

```{r, results = 'hide'}
fit_metrics %>%
  filter(!(ModelID %in% c(7,8,9))) %>%
  gf_point(
    RMSE_CV ~ ModelID,
    data  = results,
    xlab  = "Polynomial Order",
    ylab  = "Leave-One-Out Root Mean Squared Error (lower is better)")
```

And lo and behold, the 3rd order polynomial does the best! (We may have gotten a bit lucky here; this won't always return the "best" model by any means, but it's about as good a method as we have)

We don't actually have a test set in this case, but if we did we could evaluate prediction error for the final model on the test set like this.

```{r, eval = FALSE}
final_model <- lm(formula(m3), data = FakePolyData)
mean_squared_error(fitted_model = final_model, dataset = TestSet)
```

