How do I calculate the percent overlap of my 95% confidence interval ellipses in my NMDS plots? I used the 'vegan' package and my code is posted below.
The data sheet columns are by species and the rows are by plot # 1-16 which are categorized into four groups as seen in the NMDS attached below. Here is a few rows and columns as an example. Ten column examples are here: https://www.dropbox.com/s/gf5llxm986i3kld/mydata.txt?dl=0 I would love to sort by "Treatment" so that I have four ellipses
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
The code used to produce the NMDS plots is:
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"))
`
Below is a resultant NMDS plots.
