首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >R语言实现PCOA分析

R语言实现PCOA分析

作者头像
一粒沙
发布2019-12-19 11:30:33
10.5K1
发布2019-12-19 11:30:33
举报
文章被收录于专栏:R语言交流中心R语言交流中心

大家对主成分分析(principal components analysis, PCA) 都很熟悉,但是今天我们来介绍下主坐标分析(principal coordinate analysis, PCoA)。那么这两个差了个o字母具体有什么区别?首先PCA是常用的降维算法;利用线性变换,将数据变换到一个新的坐标系统中;然后再利用降维的思想,使得任何数据投影的第一大方差在第一个坐标(称为第一主成分)上,第二大方差在第二个坐标(第二主成分)上。这种降维的思想首先减少数据集的维数,同时还保持数据集的对方差贡献最大的特征,最终使数据直观呈现在二维坐标系。PCoA主要是探索数据相似度或者相异度可视化方法。可呈现研究数据相似性或差异性的可视化坐标,是一种非约束性的数据降维分析方法,可用来研究样本群落组成的相似性或相异性。其实通俗的讲,PCA主要是基于原始数据矩阵的降维;PCoA主要是基于样本的原始数据计算出来的距离矩阵的降维。如果样本数目比较多,而物种数目比较少,那肯定首选PCA;如果样本数目比较少,而物种数目比较多,那肯定首选PCoA。

接下来我们看下在R中如何去实现,首先安装ape包和vegan包,联合使用才能达到最终的目的。包的安装我们就不赘述了,其在CRAN平台,直接install.packages()。

首先是数据的导入,我们利用vegan自带的数据dune。具体的数据集的构成大家可以直接在包的信息中去看。接下来我们首先基于dune数据构造距离矩阵,需要用到的函数vegdist。

vegdist(x, method="bray",binary=FALSE, diag=FALSE, upper=FALSE,

na.rm = FALSE, ...)

其中主要的参数:

Method 支持了目前大部分的距离计算函数。"manhattan", "euclidean", "canberra","clark", "bray", "kulczynski","jaccard", "gower", "altGower","morisita", "horn", "mountford", "raup","binomial", "chao", "cao" or"mahalanobis"。

Diag 是否显示对角线的值。

Upper 是否显示对角线以上的值

library(vegan)
data(dune)
data(dune.env)
dune.dist <- vegdist(dune)#构造距离矩阵。 

接下来就是利用ape中的pcoa函数获取PCOA分析结果。当然也可以应用我们R自带的函数cmdscale。

pcoa(D, correction="none",rn=NULL)

其中主要参数:

D 不用多说就是距离矩阵

Correction 主要指需不需进行校正。

Rn 主要是指的每个样本的名称标签,可以进行自定义

接下来我们看下实例:

res=pcoa(dune.dist,correction="none")
 
head(res$value)

其中主要的值是特征值的一些相关的转换的值。

head(res$vectors)

其中主要是和PCA中主成分类似的柱坐标的值,进行了排序展示,一般选择前两个绘制二维可视散点图。

biplot(res)#可视化PCOA 的结果

为了进一步完善我们的可靠性,我们还可以利用vegan中的ANOSIM相似性分析是一种非参数检验,用来检验组间(两组或多组)差异是否显著大于组内差异,从而判断分组是否有意义。

anosim(x, grouping, permutations = 999, distance = "bray", strata = NULL, parallel = getOption("mc.cores"))
其中的参数我们不再赘述,直接看下实例:
dune.ano <- with(dune.env, anosim(dune.dist, Management))
summary(dune.ano)

plot(dune.ano)

至此,我们的PCOA的分析过程可以实现,那么如何优化我们输出的可视化图像,我们需要用到ggplot2这个包可以对我们的值进行更加友好的可视化。直接看实例:


library(ggplot2)
P= dune.ano$signif
R=round( dune.ano$statistic,3)
data=res$vectors[,1:2]
data=data.frame(data)
colnames(data)=c('x','y')
eig=as.numeric(res$value[,1])
 
type=dune.env$Management
 
data=cbind(data,type)
 
p =ggplot(data, aes(x=x, y=y, color=type))+
 
geom_point(alpha=.7, size=2) +
 
labs(x=paste("PCoA 1 (",format(100* eig[1] / sum(eig), digits=4), "%)", sep=""),
 
y=paste("PCoA 2 (", format(100*eig[2] / sum(eig), digits=4), "%)", sep=""),
 
title="PCoA")
 
p +theme_classic()+stat_ellipse(type ="t", linetype = 2)+
 
annotate("text",x=0.004,y=-0.15,parse=TRUE,size=4,label=paste('R:',R),family="serif",fontface="italic",colour="darkred")+
 
annotate("text",x=0,y=-0.17,parse=TRUE,size=4,label=paste('p:',P),family="serif",fontface="italic",colour="darkred")
本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2019-12-09,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 R语言交流中心 微信公众号,前往查看

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

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

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档