在 ggplot2 中绘制 envfit 向量（素食包）答案

【问题标题】：Plotting envfit vectors (vegan package) in ggplot2在 ggplot2 中绘制 envfit 向量（素食包）
【发布时间】：2013-02-05 15:52:37
【问题描述】：

我正在完成我在 vegan 和 ggplot2 中创建的 NMDS 图，但无法弄清楚如何将 envfit 物种加载向量添加到图中。当我尝试它说“无效的图形状态”。

下面的示例是从另一个问题 (Plotting ordiellipse function from vegan package onto NMDS plot created in ggplot2) 稍微修改的，但它准确地表达了我想要包含的示例，因为我首先使用这个问题来帮助我将 metaMDS 放入 ggplot2：

library(vegan)
library(ggplot2)
data(dune)

# calculate distance for NMDS
NMDS.log<-log(dune+1)
sol <- metaMDS(NMDS.log)

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

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

#Add species loadings
vec.sp<-envfit(sol$points, NMDS.log, perm=1000)
plot(vec.sp, p.max=0.1, col="blue")

【问题讨论】：

+1 很好的问题。请研究我的答案，看看应该如何从对象中提取分数或绘制它们。您不应该，在大多数情况下，vegan 需要执行 sol$points 之类的操作，因为我们为此提供了 scores() 方法以及 plot() 方法。

标签： r ggplot2 vegan

【解决方案1】：

（否则很好）接受的答案的问题，这解释了为什么向量在包含的图中都是相同的长度 [请注意，接受的答案现在已经被编辑以按方式缩放箭头我在下面描述，以避免遇到 Q&A] 的用户混淆，存储在 envfit() 返回的对象的 $vectors$arrows 组件中的是拟合向量的方向余弦。这些都是单位长度，因此@Didzis Elferts 图中的箭头都是相同的长度。这与plot(envfit(sol, NMDS.log)) 的输出不同，这是因为我们通过与排序配置（“轴”）的相关性来缩放矢量箭头坐标。这样，与排序配置关系较弱的物种会得到更短的箭头。缩放是通过将方向余弦乘以sqrt(r2) 来完成的，其中r2 是打印输出表中显示的值。在将向量添加到现有绘图时，vegan 还会尝试缩放向量集，以便它们填充可用绘图空间，同时保持箭头的相对长度。如何做到这一点在?envfit 的详细信息 部分进行了讨论，并且需要使用未导出的函数vegan:::ordiArrowMul(result_of_envfit)。

这是一个完整的工作示例，它使用 ggplot2 复制了 plot.envfit 的行为：

library(vegan)
library(ggplot2)
library(grid)
data(dune)

# calculate distance for NMDS
NMDS.log<-log1p(dune)
set.seed(42)
sol <- metaMDS(NMDS.log)

scrs <- as.data.frame(scores(sol, display = "sites"))
scrs <- cbind(scrs, Group = c("hi","hi","hi","md","lo","hi","hi","lo","md","md",
                              "lo","lo","hi","lo","hi","md","md","lo","hi","lo"))

set.seed(123)
vf <- envfit(sol, NMDS.log, perm = 999)

如果我们在这一点停下来看看vf：

> vf

***VECTORS

             NMDS1       NMDS2     r2 Pr(>r)    
Belper -0.78061195 -0.62501598 0.1942  0.174    
Empnig -0.01315693  0.99991344 0.2501  0.054 .  
Junbuf  0.22941001 -0.97332987 0.1397  0.293    
Junart  0.99999981 -0.00062172 0.3647  0.022 *  
Airpra -0.20995196  0.97771170 0.5376  0.002 ** 
Elepal  0.98959723  0.14386566 0.6634  0.001 ***
Rumace -0.87985767 -0.47523728 0.0948  0.429
.... <truncated>

所以r2 数据用于缩放NMDS1 和NMDS2 列中的值。最终的情节是这样产生的：

spp.scrs <- as.data.frame(scores(vf, display = "vectors"))
spp.scrs <- cbind(spp.scrs, Species = rownames(spp.scrs))

p <- ggplot(scrs) +
  geom_point(mapping = aes(x = NMDS1, y = NMDS2, colour = Group)) +
  coord_fixed() + ## need aspect ratio of 1!
  geom_segment(data = spp.scrs,
               aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
               arrow = arrow(length = unit(0.25, "cm")), colour = "grey") +
  geom_text(data = spp.scrs, aes(x = NMDS1, y = NMDS2, label = Species),
            size = 3)

这会产生：

【讨论】：

只要 coord_fixed() 保持在 ratio=1，ggplot 中的点位置是否总是与 Vegan ordiplot 的（正确）点位置匹配？
@Vinterwoo 无论有没有coord_fixed 和ratio = 1，这些点都将匹配envfit 的plot 方法中的正确点位置。需要coord_fixed 等部分来保持排序解决方案的度量属性；两个轴都应该以相同的“单位”进行缩放。
感谢加文。关于使用 envfit 向量来覆盖有关物种如何影响站点排序位置的信息 - 我（很可能是有缺陷的）理解这种信息无法从 NMDS 分析中提取出来。并且不可能获得有关排序图解释的方差的信息。似乎来自 envfit 的重要性代码，以及 spp.scores 中箭头的位置似乎提供了这种信息。如果此评论涉及原始问题的更全局性方面，而不是您出色答案的具体细节，我深表歉意。
只有当 NMDS 无法访问物种数据时，您的理解才有意义。 NMDS 本身不使用此信息；嵌入到相异矩阵 (Dij) 中，NMDS 使用它。但是，我们可以使用其他方法熟悉的想法将向量或点投影到 NMDS 解决方案中。如果您使用原始数据而不是 Dij 拟合模型，metaMDS 的 plot 方法可以添加物种点作为 NMDS 站点得分的加权平均值。 envfit 使用成熟的向量拟合方法，post hoc。
@Vinterwoo 主要需要注意的是，这些如此衍生的分数（无论您采用哪种方式）都从 NMDS 解决方案中嵌入的“物种”信息中删除了。这本身取决于所使用的相异系数。

【解决方案2】：

从添加库开始。另外库grid 是必需的。

library(ggplot2)
library(vegan)
library(grid)
data(dune)

进行 metaMDS 分析并将结果保存在数据框中。

NMDS.log<-log(dune+1)
sol <- metaMDS(NMDS.log)

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

添加物种加载并将它们保存为数据框。箭头余弦的方向存储在列表vectors 和矩阵arrows 中。要获得箭头的坐标，这些方向值应乘以存储在vectors$r 中的r2 值的平方根。更直接的方法是使用@Gavin Simpson 的回答中提供的函数scores()。然后添加包含species 名称的新列。

vec.sp<-envfit(sol$points, NMDS.log, perm=1000)
vec.sp.df<-as.data.frame(vec.sp$vectors$arrows*sqrt(vec.sp$vectors$r))
vec.sp.df$species<-rownames(vec.sp.df)

箭头添加geom_segment()，物种名称添加geom_text()。这两个任务都使用了数据框vec.sp.df。

ggplot(data = NMDS, aes(MDS1, MDS2)) + 
  geom_point(aes(data = MyMeta, color = MyMeta$amt))+
  geom_segment(data=vec.sp.df,aes(x=0,xend=MDS1,y=0,yend=MDS2),
      arrow = arrow(length = unit(0.5, "cm")),colour="grey",inherit_aes=FALSE) + 
  geom_text(data=vec.sp.df,aes(x=MDS1,y=MDS2,label=species),size=5)+
  coord_fixed()

【讨论】：

这很棒！通过一些操作，我可以让它看起来非常漂亮，奇怪的是，加载似乎比来自 envfit 的更准确。我认为那些可能已经通过线的长度显示了 p 值，但我宁愿自己报告。谢谢！ :)
嗯，这里似乎有些不对劲——向量的长度不应该相同...它们不在envfit 的输出中。而@HFBrowning 向量的长度只是表示相关的强度；它不是 p 值的直接指示。
啊，我现在明白了——你不应该深入研究对象并撕掉组件。使用访问器函数 - 在本例中为 scores(vec.sp, display = "vectors")。有关完整示例，请参阅我的答案。
是的，你是对的，这不是应该的方式。 @HFBrowning 接受您的回答后，我会立即删除我的回答。
@DidzisElferts 而不是删除它，也许只是更新它？我的意思是你做了我刚刚为我的答案复制的大部分艰苦工作，只修改了几位样式并提取plot.envfit方法使用的相同数据。

【解决方案3】：

我可以晚一点补充吗？

Envfit 提供 pvalues，有时您只想绘制重要参数（纯素食者可以在 plot 命令中使用 p.=0.05 为您做一些事情）。我很难用 ggplot2 做到这一点。这是我的解决方案，也许您找到一个更优雅的解决方案？

从Didzis上面的回答开始：

ef<-envfit(sol$points, NMDS.log, perm=1000)
ef.df<-as.data.frame(ef$vectors$arrows*sqrt(ef$vectors$r))
ef.df$species<-rownames(ef.df)

#only significant pvalues
#shortcutting ef$vectors
A <- as.list(ef$vectors)
#creating the dataframe
pvals<-as.data.frame(A$pvals)
arrows<-as.data.frame(A$arrows*sqrt(A$r))
C<-cbind(arrows, pvals)
#subset
Cred<-subset(C,pvals<0.05)
Cred <- cbind(Cred, Species = rownames(Cred))

“Cred”现在可以在上面讨论的 geom_segment-argument 中实现。

【讨论】：

【解决方案4】：

简短补充：要在ggplot2 又名“箭头长度充分利用绘图区域”中获得plot.envfit 功能的完整表示，需要应用一个因素。我不知道上面的答案是否故意遗漏了它，因为加文甚至特别提到了它？只需使用arrow_factor <- ordiArrowMul(vf) 提取所需的比例因子，然后您可以将其应用于spp.scrs 中的两个 NMDS 列，或者您可以像手动执行此操作

arrow_factor <- ordiArrowMul(vf)
spp.scrs <- as.data.frame(scores(vf, display = "vectors")) * arrow_factor
spp.scrs <- cbind(spp.scrs, Species = rownames(spp.scrs), Pvalues = vf$vectors$pvals, R_squared = vf$vectors$r)

# select significance similarly to `plot(vf, p.max = 0.01)`
spp.scrs <- subset(spp.scrs, Pvalues < 0.01)

# you can also add the arrow factor in here (don't do both!)
ggplot(scrs) +
  geom_point(mapping = aes(x = NMDS1, y = NMDS2, colour = Group)) +
  coord_fixed() + ## need aspect ratio of 1!
  geom_segment(data = spp.scrs,
               aes(x = 0, xend = NMDS1 * arrow_factor, y = 0, yend = NMDS2 * arrow_factor),
               arrow = arrow(length = unit(0.25, "cm")), colour = "grey") +
  geom_text(data = spp.scrs, aes(x = NMDS1 * arrow_factor, y = NMDS2 * arrow_factor, label = Species),
            size = 3)

【讨论】：