R语言|主成分分析的图如何读？如何用R快速实现？

原创通往学术之路通往学术之路

主成分分析（PCA）分析图到底怎么读？

主成分分析（The Principal Component Analysis，PCA），属于一种无约束排序分析。对于主成分分析，从其计算过程到出图过程，从他人论文到自己发表的论文，主成分分析基本上是好多文章的必备，相信大部分人都非常熟悉了。不过，我相信大家都应该会有这样的迷茫：主成分图到底在说什么？那些说主成分一代表了XX指标，PC2代表了YY指标，又有多少解释量；还有从夹角大小，箭头方向等等来阐述主成分分析的结果，这些信息是如何获得的呢。

相信有人也向身边的人请教过，但是回答总是不疼不痒。对于这样的结果，我们是非常拒绝的。曾经的我也问过我老师，老师的回答是：主成分分析是通过降维的方法，寻找到能代表数据变化特征的变量，这样的话我根本听不懂也想不明白，因为书上对于主成分分析的定义是：通过降维技术把多个变量转化成少数几个主成分，这些主成分是能够反应原始变量的绝大多部分信息，是原始变量的线性组合。也上过数理统计老师的课，在讲主成分分析时也非常认真听讲，但是老师仅仅讲述主成分计算的原理和过程，对于如何运用，如何解释看到的PCA图，老师几乎没有涉及。

云里雾里的说法和解释让人感觉高深莫测，也很专业，然而对实际运用一点价值都没有。个人认为，当老师在介绍某个方法时，如果不是针对特定数学专业或计算机专业的学生，最好不要过多涉及方法的底层算法逻辑，因为非数学专业的人并不需要根据方法的底层算法去设计新的方法，最重要的是他想知道通过这样的方法，获得什么样的结果，这个结果对于解释他遇到的分析对象有哪些用处，如何解释分析结果并作推广。

唠叨几句，言归正传，对于主成分分析，好在我不抛弃，不放弃的坚持探寻，终于理解了主成分结果的意义。因此今天来仔细做一个关于主成分分析的总结。

关于定义及算法并不想涉及，这里只涉及解读主成分分析图和如何用R语言实现它。

PCA图形解读

上一个vegan包中内置数据理化性质部分作的主成分分析图（数据未转化）：

下面就以这幅图来讲述它背后的意义（主要逻辑是作垂线）。

对于样点来说（虽然实际分析用不到），如果该样点处于图的原点（0，0）附近，则说明该点对应的指标数值位于相应变量的均值附近，如上图的12号和15号点，它们的所有理化数据均处于平均水平。如果该样点远离原点位置，则先看它所处的象限，然后再作它与理化箭头或物种箭头（本例中使用理化做PCA，故不包含）方向或正向/反向延长线方向的垂线，以垂线和箭头的交点到原点（0，0）的距离长短做解释，距离越长，则说明该点处对应的指标数值越大，反之越小。如下图：

图中的绿色线比紫色线长，说明24号点的S含量比16号点的要高，回到原始数据上验证一下：

确实如此，这样的情况即使是将数据标准化后也是一样的。

对于变量对主成分的贡献（因子载荷）大小来说，也是看该变量在对应主成分上的投影（变量箭头末端作与对应坐标轴的垂线，然后看垂线交点到原点的距离长短-即投影）。如果变量A比B的投影长，则说明A对相应坐标轴的贡献比B大（关系更密切），反之则小。

如下图：

上图中，S在PC1上的投影（紫色线）要高于Mn（绿色线）和N（蓝色线），说明三者对PC1的贡献排序为：S>Mn>N。而对于同向的变量来说（如下图），很明显，Ca对PC1的贡献(绿色线)明显高于S（紫色线），但Ca对PC2的贡献（绿色点线）却低于S(紫色点线)。

对来变量间的关系来说，看它们之间的夹角，夹角越小，说明它们间的相关性越强，比如下图中挑选出来的Ca和Mg，它们几乎重叠了，右边的图为采用相关系数计算后的出图，发觉它们两者相关性非常强：

因此，对于变量间的关系，个人认为如果做了主成分分析，那就没必要再加相关性分析了，实在有特殊要求时，可在材料方法部分阐述使用相关分析获得变量间的Pearson相关系数和p值即可，给更多的分析留版面。

关于主成分分析的结果，还剩最后一个关键问题：PC1和PC2代表了哪些信息？还是作垂线来比较，不过更我建议采用“视觉评判法”，比如上图中，根据视觉效果，N在PC2上的投影明显高于在PC1上的，而诸如Al、Fe也一样，那么可以尝试将PC2归为代表了土壤中N、Al和Fe等化合物形态更多的轴；其它诸如Mg、Ca等变量在PC1上的投影长于PC2，那么这些可以归类为化合物形态相对较少的轴。

其实这样归类带有主观评判的意味，具体操作过程中，可以根据背景知识来对PC1和PC2归类。实在不会，推荐的做法是尽量归类正确，然后坐等审稿专家大牛对这一归类提出修改意见即可，不过也需要遇到哪些对PCA比较了解和负责的审稿专家，遇到啥也不会的那种砖家的话，那就算了，能发文就行。虽有投机取巧之意，但也不失为一种解决办法。

最后：主成分分析不是通过降维找能代表所有数据的少数变量，而是将多维变量降维成两个（最常用）主要成分（类似潜变量的意思），而这两个成分（PC1和PC2）能最大化的代表原始数据信息。

R中PCA实现

对于R中的PCA实现，其实非常简单，下面分别使用两种方式实现：

 #加载R包，设置工作环境library(ggplot2)library(vegan)setwd("D:/Rstudio/R")library(corrplot)data("varechem")#取vegan自带环境数据进行分析mypca <- princomp(varechem,cor=TRUE) #基于原始数据的PCAenv1 <- log(varechem+1)#原始数据进行转换mypca1 <- princomp(env1,cor=TRUE)#基于转换后的PCAsummary(mypca)#查看结果1summary(mypca1)#查看结果2#出图biplot(mypca,choices = 1:2,scale = 1,col = c('red', 'blue'))abline(v=0,h=0)#左，数据未转化biplot(mypca1,choices = 1:2,scale = 1,col = c('red', 'blue'))abline(v=0,h=0)#右，数据转化

上图中左图和右图的区别初看差异不大，但如果仔细看就会有问题，比如Fe和pH,左图明显发生重叠，而右图明显Al和Fe的关系更紧密，结合本文上方提到的相关系数图会发现：数据转化后的更符合相关关系！这点点的差异说明了变量的量级变异对计算结果影响非常大，大家可以翻看上方关于理化数据部分的截图，会发现Fe和Al的中的数值差异可达100，特别是Al。因此，在进行数据分析时，最好对数据进行转化，然后再进行分析。

由于R自带的主成分分析中，并没有将位置数据和变量数据给予展示，因此推荐使用vegan包的rda()函数进行PCA分析，这里不提供解释变量，如果提供解释变量则未RDA分析，代码如下：

pca_env <- rda(env1, scale = TRUE) #没有解释变量summary(pca_env, scaling = 1) #提取信息summary(mypca1,scaling=1)#与R自带PCA做对比

结果图如下：

可以看到其解释量也是一致的，不同的是vegan包的结果中还给出了Species scores（物种得分，这里是理化变量在主成分图中的坐标位置）和Site scores（样点得分，这里是样点在主成分图中的坐标位置）。有朋友会好奇：获得的物种得分和样地得分有什么作用，先前我也看搞不懂，但现在想通了，这个一方面能方便输出到ggplot2中出更好看的图，另一方面还看到过有的大牛用PCA降维后，提取环境因子或物种因子在主成分坐标轴上的值，用以建模或做回归分析，后面再遇到这类文章时再好好总结一下。

最后，ggplot2出图部分，就留给有余力的朋友做练习了，这里给出提取相应坐标的代码：

pca_exp <- pca_env$CA$eig / sum(pca_env$CA$eig)#计算每个主成分的贡献量exp <- pca_exp[1:2]#提取PC1和PC2的值site <- scores(pca_env, choices = 1:2, scaling = 1, display = 'site')#提取样点位置env<- scores(pca_env, choices = 1:2, scaling = 2, display = 'sp')#提取环境变量位置

完整代码

#加载R包，设置工作环境library(ggplot2)library(vegan)setwd("D:/Rstudio/R")library(corrplot)data("varechem")#取vegan自带环境数据进行分析mypca <- princomp(varechem,cor=TRUE) #基于原始数据的PCAenv1 <- log(varechem+1)#原始数据进行转换mypca1 <- princomp(env1,cor=TRUE)#基于转换后的PCAsummary(mypca)#查看结果1summary(mypca1)#查看结果2#出图biplot(mypca,choices = 1:2,scale = 1,col = c('red', 'blue'))abline(v=0,h=0)#左，数据未转化biplot(mypca1,choices = 1:2,scale = 1,col = c('red', 'blue'))abline(v=0,h=0)#右，数据转化pca_env <- rda(env1, scale = TRUE) #没有解释变量summary(pca_env, scaling = 1) #提取信息summary(mypca1,scaling=1)#与R自带PCA做对比pca_exp <- pca_env$CA$eig / sum(pca_env$CA$eig)#计算每个主成分的贡献量exp <- pca_exp[1:2]#提取PC1和PC2的值site <- scores(pca_env, choices = 1:2, scaling = 1, display = 'site')#提取样点位置env<- scores(pca_env, choices = 1:2, scaling = 2, display = 'sp')#提取环境变量位置