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

```{r, warning=FALSE, message=FALSE, echo = FALSE}
library(dplyr)
library(ggplot2)
library(gridExtra)
library(purrr)
```

Most examples here are adapted from examples discussed at https://daviddalpiaz.github.io/r4sl/knn-reg.html (this is a useful companion to our text).

Suppose we have the following data:

```{r, echo = FALSE}
train_data <- data.frame(
  x = 1:5, 
  y = c(2, 1, 3, 5, 5)
)
```

```{r, fig.height=1.5}
train_data
ggplot(data = train_data, mapping = aes(x = x, y = y)) +
  geom_point()
```

### 2 Basic Approaches to Estimating $f(x)$:

1. Specify a model like $f(x) = \beta_0 + \beta_1 x$.  Estimate the *parameters* $\beta_0$ and $\beta_1$.
2. Local Approach: $f(x_0)$ should look like the data in a neighborhood of $x_0$

Models that don't specify a specific parametric form for $f$ are often called *nonparametric*.

### K Nearest Neighbors

 * The predicted value at a test point $x_0$ is the average of the $K$ training set observations that are closest to $x_0$ (the K nearest neighbors).

$$\hat{f}(x_0) = \frac{1}{K} \sum_{i \in N_0^{(k)}} y_i$$

 * Here $N_0^{(k)}$ is a set of indices for the $k$ observations that have values $x_i$ nearest to the test point $x_0$.

Using our example data above:

 * Suppose we want to make a prediction at the test point $x_0 = 3.75$
 * Set $k = 3$
 * In our training data set, what are the $k$ nearest neighbors to the test point?

\vspace{2cm}

 * What is the fitted/predicted value at $x_0$?
 
\vspace{2cm}

\newpage

The `train` and `predict` functions in the `caret` package will do these calculations for us:

```{r, fig.height = 2.4, message = FALSE}
library(caret)

# "train" the KNN model -- but note there are no parameters to estimate like in linear regression!
knn_fit <- train(
  form = y ~ x,
  data = train_data,
  method = "knn",
  trControl = trainControl(method = "none"),
  tuneGrid = data.frame(k = 3)
)

# get test set predictions.  here we have just one point in the test set
test_data <- data.frame(
  x = 3.75
)
test_data
predict(knn_fit, newdata = test_data)
```

```{r, fig.height = 3}
# to make a plot of the fitted function f hat, we can set up a function to get predicted values
predict_knn <- function(x) {
  predict(knn_fit, newdata = data.frame(x = x))
}

ggplot(data = train_data, mapping = aes(x = x, y = y)) +
#  geom_line(data = test_data_for_plot, mapping = aes(x = x, y = y_hat)) +
  stat_function(fun = predict_knn) +
  geom_point() +
  geom_point(mapping = aes(x = 3.75, y = 4.333333), color = "orange", size = 4)
```

\newpage

### Flexibility is determined by $k$

Recall we previously fit models for the relationship between nitrogen oxides concentrations and distance from Boston employment centers in 506 neighborhoods around Boston.  Below are predictions from KNN models with varying values of $k$.

```{r, fig.height = 6, warning = FALSE, message=FALSE, echo=FALSE}
library(MASS) # for Boston data

# function to make a plot
make_knn_plot <- function(k) {
  # "train" the KNN model
  knn_fit <- train(
    form = nox ~ dis,
    data = Boston,
    use.all = FALSE,
    method = "knn",
    trControl = trainControl(method = "none"),
    tuneGrid = data.frame(k = k) # note: using the value of k specific to this iteration of the loop
  )
  
  # function to get predictions
  predict_knn <- function(x) {
    predict(knn_fit, newdata = data.frame(dis = x))
  }
  
  p <- ggplot(data = Boston, mapping = aes(x = dis, y = nox)) +
    geom_point() +
    stat_function(fun = predict_knn, color = "orange") +
    ggtitle(paste0("k = ", k))
  
  return(p)
}

# plot all the plots
grid.arrange(
  make_knn_plot(1),
  make_knn_plot(5),
  make_knn_plot(10),
  make_knn_plot(50),
  make_knn_plot(100),
  make_knn_plot(250)
)
```

 * Is the bias of KNN lower for large $k$ or small $k$?

\vspace{2cm}

 * Is the variance of KNN lower for large $k$ or small $k$?

\vspace{2cm}

 * How should we choose the value of $k$?

\vspace{2cm}

\newpage

### KNN Automatically Adjusts to Different Functional Forms

I simulated three fake data sets of size $n = 100$:

 * One where the true function is linear
 * A second where the true function is quadratic
 * A third where the true function is sinusoidal

In all cases, a KNN fit with $k = 10$ nearest neighbors (shown in orange) does a good job at recovering the underlying function (shown in black).

```{r, fig.height = 6, echo = FALSE}
library(FNN)
line_reg_fun <- function(x) {
  x
}

quad_reg_fun <- function(x) {
  x ^ 2
}

sine_reg_fun <- function(x) {
  sin(x)
}

get_plot_knn_vs_true <- function(reg_fun, sample_size = 100, noise_sd = 1, k = 10) {
  x <- runif(n = sample_size, min = -5, max = 5)
  y <- rnorm(n = sample_size, mean = reg_fun(x), sd = noise_sd)
  train_data <- data.frame(x = x, y = y)
  
  test_x_grid <- data.frame(x = seq(-5, 5, by = 0.01))
  
  test_results <- test_x_grid %>%
    mutate(
      y_true_mean = reg_fun(x),
      y_knn = FNN::knn.reg(
        train = train_data %>% dplyr::select(x),
        test = test_x_grid,
        y = train_data %>% dplyr::select(y),
        k = 10)$pred
    )
  
  plot_object <- ggplot(mapping = aes(x = x)) +
    geom_point(data = train_data, mapping = aes(y = y)) +
    geom_line(data = test_results, mapping = aes(y = y_true_mean), size = 2) +
    geom_line(data = test_results, mapping = aes(y = y_knn), color = "orange", size = 1)
  
  return(plot_object)
}

set.seed(42)
line_plot <- get_plot_knn_vs_true(line_reg_fun, sample_size = 100, noise_sd = 1, k = 10)
quad_plot <- get_plot_knn_vs_true(quad_reg_fun, sample_size = 100, noise_sd = 2, k = 10)
sine_plot <- get_plot_knn_vs_true(sine_reg_fun, sample_size = 100, noise_sd = 0.5, k = 10)

grid.arrange(line_plot, quad_plot, sine_plot, ncol = 1)
```

\newpage

### KNN with multiple explanatory variables

The basic idea is the same, but now we use Euclidean distance to find the nearest neighbors.

```{r, echo = FALSE}
train_data <- data.frame(
  x1 = 1:5,
  x2 = c(2, 3, 1, 5, 4),
  y = c(2, 1, 3, 5, 5)
)
```

In our fake example, suppose we now have two x variables:

```{r}
train_data
```

We now want to make a prediction at the test set point $x_0 = (3.75, 2)$ using the $k = 3$ nearest neighbors.

We need to first find the Euclidean distance of each training set observation from the test set point:

```{r}
train_data <- train_data %>% mutate(
  distance_from_test = sqrt((x1 - 3.75)^2 + (x2 - 2)^2)
)
train_data
```

 * What are the $k = 3$ nearest neighbors to the test point?

\vspace{2cm}

 * What is the fitted/predicted value at the test point?

\newpage

### If $p > 1$, KNN performance often improves if we scale the explanatory variables

 * Divide by standard deviation, so the rescaled variable has standard deviation 1
 * Consider an example predicting based on two variables:
    * `tax`: full-value property-tax rate per \$10,000
    * `dis`: weighted mean of distances to five Boston employment centres

```{r}
# train/test split
set.seed(76520)
train_inds <- caret::createDataPartition(Boston$nox, p = 0.8)
train_boston <- Boston %>% slice(train_inds[[1]])
test_boston <- Boston %>% slice(-train_inds[[1]])

# rescale.  Here I rescale based on just training set data standard deviations,
# but actually it's ok to rescale based on combined train/test standard deviations
scale_sds <- train_boston %>%
  summarize(
    dis_sd = sd(dis),
    tax_sd = sd(tax)
  )
scale_sds

train_boston_scaled <- train_boston %>%
  mutate(
    dis = dis / scale_sds$dis_sd,
    tax = tax / scale_sds$tax_sd
  )
train_boston_scaled %>%
  summarize(
    dis_sd = sd(dis),
    tax_sd = sd(tax)
  )

test_boston_scaled <- test_boston %>%
  mutate(
    dis = dis / scale_sds$dis_sd,
    tax = tax / scale_sds$tax_sd
  )
```

\newpage

```{r}
# KNN fit from scaled data
knn_fit_scaled <- train(
  form = nox ~ dis + tax,
  data = train_boston_scaled,
  use.all = FALSE,
  method = "knn",
  trControl = trainControl(method = "cv"),
  tuneGrid = data.frame(k = 5)
)

# KNN fit from original data
knn_fit_orig <- train(
  form = nox ~ dis + tax,
  data = train_boston,
  use.all = FALSE,
  method = "knn",
  trControl = trainControl(method = "cv"),
  tuneGrid = data.frame(k = 5)
)
```

```{r}
# test set RMSE, scaled data fit
(test_boston_scaled$nox - predict(knn_fit_scaled, newdata = test_boston_scaled))^2 %>%
  mean() %>%
  sqrt()

# test set RMSE, original data fit
(test_boston$nox - predict(knn_fit_orig, newdata = test_boston))^2 %>%
  mean() %>%
  sqrt()
```

Actually, we could have had `caret` do the scaling for us by passing in a `preProcess = "scale"` argument to `train`:

```{r}
# KNN fit from scaled data, but caret does the scaling
knn_fit_scaled_by_caret <- train(
  form = nox ~ dis + tax,
  data = train_boston, # note I'm giving train my original data frame
  use.all = FALSE,
  method = "knn",
  preProcess = "scale", # this is the only new line
  trControl = trainControl(method = "none"),
  tuneGrid = data.frame(k = 5)
)

# test set RMSE, scaled data fit, scaling done by caret
# note that since caret is handling scaling, I give it my original data for prediction
(test_boston$nox - predict(knn_fit_scaled_by_caret, newdata = test_boston))^2 %>%
  mean() %>%
  sqrt()
```

\newpage

### Curse of Dimensionality: KNN performance degrades relatively quickly as the number of explanatory variables $p$ increases

 * Degradation in performance affects all methods, but affects non-parametric methods more
    * Parametric models assume a restricted parametric form for $f$ and are trying to learn only a few parameters
    * Non-parametric methods are trying to learn the functional form.  This is more difficult in higher dimensions
    * You will have a homework problem about this.

```{r, echo = FALSE}
all_x_vars <- names(Boston %>% dplyr::select(-nox, -chas))
x_vars_cor_nox <- sapply(all_x_vars, function(var_name) {
  cor(Boston[[var_name]], Boston[["nox"]])
})
vars_to_include <- data.frame(
  x_var = all_x_vars,
  cor_nox = x_vars_cor_nox,
  abs_cor_nox = abs(x_vars_cor_nox)
) %>%
  dplyr::arrange(desc(abs_cor_nox)) %>%
  dplyr::select(-abs_cor_nox) %>%
  as.data.frame()

vars_to_include <- vars_to_include[1:10, ]
```

To explore this, I pulled out the top 10 explanatory variables in the Boston data set that were most correlated with `nox`, and sorted them in decreasing order of (absolute value of) correlation:

```{r}
vars_to_include
```

```{r, echo = FALSE}
get_knn_rmse_top_p <- function(p) {
  vars_to_use <- vars_to_include$x_var[seq_len(p)]
  model_formula <- as.formula(paste0("nox ~ ", paste(vars_to_use, collapse = " + ")))
  
  knn_fit_scaled_by_caret <- train(
    form = model_formula,
    data = train_boston,
    use.all = FALSE,
    method = "knn",
    preProcess = "scale", # this is the only new line
    trControl = trainControl(method = "none"),
    tuneGrid = data.frame(k = 10)
  )
  
  # test set RMSE, scaled data fit, scaling done by caret
  # note that since caret is handling scaling, I give it my original data for prediction
  (test_boston$nox - predict(knn_fit_scaled_by_caret, newdata = test_boston))^2 %>%
    mean() %>%
    sqrt() %>%
    return()
}

get_lm_rmse_top_p <- function(p) {
  vars_to_use <- vars_to_include$x_var[seq_len(p)]
  model_formula <- as.formula(paste0("nox ~ poly(", paste(vars_to_use, collapse = ", 2) + poly("), ", 2)"))
  
  lm_fit <- train(
    form = model_formula,
    data = train_boston,
    method = "lm",
    trControl = trainControl(method = "none")
  )
  
  # test set RMSE, scaled data fit, scaling done by caret
  # note that since caret is handling scaling, I give it my original data for prediction
  (test_boston$nox - predict(lm_fit, newdata = test_boston))^2 %>%
    mean() %>%
    sqrt() %>%
    return()
}

rmse_results <- data.frame(
  p = seq_len(nrow(vars_to_include)),
  knn_rmse = sapply(seq_len(nrow(vars_to_include)), get_knn_rmse_top_p),
  lm_rmse = sapply(seq_len(nrow(vars_to_include)), get_lm_rmse_top_p)
)
```

I then fit a sequence of KNN models (k = 10) and linear models where each explanatory variable entered with a degree 2 polynomial term (no interactions).  A model with $p$ features used the $p$ features that were most correlated with `nox`.  Here are the RMSE for each of these models:

```{r, fig.height = 2.5}
rmse_results
ggplot(data = rmse_results, mapping = aes(x = p)) +
  geom_line(mapping = aes(y = knn_rmse), color = "orange") +
  geom_line(mapping = aes(y = lm_rmse))
```


```{r, eval = FALSE, echo = FALSE}
sim_knn_data <- function(n_obs = 50) {
  x1 <- seq(0, 10, length.out = n_obs)
  x2 <- runif(n = n_obs, min = 0, max = 10)
  x3 <- runif(n = n_obs, min = 0, max = 10)
  x4 <- runif(n = n_obs, min = 0, max = 10)
  x5 <- runif(n = n_obs, min = 0, max = 10)
  y <- x1 ^ 2 + rnorm(n = n_obs)
  data.frame(y, x1, x2, x3, x4, x5)
}

set.seed(42)
knn_data_train <- sim_knn_data()
knn_data_test <- sim_knn_data()

# define helper function for getting knn.reg predictions
# note: this function is highly specific to this situation and data set
get_test_rmse <- function(p = 1, k = 5) {
  x_vars <- paste0("x", seq_len(p))
  
  x_train <- knn_data_train %>%
    dplyr::select(x_vars)
  x_test <- knn_data_test %>%
    dplyr::select(x_vars)
  scale_values <- x_train %>%
    summarize_all(sd)
  x_train <- scale(x_train, scale = scale_values) %>%
    as.data.frame()
  x_test <- scale(x_test, scale = scale_values) %>%
    as.data.frame()
  
  y_train <- knn_data_train %>%
    dplyr::select(y)
  y_test <- knn_data_test %>%
    dplyr::select(y)
  
  pred_knn <- FNN::knn.reg(
    train = x_train,
    test = x_test,
    y = y_train,
    k = k)$pred
  
  lm_fit <- lm(y ~ ., data = cbind(y_train, x_train))
  pred_lm <- predict(lm_fit, newdata = x_test)
  
  data.frame(
    p = p,
    knn = rmse(actual = y_test$y, predicted = pred_knn),
    lm = rmse(actual = y_test$y, predicted = pred_lm)
  )
}

cod_results <- map_dfr(
  seq_len(5),
  get_test_rmse
)

ggplot(data = cod_results, mapping = aes(x = p)) +
  geom_line(mapping = aes(y = knn), color = "orange") +
  geom_line(mapping = aes(y = lm), color = "cornflowerblue") +
  ylab("RMSE")
```

