PCA Interpretation – Understanding Results from Principal Component Analysis

pca

As part of a University assignment, I have to conduct data pre-processing on a fairly huge, multivariate (>10) raw data set. I'm not a statistician in any sense of the word, so I'm a little confused as to what's going on. Apologies in advance for what is probably a laughably simple question – my head's spinning after looking at various answers and trying to wade through the stats-speak.

I've read that:

PCA allows me to reduce the dimensionality of my data
It does so by merging / removing attributes / dimensions that correlate a lot (and thus are a little unnecessary)
It does so by finding eigenvectors on covariance data (thanks to a nice tutorial I followed through to learn this)

Which is great.

However, I'm really struggling to see how I can apply this practically to my data. For instance (this is not the data set I'll be using, but an attempt at a decent example people can work with), if I were to have a data set with something like…

PersonID     Sex     Age Range    Hours Studied     Hours Spent on TV      Test Score     Coursework Score 
1            1       2            5                 7                      60             75
2            1       3            8                 2                      70             85 
3            2       2            6                 6                      50             77
...          ...     ...          ...               ...                    ...            ...

I'm not quite sure how I would interpret any results.

Most of the tutorials I've seen online seem to give me a very mathematical view of PCA. I've done some research into it and followed them through – but I'm still not entirely sure what this means for me, who's just trying to extract some form of meaning from this pile of data I have in front of me.

Simply performing PCA on my data (using a stats package) spits out an NxN matrix of numbers (where N is the number of original dimensions), which is entirely greek to me.

How can I do PCA and take what I get in a way I can then put into plain english in terms of the original dimensions?

Best Answer

Pages 13-20 of the tutorial you posted provide a very intuitive geometric explanation of how PCA is used for dimensionality reduction.

The 13x13 matrix you mention is probably the "loading" or "rotation" matrix (I'm guessing your original data had 13 variables?) which can be interpreted in one of two (equivalent) ways:

The (absolute values of the) columns of your loading matrix describe how much each variable proportionally "contributes" to each component.
The rotation matrix rotates your data onto the basis defined by your rotation matrix. So if you have 2-D data and multiply your data by your rotation matrix, your new X-axis will be the first principal component and the new Y-axis will be the second principal component.

EDIT: This question gets asked a lot, so I'm just going to lay out a detailed visual explanation of what is going on when we use PCA for dimensionality reduction.

Consider a sample of 50 points generated from y=x + noise. The first principal component will lie along the line y=x and the second component will lie along the line y=-x, as shown below.

enter image description here

The aspect ratio messes it up a little, but take my word for it that the components are orthogonal. Applying PCA will rotate our data so the components become the x and y axes:

enter image description here

The data before the transformation are circles, the data after are crosses. In this particular example, the data wasn't rotated so much as it was flipped across the line y=-2x, but we could have just as easily inverted the y-axis to make this truly a rotation without loss of generality as described here.

The bulk of the variance, i.e. the information in the data, is spread along the first principal component (which is represented by the x-axis after we have transformed the data). There's a little variance along the second component (now the y-axis), but we can drop this component entirely without significant loss of information. So to collapse this from two dimensions into 1, we let the projection of the data onto the first principal component completely describe our data.

enter image description here

We can partially recover our original data by rotating (ok, projecting) it back onto the original axes.

enter image description here

The dark blue points are the "recovered" data, whereas the empty points are the original data. As you can see, we have lost some of the information from the original data, specifically the variance in the direction of the second principal component. But for many purposes, this compressed description (using the projection along the first principal component) may suit our needs.

Here's the code I used to generate this example in case you want to replicate it yourself. If you reduce the variance of the noise component on the second line, the amount of data lost by the PCA transformation will decrease as well because the data will converge onto the first principal component:

set.seed(123)
y2 = x + rnorm(n,0,.2)
mydata = cbind(x,y2)
m2 = colMeans(mydata)

p2 = prcomp(mydata, center=F, scale=F)
reduced2= cbind(p2$x[,1], rep(0, nrow(p2$x)))
recovered = reduced2 %*% p2$rotation

plot(mydata, xlim=c(-1.5,1.5), ylim=c(-1.5,1.5), main='Data with principal component vectors')
arrows(x0=m2[1], y0=m2[2]
       ,x1=m2[1]+abs(p2$rotation[1,1])
       ,y1=m2[2]+abs(p2$rotation[2,1])
       , col='red')
arrows(x0=m2[1], y0=m2[2]
       ,x1=m2[1]+p2$rotation[1,2]
       ,y1=m2[2]+p2$rotation[2,2]
       , col='blue')

plot(mydata, xlim=c(-1.5,1.5), ylim=c(-1.5,1.5), main='Data after PCA transformation')
points(p2$x, col='black', pch=3)
arrows(x0=m2[1], y0=m2[2]
       ,x1=m2[1]+abs(p2$rotation[1,1])
       ,y1=m2[2]+abs(p2$rotation[2,1])
       , col='red')
arrows(x0=m2[1], y0=m2[2]
       ,x1=m2[1]+p2$rotation[1,2]
       ,y1=m2[2]+p2$rotation[2,2]
       , col='blue')
arrows(x0=mean(p2$x[,1])
      ,y0=0
      ,x1=mean(p2$x[,1])
      ,y1=1
      ,col='blue'
       )
arrows(x0=mean(p2$x[,1])
       ,y0=0
       ,x1=-1.5
       ,y1=0
       ,col='red'
)
lines(x=c(-1,1), y=c(2,-2), lty=2)


plot(p2$x, xlim=c(-1.5,1.5), ylim=c(-1.5,1.5), main='PCA dimensionality reduction')
points(reduced2, pch=20, col="blue")
for(i in 1:n){
  lines(rbind(reduced2[i,], p2$x[i,]), col='blue')
}

plot(mydata, xlim=c(-1.5,1.5), ylim=c(-1.5,1.5), main='Lossy data recovery after PCA transformation')
arrows(x0=m2[1], y0=m2[2]
       ,x1=m2[1]+abs(p2$rotation[1,1])
       ,y1=m2[2]+abs(p2$rotation[2,1])
       , col='red')
arrows(x0=m2[1], y0=m2[2]
       ,x1=m2[1]+p2$rotation[1,2]
       ,y1=m2[2]+p2$rotation[2,2]
       , col='blue')
for(i in 1:n){
  lines(rbind(recovered[i,], mydata[i,]), col='blue')
}
points(recovered, col='blue', pch=20)

Related Solutions

Solved – PCA and degrees of freedom

I think the point that I was missing is that PCA does not tell us about the bare number of degrees of freedom, but rather the variance-weighted number of degrees of freedom.

In my first example all ten variables had equal variance (barring sampling noise), and so the PCA shows each successive principal component adding a roughly equal amount of variance.

In the second example the variables I'm working with have a wide range of variances: $y_1$ has a much larger variance than $y_{10}$. The PCA therefore correctly tells me that fewer principal components are required to capture most of the variance, because the higher indexed $y$'s contribute increasingly little to the total variance of the system.

Solved – Step by step implementation of PCA in R using Lindsay Smith’s tutorial

You were very very nearly there and got caught by a subtle issue in working with matrices in R. I worked through from your final_data and got the correct results independently. Then I had a closer look at your code. To cut a long story short, where you wrote

row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))

you would have been ok if you had written

row_orig_data = t(t(feat_vec) %*% t(trans_data))

instead (because you'd zeroed out the part of trans_data that was projected on the second eigenvector). As it was you were trying to multiply a $2\times1$ matrix by a $2\times10$ matrix but R did not give you an error. The problem is that t(feat_vec[1,]) is treated as $1\times2$. Trying row_orig_data = t(as.matrix(feat_vec[1,],ncol=1,nrow=2) %*% t(trans_data)) would have given you a non-conformable arguments error. The following, possibly more along the lines of what you intended, would also have worked

row_orig_data = t(as.matrix(feat_vec[1,],ncol=1,nrow=2) %*% t(trans_data)[1,])

as it multiplies a $2\times1$ matrix by a $1\times10$ matrix (note that you could have used the original final_data matrix here). It isn't necessary to do it this way, but it's nicer mathematically because it shows that you are getting $20=2\times10$ values in row_orig_data from $12=2\times1 + 1\times10$ values on the right hand side.

I've left my original answer below, as someone might find it useful, and it does demonstrate getting the required plots. It also shows that the code can be a bit simpler by getting rid of some unnecessary transposes: $(XY)^T=Y^TX^T$ so t(t(p) %*% t(q)) = q %*% t.

Re your edit, I've added the principal component line in green to my plot below. In your question you got had the slope as $x/y$ not $y/x$.

Write

d_in_new_basis = as.matrix(final_data)

then to get your data back in its original basis you need

d_in_original_basis = d_in_new_basis %*% feat_vec

You can zero out the parts of your data that are projected along the second component using

d_in_new_basis_approx = d_in_new_basis
d_in_new_basis_approx[,2] = 0

and you can then transform as before

d_in_original_basis_approx = d_in_new_basis_approx %*% feat_vec

Plotting these on the same plot, together with the principal component line in green, shows you how the approximation worked.

plot(x=d_in_original_basis[,1]+mean(d$x),
     y=d_in_original_basis[,2]+mean(d$y),
     pch=16, xlab="x", ylab="y", xlim=c(0,3.5),ylim=c(0,3.5),
     main="black=original data\nred=original data restored using only a single eigenvector")
points(x=d_in_original_basis_approx[,1]+mean(d$x),
       y=d_in_original_basis_approx[,2]+mean(d$y),
       pch=16,col="red")
points(x=c(mean(d$x)-e$vectors[1,1]*10,mean(d$x)+e$vectors[1,1]*10), c(y=mean(d$y)-e$vectors[2,1]*10,mean(d$y)+e$vectors[2,1]*10), type="l",col="green")

enter image description here

Let's rewind to what you had. This line was ok

final_data = data.frame(t(feat_vec %*% row_data_adj))

The crucial bit here is feat_vec %*% row_data_adj which is equivalent to $Y=S^TX$ where $S$ is the matrix of eigenvectors and $X$ is your data matrix with your data in rows, and $Y$ is the data in the new basis. What this is saying is that the first row of $Y$ is the sum of (rows of $X$ weighted by the first eigenvector). And the second row of $Y$ is the sum of (rows of $X$ weighted by the second eigenvector).

Then you had

trans_data = final_data
trans_data[,2] = 0

This is ok: you're just zeroing out the parts of your data that are projected along the second component. Where it goes wrong is

row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))

Writing $\hat Y$ for the matrix of data $Y$ in the new basis, with zeros in the second row, and writing $\mathbf{e}_1$ for the first eigenvector, the business end of this code t(feat_vec[1,]) %*% t(trans_data) comes down to $\mathbf{e}_1 \hat Y$.

As explained above (this is where I realised the subtle R problem and wrote the first part of my answer), mathematically you are trying to multiply a $2\times1$ vector by a $2\times10$ matrix. This doesn't work mathematically. What you should do is take the first row of $\hat Y$ = the first row of $Y$: call this $\mathbf{y}_1$. Then multiply $\mathbf{e}_1$ and $\mathbf{y}_1$ together. The $i$th column of the result $\mathbf{e}_1\mathbf{y}_1$ is the eigenvector $\mathbf{e}_1$ weighted by the 1st coordinate only of the $i$th point in the new basis, which is what you want.

Best Answer

Related Solutions

Solved – PCA and degrees of freedom

Solved – Step by step implementation of PCA in R using Lindsay Smith’s tutorial

Related Question