2

After performing PCA I would like to project any new samples to the principal component space (I would like to see how samples cluster together). I did the PCA analysis in R:

pca <- eigen(cov(A))

I looked at the predict function and also tried new sample %*% eigenvector after scaling the new object on to the pc space. But I get an error message saying dimensions are not matching when I do the new sample %*% eigenvector cross product. Please see the dummy example below.

I did my pca analysis on the matrix A and obtained my eigenvectors:

A = matrix(c(0,1,1,1,2,0,1,1,2,1,1,0,0,2,2), nrow=5, ncol=3, byrow = TRUE)

dimnames(A) = list(c("obs1", "obs2","obs3","obs4","obs5"), 
                   c("sample1", "sample2", "sample3")
                   )

x <- (A-rowMeans(A)) / apply(A,1,sd)

pca <- eigen(cov(x))

Now I have a new set of samples that I would like to project:

B = matrix(c(0,0,1,2,2,1,0,1,2,1), nrow=5, ncol=2, byrow = TRUE) 

dimnames(B) = list(c("obs1", "obs2","obs3","obs4","obs5"),
                   c("sample4", "sample5")
                   )

I centered my new matrix according to the pca space:

y <- (B-rowMeans(A)) / apply(A,1,sd)

And then tried to multiply it with the eigenvectors:

y %*% pca$vectors
Error in y %*% pca$vectors : non-conformable arguments

Since my eigenvector is a 3x3 matrix (as I have 3 samples in matrix A and 2 samples in matrix B - with 5 same observations) I can not multiply this with a 5x2 matrix.

Any suggestions what I am doing wrong?

Karolis Koncevičius
  • 4,282
  • 7
  • 30
  • 47
Something like that
  • 121
  • 2
  • This post appears to make the same error as https://stats.stackexchange.com/questions/447714/confusion-about-covariance-matrix -- TL;DR, there are two conventions of PCA, and it's important to know which one your software is using. – Sycorax Feb 04 '20 at 14:14

1 Answers1

1

The code you wrote seems to treat samples as features. Consider this part:

> cov(x)
            sample1     sample2    sample3
sample1  0.56666667  0.08213672 -0.6488034
sample2  0.08213672  0.35119661 -0.4333333
sample3 -0.64880339 -0.43333333  1.0821367

As you can see the covariance matrix is constructed for samples, not features. Following this your pca$vectors will have 3 dimensions. But the samples from B you are trying to project have 5:

> pca$vectors
           [,1]       [,2]      [,3]
[1,] -0.5064330  0.6404626 0.5773503
[2,] -0.3014404 -0.7588151 0.5773503
[3,]  0.8078733  0.1183525 0.5773503

> B
     sample4 sample5
obs1       0       0
obs2       1       2
obs3       2       1
obs4       0       1
obs5       2       1

To fix the issue the covariance matrix should be computed on a transpose of x.

pca <- eigen(cov(t(x)))

And the obtained projections:

> t(y) %*% pca$vectors
              [,1]        [,2]       [,3]      [,4]        [,5]
sample4 -0.9933442 -0.78697392 -1.6369351 0.1304668 -0.17519903
sample5  1.5496399 -0.01146507 -0.8184675 0.0652334 -0.08759952

Which are the same with predict:

> predict(prcomp(t(x)), t(y))
               PC1        PC2       PC3
sample4  0.9933442 0.78697392 0.9869463
sample5 -1.5496399 0.01146507 0.4934732

NOTE 1: The directions of PCA are arbitrary so the signs can be inverted

NOTE 2: The "eigen" version has more PCs, but the last 3 are degenerate. This is also why the PC3 can be different. The highest number of PCs you get from a data matrix is the minimum of #columns and #rows minus 1. So in this case 2. You can also check this by looking at pca$values, and note that the eigenvalues for the last 3 vectors are 0:

> all.equal(pca$values[1], 0)
[1] "Mean relative difference: 1"
> all.equal(pca$values[2], 0)
[1] "Mean relative difference: 1"
> all.equal(pca$values[3], 0)
[1] TRUE
> all.equal(pca$values[4], 0)
[1] TRUE
> all.equal(pca$values[5], 0)
[1] TRUE
Karolis Koncevičius
  • 4,282
  • 7
  • 30
  • 47
  • By taking the transpose of x, `eigen(cov(t(x)))` are you not clustering the observations rather than the samples? – Something like that Feb 05 '20 at 20:38
  • @Somethinglikethat `cov` function computes the covariance matrix. The purpose of this matrix is not to use it for clustering (as would be the case for hierarchical clustering), but to describe relationships between observations. Out of this matrix then one selects the projections (like principal components). The projection is done in the space of observations. After that you project the samples on these new axes to get the scores for each sample. To make sure this is the case - please run `prcomp()` in R, but note that it expects samples to be in the rows of a matrix, not columns. – Karolis Koncevičius Feb 05 '20 at 21:00