---
title: "Supervised Learning"
subtitle: "Linear regression"
date: "June 21st, 2021"
---

```{r setup, echo = FALSE, warning = FALSE, message = FALSE}
library(tidyverse)
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")
```

## Simple linear regression

We assume a __linear relationship__ for $Y = f(X)$:

$$Y_{i}=\beta_{0}+\beta_{1} X_{i}+\epsilon_{i}, \quad \text { for } i=1,2, \dots, n$$

- $Y_i$ is the $i$th value for the __response__ variable
  
- $X_i$ is the $i$th value for the __predictor__ variable

--
  
- $\beta_0$ is an _unknown_, constant __intercept__: average value for $Y$ if $X = 0$
  
- $\beta_1$ is an _unknown_, constant __slope__: increase in average value for $Y$ for each one-unit increase in $X$
  
--
  
- $\epsilon_i$ is the __random__ noise: assume __independent, identically distributed__ (_iid_) from Normal distribution

$$\epsilon_i \overset{iid}{\sim}N(0, \sigma^2) \quad \text{ with constant variance } \sigma^2$$

---

## Simple linear regression estimation

We are estimating the __conditional expection (mean)__ for $Y$:

$$\mathbb{E}[Y_i| X_i] = \beta_0 + \beta_1X_i$$

- average value for $Y$ given the value for $X$

--

How do we estimate the __best fitting__ line?

--

__Ordinary least squares (OLS)__ - by minimizing the __residual sum of squares (RSS)__

$$R S S\left(\beta_{0}, \beta_{1}\right)=\sum_{i=1}^{n}\left[Y_{i}-\left(\beta_{0}+\beta_{1} X_{i}\right)\right]^{2}=\sum_{i=1}^{n}\left(Y_{i}-\beta_{0}-\beta_{1} X_{i}\right)^{2}$$

--

$$\widehat{\beta}_{1}=\frac{\sum_{i=1}^{n}\left(X_{i}-\bar{X}\right)\left(Y_{i}-\bar{Y}\right)}{\sum_{i=1}^{n}\left(X_{i}-\bar{X}\right)^{2}} \quad \text{ and } \quad \widehat{\beta}_{0}=\bar{Y}-\widehat{\beta}_{1} \bar{X}$$

where $\bar{X} = \frac{1}{n}\sum_{i=1}^n X_i$ and $\bar{Y} = \frac{1}{n}\sum_{i=1}^n Y_i$

---

## Connection to covariance and correlation

[__Covariance__](https://en.wikipedia.org/wiki/Covariance) describes the __joint variability of two variables__

$$\text{Cov}(X, Y) = \sigma_{X,Y} = \mathbb{E}[(X-\mathbb{E}[X])(Y-\mathbb{E}[Y])]$$

--

We compute the __sample covariance__ (use $n - 1$ since we are using the means and want [__unbiased estimates__](https://lazyprogrammer.me/covariance-matrix-divide-by-n-or-n-1/))

$$\hat{\sigma}_{X,Y} = \frac{1}{n-1} \sum_{i=1}^{n}\left(X_{i}-\bar{X}\right)\left(Y_{i}-\bar{Y}\right)$$

--

__Correlation__ is a _normalized_ form of covariance, ranges from -1 to 1

$$\rho_{X,Y} = \frac{\text{Cov}(X,Y)}{\sigma_X \cdot \sigma_Y}$$

__Sample correlation__ uses the sample covariance and standard deviations, e.g. $s_X^2 = \frac{1}{n-1} \sum_i (X_i - \bar{X})^2$

$$r_{X,Y} = \frac{\sum_{i=1}^{n}\left(X_{i}-\bar{X}\right)\left(Y_{i}-\bar{Y}\right)}{\sqrt{\sum_{i=1}^{n}\left(X_{i}-\bar{X}\right)^{2} \sum_{i=1}^{n}\left(Y_{i}-\bar{Y}\right)^{2}}}$$

---

## Connection to covariance and correlation

So we have the following:

$$\widehat{\beta}_{1}=\frac{\sum_{i=1}^{n}\left(X_{i}-\bar{X}\right)\left(Y_{i}-\bar{Y}\right)}{\sum_{i=1}^{n}\left(X_{i}-\bar{X}\right)^{2}} \quad \text{ compared to } \quad r_{X,Y} = \frac{\sum_{i=1}^{n}\left(X_{i}-\bar{X}\right)\left(Y_{i}-\bar{Y}\right)}{\sqrt{\sum_{i=1}^{n}\left(X_{i}-\bar{X}\right)^{2} \sum_{i=1}^{n}\left(Y_{i}-\bar{Y}\right)^{2}}}$$

--

$\Rightarrow$ Can rewrite $\hat{\beta}_1$ as:

$$\widehat{\beta}_{1} = r_{X,Y} \cdot \frac{s_Y}{s_X}$$

$\Rightarrow$ Can rewrite $r_{X,Y}$ as:

$$r_{X,Y} = \widehat{\beta}_{1} \cdot \frac{s_X}{s_Y}$$

--

Can think of $\widehat{\beta}_{1}$ weighting the ratio of variance between $X$ and $Y$...

---

## Example data: NFL teams summary

Created dataset using [`nflfastR`](https://www.nflfastr.com/) summarizing NFL team performances from 1999 to 2020

```{r, warning = FALSE, message = FALSE}
library(tidyverse)
nfl_teams_data <- read_csv("http://www.stat.cmu.edu/cmsac/sure/2021/materials/data/regression_projects/nfl_team_season_summary.csv")
nfl_teams_data
```

---

## Modeling NFL score differential

.pull-left[

Interested in modeling a team's __score differential__

```{r score-diff, echo = TRUE, eval = FALSE}
nfl_teams_data <- nfl_teams_data %>%
  mutate(score_diff = 
           points_scored - points_allowed)
nfl_teams_data %>%
  ggplot(aes(x = score_diff)) +
  geom_histogram(color = "black", 
                 fill = "darkblue",
                 alpha = 0.3) +
  theme_bw() +
  labs(x = "Score differential")
```



]
.pull-right[

```{r ref.label = 'score-diff', echo = FALSE}

```


]

---

### Relationship between score differential and passing performance

.pull-left[
- `offense_ave_epa_pass`: team's average [expected points added (EPA)]() per pass attempt

```{r epa-score-diff-plot, eval = FALSE}
nfl_plot <- nfl_teams_data %>%
  ggplot(aes(x = offense_ave_epa_pass, #<<
             y = score_diff)) +
  geom_point(alpha = 0.5) +
  theme_bw() +
  labs(x = "EPA gained per pass attempt",
       y = "Score differential")
nfl_plot
```

We fit linear regression models using `lm()`, formula is input as: `response ~ predictor`

```{r score-diff-reg}
init_lm <- lm(score_diff ~ offense_ave_epa_pass, #<<
              data = nfl_teams_data) #<<
```


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

```
]

---

## View the model `summary()`

```{r summary-score-diff-reg}
summary(init_lm)
```

---

## Inference with OLS

Reports the intercept and coefficient estimates: $\quad \hat{\beta}_0 \approx -4.89 \quad, \quad \hat{\beta}_1 \approx 566.352$ 

--

Estimates of uncertainty for $\beta$s via __standard errors__: $\quad \widehat{SE}(\hat{\beta}_0) \approx 2.561 \quad, \quad \widehat{SE}(\hat{\beta}_1) \approx 19.226$ 

--

$t$-statistics are coefficients `Estimates` / `Std. Error`, i.e., number of standard deviations from 0

  - _p-values_ (i.e., `Pr(>|t|)`): estimated probability observing value as extreme as |`t value`| __given the null hypothesis__ $\beta = 0$
  
  - p-value $<$ conventional threshold of $\alpha = 0.05$, __sufficient evidence to reject the null hypothesis that the coefficient is zero__,
  
  - Typically |`t values`| $> 2$ indicate __significant__ relationship at $\alpha = 0.05$

  - i.e., there is a __significant__ association between `offense_ave_epa_pass` and `score_diff`


---

## Be careful!

__Caveats to keep in mind regarding p-values:__

--

If the true value of a coefficient $\beta = 0$, then the p-value is sampled from a [Uniform(0,1) distribution](https://en.wikipedia.org/wiki/Continuous_uniform_distribution)

- i.e. it is just as likely to have value 0.45 as 0.16 or 0.84 or 0.9999 or 0.00001...

--

$\Rightarrow$ Hence why we typically only reject for low $\alpha$ values like 0.05

  - Controlling the Type 1 error rate at $\alpha = 0.05$, i.e., the probability of a __false positive__ mistake
  
  - 5% chance that you'll conclude there's a significant association between $x$ and $y$ __even when there is none__

--

Remember what a standard error is? $SE = \frac{\sigma}{\sqrt{n}}$

  - $\Rightarrow$ As $n$ gets large __standard error goes to zero__, and *all* predictors are eventually deemed significant

- While the p-values might be informative, we will explore other approaches to determine which subset of predictors to include (e.g., holdout performance)

---

## Back to the model summary: `Multiple R-squared`

Back to the connection between the coefficient and correlation:

$$r_{X,Y} = \widehat{\beta}_{1} \cdot \frac{s_X}{s_Y} \quad \Rightarrow \quad r^2_{X,Y} = \widehat{\beta}_{1}^2\cdot \frac{s_X^2}{s_Y^2}$$

Compute the correlation with `cor()`:

```{r example-cor}
with(nfl_teams_data, cor(offense_ave_epa_pass, score_diff)) #<<
```

--

The squared `cor` matches the reported `Multiple R-squared`

```{r example-cor-2}
with(nfl_teams_data, cor(offense_ave_epa_pass, score_diff))^2 #<<
```


---

## Back to the model summary: `Multiple R-squared`

Back to the connection between the coefficient and correlation:

$$r_{X,Y} = \widehat{\beta}_{1} \cdot \frac{s_X}{s_Y} \quad \Rightarrow \quad r^2_{X,Y} = \widehat{\beta}_{1}^2\cdot \frac{s_X^2}{s_Y^2}$$

$r^2$ (or also $R^2$) estimates the __proportion of the variance__ of $Y$ explained by $X$ 

  - More generally: variance of model predictions / variance of $Y$

```{r general-rsquared}
var(predict(init_lm)) / var(nfl_teams_data$score_diff) #<<
```


---

## Generating predictions

We can use the `predict()` function to either get the fitted values of the regression:

```{r init-predict}
train_preds <- predict(init_lm)
head(train_preds)
```

Which is equivalent to using:
```{r init-fitted}
head(init_lm$fitted.values)
```


---

## Predictions for new data

Or we can provide it `newdata` which __must contain the explanatory variables__:

.pull-left[
```{r steelers-preds, eval = FALSE}
steelers_data <- nfl_teams_data %>% 
  filter(team == "PIT", season == 2020)

new_steelers_data <- steelers_data %>%
  dplyr::select(team, offense_ave_epa_pass) %>%
  slice(rep(1, 3)) %>% #<<
  mutate(adj_factor = c(0.5, 1, 2),
         offense_ave_epa_pass = offense_ave_epa_pass * adj_factor)
new_steelers_data$pred_score_diff <- 
  predict(init_lm, newdata = new_steelers_data) #<<
nfl_plot +
  geom_point(data = new_steelers_data,
             aes(x = offense_ave_epa_pass,
                 y = pred_score_diff),
             color = "gold", size = 5)
```
]
.pull-right[
```{r ref.label='steelers-preds', echo = FALSE, fig.height=6,fig.width=6,fig.align="center"}

```


]

---

### Plot observed values against predictions

Useful diagnostic (for __any type of model__, not just linear regression!) 

.pull-left[

```{r plot-pred-obs, eval = FALSE}
nfl_teams_data %>%
  mutate(pred_vals = predict(init_lm)) %>% #<<
  ggplot(aes(x = pred_vals,
             y = score_diff)) +
  geom_point(alpha = 0.5) +
  geom_abline(slope = 1, intercept = 0, #<<
              linetype = "dashed",
              color = "red",
              size = 2) +
  theme_bw()
```

- "Perfect" model will follow __diagonal__

]
.pull-right[
```{r ref.label='plot-pred-obs', echo = FALSE, fig.height=6,fig.width=6,fig.align="center"}

```


]

---

### Plot observed values against predictions

Can augment the data with model output using the [`broom` package](https://cran.r-project.org/web/packages/broom/vignettes/broom.html)

.pull-left[

```{r plot-pred-obs-broom, eval = FALSE}
nfl_teams_data <- 
  broom::augment(init_lm, nfl_teams_data) #<<
nfl_teams_data %>%
  ggplot(aes(x = .fitted, 
             y = score_diff)) + #<<
  geom_point(alpha = 0.5) +
  geom_abline(slope = 1, intercept = 0, #<<
              linetype = "dashed",
              color = "red",
              size = 2) +
  theme_bw()
```

- Adds various columns from model fit we can use in plotting for model diagnostics

]
.pull-right[
```{r ref.label='plot-pred-obs-broom', echo = FALSE, fig.height=6,fig.width=6,fig.align="center"}

```


]

---

## Plot residuals against predicted values

.pull-left[

- Residuals = observed - predicted

- Conditional on the predicted values, the __residuals should have a mean of zero__

```{r plot-resid-fit-broom, eval = FALSE}
nfl_teams_data %>%
  ggplot(aes(x = .fitted, 
             y = .resid)) + #<<
  geom_point(alpha = 0.5) +
  geom_hline(yintercept = 0, #<<
             linetype = "dashed",
             color = "red",
             size = 2) +
  # To plot the residual mean
  geom_smooth(se = FALSE) + #<<
  theme_bw()
```

- Residuals __should NOT display any pattern__

]
.pull-right[
```{r ref.label='plot-resid-fit-broom', echo = FALSE, fig.height=6,fig.width=6,fig.align="center"}

```


]

---

## Multiple regression 

We can include as many variables as we want (assuming $n > p$!)

$$Y=\beta_{0}+\beta_{1} X_{1}+\beta_{2} X_{2}+\cdots+\beta_{p} X_{p}+\epsilon$$

OLS estimates in matrix notation ( $\boldsymbol{X}$ is a $n \times p$ matrix):
  
$$\hat{\boldsymbol{\beta}} = (\boldsymbol{X} ^T \boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{Y}$$

--

Can just add more variables to the formula in `R`

```{r}
multiple_lm <- lm(score_diff ~ offense_ave_epa_pass + defense_ave_epa_pass, #<<
                  data = nfl_teams_data)
```

- Use the `Adjusted R-squared` when including multiple variables $= 1 - \frac{(1 - R^2)(n - 1)}{(n - p - 1)}$

  - Adjusts for the number of variables in the model $p$  
  
  - Adding more variables __will always increase__ `Multiple R-squared`

---

### What about the Normal distribution assumption???

$$Y=\beta_{0}+\beta_{1} X_{1}+\beta_{2} X_{2}+\cdots+\beta_{p} X_{p}+\epsilon$$

- $\epsilon_i$ is the __random__ noise: assume __independent, identically distributed__ (_iid_) from Normal distribution

$$\epsilon_i \overset{iid}{\sim}N(0, \sigma^2) \quad \text{ with constant variance } \sigma^2$$

--

__OLS doesn't care about this assumption__, it's just estimating coefficients!

--

In order to perform inference, __we need to impose additional assumptions__

By assuming $\epsilon_i \overset{iid}{\sim}N(0, \sigma^2)$, what we really mean is:

$$Y \overset{iid}{\sim}N(\beta_{0}+\beta_{1} X_{1}+\beta_{2} X_{2}+\cdots+\beta_{p} X_{p}, \sigma^2)$$

--

So we're estimating the mean $\mu$ of this conditional distribution, but what about $\sigma^2$?
  
--

[Unbiased estimate](https://bradleyboehmke.github.io/HOML/linear-regression.html#simple-linear-regression) $\hat{\sigma}^2 = \frac{RSS}{n - (p + 1)}$, its square root is the `Residual standard error` 

  - __Degrees of freedom__:  $n - (p + 1)$, data supplies us with $n$ "degrees of freedom" and we used up $p + 1$

---

### Check the assumptions about normality with [`ggfortify`](https://cran.r-project.org/web/packages/ggfortify/vignettes/plot_lm.html)

```{r more-resid-plots, fig.height=5,fig.width=15,fig.align="center"}
library(ggfortify)
autoplot(multiple_lm, ncol = 4) + theme_bw() #<<
```

- Standardized residuals = residuals `/ sd(`residuals`)` (see also `.std.resid` from `augment`)


