来自 FactoMineR::PCA() 的 PCA 未返回与 base::prcomp() 相同的结果

PCA from FactoMineR::PCA() not returning same as base::prcomp()

进行一些 PCA 分析,并将 FactoMineR 函数 PCA 的结果与 baseprcomp 的结果进行比较时,我没有得到相同的结果。一个例子

library(ISLR)
library(FactoMineR)
data("NCI60")

df <- NCI60$data


pca_prcomp <- prcomp(df, scale. = T)
pca_facto <- FactoMineR::PCA(df, scale.unit = T, graph = F, ncp = 65)


# One column is missing

dim(pca_prcomp$x)
dim(pca_facto$ind$coord) 

# Values are similiare - but not the same

head(pca_prcomp$x[, 1:2])
head(pca_facto$ind$coord[, 1:2])


# Using scale function - does not return same values

pca_facto_scale <- PCA(scale(df), scale.unit = F, graph = F, ncp = 65)

head(pca_facto$ind$coord[, 1:2], 3)
head(pca_facto_scale$ind$coord[, 1:2], 3)

执行 PCA 的特征值和 SVD 方法之间的差异很可能存在(有关详细信息,请参阅 this great answer

来自?prcomp

The calculation is done by a singular value decomposition of the (centered and possibly scaled) data matrix, not by using ‘eigen’ on the covariance matrix. This is generally the preferred method for numerical accuracy.

来自?PCA:

Returns a list including:
eig: a matrix containing all the eigenvalues, the percentage of variance and the cumulative percentage of variance

抱歉来晚了,FactoMineR 包使用与 svd() 相同的方法,应该与 prcomp() 方法相似(但不相同),并且都列出了在 Q-mode 下,由于其数值精度,这是进行 PCA 的首选方法。但是请注意,我没有说相同,为什么? FactoMineR 使用自己的 PCA 算法计算组件数量,如下所示:

ncp <- min(ncp, nrow(X) - 1, ncol(X))

这清楚地告诉您为什么您得到的组件数量是 63 而不是 prcomp() 通常给出的 64。您的数据集是典型的基因组学数据,其中 n 行小于 p 列基因,并且上述代码显然会采用列或行,以数量较少的为准。如果您遵循 svd() 算法,它将 return 64 维而不是 63。

要进一步探索源代码,请键入 FactoMineR:::PCA

Q-mode(svdprcomp()FactoMineR::PCA())和R-mode(eigen()princomp()) 我建议访问此 answer.

旁注:对于prcomp(),您需要传递center = T 参数以便在执行PCA 之前将数据居中。另一方面,缩放会给所有基因列同等的权重。

pca_prcomp <- prcomp(df, center = T, scale. = T) # add center=T

为了缩放,prcomp() 使用 N 作为除数,而 FactoMineR::PCA() 使用 N-1 作为除数。下面的代码将证明这一点(参考上面相同的链接答案):

# this is the scaled data by scale()
df_scaled <- scale(df)

# then you need to get the standardized data matrix from the output of the FactoMineR::PCR() function, which can be done easily as follows:
df_restored <- pca_facto$svd$U %*% diag(pca_facto$svd$vs) %*% t(pca_facto$svd$V)

# the to make both FactoMineR::PCR() and scale() match up you need to do the correction
df_corrected <- df_restored * sqrt(63 / 64) # correct for sqrt(N-1/N)

head(df[, 1:5]) # glimpse the first five columns only!
head(df_scaled[, 1:5])
head(df_restored[, 1:5]) # glimpse the first five columns only!
head(df_corrected[, 1:5])
round(head(df_scaled[, 1:5]), 3) == round(head(df_corrected[, 1:5]), 3) # TRUE

R> head(df[, 1:5])
       1      2      3      4      5
V1 0.300  1.180  0.550  1.140 -0.265
V2 0.680  1.290  0.170  0.380  0.465
V3 0.940 -0.040 -0.170 -0.040 -0.605
V4 0.280 -0.310  0.680 -0.810  0.625
V5 0.485 -0.465  0.395  0.905  0.200
V6 0.310 -0.030 -0.100 -0.460 -0.205
R> head(df_scaled[, 1:5])
       1        2      3      4      5
V1 0.723  1.59461  1.315  1.345 -0.600
V2 1.584  1.73979  0.438  0.649  0.905
V3 2.173 -0.01609 -0.346  0.264 -1.301
V4 0.678 -0.37256  1.615 -0.441  1.235
V5 1.142 -0.57720  0.958  1.130  0.359
V6 0.746 -0.00289 -0.185 -0.120 -0.476
R> head(df_restored[, 1:5])
      [,1]     [,2]   [,3]   [,4]   [,5]
[1,] 0.729  1.60722  1.326  1.356 -0.605
[2,] 1.596  1.75354  0.442  0.654  0.912
[3,] 2.190 -0.01622 -0.349  0.266 -1.311
[4,] 0.683 -0.37550  1.628 -0.444  1.244
[5,] 1.151 -0.58176  0.965  1.139  0.361
[6,] 0.752 -0.00291 -0.186 -0.121 -0.480
R> head(df_corrected[, 1:5])
      [,1]     [,2]   [,3]   [,4]   [,5]
[1,] 0.723  1.59461  1.315  1.345 -0.600
[2,] 1.584  1.73979  0.438  0.649  0.905
[3,] 2.173 -0.01609 -0.346  0.264 -1.301
[4,] 0.678 -0.37256  1.615 -0.441  1.235
[5,] 1.142 -0.57720  0.958  1.130  0.359
[6,] 0.746 -0.00289 -0.185 -0.120 -0.476
R> round(head(df_scaled[, 1:5]), 3) == round(head(df_corrected[, 1:5]), 3)
      1    2    3    4    5
V1 TRUE TRUE TRUE TRUE TRUE
V2 TRUE TRUE TRUE TRUE TRUE
V3 TRUE TRUE TRUE TRUE TRUE
V4 TRUE TRUE TRUE TRUE TRUE
V5 TRUE TRUE TRUE TRUE TRUE
V6 TRUE TRUE TRUE TRUE TRUE

书摘

FactoMineR 包中还有 François Husson、Sébastien 和 Lê Jérôme Pagès 编写的名为“使用 R 的示例进行探索性多元分析”第二版的书。以下是该书第 55 页的节选,该节讨论了与您类似的基因组研究的数据集,其中 n 行 (43) 远少于 p 7407 列 chicken.csv 数据集, 你可以在他们的网站上看到更多信息,数据集本身可以从这个 link.

下载