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
方法(迭代能够处理少量缺失值的方法)。
pca
returns一个对象classpcaRes
,ggbiplot
returns出现如下错误:
#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
(标准偏差),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)
我通常使用 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
returns一个对象classpcaRes
,ggbiplot
returns出现如下错误:
#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
(标准偏差),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)