主成分方差解释率计算
通常,求得了PCA降维后的特征值,我们就可以绘图,但各个维度的方差解释率没有得到,就无法获得PC坐标的百分比。
有些工具的结果是提供了维度标准差的,如ggbiplot绘图时,直接会给你算出各个坐标的方差解释率。但我觉得这类工具绘图远不如ggplot本身,此时,就需要自己计算。
当理解了PCA的原理和含义后,就比较容易得到。网上一大堆,这里不介绍。
以ggbiplot数据为例,并将算出结果与之比较。
if(!require(devtools))
install.packages("devtools")
library(devtools)
if(!require(ggbiplot))
install_github("vqv/ggbiplot")
library(ggbiplot)
data(wine)
pca <- prcomp(wine, scale. = TRUE)
ggbiplot(pca,
# groups = wine.class,
ellipse = TRUE, circle = TRUE,
obs.scale = 1, var.scale = 1) +
scale_color_discrete(name = '') +
theme(legend.direction = 'horizontal', legend.position = 'top')
R自带函数prcomp的结果中得到PCA的5个对象结果,其中包含了标准差(sdev)和特征向量(x)。
> str(pca)
List of 5
$ sdev : num [1:13] 2.169 1.58 1.203 0.959 0.924 ...
$ rotation: num [1:13, 1:13] -0.14433 0.24519 0.00205 0.23932 -0.14199 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:13] "Alcohol" "MalicAcid" "Ash" "AlcAsh" ...
.. ..$ : chr [1:13] "PC1" "PC2" "PC3" "PC4" ...
$ center : Named num [1:13] 13 2.34 2.37 19.49 99.74 ...
..- attr(*, "names")= chr [1:13] "Alcohol" "MalicAcid" "Ash" "AlcAsh" ...
$ scale : Named num [1:13] 0.812 1.117 0.274 3.34 14.282 ...
..- attr(*, "names")= chr [1:13] "Alcohol" "MalicAcid" "Ash" "AlcAsh" ...
$ x : num [1:178, 1:13] -3.31 -2.2 -2.51 -3.75 -1.01 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:13] "PC1" "PC2" "PC3" "PC4" ...
- attr(*, "class")= chr "prcomp"
手动计算方差解释率:
> pca$sdev^2/sum(pca$sdev^2)*100
#注意平方
[1] 36.1988481 19.2074903 11.1236305 7.0690302 6.5632937 4.9358233
[7] 4.2386793 2.6807489 2.2221534 1.9300191 1.7368357 1.2982326
[13] 0.7952149
可看出,前两个主成分与图中一致。当然如果没有标准差结果,我们也可以根据特征向量计算出来:
> sdev<- apply(pca$x,2,sd)
> sdev
PC1 PC2 PC3 PC4 PC5 PC6 PC7
2.1692972 1.5801816 1.2025273 0.9586313 0.9237035 0.8010350 0.7423128
PC8 PC9 PC10 PC11 PC12 PC13
0.5903367 0.5374755 0.5009017 0.4751722 0.4108165 0.3215244
绘图示例
一个示例,可在此基础上进一步优化。如样本要再分组,可加shape。
ggplot(data=data.frame(pca$x), aes(PC1,PC2)) +
stat_ellipse(aes(fill=wine.class),type="norm",geom="polygon",alpha=0.2,color=NA)+
geom_point(size=2)+
# scale_size(guide=FALSE)+
scale_color_manual(values = col)+
geom_vline(xintercept = 0,linetype="dotted")+
geom_hline(yintercept = 0,linetype="dotted")+
labs(x=paste0("PC1", sprintf("(%0.2f%%)",100*pca$sdev[1]^2/sum(pca$sdev^2))),
y=paste0("PC2", sprintf("(%0.2f%%)",100*pca$sdev[2]^2/sum(pca$sdev^2))))