查看: 3940|回复: 9

[R语言] R统计绘图-PCA分析绘图及结果解读(误差线,多边形,双Y轴图)

  [复制链接]

版主

Rank: 10Rank: 10Rank: 10

主题
31
注册时间
2017.4.3
在线时间
210 小时

发表于 2021.12.14 18:36:19 | 显示全部楼层 |阅读模式
本帖最后由 静默人声412 于 2021.12.14 18:37 编辑

虽然PCA和RDA分析及绘图都写过教程,但是对于结果的解释都没有写的很详细,刚好最近有人询问怎样使用FactoMineR factoextra包进行PCA分析。所以使用R统计绘图-环境因子相关性热图中的不同土壤环境因子数据进行PCA绘图和结果解读推文。
一、 数据准备


  1. <i style="color: rgb(175, 175, 175); font-family: Consolas; background-color: rgb(255, 255, 255);"># 1.1 设置工作路径</i>
复制代码

图1|data.csv。环境因子分类信息可以另整理一张表。没有环境因子分类信息就整理为图2格式即可。


图2|环境因子及分组信息表,env.csv。行为样品名称,列为环境因子名称和分组信息,共有11个环境变量,3个分组信息。
二、 主成分分析(PCA)
数据都准备好后,使用FactoMineR中的PCA()函数进行PCA分析,factoextra包可用于提取PCA分析结果信息和绘图(用ggplot2绘图也可)。配套的FactoInvestigate包(http://factominer.free.fr/reporting/index.html)可以自动输出FactoMineR中PCA,MCA和CA分析的结果解释。
  1. # 2.1 调用R包
  2. #install.packages("FactoInvestigate")
  3. #install.packages("htmltools",version="0.5.2")
  4. library(FactoMineR)
  5. library(factoextra) # 用于提取PCA分析结果信息,其实不用也可以。
  6. library(FactoInvestigate)
  7. #browseVignettes("FactoMineR") # 查看帮助文档


  8. # 2.2 PCA分析
  9. #res.pca <- PCA(env[4:14], graph = TRUE,axes = c(1,2))
  10. ##graph = TRUE会默认绘制前两个轴的样本与响应变量PCA图(两张图)。
  11. #page(PCA) # 查看函数代码


  12. ## 具有附加定性变量的PCA分析:附加样本、定性和定量变量对主成分分析没有贡献,
  13. ###坐标来自于基于PCA分析结果的预测值。
  14. res.pca <- PCA(env[2:14],
  15.                scale.unit = TRUE,# 默认对数据进行scale标准化,最后响应变量间的均值为0,sd=1。
  16.                ##相同尺度的标准化避免了一些变量仅仅因为其较大的度量单位而成为主导变量。使变量具有可比性。
  17.                ncp = 11, # 结果中默认保留5,nrow(env[-ind.sup,])-1(样本数)和ncol(env[2:14][-c(quanti.sup,quali.sup)])(响应变量数)中的最小值个主成分
  18.                #ind.sup = 1, # 附加样本,1表示数据表的第一行,不是必须的数据,默认为NULL。附加样本的坐标将根据PCA分析结果进行预测。
  19.                quanti.sup = NULL, # 附加的定量响应变量数据,不是必须的,默认为NULL。基于PCA分析结果预测坐标。
  20.                quali.sup = 1:2, # 附加定性变量,分类信息最好不要用数字表示,在行数区间的数字,运行会报错。
  21.                ##1:2表示定性变量在env[2:14]中的列数。不是必须的,默认为NULL。基于PCA分析结果预测坐标,可用于对样本进行分组着色。
  22.                row.w = NULL, col.w = NULL, # 分别给样本和响应变量设定权重,默认所有样本或所有响应变量的权重都是一样的。默认相当于对变量数据进行了中心化,变量均值接近0。可以设置为1个数字,也可以设置为分别与行数或列数相等的矢量。
  23.                graph = FALSE)
复制代码
注::PCA是基于线性模型的排序方法,因此会涉及数据的中心化和标准化。数据中心化是将原数据减去其均值,中心化后数据的均值为0。样方的中心化是让每个样方的平均值为0,响应变量中心化是让每个响应变量的平均值为0。如果响应变量的量纲不同,响应变量数据就需要标准化。数据的标准化(归一化)是指将数据按照比例缩放,使之所有变量具有相同的尺度,比如常用的Z-score标准化(先对变量进行零中心化,再除以数据标准误差),或者离差标准化(数据减去最小值,再除以最大值与最小值的差值)后数据尺度为[0,1],而MaxAbsScaler绝对值最大标准化(数据除以最大绝对值,不移动中心点。不会将稀疏矩阵变得稠密)后数据尺度为[-1,1]。数据转换方法的选择,需要考虑研究在意的是数据的数值大小变化还是数据变化程度的变化,在意数值本身的变化可以使用原始数据或log转换,在意数据变化程度,则可进行标准化。当数据集中含有离群点,即异常值时,可以用z-score进行标准化,但是标准化后的数据并不理想,因为异常点的特征往往在标准化之后容易失去离群特征。此时可以用RobustScaler方法针对离群点做标准化处理。
FactoInvestigate包的Investigate()可以自动输出FactoMineR包PCA分析的结果解释,输出文件可以是word,html和PDF。默认输出离群值,主成分对样本方差的解释度检测,维度描述,样本聚类及对每个cluster贡献高的变量等信息。
  1. # 2.3 FactoInvestigate包的Investigate()输出结果报告
  2. #?Investigate # 查看函数帮助信息
  3. Investigate(res.pca,
  4.             file = "Investigate.Rmd",openFile = FALSE,
  5.             document = c("word_document"))
复制代码
注:输出的结果包括Investigate.docx或Workspace.RData和Investigate.Rmd。Investigate.docx是word版PCA分析结果报告与结果解读,写的很详细。如果因为某种原因输出报告文件报错(常见htmltools包版本或LaTeX报错),就会输出Workspace.RData和Investigate.Rmd,包含分析代码和绘图结果,两者放在在同一目录,运行Investigate.Rmd,可以复现PCA分析结果。图不是太好看,可以输出数据自己画。
  1. # 2.4 PCA分析结果的输出数据介绍
  2. str(res.pca)# 输出一个list数据集
  3. summary(res.pca,ncp=2,nbelements=Inf) # nbelements=Inf,输出所有结果


  4. ## 2.4.1 主成分(PC)特征值
  5. ## eigenvalue:特征值可以用来确定在PCA之后要保留的主成分的数量
  6. get_eigenvalue(res.pca) # 11个主成分累计代表的样本方差为100%。
  7. sum(res.pca$eig[,1]) # 所有主成分所表示的总方差为11(11个原始变量)。
  8. res.pca$eig # 包含3列数据,每个主成分的特征值,方差百分比和累积方差百分比。
复制代码


图3|主成分特征值,res.pca$eig。特征值表示每个主成分保留的原始变量的变化量。特征值从高到低排序,PC1表示变量变化最大的方向。选择代表方差>1(特征值>1表示PC比标准化数据中的一个原始变量解释的方差更大。数据未标准化或变量权重不同,不能这样进行选择),此处累计方差解释量~65%(>80%更好)的PC1和PC2进行可视化和后续分析(也可以选择累计方差解释率超过一定阈值的PC用作后续分析)。
  1. ## 2.4.2 样本PCA输出数据
  2. res.pca$ind # 包含所有样本的输出结果的矩阵列表
  3. #get_pca_ind(res.pca) # 可以使用factoextra包的函数提取样本输出数据。


  4. res.pca$ind$coord # 样本在PC的坐标=左奇异值/sqrt(eig)
  5. res.pca$ind$cos2 # 样本坐标的平方/样本Euclid(L2)范数的平方(标准化数据)。
  6. ## 未标准化数据,样本1的cos2=ind1.coord^2/每个个体与PCA重心之间的平方距离(d2)。d2=([(Var1-Var1.mean)/Var1.sd]^2+...+[(Varn-Varn.mean)/Varn.sd]^2)
  7. res.pca$ind$contrib # 所有样本对主成分的贡献度(%)=(res.pca$ind$coord^2)*100/eig^2/样本数
  8. res.pca$ind$dist # 变量标准化之后,样本L2范数(Euclid范数),即每个样本包含变量的(值平方*变量权重)的和的平方根
  9. ##可反应样本间距离。dist值差异越大,样本相似性越小。
复制代码
注:L2范数指向量中所有数据平方和的平方根。范数介绍:https://blog.csdn.net/a493823882/article/details/80569888
  1. ## 2.4.3 响应变量PCA输出数据
  2. res.pca$var # 包含所有响应变量相关结果的矩阵列表
  3. #get_pca_var(res.pca) # 可以使用factoextra包的函数提取响应变量输出数据。


  4. res.pca$var$coord # 所有响应变量(这里是环境因子)的坐标,等于右奇异值*特征根的平方根=载荷值(主成分能代表的响应变量的方差)*特征根(主成分的标准误差)。


  5. res.pca$var$cor # 样本权重一样时,每个响应变量与轴的相关性=响应变量坐标/(变量值的平方和/样本数的平方根)。绝对值越大,表明相关性越强,正负表示方向。


  6. res.pca$var$cos2 # 所有响应变量与轴相关性的平方,等于res.pca$var$cor^2。


  7. res.pca$var$contrib # 响应变量对每个主成分的方差贡献值=(var$cos2*100)/apply(res.pca$var$cos2,2,sum)
复制代码
注:PCA()也计算了响应变量的dist,计算方式为dist2.var <- as.vector(crossprod(rep(1, nrow(X)), as.matrix(X^2 * row.w)))。X为标准化后的数据矩阵,row.w=样本权重/(样本权重之和),crossprod()内积运算。
  1. # 2.4.4 附加变量预测输出结果:附加样本、定性和定量变量对主成分没有任何贡献。
  2. res.pca$quali.sup # 附加定性变量中每个分类的输出结果。最后输出的报告里会基于Wilks test分析那些定性变量能更好的解释样本之间的差异。
  3. res.pca$quali.sup$coord # 定性变量中每个分类的坐标
  4. res.pca$quali.sup$cos2 # 定性变量中每个分类被主成分代表的部分
  5. res.pca$quali.sup$v.test # 定性变量中每个分类Normal distribution准则
  6. res.pca$quali.sup$eta2 # 定性变量与每个主成分的相关性系数的平方
  7. res.pca$quali.sup$dist # 可用于比较分类变量个水平间的距离。
复制代码


图4|附加变量预测输出结果。
  1. # 2.4.5 其它输出数据
  2. library(rstatix)
  3. #env[4:11] %>% get_summary_stats(type="mean_sd")
  4. res.pca$call$centre # 响应变量的均值
  5. res.pca$call$ecart.type # 响应变量的标准误差
  6. res.pca$call$row.w.init # 样本权重
  7. res.pca$call$col.w # 响应变量权重
  8. res.pca$call$ncp # 主成分维度
  9. #env[4:11] %>% group_by(env$tillage) %>% get_summary_stats(type="mean")
  10. res.pca$call$quali.sup # 附加定性变量统计值
  11. res.pca$call$X # 原始数据
  12. res.pca$call$call # 运行代码
复制代码

三、PCA分析结果绘图
factoextra包中的绘图函数fviz_pca_biplot()可以用于绘制PCA双序图。fviz()是ggpubr包中ggscatter()函数的包装器,可以与ggpubr的函数联用。也可用提取的数据使用ggplot2绘图。这里先使用ggplot2绘图,下篇文章介绍factoextra包中的绘图函数。
3.1 提取绘图数据
FactoInvestigate包的Investigate()默认输出的图很多,但是信息没有整合,而且很多细节还需要调整。这里提取数据,自己使用ggplot2进行绘图。
  1. # 3.1.1 提取样本坐标
  2. # get_pca_ind(res.pca) # factoextra提取分析结果,可以很方便的提取prcomp()和princomp()函数的样本PCA分析结果。
  3. sam=data.frame(res.pca$ind$coord)
  4. sam
  5. sam=data.frame(sam[1:2],env[1:3]) # 添加样本分组信息
  6. sam


  7. # 3.1.2 提取特征值-每个主成分对样本方差的解释度。
  8. fviz_eig(res.pca,addlabels = TRUE) # 对每个主成分对方差的解释度绘制柱形图。
  9. #fviz_screeplot(res.pca)# 对每个主成分对方差的解释度绘制柱形图。
  10. eig.val <- get_eigenvalue(res.pca) # get_eig(res.pca)效果一样,提出的是res.pca$eig矩阵中数据。
  11. eig.val
  12. ## 提取前两个主成分是特征值,并保留两位小数点。
  13. PC1=round(eig.val[1,2],2)
  14. PC2=round(eig.val[2,2],2)
  15. PC1
  16. PC2


  17. # 3.1.3 提取响应变量坐标值
  18. var <- get_pca_var(res.pca) # 提取响应变量结果矩阵列表。
  19. res=data.frame(var$coord) # res=data.frame(res$var$coord),效果一样。
  20. res
  21. type
  22. res$type = type[match(rownames(type),rownames(res)),1] # 添加环境因子分类信息
  23. res


  24. # 3.1.4 提取附加定性变量的坐标值
  25. quali = data.frame(res.pca$quali.sup$coord)
  26. quali
复制代码


图5|样本坐标信息表,sam。
factoextra包中提取样本、变量和特征根PCA分析结果的函数,常用于prcomp()和princomp()分析结果提取。
  1. # 3.1.5 prcomp()和princomp()主成分分析
  2. prcmp = prcomp(env[4:14], scale =TRUE)
  3. princmp= princomp(env[4:14], cor = TRUE, scores = TRUE) # cor = TRUE表示对数据进行标准化,scores=TRUE表示输出所有PC结果。
  4. prcmp$sdev # 特征根=princmp$sdev
  5. prcmp$rotation # 变量载荷矩阵(列就为PC的特征向量)=princmp$loadings,可用于预测附加样本坐标。
  6. ## 载荷值*特征根得出变量在每个主成分的坐标
  7. prcmp$center# 变量均值=princmp$center
  8. prcmp$scale # 变量的standard deviations=princmp$scale
  9. prcmp$x # 样本观察值(坐标)=princmp$scores


  10. get_pca_ind(prcmp) # 提取样本坐标和代表性等信息。
  11. get_pca_var(prcmp) # 提取变量坐标和代表性等信息。
  12. get_eigenvalue(prcmp) # 提取主成分特征值、解释度和累积解释度。
复制代码
注:prcomp()和princomp()提取的结果与PCA()有些微差别,分析是对数据的处理有些微差异,prcomp()和princomp()标准误差计算用的是n(样本数)-1,而PCA()使用的n。推荐使用prcomp()和princomp()分析后,使用factoextra包函数提取坐标等信息。PCA()不会输出变量载荷矩阵。
3.2 绘图
这里将以不同颜色(tillage)的形状(depth)表示样本,箭头表示响应变量(不同分类用不同颜色表示),用名称展示分类变量tillage和depth,按照tillage计算凸点,绘制分组多边形。
  1. library(ggsci)
  2. library(scales)
  3. library(ggrepel)


  4. # 3.2.1 准备绘图颜色
  5. mypal=pal_d3("category10")(10) # 使用ggsci中d3颜色集
  6. mypal # 之后所有图都使用这些颜色
  7. show_col(mypal) # 查看颜色


  8. # 3.2.2 绘图-不绘制多边形
  9. p = ggplot(sam,aes(x=Dim.1,y=Dim.2))+
  10.   geom_point(aes(color=tillage,shape=depth),size=5)+
  11.   scale_shape_manual(values = c(15:17))+ # guide = "none",不显示图例;key_glyph="point" # 设置图例的形状
  12.   scale_color_manual(values = mypal,
  13.                      #name="Tillage", # 设置图例名称
  14.                      #guide = "none",
  15.                      limits=c("CK","NT","PT","WL","anorganic","organic","N format","pH"))+ # 设置图例显示顺序
  16.   geom_text(inherit.aes = FALSE,data=quali,
  17.              aes(x=Dim.1,y=Dim.2,label=row.names(quali)),
  18.              color="red",size=6)+
  19.   geom_segment(data=res,
  20.                aes(x=0,y=0,xend=res[,1],yend=res[,2],color=factor(type)),
  21.                size=0.75,
  22.                #arrow.fill = "white" # 与arrow()中的type="closed"联用形成空心箭头。
  23.                arrow=arrow(angle =35,length=unit(0.3,"cm")))+ # type = "closed",可以把箭头设置为实心,key_glyph="abline" #设置图例形状为斜线
  24.   geom_text(data=res,show.legend = FALSE,
  25.             aes(x=res[,1],y=res[,2],label=rownames(res),color=type),size=6,
  26.                   hjust=(1-sign(res[,1]))/2,
  27.             angle=(180/pi)*atan(res[,2]/res[,1]),
  28.                   position = "jitter")+
  29.   geom_vline(xintercept=0)+ geom_hline(yintercept=0)+
  30.   #guides(color=guide_legend(),fill=guide_legend(),colorbar=guide_colourbar()) # 根据颜色和形状系统设置图例
  31.   theme_bw()+
  32.   labs(title="PCA plot",
  33.        x=paste("PC1",PC1," %"),
  34.        y=paste("PC2",PC2," %"))+
  35.     theme(legend.position="right", # legend.position="none"图中不展示图例
  36.         legend.title = element_text(face = "bold", size = 18,colour = "black"),
  37.         legend.text = element_text(size = 12),
  38.         panel.grid=element_blank(),
  39.         axis.title = element_text(face = "bold", size = 18,colour = "black"),
  40.         axis.text = element_text(face = "bold", size = 16,color="black"))
  41. ggsave("PCA.pdf",p,device = "pdf",width = 12,height=12,family="Times") # 输出pdf图片到本地


  42. # 3.2.3 按tillage添加分组多边形
  43. ##计算凸点
  44. library(plyr)
  45. find_hull <- function(df) df[chull(df[,1], df[,2]),]
  46. hulls <- ddply(sam, "tillage", find_hull)
  47. hulls


  48. ## 绘制多边形
  49. p2 = p + geom_polygon(data=hulls,aes(x=Dim.1,y=Dim.2,fill=tillage),                       alpha=0.3,show.legend = FALSE)+
  50.   scale_fill_manual(values = mypal)
  51. p2
  52. ggsave("PCA2.pdf",p2,device = "pdf",width = 12,height=8,family="Times")
复制代码


图6|PCA双序图。图例没有区分tillage和type,pdf需要再编辑一下。每个轴样本/变量坐标的变异反映了每个PC的重要性,通过特征值来度量。特征值越大,该PC上样本的分散程度越大。因此,PC1的样本/变量坐标较第二轴的分散程度大。


图7|按分类设置多边形阴影PCA双序图。图中的信息包含的可能太多了。也可以将样本和变量分开两幅图绘制。fviz_pca_biplot()也可以绘制各种分组阴影。如图所示,NT与其它tillage沿着PC2轴存在差异;不同tillage中,沿着PC1,depth间距离存在差异,这些差异结果都要再结合多变量统计分析进行解释。可参考R统计-微生物群落结构差异分析及结果解读进行多变量样本差异分析统计分析。
四、 PCA结果解读
PCA图解释可以参照FactoInvestigate包的Investigate()输出的PCA分析结果报告进行解释。PCA分析是对高维数据进行降维的一种方法,用于从多变量数据中提取重要信息。当数据集中的变量高度相关时,PCA方法特别有用。相关性表明数据中存在冗余。由于这种冗余性,使用PCA将原始变量通过线性组合,转换为几组没有相关关系的新变量(主成分),以解释原始变量中的大部分方差。在样本数大于响应变量个数时,主成分个数小于或等于原始变量个数。每个主成分(PC)都代表样本间原始变量的一部分信息,主成分特征值(eigenvalue)表示每个轴所能代表的样本间变量信息(差异程度),等于数据在变换后的坐标上对应PC的方差(每个主成分保留的方差)。综上所述,主成分分析的主要目的是:1)识别数据集中的潜在变量,2)通过去除数据中的噪声和冗余来降低数据的维度,3)识别相关变量。
4.1 主成分(PC)
每个主成分的特征值,是数据矩阵奇异值的平方。样本和响应变量坐标都是使用样本(left singular vectors)/响应变量(right singular vectors)奇异值除以PC特征值的平方根得到的。奇异值分解降维参考:数据降维:特征值分解和奇异值分解的实战分析数据标准化之后,每个响应变量的方差为1,所以整个数据矩阵的方差为11,每个PC对应一个新变量,特征值即为PC能代表的方差,每个PC的特征值除以所有PC特征值之和(变量总方差)得到每个PC对原始变量方差的解释度。越少PC能反应原始变量的方差的比例越好越好,表示PC能很好的代表原始变量间的差异。需要将PCA分析得到的新变量(PC),用于与其它变量进行关联分析,挑选的PC需要能很好的代表原始变量的变化。1)挑选PC特征值>变量总方差/变量个数的PC;2)挑选能代表>80%原始变量差异的PC;3)PC特征值大于随机数据矩阵的平均特征值(FactoInvestigate包的Investigate()会输出模拟结果的PC解释度参考阈值)。可参考一文了解R语言数据分析 ----主成分分析
  1. # 4.1.1 特征值计算过程
  2. ## 计算变量均值,中心化
  3. X=env[4:14]
  4. centre = apply(X, 2,mean) # 求每个响应变量的均值
  5. X <- t(t(as.matrix(X)) - centre)


  6. ## 中心化数据求范数均值,标准化
  7. ## 标准误差(均方误差:[url]https://zhidao.baidu.com/question/393022200642177405.html[/url])指变量的离均差(数据-mean)平方的算数平均数(均方)的平方根(计算方式同样本标准误)。
  8. # PCA()函数中的标准化结果与sd计算结果存在略微差异
  9. #sd = apply(X, 2,sd) # sd()函数计算离均差平方的算数平均数时除以是(样本数-1),标准化后数据标准误差为1。
  10. col.w <- rep(1, ncol(X)) # 变量权重
  11. row.w <- rep(1, nrow(X)) # 样本权重
  12. ec.tab <- function(V, poids) {
  13.         ecart.type <- sqrt(as.vector(crossprod(poids/sum(poids),
  14.             as.matrix(V^2))))
  15.         ecart.type[ecart.type <= 1e-16] <- 1
  16.         return(ecart.type)
  17.     }#corssprod()矩阵内积计算,同一行数据相乘后,同一列数据相加,得到与矩阵列数相同的矢量
  18. ecart.type <- ec.tab(X, row.w) # row.w表示样本权重矢量,默认所有样本都为1。
  19. ##PCA()函数标准化过程,计算离均差平方的算数平均数时除以的是(样本数);标准化后数据使用sd()计算标准误差值接近1;
  20. ##这与vegan::rda()函数数据标准化结果存在略微不同,rda()使用scale()标准化,除以(样本数-1)。
  21. X <- t(t(X)/ecart.type)


  22. ## 奇异值分解,计算特征值
  23. ### Singular Value Decomposition of a Matrix
  24. tmp = svd.triplet(X, row.w = row.w, col.w = col.w, ncp = 5)
  25. ### 计算特征值
  26. eig <- tmp$vs^2


  27. # 4.1.2  # 绘制各主成分能代表的方差比例。
  28. sum(res.pca$eig[,1]) # PC特征值和为11。
  29. fviz_eig(res.pca,addlabels = TRUE)
  30. ggsave("Scree plot.pdf",width = 12,height = 12,family="Times")
  31. ##FactoInvestigate包的Investigate()输出结果中,
  32. ##包含模拟1552个均匀分布的相同大小数据矩阵进行PCA分析,
  33. ##最后,选取0.95分位数的PC解释度作为挑选PC的参考值。
复制代码
注:ec.tab()来自PCA()函数。R中有很多进行PCA分析的函数,有时候会发现计算出来的结果存在差异(数据处理方式不同,相关性矩阵/协方差矩阵)。比如FactoMineR::PCA()和vegan::rda()等。所以使用R中的查看函数代码功能,查看运行的函数对原始数据做了那些处理是很重要的。下面记录一下R中查看函数代码的命令。
  1. # 直接输入函数名查看
  2. PCA
  3. # page()查看
  4. page(PCA)
  5. # edit()查看,可修改
  6. edit(PCA)
  7. # methods()
  8. methods(rda)
  9. # getAnywhere()
  10. getAnywhere(rda.default)
复制代码


图8|各主成分对方差解释度。fviz_eig(res.pca)。
4.2 PCA分析-样本间
绘制的PCA图中,响应变量以箭头表示,样本和定性变量以形状点表示。样本之间的距离为欧几里德距离(Euclidean distance,关联L2范数:中心化后变量值平方和的平方根:https://zhuanlan.zhihu.com/p/26884695)。样本间连线的长度表示样本间差异,长度越短,表示样本越相似,反之,则表示样本越不相似。如图6所示,NT+10-20 cm的三个重复样本间差异,就比NT+20-30 cm的三个重复样本间的差异更大。
  1. # NT+10-20 cm重复样本和NT+20-30 cm重复样本距离比较
  2. b = rownames(env[c(env$tillage == "NT" & env$depth == "10-20 cm"),])
  3. c = rownames(env[c(env$tillage == "NT" & env$depth == "20-30 cm"),])
  4. res.pca$ind$dist[b] # 三个重复距离差异很大
  5. res.pca$ind$dist[c] # 三个重复距离相近
复制代码


图9|NT+10-20 cm的三个重复样本间距离较大。

4.3 PCA分析-样本与响应变量之间
如图6所示:原点(箭头起始点)表示响应变量平均值位置,如果样本的垂直投影点在箭头指向的方向,表示响应变量在该样本的值大于平均值,否则,就小于平均值。比如:NT+0-10 cm样本中的TK变量的值应该是大于平均值,且大于CK+20-30 cm样本中的TK值的。
  1. # NT+0-10cm样本中的TK值
  2. env[env$tillage == "NT" & env$depth == "0-10 cm","TK"] # 三个样本中都大于均值
  3. # TK均值
  4. mean(env$TK)
  5. # CK+20-30 cm样本中的TK值
  6. env[c(env$tillage == "CK" & env$depth == "20-30 cm"),"TK"] # 三个样本中都小于均值
复制代码


图10|PCA双序图中样本与响应变量之间关系解读。
4.4 PCA分析-样本与主成分之间
样本与主成分间的关系主要是1)主成分对样本的代表性及2)样本对主成分贡献度。在每个变量及样本的权重相同,变量方差都相等的情况下,主成分对样本的代表性越高,该样本对主成分的贡献度越高。代表性和贡献度数据在绘图时都可以用于设置样本点颜色,点大小,也可用于对样本进行聚类分析。
  1. # 4.4.1 主成分对样本的代表性-主成分能代表的响应变量的比例
  2. res.pca$ind$coord # 样本坐标=样本(左)奇异值*sqrt(eig)
  3. res.pca$ind$cos2 # 用于表示样本能被PC代表的程度,cos2值越大,PC越能代该样本(我认为可理解为一种相关程度R^2)。
  4. ## 每个样本在所有主成分中的cos2值的和为1
  5. rowSums(res.pca$ind$cos2) # 所有样本为1


  6. ## 基于cos2绘制PC所能代表的样本部分占该样本的比例(the quality of the individuals on the factor map)。
  7. library(corrplot)
  8. pdf("ind.heatmap.pdf",height = 12,width = 6,family="Times")
  9. corrplot(res.pca$ind$cos2,col = mypal[1:4],is.corr=FALSE)
  10. dev.off()
  11. ### fviz_cos2绘制PC1和PC2所能代表的样本部分占该样本的比例(quality)。
  12. fviz_cos2(res.pca, choice = "ind",axes = 1:2)# 绘图结果是柱形图"ind"换成“var”可绘制变量与PC相关性。


  13. # 4.4.2 每个样本对PC的贡献度
  14. res.pca$ind$contrib # 样本对主成分的贡献度(%)=res.pca$ind$cos2/apply(res.pca$ind$cos2,2,sum)*100%。


  15. ## PC中所有样本贡献度的和为100
  16. colSums(res.pca$ind$contrib) # 所有PC值为100%


  17. ## 基于contrib绘制样本对PC的贡献度
  18. ### corrplot绘制
  19. corrplot(res.pca$ind$contrib,col = mypal[1:4],is.corr=FALSE)
  20. ### fviz_contrib()绘制样本对PC1和PC2的贡献度
  21. ind.contrib = fviz_contrib(res.pca, choice = "ind", axes = 1:2)#红色虚线表示预期平均贡献度。
  22. ind.contrib
  23. ggsave(plot=ind.contrib,"ind.contrib.pdf",device = "pdf",height = 8,width = 16,family="Times")
复制代码


图11|主成分对样本的代表性ind.heatmap.pdf。cos2值越大,PC越能代该样本。图中显示PC1对L.SG.a.1、L.SG.c.1、L.SG.c.2、L.SG.c.3、L.MG.a.1、L.MG.a.2和L.MG.c.3等样品的代表性较高。


图12|样本对PC1-2的贡献度。红色虚线表示预期平均贡献度(样本贡献度相同的情况下,预期平均贡献度为(1/36*Eig1+1/36*Eig2)/(Eig1+Eig2)*100%=1/36*100%=2.78%)。样本对PC1-2贡献度高于预期贡献度,可以认为这些样本对PC1-2比较重要。
4.5 PCA分析-响应变量间
  1. # 4.5.1 直接使用cor计算响应变量相关性系数
  2. cor = data.frame(round(cor(env[4:14]),2))
  3. cor
  4. cor["Nitrate"]


  5. # 4.5.2 绘制响应变量PCA-可以判断响应变量相关关系的正负。
  6. ## factoextra包的fviz_pca_var()函数用于绘制响应变量PCA图
  7. var.p = fviz_pca_var(res.pca,
  8.              axes = c(1, 2),
  9.              repel=TRUE, # 防止变量标签重叠
  10.              alpha.var = 0.8, # alpha.var=“cos2”也可以根据cos2值大小更换透明度。
  11.              col.circle = "grey50",
  12.              col.var = "cos2", # cos2值越低的变量,颜色越黄。
  13.              gradient.cols = c(mypal[1],'#FF7F00',mypal[2]),
  14.              select.var=NULL) # select.var可根据变量名称、cos2和contrib值筛选想要展示的响应变量。
  15. ggsave(plot=var.p,"var.p.pdf",height = 8,width = 8)
复制代码


图13|响应变量PCA图。具有正相关关系的响应变量间的夹角小于90度。夹角大于90度的响应变量间关系为负相关。夹角为90度表示没有相关关系。变量与原点之间的距离表示变量在PCA中的作用,距离原点越远,越能在PCA图中得到很好的表示。所有响应变量的方差都为1,则该圆的半径为1。图中OC与OM的起点与终点一样,两者相关性为1。
4.6 PCA分析-响应变量与主成分之间
‍解释与4.4 样本与主成分之间关系差不多。可以根据主成分对变量的代表性和变量对主成分的贡献度或累计贡献度筛选主成分和变量。PCA和RDA是线性排序模型,变量坐标是变量数据对样本坐标的回归(regression relations)。因此变量坐标就是一个斜率参数,连接原点与变量点,可以得到变量在每个PC拟合多度值的变化速率。
  1. # 4.6.1 变量与主成分之间的相关性-主成分对响应变量的代表性(解释与4.4部分差不多)
  2. res.pca$var$coord # 箭头方向表示变量变化的方向。所有响应变量(这里是环境因子)的坐标,等于载荷值(主成分能代表的响应变量的方差)*特征根(主成分的标准误差)。


  3. res.pca$var$cor # 表示主成分与响应变量间的相关性。样本权重一样时,每个响应变量与轴的相关性=响应变量坐标/(变量值的平方和/样本数的平方根)。
  4. res.pca$var$cos2 # 主成分与响应变量相关性的平方,表示主成分对响应变量的代表性,等于res.pca$var$cor^2=res.pca$var$coord^2。
  5. rowSums(res.pca$var$cos2) # 对于所有变量方差为1的数据,某个变量所有主成分cos2之和为1。


  6. ## 响应变量与样本PC坐标的pearson相关性系数是响应变量的坐标。
  7. res.pca$var$coord["pH",]
  8. cor(res.pca$ind$coord[,"Dim.1"],env["pH"]) # 计算pH与样本PC1坐标的相关系数
  9. ##相关系数为-0.7542129,与cor值一致。


  10. ##可视化主成分对变量的代表性-PC代表的部分响应变量占原始变量的比例(quality)-cos2
  11. library(corrplot)
  12. pdf("var.cos2.pdf")
  13. corrplot(res.pca$var$cos2,
  14.          method="number",tl.col = "blue",
  15.          col.lim = c(0, 1),is.corr=FALSE)
  16. dev.off() # 所有主成分热图
  17. fviz_cos2(res.pca, choice = "var", axes = 1:2) # 前两个轴柱形图
  18. ##cos2越大表示主成分代表的该变量成分越多,主成分能很好的代表该变量。反之,该变量则不能很好的被主成分代表。
复制代码
注:如图13所示,AHN与pH相比,AHN与PC1的夹角(只看锐角)更小,其cos绝对值更大,表明AHN与PC1的相关性更强,相关性的正负,由箭头所在象限决定,从第一象限开始算夹角,箭头在二、三象限的角的cos值为负,一、四象限为正。看响应变量与PC2的相关性正负,则需要把第二象限作为0度,这样算起来很不方便。可以直接看响应变量箭头方向与轴方向(负→正)的一致性,以判断响应变量与PC相关性的正负。如AHN与PC1坐标从低到高的方向一致,与PC2方向相反,表明AHN与PC1正相关,与PC2负相关;pH箭头指向的方向与PC1的方向相反,与PC2一致,表明pH与PC1相关性为负,与PC2相关性为正。


图14|响应变量与样本PC坐标的pearson相关性系数是响应变量的坐标。可用于预测附加定性变量的坐标。


图15|PC对响应变量的代表性var.cos2.pdf。cos2越大表示主成分能很好的代表该变量。反之,该变量则不能很好的被主成分代表。对于所有变量方差为1的数据,某个变量所有主成分cos2之和为1。如果一个变量在PC1和PC2中的cos2和为1,则箭头终点位于圆周长上(图13)。
  1. # 4.6.2 变量对主成分的贡献度-寻找主成分主要代表的那些原始变量的方差,
  2. ##也可称为主成分可能的特征。变量在解释给定主成分的变异性方面所起的作用用百分比表示。
  3. contrib = data.frame(res.pca$var$contrib)
  4. contrib# 显示响应变量对每个主成分的方差解释度(%)。计算方式:所有样本对主成分的贡献度(%)=(res.pca$ind$coord^2)*100/eig^2/样本数
  5. ##表明响应变量与主成分的相关性系数的绝对值越大,对该主成分的贡献越大。


  6. ## 可视化所有主成分中响应变量的贡献度
  7. library(grDevices)
  8. #hcl.pals() # 查看存在的颜色板名称
  9. col1 = hcl.colors(12, "YlOrRd", rev = TRUE)
  10. #show_col(col1) # show_col()不能展示colorRampPalette()产生的颜色集
  11. pdf("var.contrib.pdf",height = 14,width = 7,family="Times")
  12. corrplot(res.pca$var$contrib,
  13.          method = c("pie"),tl.col = "black",col = col1,
  14.          col.lim = c(0, 100), #设置bar颜色区间
  15.          is.corr=FALSE)
  16. dev.off()


  17. ## 可视化对PC1和PC2合计贡献度最高的6个响应变量
  18. fviz_contrib(res.pca, choice = "var",
  19.              axes = 1:2, top = 6, # 贡献度高于预期贡献度
  20.              sort.val = c("desc"),
  21.              ggtheme = theme_bw())
  22. ggsave("PC1-2.contrib.pdf",height = 8,width = 12,family="Times")


  23. ## 绘制响应变量对PC1的贡献度和累积贡献曲线。
  24. ##超过预期贡献度或合计贡献度超过一定阈值的变量,可以认为是该主成分的维度特征,
  25. ##可以将主成分跟已知的环境条件联系起来,看是否能代表某一环境梯度。
  26. d1.contr = contrib[statnet.common::order(contrib[1],decreasing = TRUE),]
  27. rownames(d1.contr) = factor(rownames(d1.contr),levels = rownames(d1.contr))
  28. d1.contr$env = factor(rownames(d1.contr),levels = rownames(d1.contr))
  29. d1.contr # 数据表按PC1降序排列


  30. ### 使用twoord.plot()函数绘制相同X轴,双Y坐标轴图
  31. library(plotrix)
  32. pdf("twoord.pdf",height = 8,width = 14,family = "Times")
  33. twoord.plot(ly=d1.contr$Dim.1,type=c("bar","l"),
  34.             rx=factor(d1.contr$env,levels=c(d1.contr$env)),
  35.             ry=cumsum(d1.contr$Dim.1),
  36.             ylab="percentage of Variable contribution (%)",
  37.             rylab="Cumulative percentage of variable contribution (%)",
  38.             xlim=NULL,
  39.             xlab="Environmental Factor",
  40.             lylim=c(0,max(d1.contr$Dim.1)+1),
  41.             lcol=mypal[1],rcol=mypal[4],
  42.             xtickpos=1:11,
  43.             xticklab=factor(d1.contr$env,levels=c(d1.contr$env)),
  44.             halfwidth=0.2,
  45.             main="Contribution of variables to PC1",
  46.             do.first="plot_bg("white");abline(h=1/11*100,col="red",lty=2)"
  47.             )
  48. dev.off()


  49. ### lattice绘图系统doubleYScale()函数绘制双Y轴图
  50. library(latticeExtra)
  51. library(lattice)
  52. O1=barchart(Dim.1 ~ env,d1.contr,
  53.             ylab="percentage of Variable contribution (%)",
  54.             )
  55. O2=xyplot(cumsum(Dim.1) ~ env,d1.contr,
  56.           type="l",
  57.           ylab = "Cumulative percentage of variable contribution (%)")
  58. pdf("cumsum_PC1.pdf",width = 14,height = 8,family="Times")
  59. doubleYScale(O1,O2,add.axis=TRUE,under=FALSE,
  60.              text=c("Cumulative percentage of variable contribution (%)"),
  61.              add.ylab2 = TRUE,cex=1,use.style = FALSE)
  62. dev.off() # 可参考[url]https://www.jianshu.com/p/5885aaeda6c1[/url]添加参考线
复制代码


图16|响应变量对PC的贡献度var.contrib.pdf。


图17|对PC1-2的贡献度最高的6个响应变量PC1-2.contrib.pdf。红色虚线表示变量对主成分贡献度的预期平均分布。如果变量对PC的贡献度是均匀分布的,则预期平均贡献度应为1/N(响应变量数目)=(1/11*Eig1+1/11*Eig2)/(Eig1+Eig2)*100%=1/11*100%。对某个主成分贡献度高于预期平均贡献度的变量,可以认为该变量对主成分的贡献非常重要。计算两个轴中某个变量的总贡献度,及在解释被2个主成分保留的变异时的计算方式是:contrib = [(C1 * Eig1) + (C2 * Eig2)]/(Eig1 + Eig2)。Eig特征值衡量的是每个PC所保留的变量变化量(方差),C是响应变量对PC的贡献度。


图18|响应变量对PC1的贡献度和累积贡献曲线twoord.pdf。沿X轴按贡献度对变量进行降序排序。红色虚线表示变量对主成分贡献度的预期平均分布。如果变量对PC的贡献度是均匀分布的,则预期平均贡献度应为1/N(响应变量数目)=1/11*100%。对某个主成分贡献度高于预期平均贡献度的变量,可以认为该变量对主成分的贡献非常重要(TN→pH)。红色曲线为响应变量对PC1的累计贡献曲线,图中显示TN→pH对PC1的贡献度~80%。超过预期贡献度或合计贡献度超过一定阈值的变量,可以认为是该主成分的维度特征。可以将主成分跟已知的环境条件联系起来,看是否能代表某一环境梯度。与方差解释度占比高的主成分(如PC1和PC2)相关的响应变量是解释数据集变异性的最重要因素。不与任何PC相关或与最后维度相关的变量都是贡献较低的变量,可以删除以简化分析。
  1. # 4.6.3 根据贡献度对响应变量进行聚类
  2. set.seed(123)
  3. res.km <- kmeans(res.pca$var$coord, centers = 3, nstart = 25) #k值设为3
  4. grp <- as.factor(res.km$cluster)
  5. fviz_pca_var(res.pca, col.var = grp,
  6.              palette = mypal[1:3],
  7.              legend.title = "Cluster")
  8. ggsave("kmeans.cluster.pdf")
复制代码


图19|根据贡献度对响应变量进行聚类结果kmeans.cluster.pdf。
dimdesc()可用于鉴定与给定主成分具有统计学意义相关性的响应变量。
  1. # 4.6.4 主成分维度描述
  2. ##基于对各主成分的因子分析,找出该主成分最能表征的定性(样本分类信息)和定量变量。
  3. dim = dimdesc(res.pca,
  4.               axes = 1:3,# 设置进行分析的轴
  5.               proba = 0.05) # proba设定统计阈值。
  6. dim$call # 设置信息
  7. dim$Dim.1$quanti #定量响应变量与PC1的相关性r值及p值。
  8. dim$Dim.1$quali #附加定性响应变量与PC1的决定系数R^2及p值。
  9. dim$Dim.2$quali
  10. ## R^2表示有多少PC中所保留的样本间差异可由定性变量的level不同所解释。
  11. dim$Dim.1$category # 定性变量中与PC高度相关的level
  12. dim$Dim.2$category


  13. ## FactoInvestigate包的Investigate()会输出Wilks test p-value结果
  14. ##      tillage        depth
  15. ## 1.222640e-08 2.047917e-05
  16. ## 会输出在双序图中能最佳解释样本间距离的定性变量,比如这里是tillage。

复制代码


图20|PCA维度描述。结果显示PC1相关性最高的变量为TN(dim$Dim.1$quanti)。depth在PC1上分离的最好(PC1能很好的表示depth分类样本间差异,或depth能很好的解释样本沿PC1的距离)(dim$Dim.1$quali)。depth中的0-10 cm和20-30 cm与PC1高度相关(dim$Dim.1$category)。tillage能很好的解释样本沿PC2的距离(dim$Dim.2$quali)。tillage中的NT和CK与PC2高度相关(dim$Dim.2$category)。
4.6 PCA分析-附加定性/定量样本
基于PCA分析结果,附加定性变量的的坐标是该分类中样本坐标的平均值。附加定性变量可以作为分组变量,用于绘图,帮助解释数据分析结果。附加定量变量的的坐标通过计算定量变量与主成分的相关性得到。
  1. # 4.6.1 附加定性变量的的坐标是该分类中样本坐标的平均值
  2. ## 提取样本PC1、PC2坐标
  3. sam=data.frame(res.pca$ind$coord)
  4. sam
  5. sam=data.frame(sam[1:2],env[1:3]) # 添加样本分组信息
  6. sam
  7. library(rstatix)
  8. res.pca$quali.sup$coord["CK",c("Dim.1","Dim.2")] # 仅小数点取位不一致
  9. sam %>% group_by(tillage) %>% get_summary_stats(type="mean")


  10. ## 计算附加定性变量的PC1和PC2坐标-运行时已指定,就不用再自己计算了。
  11. coord.condition <- sam %>%
  12.   as_data_frame() %>%
  13.   group_by(condition) %>%
  14.   dplyr::summarise(
  15.     PC1.mean = mean(Dim.1),
  16.     PC2.mean = mean(Dim.2),
  17.     PC1.sd = sd(Dim.1),
  18.     PC2.sd = sd(Dim.2)
  19.     )
  20. coord.condition


  21. ##  绘制带误差bar的condition定性/分类变量的PCA散点图
  22. ### 对condition进行分列获得tillage和depth分组
  23. data = separate(coord.condition,
  24.                 "condition",
  25.                 into = c("tillage","depth"),
  26.                 sep ="-",
  27.                 remove = FALSE)
  28. head(data)


  29. ### 绘图,用形状区分depth,颜色区分tillage
  30. data$condition
  31. p.bar = ggplot(data)+
  32.   geom_point(aes(x=PC1.mean,y=PC2.mean,
  33.                  shape=depth,color=condition)
  34.              ,size=5)+
  35.   scale_shape_manual(values = c(16,17,15),guide="none")+
  36.   ## 为了后面合并图例,这里使用condition定义颜色,
  37.   ###然后根据tillage设置4个颜色,每个颜色重复3次的颜色集。
  38.   scale_color_manual(values = rep(mypal[1:4],each=3),
  39.                      breaks = data$condition)+
  40.   ## 使用geom_segment()函数绘制误差bar。
  41.   ### PC1 error bar
  42.   geom_segment(aes(x=PC1.mean-PC1.sd,xend=PC1.mean+PC1.sd,
  43.                    y=PC2.mean,yend=PC2.mean,
  44.                    color=condition),
  45.                size=0.75,show.legend = FALSE,
  46.                arrow=arrow(angle = 90,ends = "both",
  47.                            length = unit(0.2, "cm")))+
  48.   ### PC2 error bar
  49.   geom_segment(aes(x=PC1.mean,xend=PC1.mean,
  50.                    y=PC2.mean-PC2.sd,yend=PC2.mean+PC2.sd,
  51.                    color=condition),
  52.                size=0.75,show.legend = FALSE,
  53.                arrow=arrow(angle = 90,ends = "both",
  54.                            length = unit(0.2, "cm")))+
  55.   geom_vline(xintercept = 0,colour = "grey50",size=0.5)+
  56.   geom_hline(yintercept = 0,colour = "grey50",size=0.5,)+
  57.   ## 添加轴标签
  58.   xlab(paste("PC1 (",PC1,"%)",sep = ""))+
  59.   ylab(paste("PC2 (",PC2,"%)",sep = ""))+
  60.   theme_bw()+
  61.   ## 为bottom和left轴添加箭头
  62.   ## legend.position和legend.justification结合设置图例位置
  63.   theme(legend.position = c(1,1),
  64.         legend.justification=c(0,1),
  65.         legend.background = element_blank(),
  66.         plot.margin =margin(0.2,2,0,0.2,unit ='cm'), # 设置图边距,依次对应top, right, bottom和left。
  67.     axis.line.x.bottom = element_line(
  68.     arrow=arrow(type="closed",angle = 25,
  69.                 length = unit(0.3, "cm"))),
  70.     axis.line.y.left = element_line(
  71.     arrow=arrow(type="closed",angle = 25,
  72.                 length = unit(0.3, "cm"))))+
  73.   ## 合并颜色和形状图例
  74.   guides(color =guide_legend(title = "Treatments",
  75.                              override.aes = list(shape =c(16,17,15),
  76.                             label=TRUE)))
  77. ggsave(plot = p.bar,"p.bar.pdf",device="pdf",
  78.        height = 8,width=10,family="Times")
复制代码


图21|预测的附加变量condition的PC坐标,coord.condition。


图22|带error bar的condition定性/分类变量PCA散点图,p.bar.pdf。
  1. # 4.6.2 附加定量变量预测
  2. ## 随机模拟2个有36个样本的变量
  3. set.seed(123)
  4. var.supp=data.frame(row.names = rownames(env),
  5.                 ran1=sample(10:30,36,replace = TRUE), # 有放回的在10-30范围随机抽取36个整数。
  6.                 ran2=runif(36,20,30)) # 在20-30范围内随机产生36个小数。
  7. var.supp
  8. ## 预测定量变量坐标
  9. quanti.coord <- cor(var.supp, res.pca$ind$coord)
  10. quanti.coord
  11. quanti.cos2 <- quanti.coord^2 # cos2值
  12. quanti.cos2
  13. ## factoextra的fviz_add()将预测的定量变量加入到PCA图中
  14. p <- fviz_pca_var(res.pca)
  15. ?fviz_add # 查看帮助信息
  16. fviz_add(p, quanti.coord, color ="blue",
  17.          linetype = "dashed",
  18.          geom="arrow")
  19. ggsave("var.supp.pdf")
复制代码


图23|添加附加定性变量的PCA变量图,var.supp.pdf。
4.7 PCA分析-附加样本

基于PCA分析结果,使用stats包的predict()可以预测附加样本的坐标和cos2。
  1. # 4.7.1 随机模拟2个样本
  2. ind.supp = t(data.frame(row.names = colnames(env[4:14]),
  3.                 ran1=runif(11,10,30),
  4.                 ran2=runif(11,20,30)))
  5. ind.supp


  6. # 4.7.2 predict()预测附加样本坐标和cos2
  7. #?predict # 查看帮助信息
  8. ind.sup.coord <- predict(res.pca, newdata = ind.supp)
  9. ind.sup.coord$coord
  10. ind.sup.coord$cos2


  11. # 4.7.3 将预测的附加样本添加到PCA图中
  12. p <- fviz_pca_ind(res.pca, repel = TRUE)
  13. fviz_add(p, ind.sup.coord$coord, color ="blue")
复制代码


图24|预测的附加样本的坐标信息。预测过程经过两步:1)对附加变量进行与初始变量一致的中心化和标准化;2)标准化后每一个变量与所有PC的变量载荷值相乘,每个样本得到一个样本数*PC数矩阵,然后对每个样本中PC包含的变量值求和,得到附加样本PC坐标。
五、 PCA(因子)分析前变量筛选
根据PCA分析的目的,有时专家审稿会要求对原始变量进行Bartlett's test of sphericity(球形检验)和Kaiser-Meyer-Olkin Measure of Sampling Adequacy(KMO采样充分性检验),检验数据是否合适进行PCA(因子)分析,还要求对变量进行筛选(communality<0.5)的变量。在SPSS中做PCA分析时,在描述选项勾选KMO和Bartlett球形检验即可。在R中相先对原始变量做相关性(度量不相同)或协方差(原始变量的度量相同或标准化后数据)分析,再使用psych包中的cortest.bartlett()和KMO()进行上述检验。因子分析与排序分析的差异,数量生态学书中有讲解(看了也会忘),主要还是根据自己的需要选择分析方法。参考:https://zhuanlan.zhihu.com/p/49481213
  1. library(dplyr)
  2. # 5.1 相关性和协方差矩阵计算
  3. cor=cor(env[4:14],method="pearson")
  4. cor
  5. cov=cov(scale(env[4:14])) # 原始变量标准化之后的协方差分析与原始变量的相关性分析结果一致。
  6. cov
  7. #cov2cor(cov) # 将协方差矩阵转换为相关性矩阵
复制代码
注:协方差:衡量变量间相对于各自自身期望的变化趋势的统计指标,比如变量X大于自身期望时变量Y也大于自身期望,此时变量XY的协方差为正值,反之则为负值。相关系数:用以反映变量之间相关关系密切程度的统计指标。相关系数也可以看成是一种剔除了两个变量量纲影响、标准化后的特殊协方差,它消除了两个变量变化幅度的影响,而只是单纯反应两个变量每单位变化时的相似程度。参考:https://www.jianshu.com/p/f33f3f760de3

  1. # 5.2 Bartlett's test of sphericity
  2. #?cortest.bartlett # 查看函数帮助信息
  3. cor.bar = cortest.bartlett(cor,n = nrow(sam))# n表示样本数
  4. ##diag=TRUE,当使用协方差矩阵分析时,将矩阵的对角线值都设置为1。
  5. cov.bar = cortest.bartlett(cov,n = nrow(sam),diag=TRUE)
  6. cor.bar
  7. cov.bar
  8. ## p<0.05,表示变量间具有相关性,可以用于进行因子或主成分分析。
复制代码


图25|球形检验结果。球形检验的的null假设是相关矩阵中对角线系数为1,非对角线系数为0,既变量间是独立的。备择假设是变量间是相关的。因此p<0.05,接受备择假设,数据适合用于进行因子或主成分分析。
  1. # 5.3 Kaiser-Meyer-Olkin Measure of Sampling Adequacy
  2. ?KMO # 查看函数帮助信息和结果解释
  3. kmo = KMO(cor)
  4. kmo # KMO统计量是衡量总体和每个变量的抽样充分性(MSA)的指标>0.5
复制代码


图26|抽样充分性检验结果。因子分析中各变量的偏相关是排除所有其它变量的影响后,这些变量之间的相关性组成。KMO统计量就是描述相对于原始(零阶)相关性的偏相关关系有多小。如果变量具有公共因子,则偏相关系数应该比较小,KMO应接近1.0。当相关矩阵等于偏相关矩阵时,KMO测度等于0.5。KMO大于0.8可以认为成分或因子分析对这些变量是有用的。这通常发生在大多数原始相关性为正的时候.当大多数原始相关性为负值时,KMO值小于0.5,此时需要采取补救行动,要么删除低于0.5的变量,要么添加与KMO低于0.5变量有关的其它变量。Kaiser在1974年给出了经验原则(?KMO查看):0.9以上:适合性很好; 0.8~0.9:适合性良好;0.7~0.8:适合性中等;0.6~0.7:适合性一般;0.5~0.6:适合性不好;0.5以下:数据不适合。  ‍

EcoEvoPhylo :主要分享微生物生态和phylogenomics的数据分析教程。扫描二维码,关注EcoEvoPhylo。让我们大家一起学习,互相交流,共同进步。
学术交流QQ群 | 438942456
学术交流微信群 |  jingmorensheng412
加好友时,请备注姓名-单位-研究方向。                 

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
32
注册时间
2016.1.8
在线时间
655 小时

发表于 2021.12.15 08:48:21 | 显示全部楼层
辛苦辛苦,
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
32
注册时间
2016.1.8
在线时间
655 小时

发表于 2021.12.15 17:41:35 | 显示全部楼层
很不错哦
新的一天加油!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2021.5.17
在线时间
4 小时

发表于 2021.12.16 19:02:54 | 显示全部楼层
我以为这个地是灌水的,没想到真不错啊。
回复 支持 反对

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2017.6.6
在线时间
183 小时

灌水之王


发表于 2021.12.20 08:44:46 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2021.6.22
在线时间
64 小时

发表于 2021.12.23 13:53:34 | 显示全部楼层
感谢分享
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2021.3.1
在线时间
3 小时

发表于 2021.12.27 21:24:34 | 显示全部楼层
感谢分享,正好需要做pca分析,谢谢!
新的一天加油!
回复 支持 反对

使用道具 举报

帝王蝶

Rank: 4

主题
3
注册时间
2016.9.2
在线时间
64 小时

发表于 2022.1.4 22:38:32 | 显示全部楼层
非常详细
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2021.6.22
在线时间
64 小时

发表于 2022.1.6 16:54:01 | 显示全部楼层
学习
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2017.11.3
在线时间
351 小时

发表于 2022.1.11 15:12:04 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表