前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >PCA分析的方方面面

PCA分析的方方面面

作者头像
邓飞
发布2021-03-30 10:55:48
1.6K0
发布2021-03-30 10:55:48
举报

今天做PCA分析时,发现了这个ggord包,绘图很给力。利用官方的文档,学习一下。这其中,我又将PCA相关的分析方法和作图汇总了一下:

分析函数

  • eigen
  • prcomp
  • princomp
  • FactoMineR::PCA
  • ade4::dudi.pca

还有其它分析,比如冗余分析,MCA,ACM,DPCoA,metaMDS,对应分析等等。

另外,PCA不同风格的作图,应有尽有,内容丰富,干货满满!

1. 使用iris数据

代码语言:javascript
复制
library(ggplot2)
library(tidyverse)
代码语言:javascript
复制
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
代码语言:javascript
复制
## √ tibble  3.0.4     √ dplyr   1.0.4
## √ tidyr   1.1.2     √ stringr 1.4.0
## √ readr   1.4.0     √ forcats 0.5.1
## √ purrr   0.3.4
代码语言:javascript
复制
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
代码语言:javascript
复制
data(iris)
head(iris)
代码语言:javascript
复制
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

使用eigen手动分析PCA:

代码语言:javascript
复制
cor_mat = cor(iris[,1:4])
re0 = eigen(cor_mat)

特征值

代码语言:javascript
复制
re0$values
代码语言:javascript
复制
## [1] 2.91849782 0.91403047 0.14675688 0.02071484

特征向量

代码语言:javascript
复制
re0$vectors
代码语言:javascript
复制
##            [,1]        [,2]       [,3]       [,4]
## [1,]  0.5210659 -0.37741762  0.7195664  0.2612863
## [2,] -0.2693474 -0.92329566 -0.2443818 -0.1235096
## [3,]  0.5804131 -0.02449161 -0.1421264 -0.8014492
## [4,]  0.5648565 -0.06694199 -0.6342727  0.5235971

PCA得分

代码语言:javascript
复制
score1 = scale(iris[,1:4]) %*% re0$vectors %>% as.data.frame()
score1$Species = iris$Species
head(score1)
代码语言:javascript
复制
##          V1         V2          V3           V4 Species
## 1 -2.257141 -0.4784238  0.12727962  0.024087508  setosa
## 2 -2.074013  0.6718827  0.23382552  0.102662845  setosa
## 3 -2.356335  0.3407664 -0.04405390  0.028282305  setosa
## 4 -2.291707  0.5953999 -0.09098530 -0.065735340  setosa
## 5 -2.381863 -0.6446757 -0.01568565 -0.035802870  setosa
## 6 -2.068701 -1.4842053 -0.02687825  0.006586116  setosa

2. 使用prcomp做聚类分析

这里,对数据进行标准化之后,在进行PCA分析。

代码语言:javascript
复制
ord <- prcomp(iris[, 1:4],scale. = T,center = T)
summary(ord)
代码语言:javascript
复制
## Importance of components:
##                           PC1    PC2     PC3     PC4
## Standard deviation     1.7084 0.9560 0.38309 0.14393
## Proportion of Variance 0.7296 0.2285 0.03669 0.00518
## Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

3. 绘制第一张图

代码语言:javascript
复制
library(ggord)
p <- ggord(ord, iris$Species)
p

4. 设置不同风格的图形

代码语言:javascript
复制
p <- ggord(ord, iris$Species, cols = c('purple', 'orange', 'blue'))
p

不同品种,用不同的图形,圆形,十字型,三角形。

代码语言:javascript
复制
p + scale_shape_manual('Groups', values = c(1, 2, 3))
代码语言:javascript
复制
## Scale for 'shape' is already present. Adding another scale for 'shape', which
## will replace the existing scale.

设置传统的风格

代码语言:javascript
复制
p + theme_classic()

设置图例,放在上面

代码语言:javascript
复制
p + theme(legend.position = 'top')

通过vec_lab设置简写

代码语言:javascript
复制
new_lab <- list(Sepal.Length = 'SL', Sepal.Width = 'SW', Petal.Width = 'PW',
                Petal.Length = 'PL')
p <- ggord(ord, iris$Species, vec_lab = new_lab)
p

将不同的品种,分开显示

代码语言:javascript
复制
p <- ggord(ord, iris$Species, facet = TRUE, nfac = 3)
p

5. 使用princomp函数,进行PCA的分析

代码语言:javascript
复制
dd = scale(iris[,1:4])
ord <- princomp(dd)
summary(ord)
代码语言:javascript
复制
## Importance of components:
##                           Comp.1    Comp.2     Comp.3      Comp.4
## Standard deviation     1.7026571 0.9528572 0.38180950 0.143445939
## Proportion of Variance 0.7296245 0.2285076 0.03668922 0.005178709
## Cumulative Proportion  0.7296245 0.9581321 0.99482129 1.000000000

注意,

代码语言:javascript
复制
ggord(ord, iris$Species)

6. R中两个函数prcompprincomp的区别

  • prcomp函数,可以接受原始数据,在函数中定义scale .= TRUE,center = TRUE
  • princomp函数,需要使用标准化后的数据,即dd = scale(iris[,1:4],使用dd作为对象
  • 两者PCA结果是完全一致的,不过PC2的得分,正负是相反的,只是作图有区别,结果一致。

两个函数的比较

  • prcomp函数的用法:
代码语言:javascript
复制
library(tidyverse)
re1 = prcomp(iris[,1:4],center = T, scale. = T)
summary(re1)
代码语言:javascript
复制
## Importance of components:
##                           PC1    PC2     PC3     PC4
## Standard deviation     1.7084 0.9560 0.38309 0.14393
## Proportion of Variance 0.7296 0.2285 0.03669 0.00518
## Cumulative Proportion  0.7296 0.9581 0.99482 1.00000
代码语言:javascript
复制
re1$x %>% head
代码语言:javascript
复制
##            PC1        PC2         PC3          PC4
## [1,] -2.257141 -0.4784238  0.12727962  0.024087508
## [2,] -2.074013  0.6718827  0.23382552  0.102662845
## [3,] -2.356335  0.3407664 -0.04405390  0.028282305
## [4,] -2.291707  0.5953999 -0.09098530 -0.065735340
## [5,] -2.381863 -0.6446757 -0.01568565 -0.035802870
## [6,] -2.068701 -1.4842053 -0.02687825  0.006586116
  • princomp函数的用法:
代码语言:javascript
复制
re2 = princomp(iris[,1:4] %>% scale)
summary(re2)
代码语言:javascript
复制
## Importance of components:
##                           Comp.1    Comp.2     Comp.3      Comp.4
## Standard deviation     1.7026571 0.9528572 0.38180950 0.143445939
## Proportion of Variance 0.7296245 0.2285076 0.03668922 0.005178709
## Cumulative Proportion  0.7296245 0.9581321 0.99482129 1.000000000
代码语言:javascript
复制
re2$scores %>% head
代码语言:javascript
复制
##         Comp.1     Comp.2      Comp.3       Comp.4
## [1,] -2.257141  0.4784238  0.12727962  0.024087508
## [2,] -2.074013 -0.6718827  0.23382552  0.102662845
## [3,] -2.356335 -0.3407664 -0.04405390  0.028282305
## [4,] -2.291707 -0.5953999 -0.09098530 -0.065735340
## [5,] -2.381863  0.6446757 -0.01568565 -0.035802870
## [6,] -2.068701  1.4842053 -0.02687825  0.006586116

可以看到,两者标准差,解释百分比完全一样,PC1一样,PC2正负不一样,PC3一样。

作图的话,两者是上下颠倒的。

代码语言:javascript
复制
ggord(re1,iris$Species)

代码语言:javascript
复制
ggord(re2,iris$Species)

7. 使用FactoMineR包中的PCA函数进行分析

可以看到,结果一致。

代码语言:javascript
复制
library(FactoMineR)
re3 <- PCA(iris[, 1:4], graph = FALSE,scale.unit = T)
summary(re3)
代码语言:javascript
复制
## 
## Call:
## PCA(X = iris[, 1:4], scale.unit = T, graph = FALSE) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4
## Variance               2.918   0.914   0.147   0.021
## % of var.             72.962  22.851   3.669   0.518
## Cumulative % of var.  72.962  95.813  99.482 100.000
## 
## Individuals (the 10 first)
##                  Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3
## 1            |  2.319 | -2.265  1.172  0.954 |  0.480  0.168  0.043 | -0.128
## 2            |  2.202 | -2.081  0.989  0.893 | -0.674  0.331  0.094 | -0.235
## 3            |  2.389 | -2.364  1.277  0.979 | -0.342  0.085  0.020 |  0.044
## 4            |  2.378 | -2.299  1.208  0.935 | -0.597  0.260  0.063 |  0.091
## 5            |  2.476 | -2.390  1.305  0.932 |  0.647  0.305  0.068 |  0.016
## 6            |  2.555 | -2.076  0.984  0.660 |  1.489  1.617  0.340 |  0.027
## 7            |  2.468 | -2.444  1.364  0.981 |  0.048  0.002  0.000 |  0.335
## 8            |  2.246 | -2.233  1.139  0.988 |  0.223  0.036  0.010 | -0.089
## 9            |  2.592 | -2.335  1.245  0.812 | -1.115  0.907  0.185 |  0.145
## 10           |  2.249 | -2.184  1.090  0.943 | -0.469  0.160  0.043 | -0.254
##                 ctr   cos2  
## 1             0.074  0.003 |
## 2             0.250  0.011 |
## 3             0.009  0.000 |
## 4             0.038  0.001 |
## 5             0.001  0.000 |
## 6             0.003  0.000 |
## 7             0.511  0.018 |
## 8             0.036  0.002 |
## 9             0.096  0.003 |
## 10            0.293  0.013 |
## 
## Variables
##                 Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
## Sepal.Length |  0.890 27.151  0.792 |  0.361 14.244  0.130 | -0.276 51.778
## Sepal.Width  | -0.460  7.255  0.212 |  0.883 85.247  0.779 |  0.094  5.972
## Petal.Length |  0.992 33.688  0.983 |  0.023  0.060  0.001 |  0.054  2.020
## Petal.Width  |  0.965 31.906  0.931 |  0.064  0.448  0.004 |  0.243 40.230
##                cos2  
## Sepal.Length  0.076 |
## Sepal.Width   0.009 |
## Petal.Length  0.003 |
## Petal.Width   0.059 |
代码语言:javascript
复制
ggord(re3, iris$Species)

8. 使用ade4包中的dudi.pca函数进行分析

可以看到,分析结果一致。

代码语言:javascript
复制
library(ade4)
代码语言:javascript
复制
## Warning: package 'ade4' was built under R version 4.0.4
代码语言:javascript
复制
## 
## Attaching package: 'ade4'
代码语言:javascript
复制
## The following object is masked from 'package:FactoMineR':
## 
##     reconst
代码语言:javascript
复制
re4 <- dudi.pca(iris[, 1:4], scannf = FALSE, nf = 4)
summary(re4)
代码语言:javascript
复制
## Class: pca dudi
## Call: dudi.pca(df = iris[, 1:4], scannf = FALSE, nf = 4)
## 
## Total inertia: 4
## 
## Eigenvalues:
##     Ax1     Ax2     Ax3     Ax4 
## 2.91850 0.91403 0.14676 0.02071 
## 
## Projected inertia (%):
##     Ax1     Ax2     Ax3     Ax4 
## 72.9624 22.8508  3.6689  0.5179 
## 
## Cumulative projected inertia (%):
##     Ax1   Ax1:2   Ax1:3   Ax1:4 
##   72.96   95.81   99.48  100.00
代码语言:javascript
复制
ggord(re4, iris$Species)

9. 多重对应分析:MCA

multiple correspondence analysis with the tea dataset, MCA

代码语言:javascript
复制
data(tea, package = 'FactoMineR')
tea <- tea[, c('Tea', 'sugar', 'price', 'age_Q', 'sex')]
ord <- MCA(tea[, -1], graph = FALSE)
ggord(ord, tea$Tea, parse = FALSE) # use parse = FALSE for labels with non alphanumeric characters

10. 多重对应分析:MCA

代码语言:javascript
复制
library(MASS)
代码语言:javascript
复制
## 
## Attaching package: 'MASS'
代码语言:javascript
复制
## The following object is masked from 'package:dplyr':
## 
##     select
代码语言:javascript
复制
ord <- mca(tea[, -1])
ggord(ord, tea$Tea, parse = FALSE) # use parse = FALSE for labels with non alphanumeric characters

11. ACM分析

代码语言:javascript
复制
ord <- dudi.acm(tea[, -1], scannf = FALSE)
ggord(ord, tea$Tea, parse = FALSE) # use parse = FALSE for labels with non alphanumeric characters

12. metaMDS分析

代码语言:javascript
复制
library(vegan)
代码语言:javascript
复制
## Warning: package 'vegan' was built under R version 4.0.4
代码语言:javascript
复制
## Loading required package: permute
代码语言:javascript
复制
## Warning: package 'permute' was built under R version 4.0.4
代码语言:javascript
复制
## Loading required package: lattice
代码语言:javascript
复制
## This is vegan 2.5-7
代码语言:javascript
复制
ord <- metaMDS(iris[, 1:4])
代码语言:javascript
复制
## Run 0 stress 0.03775523 
## Run 1 stress 0.05277761 
## Run 2 stress 0.04367524 
## Run 3 stress 0.04709622 
## Run 4 stress 0.05433734 
## Run 5 stress 0.03775525 
## ... Procrustes: rmse 1.075667e-05  max resid 4.153723e-05 
## ... Similar to previous best
## Run 6 stress 0.05327828 
## Run 7 stress 0.0377552 
## ... New best solution
## ... Procrustes: rmse 1.336304e-05  max resid 0.0001216059 
## ... Similar to previous best
## Run 8 stress 0.04355786 
## Run 9 stress 0.04713709 
## Run 10 stress 0.03775521 
## ... Procrustes: rmse 8.440975e-06  max resid 3.786993e-05 
## ... Similar to previous best
## Run 11 stress 0.04355785 
## Run 12 stress 0.04367519 
## Run 13 stress 0.05299973 
## Run 14 stress 0.04713728 
## Run 15 stress 0.05093498 
## Run 16 stress 0.03775523 
## ... Procrustes: rmse 1.225092e-05  max resid 5.498372e-05 
## ... Similar to previous best
## Run 17 stress 0.03775526 
## ... Procrustes: rmse 2.028828e-05  max resid 8.412387e-05 
## ... Similar to previous best
## Run 18 stress 0.0377552 
## ... Procrustes: rmse 1.335765e-06  max resid 5.738162e-06 
## ... Similar to previous best
## Run 19 stress 0.04709624 
## Run 20 stress 0.04355785 
## *** Solution reached
代码语言:javascript
复制
ggord(ord, iris$Species)

13. linear discriminant analysis

代码语言:javascript
复制
ord <- lda(Species ~ ., iris, prior = rep(1, 3)/3)

ggord(ord, iris$Species)

14. correspondence analysis

代码语言:javascript
复制
ord <- dudi.coa(iris[, 1:4], scannf = FALSE, nf = 4)

ggord(ord, iris$Species)

15. correspondence analysis

代码语言:javascript
复制
library(ca)
代码语言:javascript
复制
## Warning: package 'ca' was built under R version 4.0.4
代码语言:javascript
复制
ord <- ca(iris[, 1:4])
ggord(ord, iris$Species)

16. double principle coordinate analysis (DPCoA)

代码语言:javascript
复制
library(ade4)
data(ecomor)
grp <- rep(c("Bu", "Ca", "Ch", "Pr"), each = 4)    # sample groups
dtaxo <- dist.taxo(ecomor$taxo)                    # taxonomic distance between species
ord <- dpcoa(data.frame(t(ecomor$habitat)), dtaxo, scan = FALSE, nf = 2)
ggord(ord, grp_in = grp, ellipse = FALSE, arrow = 0.2, txt = 3)

18. distance-based redundancy analysis

代码语言:javascript
复制
data(varespec)
data(varechem)

ord <- dbrda(varespec ~ N + P + K + Condition(Al), varechem, dist = "bray")

ggord(ord)

19. redundancy analysis

代码语言:javascript
复制
ord <- rda(varespec, varechem)
ggord(ord)

20. distance-based redundancy analysis

代码语言:javascript
复制
ord <- capscale(varespec ~ N + P + K + Condition(Al), varechem, dist = "bray")

ggord(ord)

21. canonical correspondence analysis

cca from vegan

代码语言:javascript
复制
ord <- cca(varespec, varechem)

ggord(ord)

增加标签。

代码语言:javascript
复制
ggord(ord, ptslab = TRUE, size = NA, addsize = 5, parse = TRUE)
代码语言:javascript
复制
## Warning: Removed 24 rows containing missing values (geom_point).

感觉不错,转发一波吧!!!

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2021-03-20,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 育种数据分析之放飞自我 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 1. 使用iris数据
  • 2. 使用prcomp做聚类分析
    • 3. 绘制第一张图
      • 4. 设置不同风格的图形
      • 5. 使用princomp函数,进行PCA的分析
      • 6. R中两个函数prcomp和princomp的区别
      • 7. 使用FactoMineR包中的PCA函数进行分析
      • 8. 使用ade4包中的dudi.pca函数进行分析
      • 9. 多重对应分析:MCA
      • 10. 多重对应分析:MCA
      • 11. ACM分析
      • 12. metaMDS分析
      • 13. linear discriminant analysis
      • 14. correspondence analysis
      • 15. correspondence analysis
      • 16. double principle coordinate analysis (DPCoA)
      • 18. distance-based redundancy analysis
      • 19. redundancy analysis
      • 20. distance-based redundancy analysis
      • 21. canonical correspondence analysis
      • cca from vegan
      领券
      问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档