R:如何将 ggbiplot 与 pcaRes 对象一起使用?绘制具有缺失值的数据的 PCA 结果

R: How to use ggbiplot with pcaRes object? plot PCA results of data with missing values

我通常使用 prcomp 函数执行主成分分析,并使用 ggbiplot 以奇特的方式绘制结果(或者仅使用 ggplot2 提取 pca.obj$x ).

像这样:

#install_github("vqv/ggbiplot")
library(ggbiplot)
data(iris)
pca.obj <- prcomp(iris[,1:4], center=TRUE, scale.=TRUE)
P <- ggbiplot(pca.obj,
         obs.scale = 1, 
         var.scale=1,
         ellipse=T,
         circle=F,
         varname.size=3,
         var.axes=T,
         groups=iris$Species, #no need for coloring, I'm making the points invisible
         alpha=0) #invisible points, I add them below
P$layers <- c(geom_point(aes(color=iris$Species), cex=5), P$layers) #add geom_point in a layer underneath (only way I have to change the size of the points in ggbiplot)
png(filename="test.png", height=600, width=600)
print(#or ggsave()
    P
)
dev.off()

然而,现在我面对的数据有一定数量的 NA,我正在使用 pcaMethods 包中的 pca 包装函数,应用 nipals 方法(迭代能够处理少量缺失值的方法)。

pcareturns一个对象classpcaResggbiplotreturns出现如下错误:

#introduce NAs
iris$Sepal.Length[sample(1:150, 5)] <- NA
iris$Sepal.Width[sample(1:150, 5)] <- NA
iris$Petal.Length[sample(1:150, 5)] <- NA
iris$Petal.Width[sample(1:150, 5)] <- NA
#pca.obj2 <- prcomp(iris[,1:4], center=TRUE, scale.=TRUE) #cannot use prcomp with NAs
#source("https://bioconductor.org/biocLite.R")
#biocLite("pcaMethods")
library(pcaMethods)
pca.obj2 <- pca(iris[,1:4], method="nipals", nPcs=3, center=TRUE, scale.=TRUE)
class(pca.obj2)
ggbiplot(pca.obj2)

Error in ggbiplot(pca.obj2) : Expected a object of class prcomp, princomp, PCA, or lda

我的问题是:

如何将 ggbiplot 应用于 pcaRes 对象?

如何将此对象转换为 prcomp 对象?

我可以用另一个函数而不是接受 pcaRes 对象的 ggbiplot 获得相同类型的图吗?

我是否应该只用变量的平均值替换 NA 值并像往常一样应用 prcomp 函数?

非常感谢!

首先,找到一个可以处理 NA 的 PCA 包真是太棒了。

由于ggbiplot不会接受pcaRes对象,我们可以利用pcaRes得到的数据潜入原来的prcomp对象中

显然您的真实数据已经包含 NA 值,因此我们将从该数据集开始并将它们换成一些虚拟值,以便我们 运行 第一个 prcomp pca.

iris_na<-iris

iris_na$Sepal.Length[sample(1:150, 5)] <- NA
iris_na$Sepal.Width[sample(1:150, 5)] <- NA
iris_na$Petal.Length[sample(1:150, 5)] <- NA
iris_na$Petal.Width[sample(1:150, 5)] <- NA

iris_dummy<-iris_na

iris_dummy[is.na(iris_dummy)]<-7777 #swap out your NAs with a dummy number so prcomp will run

然后我们 运行 第一个 pca 和你一样:

pca.obj <- prcomp(iris_dummy[,1:4], center=TRUE, scale.=TRUE)

这个对象有 5 个组成部分,x(分数),rotation(载荷),sdev(标准偏差),centerscale.虽然我怀疑 ggbiplot 只有分数和负载被使用,但我们会把它们全部换掉以确保安全。

查看分数组件 pca.obj$x 表明 prcomp 函数中计算了四个主要组件。

head(pca.obj$x)

#           PC1        PC2         PC3         PC4
#[1,] -2.656740  0.3176722  0.03763067 -0.04122948
#[2,] -2.688275 -0.1821744  0.19912795  0.07297624
#[3,] -2.862673 -0.1447518 -0.02134749 -0.02462359
#[4,] -2.718294 -0.3189371 -0.03318459 -0.11675762
#[5,] -2.700864  0.3274887 -0.07503096 -0.11347939
#[6,] -2.252918  0.7436711 -0.14611455 -0.08218007

因此,当我们 运行 具有 pcaRes 的下一个 pca 时,我们确保指定使用 nPcs 参数计算 4 个主成分。这里我们使用的是真实数据,其中包含 NAs.

pca.obj2 <- pca(iris_na[,1:4], method="nipals", nPcs=4, center=TRUE, scale.=TRUE)

然后只需将 pcaRes 值换成 prcomp 值并将其传递给 ggbiplot

pca.obj$x<-pca.obj2@scores 

pca.obj$rotation<-pca.obj2@loadings 

pca.obj$sdev<-pca.obj2@sDev

pca.obj$center<-pca.obj2@center

pca.obj$scale<-pca.obj2@scale

P2 <- ggbiplot(pca.obj,
              obs.scale = 1, 
              var.scale=1,
              ellipse=T,
              circle=F,
              varname.size=3,
              var.axes=T,
              groups=iris$Species, 
              alpha=0) 
P2$layers <- c(geom_point(aes(color=iris$Species), cex=5), P2$layers)