如何把环境变量加入R语言的NMDS分析中

非度量多维缩放 (NMDS) 是可用于在二维中显示多维数据的多种排序图类型之一。NMDS 图基本上突出了复杂多维数据样本之间的相似性。图上的两个点（样本）越接近，这些样本在基础数据方面就越相似。就我而言，我通常会查看样本中微生物的物种丰度。因此，两个样本越接近，它们的微生物群落就越相似。

生成 NMDS 图的方法是非度量的，这意味着该算法不假设任何特定的数据分布，当您拥有本质上奇怪且不符合钟形曲线样式的数据时，这很方便分配。该算法使用基于等级的方法，该方法基本上根据特定距离度量（即 Bray Curtis）根据彼此的接近程度对样本进行排名，并在二维空间中绘制这些关系。更多关于这里方法的细节。虽然这是一种针对不合格数据的稳健方法，但这意味着在 NMDS 排序图中，轴和坐标系是任意的，并且与任何有意义的值没有直接关系。

重要说明： NMDS 只是一种可视化技术，而不是对样本分离或相关性的统计评估。为此，您应该运行统计检验，例如分类变量的ANOSIM 检验或连续变量的Mantel 检验。如果您希望轴有意义，或者如果您想讨论变体分区，其他排序技术（如PCA、CA、CCA、RDA 等）可能对您更有用。

好了，说了这么多。有多种方法可以在 NMDS 图上叠加额外信息，以帮助可视化影响数据的潜在趋势。就我而言，潜在的环境变量与微生物群落结构共同变化，从而可能推动群落的变化。来自 R 包vegan的一系列函数旨在将额外信息覆盖在 NMDS 图等排序图上。这些函数包括envfit、ordiplot、ordiellipse和ordisurf。我将首先解释如何在这里使用 envfit 的功能。

ENVFIT

Envfit 将环境向量或因素拟合到一个排序上。envfit的默认绘图，如metaMDS，在美学上不是很讨人喜欢，所以我将展示如何使用ggplot2绘制两者。

首先加载您的库并读入您的数据：

library(vegan) 

library(ggplot2) 

df = read.csv("your_data_frame.csv", header = TRUE)

我的数据包含样本，因为行和列包含环境变量或物种丰度数据。在这个例子中，我有分类和连续环境变量。

接下来，对您的数据进行子集化，以便您拥有一个仅包含环境变量(env)的数据框和一个仅包含物种丰度数据(com)的数据框。在下面的代码中，我命名了包含相应信息的数据中的列范围。

com = df[,9:32] 

env = df[,2:8]

执行 NMDS 排序，如本教程中所述

#convert com to a matrix

m_com = as.matrix(com) 

#nmds code 

nmds = metaMDS(m_com, distance = "bray") 

nmds

现在我们使用环境数据框env运行 envfit 函数。

en = envfit(nmds, env, permutations = 999, na.rm = TRUE)

en

第一个参数是我们刚刚执行的 NMDS 排序中的 metaMDS 对象。接下来是env，我们的环境数据框。然后我们声明我们想要 999 个排列，并删除任何缺少数据的行。

当我们调用 envfit 函数的结果时，我们得到的结果被分成Vectors和Factors，分别用于您的连续变量和分类变量。如果您的数据不包含两者，那么您将只能看到向量或因子。对于向量和因子，显示了一个表格，将变量作为行给出，然后给出它们各自在 NMDS1 和 NMDS2 轴上的 NMDS 排列坐标。如果排列 > 0，则使用环境变量的排列来评估拟合向量或因子的重要性。

我解释这一点的方式是向量在你的情节中充当一种“伪轴”。如果你沿着你的图朝那个方向移动，你会看到你的样本通常相对于那个变量增加。这当然不是一个完美的关系，这就是 r2 和显着性值来解释关联强度的地方。更长的箭头，意味着更强的关联。因为这种关系并不总是线性的，envfit 的开发人员建议您也可以使用其他方法将重要变量映射到您的数据上，包括点大小或 ordisurf 函数。

属于某个类别的所有样本的质心值由因子坐标绘制，这使您可以查看来自特定类别的样本分组在一起（或不分组）的程度。

我们可以使用基数 R 绘制 NMDS 和 envfit 结果，以查看因子和向量是如何绘制的：

plot(nmds) 

plot(en)

envfit 向量和因子（蓝色）覆盖在原始 NMDS 图上，样本为黑点。这里的基础 R 情节真的很难阅读，容易人满为患，而且难以定制。我建议使用ggplot2来制作更好看的图。

为了使用ggplot2进行绘图，您需要从 nmds 和 envfit 结果中提取适当的信息。

对于 NMDS 输出，使用以下代码提取 NMDS 排序空间中的样本坐标。然后从您的原始数据中添加包含您想要包含在绘图中的信息的列。在这种情况下，我包括“季节”列

data.scores = as.data.frame(scores(nmds)) 

data.scores$season = df$Season

vegan 包的最新版本 (>2.6-2) 将 score(nmds) 对象的格式更改为列表，因此上面的代码将抛出一个错误，说 arguments 暗示不同的行数。如果您收到此错误，请尝试使用以下代码来提取您的网站分数：

#extract NMDS scores (x and y coordinates) for sites from newer versions of vegan package 

data.scores = as.data.frame(scores(nmds)$sites) 

#add 'season' column as before 

data.scores$season = df$Season

从 envfit 结果中提取所需的信息有点复杂。envfit 输出包含有关每个变量的段长度的信息。这些片段被缩放到 r2 值，因此具有较长片段的环境变量与数据的相关性比具有较短片段的那些更强烈。您可以使用分数提取此信息。然后这些长度被进一步缩放以适合绘图。这是通过特定于分析的乘数完成的，可以使用命令ordiArrowMul(en)访问。下面我将分数乘以这个乘数，以使坐标保持正确的比例。

因为我的数据包含连续的和分类的环境变量，所以我分别使用“向量”和“因子”选项从两者中提取信息。

en_coord_cont = as.data.frame(scores(en, "vectors")) * ordiArrowMul(en) en_coord_cat = as.data.frame(scores(en, "factors")) * ordiArrowMul(en)

现在数据采用适当的格式，可以使用 ggplot2 在 R 中绘制。

首先，我们将制作一个没有 envfit 变量的 nmds 图：

gg = 

ggplot(data = data.scores, aes(x = NMDS1, y = NMDS2)) + 

geom_point(data = data.scores, aes(colour = season), size = 3, alpha = 0.5) +     

scale_colour_manual(values = c("orange", "steelblue")) + 

theme(axis.title = element_text(size = 10, face = "bold", colour = "grey30"),     panel.background = element_blank(), panel.border = element_rect(fill = NA, colour = "grey30"), axis.ticks = element_blank(), axis.text = element_blank(), legend.key = element_blank(), legend.title = element_text(size = 10, face = "bold", colour = "grey30"), legend.text = element_text(size = 9, colour = "grey30")) + 

labs(colour = "Season")      

gg

下面是 NMDS 图和 envfit 数据的代码：

gg = 

ggplot(data = data.scores, aes(x = NMDS1, y = NMDS2)) +

geom_point(data = data.scores, aes(colour = season), size = 3, alpha = 0.5) +   scale_colour_manual(values = c("orange", "steelblue")) +

geom_segment(aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2)，data = en_coord_cont, size =1, alpha = 0.5, colour = "grey30") + 

geom_point(data = en_coord_cat, aes(x = NMDS1, y = NMDS2),shape = "diamond", size = 4, alpha = 0.6, colour = "navy") +      

geom_text(data = en_coord_cat, aes(x = NMDS1, y = NMDS2+0.04),label = row.names(en_coord_cat), colour = "navy", fontface = "bold") + 

geom_text(data = en_coord_cont, aes(x = NMDS1, y = NMDS2), colour = "grey30",fontface = "bold", label = row.names(en_coord_cont)) + 

theme(axis.title = element_text(size = 10, face = "bold", colour = "grey30"), panel.background = element_blank(), panel.border = element_rect(fill = NA, colour = "grey30"), axis.ticks = element_blank(), axis.text = element_blank(), legend.key = element_blank(), legend.title = element_text(size = 10, face = "bold", colour = "grey30"), legend.text = element_text(size = 9, colour = "grey30")) +

labs(colour = "Season")   

gg

看起来比原始的基本 R 示例要好得多。代码有点复杂，因为绘图现在有很多方面，包括点、线段和文本标签。此图与本教程系列中的其他图的主要区别在于，我们在这里从 3 个不同的数据帧中绘制信息，因此每次向图中添加几何图层时都必须指定数据源。

资料来源 https://jkzorz.github.io/2020/04/04/NMDS-extras.html