## Principal Component Analysis 

Today's tutorial is adapted from "An Introduction to Statistical Learning with applications in R" by James, Witten, Hastie, and Tibshirani. See class slides for an overview of PCA.

We are going to use familiar data that we have explored with other unsupervised methods: the USArrests data. 

```{r}
library(tidyverse)
library(gridExtra)

data("USArrests")
head(USArrests, 10)
```

First, we compute the variance of each variable. We do this to inspect the scales of each variable and decide whether to standardize the data.

```{r}
apply(USArrests, 2, var)

#we see a lot of variance! So let's standardize
scaled_df = apply(USArrests, 2, scale)
head(scaled_df)
```

Recall that the goal of PCA is to capture variability in the data using smaller set of variables. We might utilize PCA when:

- we have a high number of variables, and investigating pairwise scatterplots across 400 variables (for example) is excessive, and costly.
- we need lower dimensional representation of the data. This is particularly true when creating predictive pipelines or reducing computational costs. 

To begin, we will calculate the components without using a PCA package in R. Below, we first calculate the covariance matrix of the data.  Next, we calculate eigenvalues & eigenvectors. There is a bunch of matrix alegbra under the hood of PCA!

```{r}
arrests.cov = cov(scaled_df)
arrests.eigen = eigen(arrests.cov)
str(arrests.eigen)

# extract loadings
phi = arrests.eigen$vectors[,1:2]
phi = -phi
row.names(phi) = c("Murder", "Assault", "UrbanPop", "Rape")
colnames(phi) = c("PC1", "PC2")
print(phi)
```

Recall from lecture: 

- each component defines a direction in the feature space
- because eigenvectors are orthogonal to each other, these are uncorrelated 
- the table (e.g., `print phi`) shows us that PC1 is picking up a lot on crime related measures while PC2 picks up on population.

Next we can project the $n$ data points ($X_1$ to $X_n$) onto the first eigenvector, these projected values are called the principal component scores for each observation.

```{r}
# calculate PC scores (%*% is for matrix multiplication)
PC1 = as.matrix(scaled_df) %*% phi[,1]
PC2 = as.matrix(scaled_df) %*% phi[,2] 

# create a data frame with PC scores & state names
PC = data.frame(state= row.names(USArrests), PC1, PC2)
head(PC)
```

Now that we’ve calculated the first and second principal components for each state, we can plot them against each other and create a two-dimensional view of the data. 

As we saw earlier from our tabular results, the first principal component (x-axis) roughly corresponds to the rate of serious crimes. States such as California and Florida have high rates of serious crimes, while places like North Dakota and Vermont have far lower rates. The second component (y-axis) is roughly capturing as urbanization. States that are close to the origin, such as Indiana or Virginia, are close to the average in both categories.

```{r}
# create a two dimensional plot
ggplot(PC, aes(PC1, PC2)) + 
    modelr::geom_ref_line(h=0) + 
    modelr::geom_ref_line(v=0) + 
    geom_text(aes(label=state), size=3) + 
    xlab("First Principal Component") + ylab("Second Principal Component") + 
    ggtitle("First Two Principal Components of USArrests Data")
```


Here we have only looked at 2 of the 4 PCs, how do we know how much variance these capture?  Similar to our analysis using k-means, we want to look at the porportion of variance explained by each PC:

```{r}
PVE <- arrests.eigen$values / sum(arrests.eigen$values)
round(PVE, 2)
```

First PC therefore captures 62% of the variability, the second captures about 25%. We might also want to plot the PVE and cumulative PVE, this is called a scree plot.

```{r}
PVEplot = qplot(c(1:4), PVE) + 
    geom_line() + 
    xlab("Principal Component") + ylab('PVE') + 
    ggtitle("Scree Plot") + ylim(0,1) + 
    theme_minimal()

cumulativePVE <- qplot(c(1:4), cumsum(PVE)) + 
  geom_line() + 
  xlab("Principal Component") + ylab(NULL) + 
  ggtitle("Cumulative Scree Plot") + ylim(0,1) +
  theme_minimal()

grid.arrange(PVEplot, cumulativePVE, ncol = 2)
```

Since PCA is all about reduction we want the best trade off, i.e. the lowest number of components that capture the most variance. Thus given the plots above you might choose 3 or 4, for an interpretable number that also captures a high level of variance. 

# R has packages with built-in functions for the entire thing! 

## With your partner [Practice]

Below, explore the same analysis using the stats packages' `prcomp`.

```{r}
pca_result <- prcomp(USArrests, scale = TRUE)
pca_result$scale
```

Examine the output of this function. The center and scale components correspond to the means and standard deviations of the variables that were used for scaling prior to implementing PCA. The rotation matrix provides the principal component loadings; each column of `pca_result$rotation` contains the corresponding principal component loading vector.

Also, notice that PCA1 and PCA2 are opposite signs from what we computated earlier. Recall that by default, eigenvectors in R point in the negative direction. We can adjust this with a simple change. Finally, we can check out the results using a biplot.

```{r}
pca_result$rotation <- -pca_result$rotation
pca_result$rotation
pca_result$x <- - pca_result$x
biplot(pca_result, scale = 0)
```

Biplot discussion: The biplot shows the first two principal components for the USArrests data. The state names represent the scores for the first two principal components. The arrows indicate the first two principal component loading vectors. The loading for Rape on the first component is .54 and its loading on the second principal component is 9.17 (the word rape is centered at the point 0.54, 0.17). This figure is known as a biplot because it displays both the principal component scores and the principal component loadings. 

# On your own

Repeat the analysis using unscaled data and create a biplot. Compare this to your analysis of scaled data. Which variable, using unscaled data, has the highest loading? Why?

Explore `decathalon` data using PCA. 

```{r}
library("factoextra")
data(decathlon2)
```

Steps:

1. subsetting
2. standardizing 
3. Calculate PCA (can use factoextra function or without this package)
4. Inspect proportion of variance explained (fviz_eig or without this package)
5. For fun: it can also be interesting to explore 'variable contributions'. The function `fviz_contrib()` [factoextra package] can be used to draw a bar plot of variable contributions.
