而不是正常的绘图功能我在用的ggplot2
创建NMDS地块。 我想在使用该函数的曲线图NMDS显示组ordiellipse()
从所述vegan
包。
实施例的数据:
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))
我如何添加ordiellipse与创建的NMDS情节ggplot2
?
下面Didzis Elferts'回答的伟大工程。 谢谢! 不过,我现在在绘制以下ordiellipse与创建的NMDS情节感兴趣ggplot2
:
ordiellipse(sol, MyMeta$amt, display = "sites", kind = "se", conf = 0.95, label = T)
不幸的是,我没有足够的了解的如何veganCovEllipse
功能的工作原理,以便能够调整脚本自己。
首先,我添加栏目组到您的NMDS数据帧。
NMDS = data.frame(MDS1 = sol$points[,1], MDS2 = sol$points[,2],group=MyMeta$amt)
第二个数据帧包含平均MDS1和MDS2值为每个组,它将被用来显示在地块组名称
NMDS.mean=aggregate(NMDS[,1:2],list(group=group),mean)
数据帧df_ell
中包含的值,以显示椭圆。 据计算与功能veganCovEllipse
被隐藏在vegan
包。 该函数被施加到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)
思路椭圆作图从另一个计算器采用的问题 。
UPDATE - 解决方案,在这两种情况下工作
首先,使与组列NMDS数据帧。
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
包含的值来显示椭圆。 它与功能的重新计算veganCovEllipse
这是隐藏在vegan
包。 该函数被施加到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()
时,这种解决方案将与您提供此功能不同的参数运行。
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)
文章来源: Plotting ordiellipse function from vegan package onto NMDS plot created in ggplot2