将纯素包中的ordiellipse函数绘制到ggplot2中创建的NMDS图上

ilmyapht 于 2023-02-06 发布在其他

关注(0)|答案(2)|浏览(124)

我使用ggplot2创建NMDS图，而不是普通图函数。我希望使用vegan包中的ordiellipse()函数在NMDS图中显示组。
示例数据：

library(vegan)
library(ggplot2)
data(dune)
# calculate distance for NMDS
sol <- metaMDS(dune)
# Create meta data for grouping
MyMeta = data.frame(
  sites = c(2,13,4,16,6,1,8,5,17,15,10,11,9,18,3,20,14,19,12,7),
  amt = c("hi", "hi", "hi", "md", "lo", "hi", "hi", "lo", "md", "md", "lo", 
          "lo", "hi", "lo", "hi", "md", "md", "lo", "hi", "lo"),
  row.names = "sites")
# plot NMDS using basic plot function and color points by "amt" from MyMeta
plot(sol$points, col = MyMeta$amt)
# draw dispersion ellipses around data points
ordiellipse(sol, MyMeta$amt, display = "sites", kind = "sd", label = T)

# same in ggplot2
NMDS = data.frame(MDS1 = sol$points[,1], MDS2 = sol$points[,2])
ggplot(data = NMDS, aes(MDS1, MDS2)) + 
  geom_point(aes(data = MyMeta, color = MyMeta$amt))

如何向使用ggplot2创建的NMDS图添加ordiellipse？
下面Didzis Elferts的回答非常好。谢谢!但是，我现在有兴趣将下面的ordiellipse绘制到使用ggplot2创建的NMDS图中：
ordiellipse(sol, MyMeta$amt, display = "sites", kind = "se", conf = 0.95, label = T)
不幸的是，我对veganCovEllipse函数的工作原理了解不够，无法自己调整脚本。

来源：https://stackoverflow.com/questions/13794419/plotting-ordiellipse-function-from-vegan-package-onto-nmds-plot-created-in-ggplo

2条答案

按热度按时间

o8x7eapl1#

首先，我向NMDS数据框中添加了列组。

NMDS = data.frame(MDS1 = sol$points[,1], MDS2 = sol$points[,2],group=MyMeta$amt)

第二个数据框包含每组的平均MDS 1和MDS 2值，将用于在图上显示组名称

NMDS.mean=aggregate(NMDS[,1:2],list(group=group),mean)

数据框df_ell包含显示椭圆的值，使用vegan包中隐藏的函数veganCovEllipse计算，该函数应用于NMDS（组）的各级，也使用函数cov.wt计算协方差矩阵。

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))
  }

现在，使用用于绘制组名称的函数geom_path()和annotate()绘制椭圆。

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)

椭圆绘制的思想借鉴了另一个堆栈overflow question。

UPDATE -适用于两种情况的解决方案

首先，制作具有分组列的NMDS Dataframe 。

NMDS = data.frame(MDS1 = sol$points[,1], MDS2 = sol$points[,2],group=MyMeta$amt)

接下来，将函数ordiellipse()的结果保存为某个对象。

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

数据框df_ell包含显示椭圆的值，使用vegan包中隐藏的函数veganCovEllipse重新计算，该函数应用于NMDS（组）的每一级，现在使用存储在ord对象中的参数-每一级的cov、center和scale。

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))
}

绘图的方法和上一个例子一样，至于计算ordiellipse()的elipses对象的坐标，这个解决方案将使用你为这个函数提供的不同参数。

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)

赞(0）回复(0）举报 2023-02-06

roqulrg32#

两项更新：

NMDS.mean=aggregate(NMDS[,1:2],list(group=group),mean)
应更新为

NMDS.mean=aggregate(NMDS[,1:2],list(group=NMDS$group),"mean")

NMDS$group不是因子，因此在NMDS$group的水平上循环不起作用。Df_ell返回零个变量中的零个观测。

NMDS$group <- as.factor(NMDS$group)

解决了这个问题。

赞(0）回复(0）举报 2023-02-06

我来回答

将纯素包中的ordiellipse函数绘制到ggplot2中创建的NMDS图上

2条答案

UPDATE -适用于两种情况的解决方案

相关问题

热门标签

最新问答