---
title: Principal Components Analysis II
author: 36-467/667
date: 11 September 2018
output: slidy_presentation
bibliography: locusts.bib
---

\[
\newcommand{\X}{\mathbf{x}}
\newcommand{\w}{\mathbf{w}}
\newcommand{\V}{\mathbf{v}}
\newcommand{\S}{\mathbf{s}}
\newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]}
\newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]}
\newcommand{\SampleVar}[1]{\widehat{\mathrm{Var}}\left[ #1 \right]}
\newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]}
\newcommand{\TrueRegFunc}{\mu}
\newcommand{\EstRegFunc}{\widehat{\TrueRegFunc}}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator*{\argmin}{argmin}
\DeclareMathOperator{\dof}{DoF}
\DeclareMathOperator{\det}{det}
\newcommand{\TrueNoise}{\epsilon}
\newcommand{\EstNoise}{\widehat{\TrueNoise}}
\]

```{r, include=FALSE}
# General set-up options
library(knitr)
opts_chunk$set(size="small",background="white", highlight=FALSE,
               cache=TRUE, autodep=TRUE,
               tidy=TRUE, warning=FALSE, message=FALSE,
               echo=FALSE)
```


## In our last episode...

- We have multivariate data $\X$ (dimension = $[n\times p]$)
- We want the best $q$-dimensional _linear_ approximation
- Solution: **Principal components analysis**
    + Take the $q$ leading eigenvectors of $\V \equiv \frac{1}{n}\X^T \X =$ sample/empirical covariance matrix
	+ These eigenvectors = the principal components
	+ Project the data on to the principal components
- Assemble the eigenvectors into $\w$ ($[p\times q]$)
    + **Scores** for the data are $\X\w \equiv \S$ ($[n\times q]$)
	+ **Approximations** are $(\X\w) \w^T = \S \w^T$ ($[n\times p]$)

## Some properties of the PCs

- The principal components are orthonormal
    + $\vec{w}_i \cdot \vec{w}_i = 1$
    + $\vec{w}_i \cdot \vec{w}_j = 0$ (unless $i=j$)
    + $\w^T\w = \mathbf{I}$
- PC1 is the direction of maximum variance through the data
    + That variance is $\lambda_1$, biggest eigenvalue of $\V$
- PC $i+1$ is the direction of maximum variance $\perp$ PC1, PC2, $\ldots$ PC $i$
    + That variance is $\lambda_{i+1}$

## Some properties of the eigenvalues

- All eigenvalues $\geq 0$
- In general, $p$ non-zero eigenvalues
- If the data are _exactly_ in a $q$-dimensional subspace, then exactly $q$ non-zero eigenvalues
- If $n < p$, at most $n$ non-zero eigenvalues
   + Two points define a line, three define a plane, ...

## Some properties of PC scores

- Average score on each PC $=0$ (b/c we centered the data)
- Variance of score on PC $i$ $=\lambda_i$ (by construction)
- Covariance of score on PC $i$ with score on PC $j$ $=0$

\begin{eqnarray}
\Var{\text{scores}} & = & \frac{1}{n} \S^T \S\\
& = & \frac{1}{n} (\X\w)^T(\X\w)\\
& = & \frac{1}{n}\w^T \X^T \X \w\\
& = & \w^T \V\w = \w^T \mathbf{\Lambda} \w\\
& = & \mathbf{\Lambda} \w^T\w\\
& = & \mathbf{\Lambda}
\end{eqnarray}

## Some properties of PCA as a whole

- If we use all $p$ principal components, we have the **eigendecomposition** of $\V$:
\[
\V = \w \mathbf{\Lambda} \mathbf{w}^T
\]
$\mathbf{\Lambda}=$ diagonal matrix of eigenvalues $\lambda_1, \ldots \lambda_p$
- If we use all $p$ principal components,
\[
\X = \S\w^T
\]
- If we use only the top $q$ PCs, we get:
    + the best rank-$q$ approximation to $\V$
    + the best dimension-$q$ approximation to $\X$



## Another way to think about PCA

- The original coordinates are correlated
- There is always another coordinate system with _uncorrelated_ coordinates
- We're rotating to that coordinate system
    + Rotating to new coordinates $\Rightarrow$ multiplying by an orthogonal matrix
    + That matrix is $\mathbf{w}$
    + The new coordinates are the scores


## PCA can be used for any multivariate data

- Nothing in the math of PCA cares about where the data came from
    + $n$ measurements on $p$ variables is all that matters
- Areas of application:
    + Reducing multiple measurements
    + Dealing with collinearity in regression
    + Recommendation engines (e.g. Netflix)
    + Gene expression levels in molecular biology
    + [Fashion](http://graceavery.com/principalcomponentanalysisandfashion/)
- ... and, of course, spatio-temporal data

## PCA with spatial data

- $n$ locations for $p$ variables
    + We saw this!
    + Each PC is a $p$-dimensional vector
    + _Scores_ are distributed over space
- vs. $n$ variables at $p$ locations
    + Each PC is a spatial pattern
    + One score for each original variable

## Recall the states... {.smaller}

```{r, echo=TRUE}
state.pca <- prcomp(state.x77,scale.=TRUE)
signif(state.pca$rotation[,1:2], 2)
```

states are locations, PCs are patterns of variables

## Each score is spatially distributed

```{r}
plot.states_scaled <- function(sizes,min.size=0.4,max.size=2,text.color="black",
                               ...) {
  plot(state.center,type="n",...)
  out.range = max.size - min.size
  in.range = max(sizes)-min(sizes)
  scaled.sizes = out.range*((sizes-min(sizes))/in.range)
  text(state.center,state.abb,cex=scaled.sizes + min.size, col=text.color)
  invisible(scaled.sizes)
}

state.confed <- rep(0,times=50)
names(state.confed) <- rownames(state.x77)
state.confed[c("South Carolina", "Mississippi", "Florida", "Alabama",
               "Georgia", "Louisiana", "Texas", "Virginia", "Arkansas",
               "Tennessee", "North Carolina")] <- 1
state.slave <- state.confed
state.slave[c("Kentucky","Missouri","Maryland","Delaware","West Virginia")] <- 1

state.colors <- rep("blue", 50)
state.colors[state.confed==1] <- "grey"
state.colors[(state.slave==1) & (state.confed==0)] <- "lightblue"

plot.states_scaled(state.pca$x[,1],min.size=0.3,max.size=1.5,
                   xlab="longitude",ylab="latitude",
                   text.color=state.colors)
```

## Try it the other way {.smaller}

Turn the data on its side

```{r, echo=TRUE}
state.vars.pca <- prcomp(t(scale(state.x77))) # What's t()?
length(state.vars.pca$sdev) # Why 8?
head(signif(state.vars.pca$rotation[,1:2]),2)
signif(state.vars.pca$x[,1],2)
```

## The states turned on their sides...

```{r}
# Why do I stick a minus sign in here?
plot.states_scaled(-state.vars.pca$rotation[,1], min.size=0.3, max.size=1.5,
                   xlab="longitude", ylab="latitude",
                   main="PC1 of spatial variation for the state.x77 variables",
                   text.color=state.colors)
```

## Not exactly the same

```{r}
plot(x=state.pca$x[,1], y=state.vars.pca$rotation[,1],
     xlab="States' scores on PC1 when states are observations",
     ylab="PC1 when states are variables")
```

... but pretty close.  No coincidence!  (See end)

## 2nd principal component

```{r}
# Again, why the minus sign?
plot.states_scaled(-state.vars.pca$rotation[,2], min.size=0.3, max.size=1.5,
                   xlab="longitude", ylab="latitude",
                   main="PC2 of spatial variation for the state.x77 variables",
                   text.color=state.colors)
```


## A famous example

- For each of $n$ variant genes ("alleles"):
    + Measure prevalence of allele at $p$ different locations
	+ Originally by blood tests, now gene sequencing
	+ Each observation is a _fraction_
- Accumulated over many thousands of genes and of locations
- Outstanding reference: @Cavalli-Sforza-history-and-geography
    + or the good-parts version, @Cavalli-Sforza-genes-peoples-languages

##  Some maps from @Cavalli-Sforza-et-al-demic-expansions

World PC1

![](http://www.stat.cmu.edu/~cshalizi/dst/18/lectures/05/PC1_global.png)

($\approx 35\%$ of between-population variance)

## Some maps from @Cavalli-Sforza-et-al-demic-expansions

World PC2

![](http://www.stat.cmu.edu/~cshalizi/dst/18/lectures/05/PC2_global.png)

($\approx 18\%$ of between-population variance)

## Some maps from @Cavalli-Sforza-et-al-demic-expansions

World PC3

![](http://www.stat.cmu.edu/~cshalizi/dst/18/lectures/05/PC3_global.png)

($\approx 12\%$ of between-population variance)

## Some maps from @Cavalli-Sforza-et-al-demic-expansions

![](http://www.stat.cmu.edu/~cshalizi/dst/18/lectures/05/PCs1--3_global.png)

- green = PC1, blue = PC2, red = PC3
- No sharp breaks...
   + This is $\approx 65\%$ of between-population variance
   + $\approx 90\%$ of variance across people is _within_ populations

## Some maps from @Cavalli-Sforza-et-al-demic-expansions

![](http://www.stat.cmu.edu/~cshalizi/dst/18/lectures/05/PC1_Europe.png)

- PC1 for Europe and Southwest Asia
   + Agriculture starts in the Fertile Crescent
   + Farmers (not just farming) spread

## Some maps from @Cavalli-Sforza-et-al-demic-expansions {.smaller}

![](http://www.stat.cmu.edu/~cshalizi/dst/18/lectures/05/PC3_Europe.png)

- PC3 for Europe and Southwest Asia
    + _Something_ centered on the steppes north of the Caucasus mountains and the Black Sea
    + Archaeology: Domestication of the horse, chariots
    + Linguistics: origins of the Indo-European languages
    + Most likely: Indo-European speaking barbarians expanding out from the steppes [@Anthony-on-indo-europeans]
    + Similar patterns for Asia



## PCA with multiple time series

- Obvious approach: $n$ time points for $p$ variables
    + Each PC is a $p$-dimensional vector $\Rightarrow$ a pattern _across variables_
- Flip perspective: $n$ variables for $p$ time-points
    + Each PC is a $n$-dimensional vector $\Rightarrow$ a pattern _across time_


## Irish wind data

```{r, include=FALSE}
# Following the example in help(wind)
library(gstat)
data(wind)
library(sp) # For some spatial manipiulations
wind.loc$y = as.numeric(char2dms(as.character(wind.loc[["Latitude"]])))
wind.loc$x = as.numeric(char2dms(as.character(wind.loc[["Longitude"]])))
coordinates(wind.loc) = ~x+y
wind$time = ISOdate(wind$year+1900, wind$month, wind$day)
```

## Irish wind data

```{r}
head(wind)
```

## Irish wind data

```{r}
require(mapdata)
map("worldHires", xlim = c(-11,-5.4), ylim = c(51,55.5))
plot(wind.loc, add=TRUE, pch=16)
text(coordinates(wind.loc), pos=1, label=wind.loc$Station)
```

## Irish wind data --- one time series

```{r}
plot(DUB~time, wind, type= 'l', ylab = "windspeed (knots)", main = "Dublin")
```

## Irish wind data --- all the time series

```{r}
# Need to convert the time column...
require(lubridate)
matplot(decimal_date(wind$time), wind[,4:15], type="l",
        ylab="windspeed (knots)", xlab="time")
```

## PCA: $n = `r nrow(wind)`$, $p=12$

```{r, echo=TRUE}
wind.pca.1 <- prcomp(wind[,4:15])
wind.pca.1$sdev
```

## PC1: The eigenvector {.smaller}

```{r, echo=TRUE}
plot(-wind.pca.1$rotation[,1], ylim=c(0, 1))
text(1:12, -wind.pca.1$rotation[,1], pos=3,
     labels=colnames(wind)[4:15])
```

A pattern over _space_

## PC1: The eigenvector {.smaller}

```{r}
map("worldHires", xlim = c(-11,-5.4), ylim = c(51,55.5))
plot(wind.loc, add=TRUE, pch=16, cex=4,
     col=heat.colors(6)[cut(-wind.pca.1$rotation[,1], breaks=6)])
legend("topleft", legend=levels(cut(-wind.pca.1$rotation[,1], breaks=6)),
       col=heat.colors(6), pch=16, cex=0.5, pt.cex=1)
```

A function of space

## PC1: The scores

```{r}
plot(decimal_date(wind$time), wind.pca.1$x[,1],
     xlab="time", ylab="windspeed (knots)",
     type="l")
```

A function of _time_

## Break for in-class exercise:

```{r}
par(mfrow=c(1,2))
map("worldHires", xlim = c(-11,-5.4), ylim = c(51,55.5))
plot(wind.loc, add=TRUE, pch=16, cex=4,
     col=heat.colors(6)[cut(-wind.pca.1$rotation[,1], breaks=6)],
     main="PC1 vector")
legend("topleft", legend=levels(cut(-wind.pca.1$rotation[,1], breaks=6)),
       col=heat.colors(6), pch=16, cex=0.5, pt.cex=1)
plot(decimal_date(wind$time), wind.pca.1$x[,1],
     xlim=c(min(decimal_date(wind$time)), 1968),
     xlab="time", ylab="windspeed (knots)",
     type="l", main="Scores on PC1")
```

_Describe the first component_

## PCA with spatio-temporal data

- ... we just did this!
- $n$ spatial locations, each with a time series of length $p$
    + PCs are spatial patterns
- $n$ time points for $p$ locations
    + PCs are time series / temporal patterns
    + Trends, or components of trends?

##  Interpreting PCA results

- PCs are linear combinations of the original coordinates
    + $\Rightarrow$ PCs change as you add or remove coordinates
	+ Put in 1000 measures of education and PC1 is education...
- Very tempting to **reify** the PCs
    + i.e., to "make them a thing"
	+ sometimes totally appropriate...
	+ sometimes not at all appropriate...
	+ Be very careful when the _only_ evidence is the PCA
	+ Smoothing artifacts can be deadly [@Novembre-Stephens-interpreting-pca]

## PCA is exploratory analysis, not statistical inference

- We assumed no model
    + Pro: That's the best linear approximation, no matter what
	+ Con: doesn't tell us how much our results are driven by noise
- Prediction: PCA predicts nothing
- Inference: _If_ $\V \rightarrow \Var{X}$ _then_ then PCs $\rightarrow$ eigenvectors of $\Var{X}$
    + But PCA doesn't **need** this assumption
    + Doesn't tell us about uncertainty
	+ Will see ways to tackle this by simulation

## Some alternatives to PCA

- **Independent component analysis** (ICA)
    + PCA analyzes data into _uncorrelated_ components
    + But uncorrelated $\neq$ independent (unless everything's Gaussian)
	+ **ICA** analyzes data into statistically-independent components
	+ Measure dependence between components, minimize
	+ Good overview: @Stone-on-ICA
- **Nonlinear** approximation
    + PCA finds low-dimensional _linear_ approximation
	+ What of things are curved?
	+ Locally linear embedding, spectral component analysis, ...
	+ Could do worse  overview: @CRS-ADAfaEPoV

## Summing up

- PCA rotates to new, uncorrelated coordinates
- Using the first $q$ PCs gives the best $q$-dimensional approximation to the data
- These are the $q$ directions of largest variance
- We can make either the basis vectors or the scores into spatio-temporal patterns
- Interpretation needs domain knowledge, and some caution
- No inference or prediction: that comes next


## Orthogonal matrices

- Matrix $\mathbf{o}$ is **orthogonal** iff $\mathbf{o}^T = \mathbf{o}^{-1}$
- $\Leftrightarrow$ the columns of $\mathbf{o}$ are orthonormal vectors
    + Why we say "orthogonal matrix" rather than "orthonormal matrix" is lost in the mists of 19th century German mathematics
- Every rotation (around the origin) corresponds to an orthogonal matrix
- Are there orthogonal matrices which aren't rotations?

## PCA of $\X$ vs. PCA of $\X^T$

- Starting from $\X$:
    + Eigenvectors $\vec{w}_i$, eigenvalues $\lambda_i$
	+ $\X = \S \w^T$
	+ $n^{-1} \X^T \X = \w \mathbf{\Lambda} \w^T$
- Starting from $\X^T$
    + Eigenvectors $\vec{u}_i$, eigenvalues $\psi_i$
	+ $n^{-1} \X \X^T = \mathbf{u}\mathbf{\Psi}\mathbf{u}^T$

## PCA of $\X$ vs. PCA of $\X^T$

\begin{eqnarray}
\mathbf{u}\mathbf{\Psi}\mathbf{u}^T & = & n^{-1} \X \X^T\\
\mathbf{u}\mathbf{\Psi}\mathbf{u}^T & = & n^{-1} \S \w^T \w \S^T\\
\mathbf{u}\mathbf{\Psi}^{1/2} \mathbf{\Psi}^{1/2}\mathbf{u}^T & = & n^{-1/2} \S \w^T \w \S^T n^{-1/2}\\
(\mathbf{u}\mathbf{\Psi}^{1/2}) (\mathbf{u}\mathbf{\Psi}^{1/2})^T & = & n^{-1/2} \S \S^T n^{-1/2}\\
\mathbf{u} & = & n^{-1/2} \mathbf{\Psi}^{-1/2} \S
\end{eqnarray}

New PC1 $\propto$ old scores on PC1, etc.


## No, really, PCA doesn't do statistical inference

- PCA is _close_ to a statistical model called factor analysis
    + but PCs $\neq$ factors
	+ Will see factor models late in the course

## Other alternatives to PCA

	
- **Slow feature analysis**
    + PCs of multiple time series are trend-ish
	+ Trends ought to change slowly
	+ SFA finds components with high correlation over time
- **Forecastable** component analysis [@GMG-ForeCA]
    + we'll come back to after looking at prediction




## References




