# How to perform dimensionality reduction with PCA in R

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?

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(res$scale != FALSE){
trunc <- scale(trunc, center = FALSE , scale=1/res$scale) } if(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, RAN,,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. • Marc, you don't need to record center and scale explicitly, prcomp does that for you. Have a look at res$center and res$scale. IMHO it is less error prone to use these (no accidental difference about centering or not / scaling or not between the explicit call to scale and the prcomp call). – cbeleites unhappy with SX Apr 28 '13 at 9:58 • This answer needs expansion because it does not answer the OP's questions about what to do after calculating the PCA or how do I reduce my dataset exactly? Given that the OP has conducted PCA on his sample, his question is what to do with it, and what actually goes on with this subsamples; not how to do PCA. We might as well propose doing E <- eigen(cov(Sample)); A<- scale(scale=F, Sample) %*% E$vectors to get yet another way to get the scores also (which is actually what princomp does stats:::princomp.default). – usεr11852 Apr 28 '13 at 14:12
• @user11852 - the question specifically makes reference to reduction of the data set (i.e. the truncation that I have demonstrated here). I'll leave it up to him to decide whether or not this was what he was looking for. – Marc in the box Apr 28 '13 at 20:24
• @Marc, thanks for the response. I think I might need to step back and re-read everything again, because I am stuck on how any of the answer above deals with dimensionality reduction. Because as you show, dim(trunc) = dim(Xt). What was the benefit of it, the dimenions did not get reduced. – B_Miner Jul 18 '13 at 18:04
• @B_Miner - Just keep in mind that truncation is used to focus on the main patterns in data and to filter out small scale patterns and noise. The truncated data is not smaller in terms of it's dimensions, but "cleaner". However, truncation does reduce the amount of data in that the entire matrix can be reconstructed with just a few vectors. A nice example is in the use of PCA for image compression, where a smaller number of PCs can be used to reconstruct the image. This smaller subset of vectors takes up less memory, but the reconstruction will have some loss in small scale detail. – Marc in the box Jul 18 '13 at 19:51

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!

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").

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. • As you mention Principal Component Regression, in R that is provided by package pls. As for the number of components to retain, I don't see any real advantage of deciding % variance over no. of components (maybe that is because I work with data that have very different levels of noise. As @Marc-in-the-box mentions, there are lots of different approaches to determine an appropriate no. of PCs, and the strategy will (should) depend both on the type data and on the type of data analysis that is to follow. – cbeleites unhappy with SX Apr 28 '13 at 10:04 • I said naive way; no scree plots, no average eigenvalue, no log-eigenvalue diagram, no partial correlation test. No probabilistic framework for your model. Because it is not what the OP asks. (I didn't mention the pls package by the way, princomp {stats}) There is a clear advantage of using % of cum. variance over a single arbitrary number: You know the quality of your reconstruction. If you just want to do dim. reduction to$K$dimensions ($K$<$D$,$D\$ being your original sample dimensionality) that's fine too but you don't really know the quality of your dimensionality reduction. – usεr11852 Apr 28 '13 at 13:28