---
title: "Principal component analysis (PCA)"
author: "Trevor Hastie and Robert Tibshirani"
output:
  pdf_document: default
---

# Unsupervised Learning

**Here, I am adapting part of the lab associated with Chapter 12 of the textbook.** 

## Principal Components Analysis

In this lab, we perform PCA on the `USArrests` data set, which is part of the base `R` package.
The rows of the data set contain the 50 states, in alphabetical order.

```{r chunk1}
states <- row.names(USArrests)
states
```

The columns of the data set contain the four variables.

```{r chunk2}
names(USArrests)
```

We first briefly examine the data. We notice that the variables have vastly different means.

```{r chunk3}
apply(USArrests, 2, mean)
```

Note that the `apply()` function allows us to apply a function---in this case, the `mean()` function---to each row or column of the data set. The second input here denotes whether we wish to compute the mean of the rows, $1$, or the columns, $2$. We see that there are on average three times as many rapes as murders, and more than eight times as many assaults as rapes.
We can also examine the variances of the four variables using the `apply()` function.

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

Not surprisingly, the variables also have vastly different variances:
 the `UrbanPop` variable measures the percentage of the population in each state living in an urban area, which is not a comparable number to the number of rapes
in each state per 100,000 individuals.
If we failed to scale the variables before performing PCA, then most of the principal components that we observed would be driven by the `Assault` variable, since it has by far the largest mean and variance.
Thus, it is important to standardize the variables to have mean zero and standard deviation one before performing PCA.

We now perform principal components analysis using the `prcomp()` function, which is one of several functions in `R` that perform PCA.

```{r chunk5}
pr.out <- prcomp(USArrests, scale = TRUE)
pr.out
```

By default, the `prcomp()` function centers the variables to have
mean zero. By using the option `scale = TRUE`, we scale the
variables to have standard deviation one. The output from
`prcomp()` contains a number of useful quantities.

```{r chunk6}
names(pr.out)
```

The `center` and `scale` components correspond to the means and standard deviations of the variables that were used for scaling prior to implementing PCA.

```{r chunk7}
pr.out$center
pr.out$scale
```

The `rotation` matrix provides the principal component loadings;
each column of `pr.out$rotation` contains the corresponding
principal component loading vector.

**This function names it
  the rotation matrix, because when we matrix-multiply the $\bf X$
  matrix by `pr.out$rotation`, it gives us the coordinates of the
  data in the rotated coordinate system. These coordinates are the
  principal component scores.**


```{r chunk8}
pr.out$rotation
```

We see that there are four distinct principal components. This is to
be expected because there are in general $\min(n-1,p)$ informative
principal components in a data set with $n$ observations and $p$
variables.

Using the `prcomp()` function, we do not need to explicitly multiply the data by the principal component loading vectors  in order to obtain the principal component score vectors. Rather the $50 \times 4$ matrix `x` has as its columns the principal component score vectors. That is, the $k$th column is the $k$th principal component score vector.

```{r chunk9}
dim(pr.out$x)
```

We can plot the first two principal components as follows:

```{r chunk10}
biplot(pr.out, scale = 0, cex=0.5)
```

The `scale = 0` argument to `biplot()` ensures that the arrows are scaled to represent the loadings; other values for `scale` give slightly different biplots with different interpretations.

Notice that this figure is a mirror image of **the left figure on slide 19**. Recall that the principal components are only unique up to a sign change, so we can reproduce the **figure from the slides** by making a few small changes:

```{r chunk11}
pr.out$rotation = -pr.out$rotation
pr.out$x = -pr.out$x
biplot(pr.out, scale = 0, cex=0.5)
```



The `prcomp()` function also outputs the standard deviation of each principal component. For instance, on the `USArrests` data set, we can access these standard deviations as follows:

```{r chunk12}
pr.out$sdev
```

The variance explained by each principal component is obtained by squaring these:

```{r chunk13}
pr.var <- pr.out$sdev^2
pr.var
```

To compute the proportion of variance explained by each principal component, we simply divide the variance explained by each principal component by the total variance explained by all four principal components:

```{r chunk14}
pve <- pr.var / sum(pr.var)
pve
```

We see that the first principal component explains $62.0 \%$ of the variance in the data, the next principal component explains $24.7 \%$ of the variance, and so forth.
We can plot the PVE explained by each component, as well as the cumulative PVE, as follows:

```{r chunk15}
par(mfrow = c(1, 2))
plot(pve, xlab = "Principal Component",
    ylab = "Proportion of Variance Explained", ylim = c(0, 1),
    type = "b", col="peru", pch=20)
plot(cumsum(pve), xlab = "Principal Component",
    ylab = "Cumulative Proportion of Variance Explained",
    ylim = c(0, 1), type = "b", col="peru", pch=20)
```

Note that the function `cumsum()` computes the cumulative sum of the elements of  a numeric vector.