--- title: Principal Components Analysis author: 36-740/620 date: Lecture 3 (1 September 2026) output: slidy_presentation: math_method: engine: mathjax url: https://bactra.org/mathjax/tex-svg.js bibliography: locusts.bib --- ```{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) ``` \[ \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}} \] ## In our last episode... > - Data is a vector >- Break the vector up into additive components > - Each component is a vector, which is a pattern > - Which components? > - Before: components were eigenvectors of the smoother matrix > + Made sense to work with them because we were smoothing > - Can we get the data to tell us what the "right" components are? ## What the data looks like {.smaller} - Data is $\mathbf{x} = n\times p$ matrix + Could be $p$-dimensional multivariate data at $n$ time points + Or multivariate data at $n$ locations + Or scalar data at $n$ time points and $p$ locations - Write $\vec{x}_i$ for row $i$ ($1\times p$ matrix) - First, center the data, just to reduce book-keeping + i.e., subtract the mean of each column + Optional: scale each column to equal variance ## Finding a "principal" component - We don't want to keep track of $p$ dimensions - We want _one_ dimension - We also don't want to distort the data too much - Pick a direction in the $p$-dimensional space, and a length-1 vector $\vec{w}$ - What's the best $\vec{w}$? ## Projections - $\vec{x}_i \cdot \vec{w} =$ length of $\vec{x}_i$'s projection on to the direction of $\vec{w}$ - $(\vec{x}_i \cdot \vec{w})\vec{w} =$ the actual projected vector ## Projections ```{r} # Make an arc demo.theta <- runif(10,min=0,max=pi/2) demo.x <- cbind(cos(demo.theta),sin(demo.theta)) # Center the coordinates demo.x <- scale(demo.x, center=TRUE, scale=FALSE) ``` ```{r, basic.plot} # Square plot! # Otherwise the aspect ratio's all wrong par(pty="s") # Plot the points plot(demo.x,xlab=expression(x^1),ylab=expression(x^2), xlim=c(-1,1), ylim=c(-1,1), main="Projections of data on to a vector") # Plot the origin points(x=0, y=0, pch=3, col="grey", cex=2) # pch=3 is a plus sign # Pick a direction (not a very good one), by a unit vector demo.w <- c(cos(-3*pi/8), sin(-3*pi/8)) # Draw an arrow for the unit vector arrows(0,0,demo.w[1],demo.w[2],col="blue") text(demo.w[1],demo.w[2],pos=4,labels=expression(w)) # Draw a dashed line for the whole line abline(0,b=demo.w[2]/demo.w[1],col="blue",lty="dashed") # Get the length of the projection of each data point on to the line projection.lengths <- demo.x %*% demo.w # Get the actual projections projections <- projection.lengths %*% demo.w # Draw those as points along the line points(projections,pch=16,col="blue") # Draw lines from each data point to its projection segments(x0=demo.x[,1], y0=demo.x[,2],x1=projections[,1],y1=projections[,2], col="grey") ``` ## Projections - $\vec{x}_i \cdot \vec{w} =$ length of $\vec{x}_i$'s projection on to the direction of $\vec{w}$ - $(\vec{x}_i \cdot \vec{w})\vec{w} =$ the actual projected vector - $\vec{x}_i - (\vec{x}_i \cdot \vec{w})\vec{w} =$ vector difference between data and projected vectors - $\|\vec{x}_i - (\vec{x}_i \cdot \vec{w})\vec{w}\|^2=$ squared error from replacing data vector with projected vector ## How well does the projection approximate the original? Do it for one vector first: \begin{eqnarray} {\|\vec{x_i} - (\vec{w}\cdot\vec{x_i})\vec{w}\|}^2 & =& \left(\vec{x_i} - (\vec{w}\cdot\vec{x_i})\vec{w}\right)\cdot\left(\vec{x_i} - (\vec{w}\cdot\vec{x_i})\vec{w}\right)\\ & = & \vec{x_i}\cdot\vec{x_i} -\vec{x_i}\cdot (\vec{w}\cdot\vec{x_i})\vec{w}\\ \nonumber & & - (\vec{w}\cdot\vec{x_i})\vec{w}\cdot\vec{x_i} + (\vec{w}\cdot\vec{x_i})\vec{w}\cdot(\vec{w}\cdot\vec{x_i})\vec{w}\\ & = & {\|\vec{x_i}\|}^2 -2(\vec{w}\cdot\vec{x_i})^2 + (\vec{w}\cdot\vec{x_i})^2\vec{w}\cdot\vec{w}\\ & = & \|\vec{x_i}\|^2 - (\vec{w}\cdot\vec{x_i})^2 \end{eqnarray} (This is the Pythagorean theorem) ## How well does the projection approximate the original? Average across all the data vectors: \begin{eqnarray} MSE(\vec{w}) & = & \frac{1}{n}\sum_{i=1}^{n}{ {\|\vec{x_i} - (\vec{w}\cdot\vec{x_i})\vec{w}\|}^2}\\ & = & \frac{1}{n}\sum_{i=1}^{n}{\left( \|\vec{x_i}\|^2 -{(\vec{w}\cdot\vec{x_i})}^2 \right)}\\ & = & \frac{1}{n}\sum_{i=1}^{n}{\|\vec{x_i}\|^2} -\frac{1}{n}\sum_{i=1}^{n}{(\vec{w}\cdot\vec{x_i})^2} \end{eqnarray} - First bit doesn't depend on $\vec{w}$, so doesn't matter for minimizing - So we want to maximize \[ L(\vec{w}) = \frac{1}{n}\sum_{i=1}^{n}{{(\vec{w}\cdot\vec{x_i})}^2} \] ## Minimizing MSE is maximizing variance \begin{eqnarray} L(w) & = & \frac{1}{n}\sum_{i=1}^{n}{{(\vec{w}\cdot\vec{x_i})}^2}\\ & = & {\left(\frac{1}{n}\sum_{i=1}^{n}{\vec{x_i}\cdot\vec{w}}\right)}^2 + \SampleVar{\vec{w}\cdot\vec{x_i}} \end{eqnarray} (Because: $\Expect{Z^2} = (\Expect{Z})^2 + \Var{Z}$) Centering simplifies our book-keeping (as promised): \[ \frac{1}{n}\sum_{i=1}^{n}{\vec{x_i} \cdot \vec{w}} = 0 \] $\therefore$ \[ L(\vec{w}) = \SampleVar{\vec{w}\cdot\vec{x_i}} \] ## Minimizing MSE is maximizing variance > The direction which gives us the best approximation of the data is the direction with the greatest variance ## OK, how do we find this magic direction? Matrix form: all the _lengths_ of projections is $\mathbf{x}\mathbf{w}$ $[n\times 1]$ \begin{eqnarray} \SampleVar{\vec{w}\cdot\vec{x_i}} & = & \frac{1}{n}\sum_{i}{{\left(\vec{x_i} \cdot \vec{w}\right)}^2}\\ & = & \frac{1}{n}{\left(\X \w\right)}^{T} \left(\X \w\right)\\ & = & \frac{1}{n} \w^T \X^T \X \w\\ & = & \w^T \frac{\X^T \X}{n} \w\\ \end{eqnarray} - Fact: $\V \equiv \frac{\X^T \X}{n} =$ sample covariance matrix of the vectors + Because we've centered; $v_{jk} = n^{-1}\sum_{i=1}^{n}{x_{ij} x_{ik}} =$ sample covariance of coordinate $x^{(j)}$ with coordinate $x^{(k)}$ + Same reason $\X^T \X$ shows up in formulas for linear regression ## OK, how do we find this magic direction? > - We need to maximize \begin{equation} \SampleVar{\vec{w}\cdot\vec{x_i}} = \w^T \V \w \end{equation} > - Constraint: $\vec{w}$ has length 1 $\Leftrightarrow$ $\w^T \w = 1$ > - Introduce a [**Lagrange multiplier**]{#lagrange} $\lambda$ \begin{eqnarray} \mathcal{L}(\w,\lambda) & \equiv & \w^T\V\w - \lambda(\w^T \w -1)\\ \frac{\partial \mathcal{L}}{\partial \lambda} & = & \w^T \w -1\\ \frac{\partial \mathcal{L}}{\partial \w} & = & 2\V\w - 2\lambda\w \end{eqnarray} > - Set derivatives to zero: \begin{eqnarray} \w^T \w & = & 1\\ \V \w & = & \lambda \w \end{eqnarray} ## The magic direction is an eigenvector \begin{eqnarray} \w^T \w & = & 1\\ \V \w & = & \lambda \w \end{eqnarray} > THIS IS AN EIGENVALUE/EIGENVECTOR EQUATION! The _value_ of the the solution is \[ \SampleVar{\vec{w}\cdot\vec{x_i}} = \w^T \V \w = \w^T \lambda \w = \lambda \] so the maximum is the _leading_ eigenvector of $\V$ ## About the sample covariance matrix $\V$ is a special matrix: [symmetric and non-negative definite]{#variance-is-symmetric-and-nnd} - Eigenvalues are all real (b/c symmetric) - If $\lambda_i \neq \lambda_j$, then $v_i \perp v_j$ (b/c symmetric) - Eigenvectors form a basis, and are orthonormal + Or can be chosen to be orthonormal - All $\lambda_i \geq 0$ (b/c non-negative definite) \begin{eqnarray} \text{Lead eigenvector of}\ \V & = & 1^{\mathrm{st}}\ \text{principal component}\\ & = & \text{Direction of maximum variance}\\ & = & \text{Best 1D approximation to the data} \end{eqnarray} ## ```{r, include=FALSE} # Not shown, in the interest of time, but you should look at these, # person who is examining the Rmd file head(demo.x) var(demo.x) eigen(var(demo.x)) ``` ```{r} <> ``` ## ```{r} <> # Find the variance-covariance matrix of the data v <- var(demo.x) # First principal component is the leading eigenvector pc.1 <- eigen(v)$vectors[,1] # Draw the arrow and the dashed line, but in a different color arrows(0,0,pc.1[1],pc.1[2],col="red") abline(0,b=pc.1[2]/pc.1[1],col="red",lty="dashed") text(pc.1[1], pc.1[2], pos=4, labels="PC1") # Get the projection lengths and actual projections projection.lengths <- demo.x %*% pc.1 projections <- projection.lengths %*% pc.1 # Draw the points points(projections, pch=16, col="red") segments(x0=demo.x[,1], y0=demo.x[,2], x1=projections[,1], y1=projections[,2], col="pink") legend("topright", legend=c("Data", "Projection on arbitrary w", "Projection on PC1", "Error with arbitrary w", "Error with PC1"), pch=c(1, 16, 16, NA, NA), lty=c("blank", "blank", "blank", "solid", "solid"), col=c("black", "blue", "red", "grey", "pink"), cex=0.75) ``` ## Multiple principle components > - What about approximating by a plane, hyper-plane, hyper-hyper-plane, etc.? > - Intuition: take the direction of maximum variance $\perp$ the first principal component > - Then direction of maximum variance $\perp$ the first two principal components > - These are the eigenvectors of $\V$, in order of decreasing $\lambda$ > - The [gory details can be found at the end of these slides]{#gorydetails} ## 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$ - If $n \geq p$, then _generally_ $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 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, \begin{eqnarray*} \S & = & \X \w\\ \X & = & \S\w^T \end{eqnarray*} - 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]{#orthogonal} + 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 (original idea) + [Dealing with collinearity or high-dimensional covariates in regression]{#pcreg} + [Recommendation engines](http://www.stat.cmu.edu/~cshalizi/dm/22/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 [cf. "eigendresses"](https://www.stat.cmu.edu/~cshalizi/dm/22/hw/10/hw-10.pdf)) - ... and, of course, spatio-temporal data ## PCA with spatial data - $n$ locations for $p$ variables + Each PC is a $p$-dimensional vector (in the space of variables) + _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]{#exchanging-rows-and-columns} for examples of doing this with the states - These are _almost_ interchangeable + Again, see [backups]{#exchanging-rows-and-columns} for the math ## In R - `prcomp` is the usual built-in PCA command - Slightly complicated object returned - Understand it with a (spatial) data example ## USA, $\approx 1977$ Dataset pre-loaded in R: ```{r, echo=TRUE} head(state.x77) ``` ## Principal components of the USA, $\approx 1977$ ```{r, echo=TRUE} state.pca <- prcomp(state.x77,scale.=TRUE) str(state.pca) ``` ## Principal components of the USA, $\approx 1977$ The weight/loading matrix $\w$ gets called `$rotation` (why?): ```{r, echo=TRUE} signif(state.pca$rotation[,1:2], 2) ``` Each column is an eigenvector of $\V$ ## Principal components of the USA, $\approx 1977$ ```{r, echo=TRUE} signif(state.pca$sdev, 2) ``` Standard deviations along each principal component $=\sqrt{\lambda_i}$ If we keep $k$ components, \[ R^2 = \frac{\sum_{i=1}^{k}{\lambda_i}}{\sum_{j=1}^{p}{\lambda_j}} \] (Denominator $=\tr{\V}$ [why?]) ## Principal components of the USA, $\approx 1977$ ```{r, echo=TRUE} signif(state.pca$x[1:10, 1:2], 2) ``` Columns are $\vec{x}_i \cdot \vec{w}_1$ and $\vec{x}_i \cdot \vec{w}_2$ ## PC1 is kinda southern ```{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) } plot.states_scaled(state.pca$x[,1],min.size=0.3,max.size=1.5, xlab="longitude",ylab="latitude") ``` - size of state abbreviation $\propto$ projection on to PC1 - coordinates = state capitols, except for AK and HI ## PC1 is kinda the legacy of slavery ```{r} 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) ``` - Correlation of PC1 with having been a slave state in 1861 is `r signif(cor(state.pca$x[,1],state.slave),2)` - Correlation of PC1 with having been in the Confederacy is `r signif(cor(state.confed, state.pca$x[,1]), 2)` ## A famous spatial 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) - [Similar analyses on smaller geographic scales say more about history]{#more-maps} ## 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_ + E.g., across space - 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$ time points for $p$ locations + PCs are spatial patterns + We get a time series of scores for each component - $n$ spatial locations, each with a time series of length $p$ + PCs are time series / temporal patterns + Trends, or components of trends? + We get a score for each location for each component = how much that location participated in that component (over all time) ## 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}$ and $\lambda$s $\rightarrow$ eigenvalues 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, appendix - [Yet more]{#alternatives} alternatives... ## 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 ## Backup: Historical asides {.smaller #history} - PCA was invented by @Pearson-invents-PCA + In the minimize-distance-to-a-plane formulation - PCA was (independently) re-invented by @Hotelling-re-invents-PCA-1, @Hotelling-re-invents-PCA-2 + In the maximize-variance formulation; also provided the name + Also implied (incorrectly!) that $X$ needs to be Gaussian - I am not sure of the subsequent history, and in particular of - The continuous-coordinate version of PCA is the "Karhunen-Loeve transform" (KL transform, KLT), worked out by Karhunen and Loeve in the 1940s, see @Loeve-probability, ch. X, sec. 34.5 + [Hand-wavy version below]{#karhunen-loeve} + Loeve didn't know about either Pearson or Hotelling; I don't know about Karhunen - In many branches of the physical sciences, PCA/KLT is also called "empirical orthogonal functions" - On Cavalli-Sforza's life and work, see @Stone-and-Lurquin-on-Cavalli-Sforza ## Backup: Karhunen-Loeve in one (hand-wavy) slide {.smaller #karhunen-loeve} - **Given**: $X(t)$ defined on a finite domain (e.g., interval) $D$, $\Expect{X(t)}=0$ (centering) - **Desired** $X(t) = \sum_{i=1}^{\infty}{S_i \phi_i(t)}$ with random coefficients $S_i$ + Implies $\Expect{S_i}=0$ - **Require** orthonormality of the functions: $\int_{D}{\phi_i(t) \phi_j(t) dt} = \delta_{ij}$ - **Implies** $\Cov{X(t), X(s)} = \sum_{i=1}^{\infty}{\sum_{j=1}^{\infty}}{\phi_i(t) \phi_j(s) \Cov{S_i, S_j}}$ - **Require** minimal mean squared error, so we pick the $\phi_i$ to minimize $\Expect{\int_{D}{\left(X(t) - \sum_{i=1}^{q}{S_i \phi_i(t)}\right)^2 dt}}$ at each $q$, subject to the orthonormality constraint - **Implies** $\int_{D}{\Cov{X(t), X(s)} \phi_i(s) ds} = \lambda_i \phi_i(t)$ + An eigenproblem again, but now solved by an eigenfunction, not an eigenvector + Eigenproblems for integral operators are classic early-20th-century math [@Courant-Hilbert-methods-of-math-phys, ch. III] * $\Cov{X(t), X(s)}$ is the **kernel** of the integral operator (and yes, stuff like this is [where our "kernel" comes from](https://bactra.org/weblog/two-kinds-of-kernel-methods.html)) - **Implies** $\Cov{S_i, S_j} = \lambda_i \delta_{ij}$ + $\therefore$ $\Cov{X(t), X(s)} = \sum_{i}^{\infty}{\phi_i(t) \phi_i(s) \lambda_i}$ - We need some sort of estimate of the covariance function... ## Backup: Even more abstract K-L - **Given**: function $X$ in some space with an inner product $\langle X, Y \rangle$, with $\Expect{X} = 0$ (centering) + Associated norm $\|X\|^2 = \langle X, X \rangle$ - **Desired**: $X = \sum_{i=1}^{\infty}{S_i \phi_i}$ with random coefficients $S_i$ - **Require** orthonormality, $\langle \phi_i, \phi_j \rangle = \delta_{ij}$ - **Implies** $\Cov{X(t), X(s)} = \sum_{i}{\sum_{j}{\phi_i(t) \phi_j(s) \Cov{S_i, S_j}}}$ - **Require** minimal error, minimize $\Expect{\| X - \sum_{i=1}^{q}{S_i \phi_i}\|^2}$ over choice of $\phi_i$ - Define $\gamma_t(s) \equiv \Cov{X(t), X(s)}$ - Doing the minimization, with the Lagrange multiplier to enforce orthonormality, gives us $\langle \gamma_t, \phi_i \rangle = \lambda_i \phi_i(t)$ which is an eigenvalue equation + $\langle \gamma_t, af + bg \rangle = a \langle \gamma_t, f \rangle + b \langle \gamma_t, g \rangle$ so this is a linear operator - And we get the implication that $\Cov{S_i, S_j} = \lambda_i \delta_{ij}$ ## Backup: $\V$ is symmetric and non-negative definite {.smaller #variance-is-symmetric-and-nnd} - Symmetry: $V_{ij} = \Cov{X^{(i)}, X^{(j)}} = \Cov{X^{(j)}, X^{(i)}} = \V_{ji}$ - Non-negative-definiteness: Pick any $p$-dimensional vector $a$, so $a \cdot X$ is a scalar, therefore $\Var{a \cdot X} \geq 0$ \begin{eqnarray} \Var{a \cdot X} & = & a \cdot \Var{X} a\\ & = & a \cdot \V a \end{eqnarray} But "$\V$ is non-negative definite" _means_ "$a \cdot \V a \geq 0$ for all vectors $a$" To see the matrix trick: \begin{eqnarray} \Var{a \cdot X} & = & \Cov{a \cdot X, a \cdot X}\\ & = & \Cov{\sum_{i=1}^{p}{a_i X_i}, \sum_{j=1}^{p}{a_j X_j}}\\ & = & \sum_{i=1}^{p}{a_i \sum_{j=1}^{p}{\Cov{X_i, X_j} a_j}}\\ & = & \sum_{i=1}^{p}{a_i \sum_{j=1}^{p}{\V_{ij} a_j}}\\ & = & \sum_{i=1}^{p}{a_i (\V a)_i}\\ & = & a \cdot \V a \end{eqnarray} ## (Backup) The gory details for multiple PCs {#gorydetails} Use $k$ directions in a $p\times k$ matrix $\w$ _Require:_ $\mathbf{w}^T\mathbf{w} = \mathbf{I}$, the basis vectors are orthonormal $\X \w =$ matrix of projection lengths $[n\times k]$ $\X \w \w^T =$ matrix of projected _vectors_ $[n\times p]$ $\X - \X \w \w^T =$ matrix of vector residuals $[n\times p]$ $(\X-\X\w\w^T)(\X-\X\w\w^T)^T =$ matrix of inner products of vector residuals $[n\times n]$ $\tr{((\X-\X\w\w^T)(\X-\X\w\w^T)^T)} =$ sum of squared errors $[1\times 1]$ ## Backup: The gory details (cont'd.) \begin{eqnarray} MSE(\w) & = & \frac{1}{n} \tr{((\X-\X\w\w^T)(\X^T - \w\w^T \X^T))}\\ & = & \frac{1}{n} \tr{(\X \X^T - \X\w\w^T\X^T - \X\w\w^T\X^T + \X\w\w^T\w\w^T\X^T)}\\ & = & \frac{1}{n}\left(\tr{(\X\X^T)} - 2\tr{(\X\w\w^T\X^T)} + \tr{(\X\w\w^T\X^T)}\right)\\ & = & \frac{1}{n}\tr{(\X\X^T)} - \frac{1}{n}\tr{(\X\w\w^T\X^T)} \end{eqnarray} so maximize $\frac{1}{n}\tr{(\X\w\w^T\X^T)}$ ## Backup: The gory details (cont'd.) "trace is cyclic" so \[ \tr{(\X\w\w^T\X^T)} = \tr{(\X^T\X\w\w^T)} = \tr{(\w^T\X^T\X\w)} \] so we want to maximize \[ \tr{\left(\w^T \frac{\X^T \X}{n}\w\right)} \] under the constraint \[ \w^T \w = \mathbf{I} \] This is the same form we saw before, so it has the same sort of solution: each column of $\w$ must be an eigenvector of $\V$. ## Backup: The Lagrange multiplier trick {#lagrange} - Want to solve \[ \max_{w}{L(w)} \] with **constraint** $f(w) = c$ - Option I: Use $f(w)=c$ to eliminate one coordinate of $w$ + Then $w=g(v,c)$ for some function $g$ and _unconstrained_ $v$ + Do $\max_{v}{L(g(v,c))}$ + An unconstrained problem with one less variable + Drawback: algebra is hard! + Drawback (more general): who wants to study operations research? - Option II: Solve an unconstrained problem with _more_ variables + Drawback: Only a French mathematician could think this made sense ## Backup: The Lagrange multiplier trick (cont'd.) - Define a **Lagrangian** \[ \mathcal{L}(w,\lambda) \equiv L(w) - \lambda(f(w) - c) \] + $=L(w)$ when the constraint holds - $\lambda$ is the **Lagrange multiplier** which **enforces** the constraint $f(w)=c$ - Now do an unconstrained optimization over $w$ and $\lambda$: \[ \max_{w, \lambda}{\mathcal{L}(w,\lambda)} \] ## Backup: The Lagrange multiplier trick (cont'd.) \[ \max_{w, \lambda}{\mathcal{L}(w,\lambda)} \] - Take derivatives: \begin{eqnarray} \frac{\partial \mathcal{L}}{\partial \lambda} & = & -(f(w)-c)\\ \frac{\partial \mathcal{L}}{\partial w} & = & \frac{\partial L}{\partial w} - \lambda\frac{\partial f}{\partial w} \end{eqnarray} - Set to 0 at the maximum: \begin{eqnarray} f(w) & = & c\\ \frac{\partial L}{\partial w} & = & \lambda\frac{\partial f}{\partial w} \end{eqnarray} + We've automatically recovered the constraint! - One equation per unknown $\Rightarrow$ Solve for $\lambda$, $w$ together ## Backup: The Lagrange multiplier trick (cont'd.) - More equality constraints $\Rightarrow$ more Lagrange multipliers - Inequality constraints, $g(w) \leq d$, are trickier - Is the unconstrained maximum inside the **feasible set**? + Yes: problem solved + No: constrained maximum is on the boundary + Boundary is _usually_ $g(w) = d$ $\Rightarrow$ treat like an equality + There are subtleties; sometimes we need to learn some operations research ## Backup: The Lagrange multiplier trick (cont'd.) - To think through: Why not to do this instead? \[ \max_{\lambda, w}{L(w) + \lambda(f(w)-c))^2} \] ## Backup: Some more maps from @Cavalli-Sforza-et-al-demic-expansions {#more-maps} ![](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 ## Backup: 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 ## Backup: Principal Components Regression {.smaller #pcreg} - 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] ## Backup: Orthogonal matrices {#orthogonal-matrix} - 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? ## Backup: Exchanging locations and variables {.smaller #exchanging-rows-and-columns} Recall the states: ```{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 ## Backup: Exchanging locations and 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") ``` ## Backup: Exchanging locations and variables {.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) ``` ## Backup: Exchanging locations and variables 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) ## Backup: Exchanging locations and variables = 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") ``` ## Backup: 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 {#alternatives} - **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