后续行动:将 vegan 包中的 ordiellipse 函数绘制到在 ggplot2 中创建的 NMDS 图上

follow-up: Plotting ordiellipse function from vegan package onto NMDS plot created in ggplot2

我正在尝试做一些类似于旧 post 的事情:plotting - original post

对于我的分析,我感兴趣的是不同的哺乳动物宿主是否有不同的跳蚤群落。我链接到的原始 post 有 2 种不同的省略号解决方案。我的问题是,当我 运行 第一个解决方案和一般解决方案时,我得到的图看起来截然不同,但我认为它们应该非常相似。下面是我的代码。

我的问题是:我是不是做错了什么或者哪个代码生成了正确的数字?或者我应该使用更好的新代码来显示不同寄主物种的跳蚤群落差异吗?

谢谢, 阿曼达

Link 到 ARG_comm 数据 Link 到 ARG_env 数据

library(vegan)
library(BiodiversityR)
library(MASS)

comm_mat <- read.csv("d:/fleaID/ARG_comm.csv",header=TRUE)
env <- read.csv("d:/fleaID/ARG_env.csv",header=TRUE)

library(ggplot2)
sol <-metaMDS(comm_mat)
MyMeta=env

#originalresponse
NMDS = data.frame(MDS1 = sol$points[,1],     MDS2=sol$points[,2],group=MyMeta$host)
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))
}

ggplot(data = NMDS, aes(MDS1, MDS2)) + geom_point(aes(color = group)) +
  geom_path(data=df_ell, aes(x=MDS1, y=MDS2,colour=group), size=1, linetype=2)+
  annotate("text",x=NMDS.mean$MDS1,y=NMDS.mean$MDS2,label=NMDS.mean$group)

#update - can use se (standard error) or sd (standard dev)
#update
NMDS = data.frame(MDS1 = sol$points[,1], MDS2 =        sol$points[,2],group=MyMeta$host)

plot.new()
ord<-ordiellipse(sol, MyMeta$host, display = "sites", 
             kind = "se", conf = 0.95, label = T)

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(ord[[g]]$cov,ord[[g]]$center,ord[[g]]$scale)))
                            ,group=g))
}

ggplot(data = NMDS, aes(MDS1, MDS2)) + geom_point(aes(color = group)) +
geom_path(data=df_ell, aes(x=NMDS1, y=NMDS2,colour=group), size=1, linetype=2)

第一种和第二种绘图方法之间有两个主要区别。 第一种方法是使用标准差计算椭圆路径,不缩放。第二种方法是使用标准误差,并且也在缩放数据。因此,通过对 ordiellipse 函数(kind='sd',而不是 'se')进行必要的更改,并删除比例( ord[[g]]$scale) 来自 veganCovEllipse 函数。我在下面包含了这个,所以你可以自己看看。

最终,两个情节都是正确的,这取决于你想展示什么。只要您指定使用标准误差或偏差,都可以使用。至于要不要缩放,这个真的要看你的数据了。 link 可能会有帮助:Understanding `scale` in R

第一种方法:

sol <-metaMDS(comm_mat)
MyMeta=env

#originalresponse
NMDS = data.frame(MDS1 = sol$points[,1],     MDS2=sol$points[,2],group=MyMeta$host)
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))
}

ggplot(data = NMDS, aes(MDS1, MDS2)) + geom_point(aes(color = group)) +
  geom_path(data=df_ell, aes(x=MDS1, y=MDS2,colour=group), size=1, linetype=2)+
  annotate("text",x=NMDS.mean$MDS1,y=NMDS.mean$MDS2,label=NMDS.mean$group)

给出:

第二种方法:

plot.new()

ord<-ordiellipse(sol, MyMeta$host, display = "sites", 
                 kind = "sd", conf = .95, label = T)

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(ord[[g]]$cov,ord[[g]]$center)))
            ,group=g))
}

plot2<-ggplot(data = NMDS, aes(MDS1, MDS2)) + geom_point(aes(color = group)) +
  geom_path(data=df_ell, aes(x=NMDS1, y=NMDS2,colour=group), size=1, linetype=2)+
  annotate("text",x=NMDS.mean$MDS1,y=NMDS.mean$MDS2,label=NMDS.mean$group)

plot2

同时给出: