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

Question

我通常使用 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 方法(一种能够处理少量缺失值的迭代方法)。

pca 返回类 pcaRes 的对象，ggbiplot 返回以下错误:

#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 函数？

非常感谢!

Answer 1

首先，找到一个可以处理 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 (标准偏差)，center和 scale .虽然我怀疑 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)