---
title: "Dimension Reduction"
subtitle: "Principal components analysis (PCA)"
date: "July 6th, 2022"
---

```{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")
```

## What is the goal of dimension reduction?

We have $p$ variables (columns) for $n$ observations (rows) __BUT__ which variables are __interesting__?

--

Can we find a smaller number of dimensions that captures the __interesting__ structure in the data?

  - Could examine all pairwise scatterplots of each variable - tedious, manual process

  - Tuesday: clustered variables based on correlation
  
--

  - Can we find a combination of the original $p$ variables?


--

__Dimension reduction__: 

- Focus on reducing the dimensionality of the feature space (i.e., number of columns), 

- While __retaining__ most of the information / __variability__ in a lower dimensional space (i.e., reducing the number of columns)

---

## [Principal components analysis (PCA)](https://en.wikipedia.org/wiki/Principal_component_analysis)

```{r out.width='110%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://bradleyboehmke.github.io/HOML/images/pcr-steps.png")
```


---

## [Principal components analysis (PCA)](https://en.wikipedia.org/wiki/Principal_component_analysis)

- PCA explores the __covariance__ between variables, and combines variables into a smaller set of __uncorrelated__ variables called __principal components (PCs)__

- PCs are __weighted__, linear combinations of the original variables

  - Weights reveal how different variables are ___loaded___ into the PCs

- We want a __small number of PCs__ to explain most of the information / variance in the data

--

__First principal component__:

$$Z_1 = \phi_{11} X_1 + \phi_{21} X_2 + \dots + \phi_{p1} X_p$$

--

  - $\phi_{j1}$ are the weights indicating the contributions of each variable $j \in 1, \dots, p$
  
  - Weights are normalized $\sum_{j=1}^p \phi_{j1}^2 = 1$
  
  - $\phi_{1} = (\phi_{11}, \phi_{21}, \dots, \phi_{p1})$ is the __loading vector__ for PC1

--
  
  - $Z_1$ is a linear combination of the $p$ variables that has the __largest variance__

---

## [Principal components analysis (PCA)](https://en.wikipedia.org/wiki/Principal_component_analysis)

__Second principal component__:

$$Z_2 = \phi_{12} X_1 + \phi_{22} X_2 + \dots + \phi_{p2} X_p$$

  - $\phi_{j2}$ are the weights indicating the contributions of each variable $j \in 1, \dots, p$
  
  - Weights are normalized $\sum_{j=1}^p \phi_{j1}^2 = 1$
  
  - $\phi_{2} = (\phi_{12}, \phi_{22}, \dots, \phi_{p2})$ is the __loading vector__ for PC2
  
  - $Z_2$ is a linear combination of the $p$ variables that has the __largest variance__
  
    - __Subject to constraint it is uncorrelated with $Z_1$__ 
    
--

We repeat this process to create $p$ principal components

---

## [Visualizing PCA](https://www.stevejburr.com/post/scatter-plots-and-best-fit-lines/) in two dimensions

```{r out.width='50%', echo = FALSE, fig.align='center'}
# Ensure reproducibility by setting random number seed
set.seed(123) 
plot_data <- tibble("x" = rnorm(50, mean = 100, sd = 20)) %>%
  mutate(y =  0.8 * x + rnorm(50, mean = 0, sd = 10))
basic_scatter <- ggplot(plot_data) +
  geom_point(aes(x, y), color = "black")+
  coord_equal() +
  theme_bw()
basic_scatter
```

---

## [Visualizing PCA](https://www.stevejburr.com/post/scatter-plots-and-best-fit-lines/) in two dimensions

```{r out.width='50%', echo = FALSE, fig.align='center'}
#fit the model
line1 <- lm(y ~ x, plot_data)$coef
#extract the slope from the fitted model
line1.slope <- line1[2]
#extract the intercept from the fitted model
line1.intercept <- line1[1]
basic_scatter_yfit <- basic_scatter +
  geom_abline(aes(slope = line1.slope, intercept = line1.intercept),
              colour = "darkred") +
  annotate("text", x = 75, y = 120, label = "y ~ x", color = "darkred",
           size = 9)
basic_scatter_yfit
```

---

## [Visualizing PCA](https://www.stevejburr.com/post/scatter-plots-and-best-fit-lines/) in two dimensions

```{r out.width='50%', echo = FALSE, fig.align='center'}
#fit the model
line2 <- lm(x ~ y, plot_data)$coef
#extract the slope from the fitted model
line2.slope <- 1 / line2[2]
#extract the intercept from the fitted model
line2.intercept <- -(line2[1] / line2[2])
basic_scatter_xyfit <- basic_scatter_yfit +
  geom_abline(aes(slope = line2.slope, intercept = line2.intercept),
              colour = "blue") +
  annotate("text", x = 125, y = 55, label = "x ~ y", color = "blue",
           size = 9)
basic_scatter_xyfit
```

---

## [Visualizing PCA](https://www.stevejburr.com/post/scatter-plots-and-best-fit-lines/) in two dimensions

```{r out.width='50%', echo = FALSE, fig.align='center'}
pca <- prcomp(cbind(plot_data$x, plot_data$y))$rotation
pca.slope <- pca[2,1] / pca[1,1]
pca.intercept <- mean(plot_data$y) - (pca.slope * mean(plot_data$x))

basic_scatter_xy_pca_fit <- basic_scatter_xyfit +
  geom_abline(aes(slope = pca.slope, intercept = pca.intercept),
              colour = "darkorange") +
  annotate("text", x = 75, y = 90, label = "PCA", color = "darkorange",
           size = 9)
basic_scatter_xy_pca_fit
```

---

## [Visualizing PCA](https://www.stevejburr.com/post/scatter-plots-and-best-fit-lines/) in two dimensions

```{r fig.width=14, fig.height=6, echo = FALSE, fig.align='center'}
plot_data %>%
  #calculate the positions using the line equations:
  mutate(yhat_line1=(x*line1.slope+line1.intercept),
         xhat_line1=x,
         yhat_line2=y,
         xhat_line2=(y-line2.intercept)/line2.slope,
         #https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line
         a=pca.slope,
         b=-1,
         c=pca.intercept,
         xhat_line3=(b*(b*x-a*y)-(a*c))/((a*a)+(b*b)),
         yhat_line3=(a*(-b*x+a*y)-(b*c))/((a*a)+(b*b)),
         #add the slopes/intercepts to this data frame:
         slope_line1=line1.slope,
         slope_line2=line2.slope,
         slope_line3=pca.slope,
         intercept_line1=line1.intercept,
         intercept_line2=line2.intercept,
         intercept_line3=pca.intercept
         )%>% 
  #drop intermediate variables
  select(-c(a,b,c)) %>%
  pivot_longer(yhat_line1:intercept_line3, names_to = "key", values_to = "value") %>%
  #transpose to a long form
  #gather(key="key",value="value",-c(x,y)) %>% 
  # have "yhat_line1", want two colums of "yhat" "line1"
  separate(key,c("type", "line"), "_") %>% 
  #then transpose to be fatter, so we have cols for xhat, yhat etc
  pivot_wider(names_from = "type", values_from = "value") %>%
  #spread(key="type",value="value") %>%
  #relable the lines with more description names, and order the factor for plotting:
  mutate(line=case_when(
           line=="line1" ~ "y ~ x",
           line=="line2" ~ "x ~ y",
           TRUE ~ "PCA"
         ),
         line = fct_relevel(line, "y ~ x", "x ~ y", "PCA")) %>% 
  ggplot() +
  geom_point(aes(x = x, y = y, color = line)) +
  geom_abline(aes(slope = slope, intercept = intercept, color = line)) +
  geom_segment(aes(x = x, y = y, xend = xhat, yend = yhat, color = line)) +
  facet_wrap(~ line, ncol = 3) +
  scale_color_manual(values = c("darkred", "blue", "darkorange")) +
  theme_bw() +
  theme(strip.background = element_blank(),
        legend.position = "none",
        strip.text = element_text(size = 16))
```

---

## Searching for variance in orthogonal directions

```{r out.width='60%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://bradleyboehmke.github.io/HOML/15-pca_files/figure-html/create-pca-image-1.png")
```


---


## PCA: [__singular value decomposition (SVD)__](https://en.wikipedia.org/wiki/Singular_value_decomposition)

$$
X = U D V^T
$$

- Matrices $U$ and $V$ contain the left and right (respectively) __singular vectors of scaled matrix $X$__

- $D$ is the diagonal matrix of the __singular values__

--

- SVD simplifies matrix-vector multiplication as __rotate, scale, and rotate again__


--
$V$ is called the __loading matrix__ for $X$ with $\phi_{j}$ as columns, 

  - $Z = X  V$ is the PC matrix

--

BONUS __eigenvalue decomposition__ (aka spectral decomposition)

- $V$ are the __eigenvectors__ of $X^TX$ (covariance matrix, $^T$ means _transpose_)

- $U$ are the __eigenvectors__ of $XX^T$

- The singular values (diagonal of $D$) are square roots of the __eigenvalues__ of $X^TX$ or $XX^T$

- Meaning that $Z = UD$

---

## Eigenvalues solve time travel?


```{r out.width='70%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://thumbs.gfycat.com/RealisticFragrantHerculesbeetle-size_restricted.gif")
```


---

## Probably not... but they guide dimension reduction

We want to choose $p^* < p$ such that we are explaining variation in the data

--

Eigenvalues $\lambda_j$ for $j \in 1, \dots, p$ indicate __the variance explained by each component__

  - $\sum_j^p \lambda_j = p$, meaning $\lambda_j \geq 1$ indicates $\text{PC}j$ contains at least one variable's worth in variability
  
  - $\lambda_j / p$ equals proportion of variance explained by $\text{PC}j$
  
  - Arranged in descending order so that $\lambda_1$ is largest eigenvalue and corresponds to PC1
  
--
  
  - Can compute the cumulative proportion of variance explained (CVE) with $p^*$ components:
  
$$\text{CVE}_{p^*} = \frac{\sum_j^{p*} \lambda_j}{p}$$

Can use [__scree plot__](https://en.wikipedia.org/wiki/Scree_plot) to plot eigenvalues and guide choice for $p^* <p$ by looking for "elbow" (rapid to slow change)

---

## Example data: NFL teams summary

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

```{r, warning = FALSE, message = FALSE}
library(tidyverse)
nfl_teams_data <- read_csv("http://www.stat.cmu.edu/cmsac/sure/2022/materials/data/sports/regression_examples/nfl_team_season_summary.csv")
nfl_model_data <- nfl_teams_data %>%
  mutate(score_diff = points_scored - points_allowed) %>%
  # Only use rows with air yards
  filter(season >= 2006) %>%
  dplyr::select(-wins, -losses, -ties, -points_scored, -points_allowed, -season, -team)
```

---

## NFL PCA example

Use the `prcomp` function (uses SVD) for PCA on __centered__ and __scaled__ data

```{r}
model_x <- as.matrix(dplyr::select(nfl_model_data, -score_diff))
pca_nfl <- prcomp(model_x, center = TRUE, scale = TRUE) #<<
summary(pca_nfl)
```


---

## Proportion of variance explained

`prcomp$sdev` corresponds to the singular values, i.e., $\sqrt{\lambda_j}$, what is `pca_nfl$sdev^2 / ncol(model_x)`?

--

Can use the `broom` package easily tidy `prcomp` summary for plotting

.pull-left[
```{r plot-pve, eval = FALSE}
library(broom)
pca_nfl %>%
  tidy(matrix = "eigenvalues") %>% #<< 
  ggplot(aes(x = PC, y = percent)) +
  geom_line() + geom_point() +
  geom_hline(yintercept = 1 / ncol(model_x), #<<
             color = "darkred", 
             linetype = "dashed") +
  theme_bw()
```

- Add reference line at $1/p$, _why_?

]
.pull-right[
```{r ref.label='plot-pve', out.width='80%', echo = FALSE, fig.align='center'}

```

]

---

## Display data in lower dimensions

`prcomp$x` corresponds to the matrix of __principal component scores__, i.e., $Z = XV$

.pull-left[

Can `augment` dataset with PC scores for plotting

  - Add team and season for context

```{r pc-scores, eval = FALSE}
pca_nfl %>%
  augment(nfl_model_data) %>% #<<
  bind_cols({
    nfl_teams_data %>% 
      filter(season >= 2006) %>%
      dplyr::select(season, team)
  }) %>%
  unite("team_id", team:season, sep = "-", #<<
        remove = FALSE) %>%
  ggplot(aes(x = .fittedPC1, y = .fittedPC2, 
             color = season)) +
  geom_text(aes(label = team_id), alpha = 0.9) +
  scale_color_gradient(low = "purple", high = "green") +
  theme_bw() + theme(legend.position = "bottom")
```


]
.pull-right[

```{r, ref.label='pc-scores',  out.width='90%', echo = FALSE, fig.align='center'}

```

]

---

## What are the [loadings of these dimensions](https://clauswilke.com/blog/2020/09/07/pca-tidyverse-style/)?

`prcomp$rotation` corresponds to the __loading matrix__, i.e., $V$

.pull-left[

```{r pc-loadings, eval = FALSE}
arrow_style <- arrow(
  angle = 20, ends = "first", type = "closed", 
  length = grid::unit(8, "pt")
)
library(ggrepel)
pca_nfl %>%
  tidy(matrix = "rotation") %>%
  pivot_wider(names_from = "PC", names_prefix = "PC", 
              values_from = "value") %>%
  mutate(stat_type = ifelse(str_detect(column, "offense"),
                            "offense", "defense")) %>%
  ggplot(aes(PC1, PC2)) +
  geom_segment(xend = 0, yend = 0, arrow = arrow_style) +
  geom_text_repel(aes(label = column, color = stat_type),
                  size = 3) +
  scale_color_manual(values = c("darkred", "darkblue")) +
  theme_bw() +
  theme(legend.position = "bottom")

```


]
.pull-right[

```{r, ref.label='pc-loadings',  out.width='90%', echo = FALSE, fig.align='center'}

```

]

---

## PCA analysis with `factoextra`

Visualize the proportion of variance explained by each PC with [`factoextra`](http://www.sthda.com/english/wiki/factoextra-r-package-easy-multivariate-data-analyses-and-elegant-visualization)

```{r, fig.align='center', fig.height=5}
library(factoextra)
fviz_eig(pca_nfl)
```

---

## PCA analysis with `factoextra`

Display observations with first two PC

```{r, fig.align='center', fig.height=5}
fviz_pca_ind(pca_nfl)
```


---

## PCA analysis with `factoextra`

Projection of variables - angles are interpreted as correlations, where negative correlated values point to opposite sides of graph

```{r, fig.align='center', fig.height=5}
fviz_pca_var(pca_nfl)
```

---

## PCA analysis with `factoextra`

__Biplot__ displays both the space of observations and the space of variables 

  - Arrows represent the directions of the original variables

```{r, fig.align='center', fig.height=5}
fviz_pca_biplot(pca_nfl)
```


