--- 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