R、冗余分析（RDA）、ggplot2、置信椭圆

R、冗余分析（RDA）、ggplot2、置信椭圆

 在生态环境领域中（实际中，其他专业也用到），冗余分析（RDA）是我们常用的分析方法，分析目的为“解释变量”对“响应变量”的影响情况。类似RDA的方法，还有CCA。这里以RDA为例→数据处理、分析过后，我们需要对结果进行可视化，R语言ggplt2程序包无疑是可视化神器，然而，怎样利用ggplot2对RDA结果进行可视化，需要我们对RDA结果进行了解，提取需要展示的元素。许多论文中在排序分析的图中添加了置信椭圆，ggplot2自带函数目前支持每个分组样本量≥4时画椭圆，这里我们使用gglayer包的geom_ord_ellipse函数给样方添加置信椭圆（支持样本量≥3）。

 library(vegan)
 library(plyr)
  library(gglayer)
 library(ggplot2)
 library(ggrepel)
fc=read.csv("D:\\R\\factor.csv",header = T,row.names = 1)#读取解释变量数据
sp=read.csv("D:\\R\\sp.csv",header = T,row.names = 1)#读取响应变量数据
spp=decostand(sp,method = "hellinger")#对响应变量做转化
fcc=log10(fc)#对解释变量做转化
uu=rda(spp~.,fcc)#RDA分析
ii=summary(uu)  #查看分析结果
sp=as.data.frame(ii$species[,1:2])*2#可根据出图结果，对画图数据做一定的放大或缩小，下同
st=as.data.frame(ii$sites[,1:2])
yz=as.data.frame(ii$biplot[,1:2])
grp=as.data.frame(c(rep("a",3),rep("b",3),rep("c",4),rep("d",6)))#根据样方类型分组，“a”有3个样本，“b”有3个样本……，共16个。注意样本的顺序和个数！
colnames(grp)="group"
ggplot() +
  geom_point(data = st,aes(RDA1,RDA2,shape=grp$group,fill=grp$group),size=4)+
  scale_shape_manual(values = c(21:25))+
  geom_ord_ellipse(aes(st$RDA1,st$RDA2,color=grp$group,group=grp$group),###注意，是在这里添加椭圆
                   ellipse_pro = 0.68,linetype=3,size=1)+###注意，是在这里添加椭圆
  geom_segment(data = sp,aes(x = 0, y = 0, xend = RDA1, yend = RDA2), 
               arrow = arrow(angle=22.5,length = unit(0.35,"cm"),
                             type = "closed"),linetype=1, size=0.6,colour = "red")+
  geom_text_repel(data = sp,aes(RDA1,RDA2,label=row.names(sp)))+
  geom_segment(data = yz,aes(x = 0, y = 0, xend = RDA1, yend = RDA2), 
               arrow = arrow(angle=22.5,length = unit(0.35,"cm"),
                             type = "closed"),linetype=1, size=0.6,colour = "blue")+
  geom_text_repel(data = yz,aes(RDA1,RDA2,label=row.names(yz)))+
  labs(x=paste("RDA 1 (", format(100 *ii$cont[[1]][2,1], digits=4), "%)", sep=""),
       y=paste("RDA 2 (", format(100 *ii$cont[[1]][2,2], digits=4), "%)", sep=""))+
  geom_hline(yintercept=0,linetype=3,size=1) + 
  geom_vline(xintercept=0,linetype=3,size=1)+
  guides(shape=guide_legend(title=NULL),color=guide_legend(title=NULL),
         fill=guide_legend(title=NULL))+
  theme_bw()+theme(panel.grid=element_blank())

  RDA探索完毕，上面的图比较骚气，但，不能用，因为需要检验模型是否显著、解释变量是否存在共线性，在论文写作中应删掉某些图层，使图更加美观、更简洁。这里只是简单的举个例，在进行约束排序分析之前，我们先要检查数据，是符合RDA还是CCA，网上有很多例子，这里不再赘述。当RDA或者CCA不能很好的解决我们的问题时，我们需要结合其它分析方法，如gbm、RF能求出相对重要性等，也可利用SEM去探索。原则上，不同的数据需要用不同的方法探索，选择一种最理想的结果。
 刚接触R的朋友，可能会因为R的“难”入门而选择较容易的CANOCO软件，当然后者也是生态环境领域的数据分析神器，但是，当我们想要个性化分析、个性化出图时，目前后者无法满足，R，应是首选。实验设计、实验过程、数据处理及分析、绘图、写作等是科研的必然过程，为促进相互进步、资源共享，我们创建了学术交流QQ群：335774366。欢迎有兴趣的朋友加入→指导。
声明：以上内容仅为作者个人理解，有不对的地方，欢迎指正。