--- title: Principal Components Analysis II author: 36-467/667 date: Lecture 6 (17 September 2020) 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]$) + There are $n$ scores on PC $i$, $\vec{x}_1 \cdot \vec{w}_i, \vec{x}_2 \cdot \vec{w}_2, \ldots \vec{x}_n \cdot \vec{w}_i$ + The matrix $\S$ assembles all the scores - **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$) + Or in matrix form: $\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 (but not always) - 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 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$ ## 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 ~\text{ by definition of} ~ \V\\ & = & \w^T ( \w \mathbf{\Lambda} \mathbf{w}^T) \w ~\text{by eigendecomposition}\\ & = & (\w^T \w) \mathbf{\Lambda} (\w^T\w)\\ & = & \mathbf{\Lambda} \end{eqnarray} ## 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 or high-dimensional covariates in regression (see backup) + [Recommendation engines](http://www.stat.cmu.edu/~cshalizi/dm/20/lectures/22/lecture-22.html) (e.g. Netflix [@Statistical-significance-of-netflix]) * $n$ users ratings or engagement with $p$ items of content, used to predict which users will like/engage with which items * Your feed on Facebook, Twitter, etc. is a recommendation engine that's trying to maximize your engagement + Gene expression levels in molecular biology [@Wall-Rechststeiner-Rocha-SVD-and-PCA] + [Fashion](http://graceavery.com/principalcomponentanalysisandfashion/) - ... and, of course, spatio-temporal data ## PCA with spatial data - $n$ locations for $p$ variables + We saw this [last time](http://www.stat.cmu.edu/~cshalizi/dst/20/lectures/05/lecture-05.html) with the states + Each PC is a $p$-dimensional vector (in "feature space") + _Scores_ are distributed over physical space * That is, each location has a score (on each PC) - vs. $n$ variables at $p$ locations + Each PC is a spatial pattern + One score for each original variable + See backup for examples of doing this with the states from last time - These are _almost_ interchangeable + Again, see backups for the math ## 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_ $\in [0,1]$ - Accumulated over many thousands of genes ($n$) and of locations ($p$) - 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 - Can do similar analyses on smaller geographic scales and learn more about history (see backups) ## PCA with multiple time series - One 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, col="darkblue") ``` ## 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_ ## Try to describe the first component here ```{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") ``` ## 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... > * Sometimes this is even what you want [@Zeller-Carmines-measurement] > - 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 [@Glymour-on-bell-curve] > + Smoothing artifacts can be deadly [@Novembre-Stephens-interpreting-pca] > * In essence: Yule and Slutsky strike again ## PCA is exploratory analysis, not statistical inference > - We assumed no model > + Pro: That's the best linear approximation _to these data_, no matter what > + Con: doesn't tell us where the data came from, what other data will look like, or how much our results are driven by noise > - Prediction: PCA predicts nothing > - Inference: _If_ $\V \rightarrow \Var{X}$ _then_ PCs $\rightarrow$ eigenvectors of $\Var{X}$ > + But PCA doesn't **need** this assumption > + Doesn't tell us about uncertainty > + We 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** tries to break data into statistically-independent additive components + Measure dependence between components, minimize + Good overview: @Stone-on-ICA - **Nonlinear** approximation + PCA finds low-dimensional _linear_ approximation + What if the real structure has curves? + 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 - PCA does no inference or prediction ## Details and asides ## Some more 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 out ## Some more 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 plains ("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 explanation: Indo-European-speaking barbarians using horses and chariots to expand out from the steppes [@Anthony-on-indo-europeans] + Similar patterns for Central and South Asia ## Principal Components Regression - Say we want to linearly regress $Y$ on $X$, with $X$ being $p$-dimensional - If $p$ is too big, for whatever reason, we could do PCA on $X$ and get a $q$-dimensional $S$ - Then we regress $Y$ on $S$ and get our predictions that way - This is **principal components regression** - Variant: $X=(Z,U)$ where $Z$ is the variable we really care about and $U$ is a big mess of controls + Do PCA on $U$, get $q$-dimensional scores $S$, and then regress $Y$ on $Z$ and $S$ + This is less common - Major application of both forms of PC regression: $n < p$, we've measured too many variables on each observation to get a fit by ordinary least squares + Also called high-dimensional regression - Another major application: dealing with collinearity or multi-collinearity + Multicollinearity means that the $p$ predictor variables really live in a lower, $q$-dimensional subspace + Could drop predictor variables until they're no longer collinear + _Or_ you could just find the scores on the $q$ non-zero principal components - No guarantee that the directions of maximum variance through the predictors are the directions along which $Y$ varies $\Rightarrow$ PC regression might not be very good - There are other ways of dealing with high-dimensional or collinear predictors... + Ridge regression, the lasso, the elastic net, etc. + But some recent work suggest PC regression is surprisingly good compared to those modern methods [@Dhillion-et-al-ridge-regression-vs-pca-regression] ## 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 + E.g. to rotate by an angle $\theta$ in two dimensions, the matrix is $\mathbf{o} = \left[ \begin{array}{cc} \cos{\theta} & \sin{\theta} \\ -\sin{\theta} & \cos{\theta} \end{array}\right]$ * If that looks funny compared to what you saw in an earlier class, remember we're writing vectors as $[1 \times p]$ matrices, so the rotation works as $\X\mathbf{o}$ + You can check that $\mathbf{o}^T \mathbf{o} = \left[\begin{array}{cc} 1 & 0 \\ 0 & 1\end{array}\right] = \mathbf{I}$ (because $\cos^2{\theta} + \sin^2{\theta} = 1$) - Are there orthogonal matrices which aren't rotations? ## 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, main="Size scaled to score on PC1") ``` ## 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]),4) 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) ``` - This is the same as the other map (up to rounding errors) - This is no coincidence (see below) ## PCA of $\X$ vs. PCA of $\X^T$ {.smaller} - Starting from $\X$ (an $n\times p$ matrix) + Eigenvectors $\vec{w}_i$, eigenvalues $\lambda_i$ + $\X = \S \w^T$ + Eigendecomposition of the variance matrix is $n^{-1} \X^T \X = \w \mathbf{\Lambda} \w^T$ (a $p\times p$ matrix) - Starting from $\X^T$ (a $p\times n$ matrix) + Eigenvectors $\vec{u}_i$, eigenvalues $\psi_i$. diagonal matrix $\mathbf{\Psi}$ + Eigendecomposition of this variance is $p^{-1} \X \X^T = \mathbf{u}\mathbf{\Psi}\mathbf{u}^T$ (an $n\times n$ matrix) \begin{eqnarray} \mathbf{u}\mathbf{\Psi}\mathbf{u}^T & = & p^{-1} \X \X^T\\ \mathbf{u}\mathbf{\Psi}\mathbf{u}^T & = & p^{-1} \S \w^T (\S \w^T)^T\\ \mathbf{u}\mathbf{\Psi}\mathbf{u}^T & = & p^{-1} \S \w^T \w \S^T\\ \mathbf{u}\mathbf{\Psi}^{1/2} \mathbf{\Psi}^{1/2}\mathbf{u}^T & = & p^{-1/2} \S \w^T \w \S^T p^{-1/2}\\ (\mathbf{u}\mathbf{\Psi}^{1/2}) (\mathbf{u}\mathbf{\Psi}^{1/2})^T & = & p^{-1/2} \S \S^T p^{-1/2}\\ \mathbf{u} & = & p^{-1/2} \mathbf{\Psi}^{-1/2} \S \end{eqnarray} New PC1 vector $\propto$ old scores on PC1, etc. ```{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") ``` ## No, really, PCA doesn't do statistical inference {.smaller} - PCA _resembles_ a statistical model called factor analysis - The factor model is $\vec{X} = \mathbf{\omega} \vec{F} + \vec{\epsilon}$ + $\vec{X}$ is the $p$-dimensional random vector we observe + $\vec{F}$ is the $q$-dimensional random vector of hidden ("latent") **factors** or **factor scores** * Usually assume $q < p$ * Usually assume $\Var{\vec{F}} = \mathbf{I}$ + $\mathbf{\omega}$ is a $[p\times q]$ matrix of **loadings** + $\vec{\epsilon}$ is random noise uncorrelated with $\vec{F}$ - This is a statistical model which can generate new data - It's can also make predictions: + Distribution of new data points + If $\vec{X}$ has some missing values, can still estimate $\vec{F}$ and then use that to predict the unobserved entries in $\vec{X}$ * [This is really what recommendation engines do](http://www.stat.cmu.edu/~cshalizi/dm/20/lectures/22/lecture-22.html) - This looks similar, but PCs $\neq$ factors + $\mathbf{\omega} \neq \mathbf{w}$ + Scores on PCs aren't even _estimates_ of factor values - We'll come back to factor models later, but you could do worse than to read the factor-analysis chapter of @CRS-ADAfaEPoV ## 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] + Finds _highly predictable_ components ## References