How to perform dimensionality reduction with PCA in R

Question

I have a big dataset and I want to perform a dimensionality reduction.

Now everywhere I read that I can use PCA for this. However, I still don't seem to get what to do after calculating/performing the PCA. In R this is easily done with the command princomp.

But what to do after calculating the PCA? If I decided I want to use the first $100$ principal components, how do I reduce my dataset exactly?

Answer 1

I believe what you are getting at in your question concerns data truncation using a smaller number of principal components (PC). For such operations, I think the function prcompis more illustrative in that it is easier to visualize the matrix multiplication used in reconstruction.

First, give a synthetic dataset, Xt, you perform the PCA (typically you would center samples in order to describe PC's relating to a covariance matrix:

#Generate data
m=50
n=100
frac.gaps <- 0.5 # the fraction of data with NaNs
N.S.ratio <- 0.25 # the Noise to Signal ratio for adding noise to data

x <- (seq(m)*2*pi)/m
t <- (seq(n)*2*pi)/n

#True field
Xt <- 
 outer(sin(x), sin(t)) + 
 outer(sin(2.1*x), sin(2.1*t)) + 
 outer(sin(3.1*x), sin(3.1*t)) +
 outer(tanh(x), cos(t)) + 
 outer(tanh(2*x), cos(2.1*t)) + 
 outer(tanh(4*x), cos(0.1*t)) + 
 outer(tanh(2.4*x), cos(1.1*t)) + 
 tanh(outer(x, t, FUN="+")) + 
 tanh(outer(x, 2*t, FUN="+"))

Xt <- t(Xt)

#PCA
res <- prcomp(Xt, center = TRUE, scale = FALSE)
names(res)

In the results or prcomp, you can see the PC's (res$x), the eigenvalues (res$sdev) giving information on the magnitude of each PC, and the loadings (res$rotation).

res$sdev
length(res$sdev)
res$rotation
dim(res$rotation)
res$x
dim(res$x)

By squaring the eigenvalues, you get the variance explained by each PC:

plot(cumsum(res$sdev^2/sum(res$sdev^2))) #cumulative explained variance

Finally, you can create a truncated version of your data by using only the leading (important) PCs:

pc.use <- 3 # explains 93% of variance
trunc <- res$x[,1:pc.use] %*% t(res$rotation[,1:pc.use])

#and add the center (and re-scale) back to data
if(all(res$scale != FALSE)){
	trunc <- scale(trunc, center = FALSE , scale=1/res$scale)
}
if(all(res$center != FALSE)){
    trunc <- scale(trunc, center = -1 * res$center, scale=FALSE)
}
dim(trunc); dim(Xt)

You can see that the result is a slightly smoother data matrix, with small scale features filtered out:

RAN <- range(cbind(Xt, trunc))
BREAKS <- seq(RAN[1], RAN[2],,100)
COLS <- rainbow(length(BREAKS)-1)
par(mfcol=c(1,2), mar=c(1,1,2,1))
image(Xt, main="Original matrix", xlab="", ylab="", xaxt="n", yaxt="n", breaks=BREAKS, col=COLS)
box()
image(trunc, main="Truncated matrix (3 PCs)", xlab="", ylab="", xaxt="n", yaxt="n", breaks=BREAKS, col=COLS)
box()

And here is a very basic approach that you can do outside of the prcomp function:

#alternate approach
Xt.cen <- scale(Xt, center=TRUE, scale=FALSE)
C <- cov(Xt.cen, use="pair")
E <- svd(C)
A <- Xt.cen %*% E$u

#To remove units from principal components (A)
#function for the exponent of a matrix
"%^%" <- function(S, power)
     with(eigen(S), vectors %*% (values^power * t(vectors)))
Asc <- A %*% (diag(E$d) %^% -0.5) # scaled principal components

#Relationship between eigenvalues from both approaches
plot(res$sdev^2, E$d) #PCA via a covariance matrix - the eigenvalues now hold variance, not stdev
abline(0,1) # same results

Now, deciding which PCs to retain is a separate question - one that I was interested in a while back. Hope that helps.

Answer 2

These other answers are very good and detailed, but I'm wondering if you're actually asking a much more basic question: what do you do once you have your PCs?

Each PC simply becomes a new variable. Say PC1 accounts for 60 % of the total variation and PC2 accounts for 30 %. As that's 90 % of the total variation accounted for, you could simply take these two new variables (PCs) as a simplified version of your original variables. This means fitting them to models, if that's what you're interested in. When it comes time to interpret your results, you do so in the context of the original variables that are correlated with each PC.

Sorry if I've underestimated the scope of the question!

Answer 3

I believe your original question stems from being a bit uncertain about what PCA is doing. Principal Component Analysis allows you to identify the principal mode of variation in your sample. Those modes are emperically calculated as the eigenvectors of your sample's covariance matrix (the "loadings"). Subsequently those vectors serve as the new "coordinate system" of your sample as you project your original sample in the space they define (the "scores"). The proportion of variation associated with the $i$-th eigenvector/ mode of variation/ loading/ principal component equals $\frac{\lambda_i}{\Sigma_{k=1}^{p} \lambda_k}$ where $p$ is your sample's original dimensionality ($p =784$ in your case). [Remember because your covariance matrix is non-negative Definite you'll have no negative eigenvalues $\lambda$.] Now, by definition the eigenvectors are orthogonal to each other. That means that their respective projections are also orthogonal and where originally you had a big possibly correlated sample of variables, now you have a (hopefully significantly) smaller linearly independent sample (the "scores").

Practically with PCA you are using the projections of the PCs (the "scores") as surrogate data for your original sample. You do all your analysis on the scores, and afterwards you reconstruct your original sample back using the PCs to find out out what happened on your original space (that's basically Principal Component Regression). Clearly, if you are able to meaningful interpreter your eigenvectors ("loadings") then you are in an even better position: You can describe what happens to your sample in the mode of variation presented by that loading by doing inference on that loading directly and not care about reconstruction at all. :)

In general what do you "after calculating the PCA" depends on the target of your analysis. PCA just gives you a linearly independent sub-sample of your data that is the optimal under an RSS reconstruction criterion. You might use it for classification, or regression, or both, or as I mentioned you might want to recognise meaningful orthogonal modes of variations in your sample.

A comment : I think the best naive way to decide the number of components to retain is to base your estimate on some threshold of sample variation you would like to retain in your reduced dimensionality sample rather than just some arbitrary number eg. 3, 100, 200. As user4959 explained you can check that cumulative variation by checking the relevant field of the list under the $loadingsfield in the list object produced by princomp.

Answer 4

After doing the PCA then you may select the first two components and plot.. You can see the variation of the components using a scree plot in R. Also using summary function with loadings=T you can fins the variation of features with the components.

You can also look at this http://www.statmethods.net/advstats/factor.html and http://statmath.wu.ac.at/~hornik/QFS1/principal_component-vignette.pdf

Try to think what you want. You can interpret lots of things from PCA analysis.

Best Abhik