主成分方差解释率计算

通常，求得了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))))

https://www.jianshu.com/p/39d22980dd61

巴特西

PCA方差解释比例求解与绘图？

主成分方差解释率计算

绘图示例

最新文章

热门文章