3

如何计算 NMDS 图中 95% 置信区间椭圆的重叠百分比？我使用了“vegan”包，我的代码发布在下面。

数据表的列按物种，行按地块#1-16 分为四组，如下面的 NMDS 所示。这里以几行和几列为例。十列示例在这里：https ://www.dropbox.com/s/gf5llxm986i3kld/mydata.txt?dl=0 我很想按“治疗”排序，这样我就有四个椭圆

  Plot   BHW   YHW   Clown   Tobacco   Unwrasse   DuskDam   YTD   Beau   BridG   GoldG   
 ------ ----- ----- ------- --------- ---------- --------- ----- ------ ------- -------- 
     1    20     0       0         0          0         0     0      0       2        0  
     2    15     0       0         0          0         0     0      0       4        3  
     3    10     0       0         0          0         0     0      0       3        3  
     4     6     0       0         0          0         0     0      0       3        2  
     5    33     0       0         0          0         0     0      0       0        7  
     6     7     0       0         0          0         0     0      0       1        0  
     7    27     1       0         0          0         0     0      0       2        3  
     8    22     0       0         0          0         0     0      0       0        3

用于生成 NMDS 图的代码是：

abund_fish_06_01to_06_03<-abund_fish[abund_fish$Date>="2016-06-01" &abund_fish$Date<="2016-06-03",]
    abund_fish_06_01to_06_03_rn<- abund_fish_06_01to_06_03[,2:62] #to remove date and plot

sol <- metaMDS(abund_fish_06_01to_06_03_rn, distance = "bray", k = 3, trymax = 100) 
NMDS = data.frame(MDS1 = sol$points[,1], MDS2 = 
sol$points[,2],group=MyMeta$Treatment)
NMDS.mean=aggregate(NMDS[,1:2],list(group=NMDS$group),mean)
veganCovEllipse<-function (cov, center = c(0, 0), scale = 1, npoints = 100) 
{
theta <- (0:npoints) * 2 * pi/npoints
Circle <- cbind(cos(theta), sin(theta))
t(center + scale * t(Circle %*% chol(cov)))
}

df_ell <- data.frame()
for(g in levels(NMDS$group)){
df_ell <- rbind(df_ell, cbind(as.data.frame(with(NMDS[NMDS$group==g,],veganCovEllipse(cov.wt(cbind(MDS1,MDS2),wt=rep(1/length(MDS1),length(MDS1)))$cov,center=c(mean(MDS1),mean(MDS2)))))
                            ,group=g))}
nmds6163bw<-ggplot(data = NMDS, aes(MDS1, MDS2)) + geom_point(aes(color = 
group, shape = group), stroke =1.25,size=2.5) +
theme_bw()+
theme(axis.text=element_text(size=14),
    axis.title=element_text(size=16,face="bold")
 )+
theme(legend.text = element_text(size = 14)
)+
theme(legend.title = element_text(size = 14)
)+
geom_path(data=df_ell, aes(x=MDS1, y=MDS2, color = group, linetype=group), size =1)+
scale_shape_manual(name = "Treatment", labels = c("10m Control", "20m Control", "March Outplant", "June Outplant"), values = c(1, 19, 15,0)) +
scale_colour_manual(name = "Treatment", labels = c("10m Control", "20m Control", "March Outplant", "June Outplant"), values = c("grey50", "grey50", "black", " black "))+
scale_linetype_manual(name = "Treatment", labels = c("10m Control", "20m Control", "March Outplant", "June Outplant"), values = c(2, 1, 1, 2))
nmds6163bwnl<-nmds6163bw+theme(legend.position="none") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black"))

`

下面是生成的 NMDS 图。

NMDS 图

4

0 回答 0