Predict only the first N principal components in a PCA analysis

Question

I'm using R to analyze a very large dataset. I conduct a PCA on one dataset,

PCA <- prcomp(formula = ~., data = train, scale = T, na.action=na.exclude)

and then want to apply the PCA on another dataset,

test_rot <- data.frame(predict(PCA,test,na.action=na.omit))

This works however it is fairly memory intensive and I'm really only interested in the first N principal components (say: 50). Is there a way to only predict the first 50 PCs in that second line of code?

Answer 1

You can just use the loadings matrix, PCA$rotation, to extract a matrix with the number of components (columns) you want to use. Rather than use predict(), just multiply your test data matrix by this loading matrix to get the transformed data returned by predict.

Answer 2

I have elaborated an example (below) of how to do this manually. You could technically do something similar using prcomp by setting your tolerances high enough to retain only the dominant PCs. In the example, the projection uses only a subset of principal components loadings defined by the keep vector.

The suggestion of @cbeleites is an important one, but for prediction of PCs it's important that the centering and scaling values come from the original data set - not the new data set. This is because you are projecting your data onto the weightings provided by your original PCA.

Example

library(MASS)

################
###train data###
################
data(iris)
dat <- as.matrix(iris[,-5])

##############
###new data###
##############
#here is some new data that has the same covariances bewteen variables ar the iris data set
require(MASS)
setosa.mean <- apply(iris[iris$Species=="setosa",-5], 2, mean)
setosa.cov <- cov(iris[iris$Species=="setosa",-5])

versicolor.mean <- apply(iris[iris$Species=="versicolor",-5], 2, mean)
versicolor.cov <- cov(iris[iris$Species=="versicolor",-5])

virginica.mean <- apply(iris[iris$Species=="virginica",-5], 2, mean)
virginica.cov <- cov(iris[iris$Species=="virginica",-5])

set.seed(1)
n <- 30
newdat <- as.data.frame(
  rbind(
    mvrnorm(n, setosa.mean, setosa.cov),
    mvrnorm(n, versicolor.mean, versicolor.cov),
    mvrnorm(n, virginica.mean, virginica.cov)
  )
)
newdat <- cbind(newdat, SPP=rep(c("setosa", "versicolor", "virginica"), each=n))
head(newdat)



##############
###full pca###
##############
#scale data
dat.sc <- scale(dat, center=TRUE, scale=FALSE)

#PCA
pca <- svd(dat.sc)
pcs <- dat.sc %*% pca$v # results are similar to pca$u but includes the units contained in pca$d
plot(scale(pca$u), scale(pcs)); abline(0,1,col=2) # demonstration of the similarity of pca$u and pcs

#prediction
keep <- 1:2 # there are many ways to define which PCs you want to use for prediction (e.g. "tol" setting in prcomp)

newdat.sc <- scale(newdat[,1:4], center=attr(dat.sc, "scaled:center"), scale=FALSE) # you must use the variable centering and scaling from the original data
pred <- newdat.sc %*% pca$v[,keep]


###Plot PCs 1 & 2 for orig. and new data
COL <- 2:4
plot(pcs, col=COL[iris$Species], cex=1, xlab=paste0("PC 1", " (", round(pca$d[1]^2/sum(pca$d^2)*100,0), "%)"), ylab=paste0("PC 2", " (", round(pca$d[2]^2/sum(pca$d^2)*100,0), "%)"))
points(pred, col=COL[match(newdat$SPP, levels(iris$Species))], pch=16)
legend("topright", legend=levels(iris$Species), col=COL, pch=16)
legend("topleft", legend=c("Original data", "New data"), col=1, pch=c(1,16))

Here are the predictions using prcomp for comparison

#compare to predict.prcomp
pca2 <- prcomp(dat, center=TRUE, scale=FALSE)
pred2 <- predict(pca2, as.data.frame(newdat))

plot(pred[,1], pred2[,1]) # 1st PC
points(pred[,2], pred2[,2], col=2) # 2nd PC
abline(0,1)