---
title: "Supervised Learning"
subtitle: "From nonparametric regression to GAMs"
date: "July 14th, 2022"
---

```{r setup, echo = FALSE}
knitr::opts_chunk$set(fig.retina = 3, warning = FALSE, message = FALSE)
```

```{r xaringan-themer, include=FALSE, warning=FALSE}
library(xaringanthemer)
style_mono_accent(base_color = "#cc002b")
```

## Model flexibility vs interpretability

[Figure 2.7, Introduction to Statistical Learning with Applications in R (ISLR)](https://www.statlearning.com/)

```{r out.width='50%', echo = FALSE, fig.align='center'}
knitr::include_graphics("http://www.stat.cmu.edu/~pfreeman/flexibility.png")
```

__Tradeoff__ between model's _flexibility_ (i.e. how "curvy" it is) and how __interpretable__ it is

- Simpler, parametric form of the model $\Rightarrow$ the easier it is to interpret

---

## Model flexibility vs interpretability

```{r out.width='50%', echo = FALSE, fig.align='center'}
knitr::include_graphics("http://www.stat.cmu.edu/~pfreeman/flexibility.png")
```

- __Parametric__ models, for which we can write down a mathematical expression for $f(X)$ __before observing the data__, _a priori_ (e.g. linear regression), __are inherently less flexible__


--

- __Nonparametric__ models, in which $f(X)$ is __estimated from the data__ (e.g. kernel regression)

---

## K Nearest Neighbors (KNN)

- Find the $k$ data points __closest__ to an observation $x$, use these to predit

  - Need to use some measure of distance, e.g., Euclidean distance

--

- KNN is data-driven, but we can actually write down the model _a priori_

--

- Regression:
$$
{\hat Y} \vert X = \frac{1}{k} \sum_{i=1}^k Y_i \,,
$$
- Classification:
$$
\hat{P}[Y = j \vert X] = \frac{1}{k} \sum_{i=1}^k 1(Y_i = j) \,,
$$

  - $1(\cdot)$ is the indicator function: returns 1 if TRUE, and 0 otherwise. 

  - Summation yields the proportion of neighbors that are of class $j$

---

## Finding the optimal number of neighbors $k$

__The number of neighbors $k$ is a tuning parameter__ (like $\lambda$ is for ridge / lasso)

--

Determining the optimal value of $k$ requires balancing bias and variance:

- If $k$ is too small, the resulting model is *too flexible*,

  - low bias (it is right on average...if we apply KNN to an infinite number of datasets sampled from the same parent population) 
  
  - high variance (the predictions have a large spread in values when we apply KNN to our infinite data). See the panels to the left on the next slide.

--

- If $k$ is too large, the resulting model is *not flexible enough*, 

  - high bias (wrong on average) and 
  
  - low variance (nearly same predictions, every time). See the panels to the right on the next slide.

---

## Finding the optimal number of neighbors $k$

```{r out.width='40%', echo = FALSE, fig.align='center'}
knitr::include_graphics("http://www.stat.cmu.edu/~pfreeman/Fig_3.16.png")
```

```{r out.width='40%', echo = FALSE, fig.align='center'}
knitr::include_graphics("http://www.stat.cmu.edu/~pfreeman/Fig_2.16.png")
```


(Figures 3.16 [top] and 2.16 [bottom], *Introduction to Statistical Learning* by James et al.)


---

## What does KNN remind you of?...

```{r out.width='40%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://media1.giphy.com/media/12Gyz2J1b9SjD2/200w.gif")
```

---

## Kernels

A kernel $K(x)$ is a weighting function used in estimators, and technically has only one required property:

- $K(x) \geq 0$ for all $x$

However, in the manner that kernels are used in statistics, there are two other properties that are usually satisfied:

- $\int_{-\infty}^\infty K(x) dx = 1$; and 

- $K(-x) = K(x)$ for all $x$.

In short: __a kernel is a symmetric PDF!__

---

## Kernel density estimation

__Goal__: estimate the PDF $f(x)$ for all possible values (assuming it is continuous / smooth)

--

$$
\text{Kernel density estimate: } \hat{f}(x) = \frac{1}{n} \sum_{i=1}^n \frac{1}{h} K_h(x - x_i)
$$

--
- $n =$ sample size, $x =$ new point to estimate $f(x)$ (does NOT have to be in dataset!)


--
- $h =$ __bandwidth__, analogous to histogram bin width, ensures $\hat{f}(x)$ integrates to 1

- $x_i =$ $i$th observation in dataset


--
- $K_h(x - x_i)$ is the __Kernel__ function, creates __weight__ given distance of $i$th observation from new point 
  - as $|x - x_i| \rightarrow \infty$ then $K_h(x - x_i) \rightarrow 0$, i.e. further apart $i$th row is from $x$, smaller the weight
  
  - as __bandwidth__ $h \uparrow$ weights are more evenly spread out (as $h \downarrow$ more concentrated around $x$) 

  - typically use [__Gaussian__ / Normal](https://en.wikipedia.org/wiki/Normal_distribution) kernel: $\propto e^{-(x - x_i)^2 / 2h^2}$
  
  - $K_h(x - x_i)$ is large when $x_i$ is close to $x$
  

---

## Commonly Used Kernels

```{r out.width='40%', echo = FALSE, fig.align='center'}
knitr::include_graphics("http://www.stat.cmu.edu/~pfreeman/kernels.png")
```


A general rule of thumb: the choice of kernel will have little effect on estimation, particularly if the sample size is large! The Gaussian kernel (i.e., a normal PDF) is by far the most common choice, and is the default for `R` functions that utilize kernels.


---

## Kernel regression

We can apply kernels in the regression setting as well as in the density estimation setting!

The classic kernel regression estimator is the __Nadaraya-Watson__ estimator:

$$\hat{y}_h(x) = \sum_{i=1}^n w_i(x) Y_i \,,$$

where
$$w_i(x) = \frac{K\left(\frac{x-X_i}{h}\right)}{\sum_{j=1}^n K\left(\frac{x-X_j}{h}\right)} \,.$$

Regression estimate is the average of all the *weighted* observed response values; 

- Farther $x$ is from observation $\Rightarrow$ less weight that observation has in determining the regression estimate at $x$

---

## Kernel regression

__Nadaraya-Watson kernel regression__

- given training data with explanatory variable $x$ and continuous response $y$ 

- _bandwidth_ $h > 0$

- and a new point $(x_{new}, y_{new})$: 

$$\hat{y}_{new} = \sum_{i=1}^n w_i(x_{new})\cdot y_i \,,$$

where

$$w_i(x) = \frac{K_h\left(|x_{new}-x_i|\right)}{\sum_{j=1}^n K_h\left(|x_{new}-x_j|\right)} \text{ with } K_h(x) = K(\frac{x}{h})$$

Example of a __linear smoother__

- class of models where predictions are _weighted_ sums of the response variable

---

## Local regression

We can fit a linear model __at each point $x_{new}$__ with weights given by kernel function centered on $x_{new}$

- we can additionally combine this with _polynomial regression_

--

Local regression of the $k^{th}$ order with kernel function $K$ solves the following:

$$\hat{\beta}(x_{new}) = \underset{\beta}{\text{arg min}}\Big\{ \sum_i K_h(|x_{new} - x_i|) \cdot (y_i - \sum_{j=0}^k x_i^k \cdot \beta_k )^2 \Big\}$$

--
__Yes, this means every single observation has its own set of coefficients__


--

Predicted value is then:

$$\hat{y}_{new} = \sum_{j=0}^k x_{new}^k \cdot \hat{\beta}_k(x_{new})$$


--

Smoother predictions than kernel regression but comes at __higher computational cost__ 

- __LOESS__ replaces kernel with k nearest neighbors
  
  - faster than local regression but discontinuities when neighbors change


---

## Smoothing splines


Use __smooth function__ $s(x)$ to predict $y$, control smoothness directly by minimizing the __spline objective function__:

$$\sum_{i=1}^n (y_i - s(x_i))^2 + \lambda \int(s''(x))^2dx$$


--

$$= \text{fit data} + \text{impose smoothness}$$

--

$$\Rightarrow \text{model fit} = \text{likelihood} - \lambda \cdot \text{wiggliness}$$
Estimate the __smoothing spline__ $\hat{s}(x)$ that __balances the tradeoff between the model fit and wiggliness__


--
```{r out.width='70%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://github.com/noamross/gams-in-r-course/blob/master/images/diffsmooth-1.png?raw=true")
```

---

## Basis functions

Splines are _piecewise cubic polynomials_ with __knots__ (boundary points for functions) at every data point


--

Practical alternative: linear combination of set of __basis functions__


--
__Cubic polynomial example__: define four basis functions:

- $B_1(x) = 1$, $B_2(x) = x$, $B_3(x) = x^2$, $B_4(x) = x^3$

where the regression function $r(x)$ is written as:

$$r(x) = \sum_j^4 \beta_j B_j(x)$$


--
- __linear in the transformed variables__ $B_1(x), \dots, B_4(x)$ but it is __nonlinear in $x$__


--
We extend this idea for splines _piecewise_ using indicator functions so the spline is a weighted sum:

$$s(x) = \sum_j^m \beta_j B_j(x)$$

---

## Number of basis functions is another tuning parameter

```{r out.width='80%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://github.com/noamross/gams-in-r-course/blob/master/images/diffbasis-1.png?raw=true")
```



---

## Generalized additive models (GAMs)

GAMs were created by [Trevor Hastie and Rob Tibshirani in 1986](https://projecteuclid.org/euclid.ss/1177013604) with intuitive construction:

- relationships between individual explanatory variables and the response variable are smooth (either linear or nonlinear via basis functions)

- estimate the smooth relationships __simultaneously__ to predict the response by just adding them up

__Generalized__ like GLMs where $g()$ is the link function for the expected value of the response $E(Y)$ and __additive__ over the $p$ variables:

$$g(E(Y)) = \beta_0 + s_1(x_1) + s_2(x_2) + \dots + s_p(x_p)$$

```{r out.width='70%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://multithreaded.stitchfix.com/assets/images/blog/fig1.svg")
```


--

- can be a convenient balance between flexibility and interpretability

- you can combine linear and nonlinear terms!

---

## Example: predicting MLB HR probability

Used the [`baseballr`](http://billpetti.github.io/baseballr/) package to scrape all batted-balls from 2022 season for EDA project:

```{r load-data, warning = FALSE, message = FALSE}
library(tidyverse)
batted_ball_data <- read_csv("http://www.stat.cmu.edu/cmsac/sure/2022/materials/data/sports/eda_projects/mlb_batted_balls_2022.csv") %>%
  mutate(is_hr = as.numeric(events == "home_run")) %>%
  filter(!is.na(launch_angle), !is.na(launch_speed),
         !is.na(is_hr))
head(batted_ball_data)
table(batted_ball_data$is_hr)
```

---

## Predict HRs with launch angle and exit velocity?

.pull-left[
```{r plot-baseball, eval = FALSE}
batted_ball_data %>%
  ggplot(aes(x = launch_speed, 
             y = launch_angle,
             color = as.factor(is_hr))) +
  geom_point(alpha = 0.5) +
  ggthemes::scale_color_colorblind(labels = c("No", "Yes")) +
  labs(x = "Exit velocity", 
       y = "Launch angle", 
       color = "HR?") +
  theme_bw() +
  theme(legend.position = "bottom")
```

- HRs are relatively rare and confined to one area of this plot

]
.pull-right[
```{r ref.label='plot-baseball', echo = FALSE}

```

]

---

## Fitting GAMs with [`mgcv`](https://cran.r-project.org/web/packages/mgcv/index.html)

First set-up training data

```{r init-train}
set.seed(2004)
batted_ball_data <- batted_ball_data %>%
  mutate(is_train = sample(rep(0:1, length.out = nrow(batted_ball_data))))
```

Next fit the initial function using smooth functions via `s()`:

```{r init-gam}
library(mgcv)
init_logit_gam <- gam(is_hr ~ s(launch_speed) + s(launch_angle), #<<
                      data = filter(batted_ball_data, is_train == 1), 
                      family = binomial, method = "REML") #<<
```

 - Use [REML](https://en.wikipedia.org/wiki/Restricted_maximum_likelihood) instead of the default for more stable solution

---

### GAM summary

```{r gam-summary}
summary(init_logit_gam)
```


---

### Visualizing partial response functions with [gratia](https://gavinsimpson.github.io/gratia/index.html)

Displays the partial effect of each term in the model $\Rightarrow$ add up to the overall prediction 

```{r, fig.align = 'center', fig.height=5}
library(gratia)
draw(init_logit_gam)
```

---

## Convert to probability scale with `plogis` function


```{r, fig.align = 'center', fig.height=5}
draw(init_logit_gam, fun = plogis) #<<
```

--

- centered on average value of 0.5 because it's the partial effect without the intercept

---

## Include intercept in plot...

```{r, fig.align = 'center', fig.height=5}
draw(init_logit_gam, fun = plogis, constant = coef(init_logit_gam)[1]) #<<
```

Intercept reflects relatively rare occurence of HRs!

---

## Model checking for number of basis functions

Use `gam.check()` to see if we need more basis functions based on an approximate test

```{r, fig.show='hide'}
gam.check(init_logit_gam)
```


---

## Check the predictions?

```{r}
batted_ball_data <- batted_ball_data %>%
  mutate(init_gam_hr_prob = 
           as.numeric(predict(init_logit_gam,
                              newdata = batted_ball_data,
                              type = "response")),
         init_gam_hr_class = as.numeric(init_gam_hr_prob >= 0.5))
batted_ball_data %>%
  group_by(is_train) %>%
  summarize(correct = mean(is_hr == init_gam_hr_class))
```

---

## What about the linear model?

```{r}
init_linear_logit <- glm(is_hr ~ launch_speed + launch_angle, 
                         data = filter(batted_ball_data, is_train == 1), 
                         family = binomial)
batted_ball_data <- batted_ball_data %>%
  mutate(init_glm_hr_prob = predict(init_linear_logit,
                                    newdata = batted_ball_data,
                                    type = "response"),
         init_glm_hr_class = as.numeric(init_glm_hr_prob >= 0.5))
batted_ball_data %>%
  group_by(is_train) %>%
  summarize(correct = mean(is_hr == init_glm_hr_class))
```


__Very few situations in reality where linear regressions perform better than an additive model using smooth functions__ - especially since smooth functions can just capture linear models...

---

## Continuous interactions as a smooth surface

```{r multi-gam}
multi_logit_gam <- gam(is_hr ~ s(launch_speed, launch_angle), 
                       data = filter(batted_ball_data, is_train == 1), 
                       family = binomial)
```

Plot the predicted heatmap:

```{r, fig.align ='center', fig.height=4}
draw(multi_logit_gam)
```

---

## Check the predictions?

```{r}
batted_ball_data <- batted_ball_data %>%
  mutate(multi_gam_hr_prob =
           as.numeric(predict(multi_logit_gam,
                              newdata = batted_ball_data,
                              type = "response")),
         multi_gam_hr_class = as.numeric(multi_gam_hr_prob >= 0.5))
batted_ball_data %>%
  group_by(is_train) %>%
  summarize(correct = mean(is_hr == multi_gam_hr_class))
```


- This has one smoothing parameter for the 2D smooth...

---

### Separate interactions from individual terms with tensor smooths

```{r tensor-gam}
tensor_logit_gam <- gam(is_hr ~ s(launch_speed) + s(launch_angle) +
                          ti(launch_speed, launch_angle), #<<
                        data = filter(batted_ball_data, is_train == 1), family = binomial)
```

--

Check the predictions

```{r}
batted_ball_data %>%
  mutate(tensor_gam_hr_prob =
           as.numeric(predict(tensor_logit_gam, newdata = batted_ball_data, type = "response")),
         tensor_gam_hr_class = as.numeric(tensor_gam_hr_prob >= 0.5)) %>%
  group_by(is_train) %>% 
  summarize(correct = mean(is_hr == tensor_gam_hr_class))
```

--

__More complicated model but yet it does not help!__

---

## Recap and useful resources

.pull-left[

<blockquote class="twitter-tweet"><p lang="en" dir="ltr">140 char vrsn<br><br>1 GAMs are just GLMs<br>2 GAMs fit wiggly terms<br>3 use + s(foo) not foo in frmla<br>4 use method = &quot;REML&quot;<br>5 gam.check()</p>&mdash; Dr Gavin Simpson 😷🇪🇺🇩🇰 (@ucfagls) <a href="https://twitter.com/ucfagls/status/842444686513991680?ref_src=twsrc%5Etfw">March 16, 2017</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

]
.pull-right[

- [GAMs in R by Noam Ross](https://noamross.github.io/gams-in-r-course/)

- [`mgcv` course](https://eric-pedersen.github.io/mgcv-esa-workshop/)

- [Stitch Fix post: GAM: The Predictive Modeling Silver Bullet](https://multithreaded.stitchfix.com/blog/2015/07/30/gam/)

- Chapters 7 and 8 of [Advanced Data Analysis from an Elementary Point of View by Prof Cosma Shalizi](https://www.stat.cmu.edu/~cshalizi/ADAfaEPoV/ADAfaEPoV.pdf) 

  - I strongly recommend you download this book, and you will refer back to it for the rest of your life
  
]



---
class: inverse center middle

# BONUS CONTENT



---

## KNN in context

Here are two quotes from ISLR to keep in mind when thinking about KNN:

- "As a general rule, parametric methods [like linear regression] will tend to outperform non-parametric approaches [like KNN] when there is a small number of observations per predictor." This is the *curse of dimensionality*: for data-driven models, the amount of data you need to get similar model performance goes up exponentially with $p$.

--

$\Rightarrow$ KNN might not be a good model to learn when the number of predictor variables is very large.

--

- "Even in problems in which the dimension is small, we might prefer linear regression to KNN from an interpretability standpoint. If the test MSE of KNN is only slightly lower than that of linear regression, we might be willing to forego a little bit of prediction accuracy for the sake of a simple model..."

--

$\Rightarrow$ KNN is not the best model to learn if inference is the goal of an analysis.

---

## KNN: two critical points to remember

1. To determine which neighbors are the nearest neighbors, pairwise Euclidean distances are computed...so we may need to scale (or standardize) the individual predictor variables so that the distances are not skewed by that one predictor that has the largest variance.

--

2. Don't blindly compute a pairwise distance matrix! For instance, if $n$ = 100,000, then your pairwise distance matrix will have $10^{10}$ elements, each of which uses 8 bytes in memory...resulting in a memory usage of 80 GB! Your laptop cannot handle this. It can barely handle 1-2 GB at this point. If $n$ is large, you have three options:
    a. subsample your data, limiting $n$ to be $\lesssim$ 15,000-20,000;
    b. use a variant of KNN that works with sparse matrices (matrices that can be compressed since most values are zero); or
    c. make use of a "kd tree" to more effectively (but only approximately) identify nearest neighbors.
  
The [`FNN` package in `R`](https://daviddalpiaz.github.io/r4sl/knn-reg.html) has an option to search for neighbors via the use of a kd tree.


--
But instead we will use the [`caret`](http://topepo.github.io/caret/index.html) package...

---

## Example data: MLB 2021 batting statistics

Downloaded MLB 2021 batting statistics leaderboard from [Fangraphs](https://www.fangraphs.com/leaders.aspx?pos=all&stats=bat&lg=all&qual=y&type=8&season=2019&month=0&season1=2019&ind=0)

```{r load-data-mlb, warning = FALSE, message = FALSE}
library(tidyverse)
mlb_data <- read_csv("http://www.stat.cmu.edu/cmsac/sure/2022/materials/data/sports/fg_batting_2021.csv")
head(mlb_data)
```


---

## Data cleaning

- [`janitor`](http://sfirke.github.io/janitor/) package has convenient functions for data cleaning like `clean_names()`

- `parse_number()` function provides easy way to convert character to numeric columns

```{r clean-data}
library(janitor)
mlb_data_clean <- clean_names(mlb_data)
mlb_data_clean <- mlb_data_clean %>%
  mutate_at(vars(bb_percent:k_percent), parse_number)
head(mlb_data_clean)
```

---

## KNN example

`caret` is a package of functions designed to simplify training, tuning, and testing statistical learning methods

- first create partitions for training and test data using `createDataPartition()`

```{r split-data}
library(caret)
set.seed(1960)
train_i <- createDataPartition(y = mlb_data_clean$w_oba, p = 0.7, list = FALSE) %>%
  as.numeric()
train_mlb_data <- mlb_data_clean[train_i,]
test_mlb_data <- mlb_data_clean[-train_i,]
```

--

- next [`train()`](http://topepo.github.io/caret/model-training-and-tuning.html) to find the optimal `k` on the training data with cross-validation

```{r train-knn}
set.seed(1971)
init_knn_mlb_train <- train(w_oba ~ bb_percent + k_percent + iso, 
                            data = train_mlb_data, method = "knn",
                            trControl = trainControl("cv", number = 10),
                            preProcess = c("center", "scale"),
                            tuneLength = 10)
```

---

## KNN example

```{r plot-knn, fig.align ='center', fig.height=5}
ggplot(init_knn_mlb_train) + theme_bw()
```


---

## KNN example

Can manually create a __tuning grid__ to search over for the tuning parameter `k`

```{r tune-k}
set.seed(1979)
tune_knn_mlb_train <- train(w_oba ~ bb_percent + k_percent + iso, 
                            data = train_mlb_data, method = "knn",
                            trControl = trainControl("cv", number = 10),
                            preProcess = c("center", "scale"),
                            tuneGrid = expand.grid(k = 2:20)) #<<
tune_knn_mlb_train$results
```

---

## KNN example

```{r plot-knn-tune, fig.align ='center', fig.height=5}
ggplot(tune_knn_mlb_train) + theme_bw()
```

---

## KNN example

```{r test-preds}
tune_knn_mlb_train$bestTune
test_preds <- predict(tune_knn_mlb_train, test_mlb_data)
head(test_preds)
RMSE(test_preds, test_mlb_data$w_oba)
```


---

## Kernel regression with `np`

Use the `npregbw` function to tune bandwidth using [__generalized cross-validation__](https://bookdown.org/egarpor/NP-UC3M/kre-i-bwd.html#kre-i-bwd-cv)

```{r, eval = FALSE}
library(np)
mlb_bw0 <- npregbw(w_oba ~ bb_percent + k_percent + iso, 
                            data = train_mlb_data)
```

```{r, include = FALSE}
library(np)
mlb_bw0 <- npregbw(w_oba ~ bb_percent + k_percent + iso, 
                            data = train_mlb_data)
```

Generate predictions with `npreg` with provided bandwidth object

```{r}
mlb_test_npreg <- npreg(mlb_bw0, newdata = test_mlb_data) #<<
RMSE(mlb_test_npreg$mean, test_mlb_data$w_oba) 
```
