---
title: "Principal Component Analysis (Chapter 10.2 in ISLR)"
author: "Torben"
date: "May, 2019"
output: 
  html_document:
    theme: united
    toc: yes
    number_sections: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, comment = NA, warning = FALSE, fig.width = 12, fig.asp = 0.62)
```

```{r}
library(tidyverse)
theme_set(theme_bw())
```

# Principal Component Analysis

In high dimensional problems it may be hard to obtain insights from pairwise 
plots of the different variables in the data. There are ${p \choose 2} = p(p-1)/2$ combinations 
(when ignoring trivial flips of the axes). Hence, even for moderate $p$ this is an intractable
number of plots.. 

A solution is to rotate the data to obtain views of the data which contain as much information 
as possible. One criterion for "information" is *variablity*.

Recall that for two random variables $X$ and $Y$, we know that the linear combination $Z = aX+bY$ for scalars 
$a$ and $b$ has variance given by:
$$\mathbb{V}(Z) = \mathbb{V}(aX + bY) = a^2\mathbb{V}(X) + b^2\mathbb{V}(Y) + 2ab\mathbb{C}(X,Y),$$
where $\mathbb{V}(X)$ is the variance of $X$ and $\mathbb{C}(X,Y)$ the covariance of $X$ and $Y$.

We can choose $a$ and $b$ arbitrarily, but we would like to have the variance of $aX+bY$ maximal.
In order for this to be well-defined we must limit $a$ and $b$: $a^2 + b^2 = 1$, otherwise
we could make the variance arbitrarily large. The values of $a$ and $b$ defines the rotation 
of $X$ and $Y$ holding maximum variance.

This idea is the basis of PCA - rotate the data such that the resulting components have 
maximum variance. For the next components to differ from the previous components we impose
the requirement that they are pairwise __uncorrelated__ with the previous components.

## Illustrative example (in two dimensions)

The data shown below gives and example of how the *intrinsic* dimension of a data set may be smaller than the actual number of columns in the data.

![ISLR Fig 6.14](day05_fig/6_14.png)

We see that the first and second PC necessarily needs to be perpendicular to each other in this simple case to be uncorrelated. 
Below, the green line show PC1, whereas the dashed lines show the perpendicular distances from each observation to PC1.
From the definition of PC2, these distances (computed with sign) are therefore the PC2 coefficients.

![ISLR Fig 6.15](day05_fig/6_15.png)

In the plot below we see how PC1 *loads* on to the underlying variables. Because PC1 is given by $a_{1,1}$`Population` + $a_{1,2}$`Ad Spending` 
we expect this linear relationship between PC1 and the variables. However, for higher dimensional data the loads $a_{k, j}$ may be zero for 
various combinations of $k$ (PC) and $j$ (variable).

![ISLR Fig 6.16](day05_fig/6_16.png)

In the picture below we see that PC2 is less correlated to the data, which also suggest that PC1 captures most of the variability.

![ISLR Fig 6.17](day05_fig/6_17.png)

## PCA is found by eigenvalues 

The general problem of finding the PCA for $X$ being $n\times p$ data can be solved by finding the *eigenvalues*
of the covariance (or correlation) matrix. Let $\Sigma$ denote the covariance matrix of $X$. That is, $\mathbb(X) = \Sigma$.

For a linear transformation $X a$, where $a$ is a $p\times 1$ dimensional vector, the covariance is $a^\top\Sigma a$.
Hence, we want to maximize this variance with respect to $a$ while having $a^\top a = 1$ (as argued for above).
This corresponds to solving $$(\Sigma - \lambda I)a = 0,$$ where $\lambda$ is a Lagrangian multiplier.
The PCs are unique up to a flip of sign. 

For data where the variables are measured at different scales it is typically a good idea to scale the variables first.
This corresponds to using the correlation matrix instead.

## The best approximating low-dimensional linear surface

The PCA can also be viewed as the best linear surface approximating the data in lower 
dimension. 

![ISLR Fig 10.2](day05_fig/10_2a.png) ![ISLR Fig 10.2](day05_fig/10_2b.png)

# Examples

## PCA analysis of `USArrest` data

First we may make a pairs plot of the data
```{r, cache=TRUE}
pairs(USArrests)
GGally::ggpairs(USArrests)
```

We see that `Murder`, `Assault` and `Rape` are somewhat correlated and more scatter
when plotting `UrbanPop` against the three measures of crimes.

We clearly see that the observations have different scales. Hence, it is 
advisable to use the correlation matrix rather than the covariance matrix 
for the PCA rotation.

```{r}
head(USArrests)
summary(USArrests)
```

Hence, we use the `scale. = TRUE` in the call to standardise (that is use the 
correlation matrix) rather than the unit based covariances 

```{r}
USArrests_pca <- prcomp(USArrests, scale. = TRUE)
USArrests_pca
plot(USArrests_pca)
```

After the PCA rotation, we can yet again make a pairs plot
```{r, cache=TRUE}
GGally::ggpairs(USArrests_pca$x %>% as_tibble())
```

We notice the null-correlations of the PCs

Proportion of variance explained
```{r}
cumsum(USArrests_pca$sdev^2)/sum(USArrests_pca$sdev^2)
```

### Biplots

Showing both the data and variables in the same plots

```{r, fig.width = 12, fig.asp=1}
par(mfrow = c(2,2), mar = c(1,1,1,1)*4)
biplot(USArrests_pca, scale = 0, choices = c(1,2))
biplot(USArrests_pca, scale = 0, choices = c(1,3))
biplot(USArrests_pca, scale = 0, choices = c(2,3))

USArrests_pca_ <- prcomp(USArrests)
biplot(USArrests_pca_, scale = 0, choices = c(1,2))
```

Note, in ISLR they flip the sign of the PCs to obtain an easier interpretation of the 
PCs - however, this is merely for interpretation purposes.
```{r, eval = FALSE}
USArrests_pca$x <- -USArrests_pca$x
```

## Crabs

In the `MASS` package we have the data set `crabs``that contain measurements on 
a specific species of crabs (Leptograpsus).

![Crabs](day05_fig/crabs.png)

The data contains data from 50 males and 50 female crabs of each of two subspecies: 
blue and orange.

```{r}
ccol <- function(sp) ifelse(sp=="B","#0f7fa9","#fa8d0f")
cpch <- function(sx) 1 + 15*(crabs$sex=="M")

data(crabs, package = "MASS")
pairs(select(crabs,FL:BD), col = ccol(crabs$sp), pch = cpch(crabs$crabs))
```

Let us consider a PCA rotation

```{r}
(crab_pca <- prcomp(select(crabs, FL:BD), scale. = TRUE))
crabs_pca <- predict(crab_pca) %>% as_tibble() ## Rotate the data
```

Plot the results
```{r, fig.width=11, fig.asp=1}
par(mfrow = c(2,2), mar = rep(4,4))
plot(crab_pca, main = "Screeplot for correlations matrix")
plot(PC1~PC2,data=crabs_pca, col=ccol(crabs$sp), pch=cpch(crabs$sex))
plot(PC1~PC2,data=crabs_pca, col=ccol(crabs$sp), pch=cpch(crabs$sex))
plot(PC2~PC3,data=crabs_pca, col=ccol(crabs$sp), pch=cpch(crabs$sex))
```
 
## $t$-SNE - an alternative to PCA for visualisation

One of the more recent alternatives to PCA is the so called $t$-SNE method. In R this is implemented in `Rtsne` package. Different to PCA, the $t$-SNE method works by finding a 2- or 3-dimensional (highly) non-linear mapping of higher dimensional distances. There are several tuning parameters one can take into account. See Google Brain's ["How to Use t-SNE Effectively"](https://distill.pub/2016/misread-tsne/) for a great introduction.

### PCA vs. t-SNE

Although both PCA and t-SNE have their own advantages and disadvantages, some key differences between PCA and t-SNE can be noted as follows:

* $t$-SNE is computationally expensive and can take several hours on million-sample data sets where PCA will finish in seconds or minutes.
* PCA it is a mathematical technique, but t-SNE is a probabilistic one.
*    * Linear dimensionality reduction algorithms, like PCA, concentrate on placing dissimilar data points far apart in a lower dimension representation. But in order to represent high dimension data on low dimension, non-linear manifold, it is essential that similar data points must be represented close together, which is something t-SNE does not PCA.
* Sometimes in $t$-SNE different runs with the same hyper-parameters may produce different results hence multiple plots must be observed before making any assessment with $t$-SNE, while this is not the case with PCA.
* Since PCA is a linear algorithm, it will not be able to interpret the complex polynomial relationship between features while $t$-SNE is made to capture exactly that.

