企业官网建设流程全解析

简介：一份面向R语言生态学数据分析用户的NMDS排序与一般加性模型映射实战教程，适合需要处理物种组成、环境因子与样点分组数据的科研人员及学生。内容系统讲解非度量多维标尺排序的原理，以及与PCA、PCoA的差异，并通过vegan包metaMDS()逐步演示数据准备、排序执行、结果提取和ggplot2可视化，还包含使用mgcv包进行GAM建模、评估及将环境因子平滑项映射到排序图的完整代码。教程中嵌入了可复用的R脚本与分组凸包、等高线映射等进阶绘制方法。资源包共1个PDF文档，大小45.43MB，图文与代码对照清晰，便于按步骤实操学习。目前已有1204人浏览学习，适合具有一定R基础、希望系统掌握NMDS与GAM结合分析的研究者。

1. 拿到物种数据想画一张能讲故事的排序图：NMDS 加 GAM 映射为什么值得学

我最早被 NMDS 折腾，是在处理一批土壤微生物群落数据的时候。几十个样点、一百多个物种，用 PCA 排出来的图一团糟，样点全挤在原点附近，物种贡献率也解释不清。后来换成 R 语言里的非度量多维标尺排序（NMDS），图形结构一下子出来了，环境因子的梯度方向也清楚了。但排序图只是第一步，真正让分析有说服力的，是把环境变量以等值线的形式映射到排序空间里——这一步用一般加性模型（GAM）做最稳。

NMDS 解决的是「不用假设线性关系、也能还原样点间生态距离」的问题，GAM 解决的是「环境梯度到底在排序图里怎么弯曲变化」的问题。两段式工作流是生态学、环境监测、微生物群落研究里最常用也最容易被误用的组合之一。这篇文章会带你跑通完整流程，也把参数选择和踩坑点讲透。

2. 为什么生态学排序偏好 NMDS：秩次逻辑与三个选型理由

2.1 从「距离」到「排队」：NMDS 到底在算什么

传统的 PCA 保留的是样点间的实际距离数值，并且要求物种对环境的响应是线性的。但真实的群落数据很少满足这个假设——物种丰度大量为零，响应曲线多是单峰型，环境梯度一长，PCA 的排序结果就会出现马蹄形扭曲。NMDS 的思路完全不同：它先计算样点对之间的生态距离（比如 Bray-Curtis 距离），然后把距离的绝对数值扔掉，只保留大小顺序，也就是秩次。之后在低维空间里找一组坐标，使得空间距离的秩次与原始距离的秩次尽可能一致。

这个「只看排队顺序、不看具体数值」的机制，就是「非度量」三个字的含义。它带来的好处非常实际：数据不服从正态分布没关系，物种丰度里一堆零值也扛得住，环境梯度是非线性响应也不会把你带偏。所以我处理植被、底栖动物、微生物这类高零值、高偏态的丰度矩阵时，默认首选就是 NMDS。

排序质量用 stress 值衡量。stress 反映的是低维空间中距离秩次与原始距离秩次的偏差程度。经验阈值是：小于 0.05 为优秀，0.05 到 0.10 为良好，0.10 到 0.20 尚可接受，超过 0.20 就说明二维排序图失真明显，需要增加维度。R 语言里跑完 metaMDS 后，结果对象里的stress字段就是你要看的第一个数字。

2.2 metaMDS 背后的策略：随机起始与 procrustes 检验

NMDS 的求解不像 PCA 那样有解析解，它是一个迭代寻优过程：随机生成一组初始坐标，反复调整样点位置，让 stress 逐步下降。这个迭代非常依赖起始位置，不同起始位置可能收敛到不同的局部最优解。vegan 包里的metaMDS()做了一个很聪明的包装：它用trymax次随机起始重复运行，从中挑选 stress 最小的解作为最终结果，并且会对多次运行结果做 procrustes 旋转比较，评估解是否稳定。

procrustes 分析做的事情是：把一个配置旋转、平移、缩放后与另一个配置最大程度重合，然后计算残差。metaMDS()内部用它来对比多次随机起始的收敛情况，如果每次起点不同但最终排序高度一致，说明这个解是可信的。这也是为什么第 6 章里我要专门演示一次手动的 procrustes 验证——它不是锦上添花，而是 NMDS 结果可重复性的基本保障。

需要提醒的是，NMDS 不输出特征值和贡献率。论文里 PCA 图下面写「前两轴解释了 45.2% 的方差」，NMDS 图下面不这么写，只写stress = 0.08。新手最容易犯的错误就是在 NMDS 结果里去找轴上百分数，这是概念上的根本区别。如果你需要报告每个环境因子对群落变异的贡献占比，那要用envfit()或变差分解，而不是 NMDS 本身。

2.3 两个容易被忽略的函数参数：autotransform 与 trymax

metaMDS()的两个参数，我见过大量使用者没搞清楚就默认跑完了。第一个是autotransform，默认值为 TRUE。它会在计算距离前自动对丰度数据做 Wisconsin 双标准化：先把每个物种的丰度除以该物种在所有样点的总量，得到相对优势度，再把每个样点的所有值按样点总量标准化。这个变换对常见生态数据效果不错，但它是「静默地」修改了你的距离矩阵。如果你提前用 Hellinger 转换或其他方法预处理过数据，就必须把autotransform设为 FALSE，否则等于做了两遍变换，结果解释起来会很别扭。

第二个是trymax。默认值只有 20，意味着只尝试 20 次随机起始。数据量大、群落结构复杂时，20 次常常找不到低 stress 的解，运行结束会报「结果未收敛」的警告。我一般会在最终分析时把trymax提到 100 甚至 200，虽然耗时增加，但换来的是更低的 stress 和更稳定的排序。set.seed()也一定要在运行前固定，否则每次结果都有细微旋转差异，后续的 GAM 映射图和文字描述就对不上了。

3. 一般加性模型凭什么能映射环境梯度：从样条到等值线

3.1 一个公式理解 GAM 映射

NMDS 排序图本身只有样点和物种的位置，环境变量如何在这个空间中变化，需要额外建模。一般加性模型（GAM，generalized additive model）的思路很直观：把样点在排序图上的坐标当作解释变量，把某个连续环境变量当作响应变量，拟合一个平滑曲面，然后在排序图上画出这个曲面的等值线。模型写出来就是：

A1_i = β0 + f(NMDS1_i, NMDS2_i) + ε_i

这里的 f 是一个二维光滑函数，由样条基函数组合而成。举个例子，土壤厚度 A1 这个变量，在排序空间里可能从左下到右上逐渐增加，也可能是中间高、四周低的斑块状分布。线性模型只能拟合前者，二次多项式勉强拟合后者但边界容易振荡。GAM 的好处是让数据自己决定弯曲程度：样条基函数的数量决定了曲面的复杂度，光滑参数控制曲线的抖动幅度，两者都由数据驱动地优化。

用 R 语言的 mgcv 包拟合这个模型只需要一行gam()调用，输出里会给出解释方差和各项显著性。你得到的模型对象可以直接喂给 vegan 的ordisurf()，或者手动生成预测网格画等值线。这也是 GAM 映射相比其他拟合方式最舒服的地方——建模和可视化工具链是现成的。

3.2 mgcv 里 k 与 method 的选择

拟合 GAM 时最常被问的参数就是 k，也就是样条基函数的数量。k 太小，曲面过于平滑，真实的环境梯度被抹平；k 太大，模型开始拟合噪声，等值线图会出现诡异的波浪状。一个保守的起点是 k = 6 到 8，因为生态学里绝大多数环境梯度用 6 个基函数已经能描述。如果样本量只有二三十个点，k 不要超过样本量的一半，否则会出现奇异拟合。

method 参数推荐用 "REML"。mgcv 提供了 GCV 和 REML 两种平滑参数选择方法，老版本默认 GCV，新版本推荐 REML。我的习惯是统一用method = "REML"，因为 REML 对光滑参数的估计更稳定，尤其是数据量小或存在强相关性时，不容易过拟合。拟合完之后用gam.check()做诊断，重点看 k-index 是否接近 1。如果显著小于 1，说明基函数不够，需要调大 k。这是验证模型设定是否合理的标准做法，也是和「跑完就看 p 值」式分析的差距所在。

3.3 为什么不直接用 envfit 画箭头

vegan 包里的envfit()可以快速把环境变量投影到排序图，画成箭头或因子重心，很多论文里 NMDS 图就配一组箭头。那为什么还要 GAM 映射？因为envfit()做的是线性拟合（或者说是多元回归在排序空间的投影），它只能表达单调递增或递减的环境梯度。

GAM 映射的价值在于展示非线性。比如土壤湿度对群落的影响可能在中间水平达到最优，两侧群落结构都变差，这在排序空间里就是一个驼峰形等值线。envfit 箭头对这种响应会显示为长度很短、方向不明显，因为线性拟合的斜率高不起来。而 GAM 等值线能完整保留这个生态学上有意义的结构。

我的做法是两者结合使用：因子型环境变量（比如土地利用类型、处理组别）用envfit()画重心点，连续环境变量用 GAM 画等值线。前者分组解释，后者梯度解释，互为补充。如果你只看重线性趋势，envfit 就够了；但要做精细的环境解释，GAM 映射基本是必须的。

4. 用 vegan + mgcv 跑通 NMDS 与 GAM 映射：dune 数据全流程

4.1 准备数据：物种矩阵与环境矩阵的对齐约定

这一节用 vegan 包内置的dune数据做完整演示。dune是荷兰一处沙丘草地的 20 个样点、29 个物种的丰度数据，配套的dune.env提供了土壤厚度（A1）、水分等级、土地利用方式等环境变量。数据量小但结构典型，非常适合用来验证流程。

library(vegan) library(mgcv) data(dune) data(dune.env) dim(dune) head(dune.env)

我的习惯是第一步先看数据的行名是否对齐。dune和dune.env的行名顺序是一致的，但在你的实际数据里不一定，合并前用rownames()核对一遍，或者直接cbind()后一起过滤缺失值。这一步不做好，GAM 建模时坐标和环境变量就会错位，拟合出来全是噪声，还很难排查。

环境变量要分清数据类型：A1 是连续数值变量，可以直接进 GAM；Moisture 是有序因子，适合用envfit()处理。如果你自己的数据里有 pH、总氮这类数值变量，直接和排序坐标一起建数据框就行。

4.2 跑 NMDS：从 vegdist 到 metaMDS

首先是物种数据的变换。我选择 Hellinger 转换，它对高零值丰度数据友好，又不像 Wisconsin 标准化那样抹掉物种间的相对差异。

set.seed(42) dune_hel <- decostand(dune, method = "hellinger") dune_dist <- vegdist(dune_hel, method = "bray") sol <- metaMDS(dune_dist, k = 2, trymax = 200, autotransform = FALSE) sol$stress

decostand()做 Hellinger 转换，把每个样点的物种丰度除以样点总和后取平方根，降低极端高丰度物种的影响。vegdist()计算 Bray-Curtis 距离，这是群落生态学的默认距离，对丰度矩阵的零值容忍度很高。metaMDS()里autotransform = FALSE是因为我们已经手动做了一次变换，不需要它再自动处理。trymax = 200保证有足够多的随机起始次数来找低 stress 解。

跑完后查看sol$stress。这个值一般在 0.1 上下浮动，属于可接受范围。如果有警告说找不到收敛解，第一反应是把 k 从 2 改成 3，而不是继续无脑加大 trymax。维度不够时，再多随机起始也只是反复困在局部最优里。

顺便对比一下不做 Hellinger 转换、直接用默认参数的效果：

set.seed(42) sol_default <- metaMDS(dune, k = 2, autotransform = TRUE) sol_default$stress

两个方案的 stress 值会有差异。这提醒你一个重要的方法学决策：用了什么样的数据预处理，必须在论文方法部分写清楚。Hellinger + Bray-Curtis 和默认的 Wisconsin 双标准化是两套完全不同的预处理逻辑，结果不能直接互换。

4.3 跑 GAM 映射：ordisurf 与 ordispGRID 两种出图路径

拿到 NMDS 对象后，最简单的 GAM 映射出图方式是用ordisurf()。它内部会调用 mgcv 拟合光滑曲面，并把等值线直接叠加到排序图上。

ordiplot(sol, display = "sites", type = "t") ordisurf(sol, dune.env$A1, add = TRUE, knots = 6)

ordiplot()先画出样点的排序图，type = "t"表示显示样点行名。ordisurf()的add = TRUE表示在已有图上叠加等值线，knots = 6控制样条基函数的数量。绘图后你可以用summary(ordisurf(sol, dune.env$A1, knots = 6))查看模型的解释方差。

如果想把建模过程掌握在自己手里，就需要用gam()手动拟合，再自己生成网格画等值线。这样自由度更高，也能输出模型诊断结果。

site_scores <- as.data.frame(scores(sol, display = "sites")) colnames(site_scores) <- c("NMDS1", "NMDS2") data_gam <- cbind(site_scores, A1 = dune.env$A1) gam_fit <- gam(A1 ~ s(NMDS1, NMDS2, k = 6, bs = "tp"), data = data_gam, method = "REML") summary(gam_fit) gam.check(gam_fit)

scores(sol, display = "sites")提取样点在 NMDS 空间中的坐标。注意新版 vegan 的列名是 NMDS1、NMDS2，老版本可能是 MDS1、MDS2，用colnames()强制统一命名最保险。s()是 mgcv 的光滑项，bs = "tp"是薄板样条，适合二维曲面拟合。method = "REML"是平滑参数选择方法，前面说过，稳定性优于 GCV。

summary(gam_fit)里有两个关键输出：R-sq.(adj) 和 Deviance explained。前者表示排序坐标能解释土壤厚度变异的百分比，后者是广义版本的偏差解释率。如果解释率过高（比如超过 0.8），先怀疑过拟合，检查 k 是否设置过大；如果过低，可能这个环境变量与群落结构关系确实弱，或者样点数量不够建模。gam.check()的 k-index 低于 1 且 p 值小，说明基函数数量不足，需要调大 k。

手动绘制等值线需要用预测网格：

grd <- expand.grid( NMDS1 = seq(min(site_scores$NMDS1) * 1.1, max(site_scores$NMDS1) * 1.1, length = 60), NMDS2 = seq(min(site_scores$NMDS2) * 1.1, max(site_scores$NMDS2) * 1.1, length = 60) ) grd$fit <- predict(gam_fit, newdata = grd) ordiplot(sol, display = "sites", type = "n") contour(unique(grd$NMDS1), unique(grd$NMDS2), matrix(grd$fit, nrow = 60), add = TRUE, col = "grey30") points(sol, display = "sites", pch = 16)

expand.grid()生成 NMDS1 和 NMDS2 的规则网格，范围向外扩 10% 是为了让等值线不贴边。predict()用拟合好的 GAM 模型算出网格上每个点的土壤厚度预测值。contour()把预测值矩阵画成等值线，注意它要求横纵坐标是递增的唯一点序列，所以要用unique()去重，并把预测向量转成 60 行矩阵。最后用points()把样点叠回去。

4.4 输出高分辨率图：pdf 与 png 的参数和中文字体问题

分析做完，导出图是最后一道工序。我一般用 PDF 和 TIFF 双份输出：PDF 给编辑部矢量图，TIFF 给报告预览。

pdf("nmds_gam_A1.pdf", width = 7, height = 6) ordiplot(sol, display = "sites", type = "n") ordisurf(sol, dune.env$A1, add = TRUE, knots = 6) points(sol, display = "sites", pch = 16, col = "darkgreen") dev.off() tiff("nmds_gam_A1.tiff", width = 2400, height = 2000, res = 300) ordiplot(sol, display = "sites", type = "n") ordisurf(sol, dune.env$A1, add = TRUE, knots = 6) points(sol, display = "sites", pch = 16) dev.off()

pdf()的width和height单位是英寸，7×6 是常规双栏排版尺寸。tiff()的res = 300保证打印清晰度，width和height按像素写，2400×2000 对应 8×6.67 英寸。如果图上要加中文标注，pdf()里需要指定支持中文字体的family参数，否则中文会变成方块；更省事的办法是图内只用英文标注，中文说明放在图注里。

5. NMDS 与 GAM 映射的五处翻车点：现象、原因与排查

5.1 stress 明明很低，排序图却看不出结构

现象：stress = 0.06，很漂亮，但在 ordiplot 里样点挤成一团，看不出任何生态梯度。原因：stress 低只表示低维空间里的距离秩次与原始距离秩次一致，不代表样点在二维平面有分散的分布。如果你的群落组成高度相似，样点之间距离本来就小，排出来自然挤在一起。这是数据结构决定的，不是模型错误。解决方法是先检查原始距离矩阵：用summary(dune_dist)看距离值的分布范围，如果大部分距离集中在很窄的区间，说明样点间差异本来就小，考虑换距离算法，或者检查是否需要对数据做更强烈的变换（比如把 Hellinger 换成标准化后的 Bray-Curtis）。还有一个常见但容易忽略的原因：你画的图里只有样点没有物种，当样点间差异主要体现在少数物种上时，图上叠加物种点（display = "species"）能看出解释方向。

5.2 物种矩阵零值过多导致距离计算异常

现象：跑vegdist()时出现警告，提示「您有零距离的样品对」，而且算出来的距离矩阵里有 NaN。原因：某两个样点完全没有共同物种，或者某些物种只出现了一次，导致 Bray-Curtis 计算不稳定。尤其是样本量小、物种稀少的调查数据，这种现象非常常见。解决办法分两步走：先做物种筛选，过滤掉只在 1 到 2 个样点出现的物种，用dune[, colSums(dune > 0) >= 3]这类方式实现；再考虑用decostand()做 Hellinger 转换，它能有效压缩零值造成的距离失真。如果筛选后仍然有样点的物种数为 0，直接删除这个样点，因为它对排序没有任何信息贡献。

5.3 把公式方向写反，等值线变成乱码

现象：GAM 模型summary()显示解释率很高，但等值线图完全是乱的，等值线走向和样点分布毫无逻辑。原因：建模时把排序坐标和环境变量的位置搞反了。正确写法是A1 ~ s(NMDS1, NMDS2)，也就是环境变量当响应变量，排序坐标当解释变量。有人会顺手写成NMDS1 ~ s(A1, NMDS2)或者把 A1 放到了光滑项里面，拟合结果看起来也有显著性，但画出来的等值线表达的就不是「环境变量在排序空间中的梯度」了。这是方向性错误，不是参数问题。排查方法很简单：打印formula(gam_fit)，看一眼模型公式左右两边是不是符合预期。还有一个连带错误：cbind()对齐时行名顺序被打乱，人为制造了错位数据。每次合并后都跑一次all(rownames(data_gam) == rownames(dune.env))做校验。

5.4 k 值设置超过样本量，模型直接奇异

现象：gam.fit报错「singular covariance matrix」，或者模型拟合完成但gam.check()显示 k-index 极低。原因：k 值代表样条基函数的个数，它不可能大于有效的样本自由度。20 个样点，k 设成 20，每个数据点几乎都被单独建模，协方差矩阵必然奇异。解决办法是把 k 控制为样本量的三分之一到二分之一，20 个点时 k = 6 到 8 就足够；同时用gam.check()的 k-index 反向验证——k-index 显著小于 1 才需要增加 k，不要一开始就堆大数。还有一个容易被忽略的细节：ordisurf()里也有knots参数，默认值不一定适合你的数据量，手动建模并诊断后再画图，是更稳妥的路径。

5.5 metaMDS 每次运行结果都不一样，图也对不上

现象：同一份数据，跑了两次metaMDS()，排序图大方向相似但轴旋转方向不同，后面叠加的 GAM 等值线方向也跟着变。原因：NMDS 是迭代寻优，随机起始位置不同导致收敛方向不同。排序的绝对值没有唯一解，只有旋转后的一致性才是可比的。这不仅影响复现，还会让你分析写到一半发现前后两张图方向不一致。解决办法是两层的：分析最开始时执行set.seed(42)（或任一个固定数字），并在方法部分写明随机种子；投稿前用两次不同种子的结果做 procrustes 验证（见第 6 章），把比较结果写进补充材料。如果 procrustes 的相关系数低于 0.9，说明 k = 2 的排序解不够稳定，需要尝试 k = 3 或检查数据中是否存在极端异常样点。

6. 稳定性和可发表图：procrustes 验证与 α 多样性叠加

6.1 用 protest 验证两次独立排序的一致性

这一步建议在每次正式分析出图前做一遍，成本很低但能避免大量返工。方法是用两个不同随机种子分别跑 NMDS，再用procrustes()和protest()做旋转对比。

set.seed(1) sol_1 <- metaMDS(dune_dist, k = 2, trymax = 100, autotransform = FALSE) set.seed(2) sol_2 <- metaMDS(dune_dist, k = 2, trymax = 100, autotransform = FALSE) pro <- procrustes(sol_1, sol_2) set.seed(123) test <- protest(sol_1, sol_2, permutations = 999) test

protest()里的permutations = 999是置换检验次数，检验两个配置的相似性是否显著优于随机排列。输出里的 correlation 越接近 1，说明两次独立排序越一致。我的判断标准是：correlation 大于 0.95 直接用第一套结果；在 0.85 到 0.95 之间，检查是否有某个样点贡献了极大残差；低于 0.85，考虑增加维度到 k = 3 或回到数据处理阶段排查异常样点。

procrustes()对象本身也可以用plot(pro)画残差图，横轴是对比配置的坐标，纵轴是每个样点的旋转残差。残差突出的点就是从排序稳定性角度值得复审的样点，这些往往也是后来在做 GAM 映射时预测误差最大的样点。

6.2 把 α 多样性叠加到排序图，形成一张综合图

环境等值线配 α 多样性气泡，是我做群落分析最喜欢用的组合方式。它在一张图里同时回答了「环境怎么变化」和「多样性在哪里更高」。用 vegan 的diversity()计算 Shannon 指数，然后按指数大小控制样点大小叠加到已有的 NMDS 和 GAM 图上。

shannon <- diversity(dune, index = "shannon") ordiplot(sol, display = "sites", type = "n") ordisurf(sol, dune.env$A1, add = TRUE, knots = 6, col = "grey40", labcex = 0.8) points(sol, display = "sites", pch = 16, cex = 1.2 + 2 * shannon / max(shannon), col = rgb(0.1, 0.4, 0.1, 0.7))

cex参数按 Shannon 指数的缩放倍率设置，1.2 是基础大小，最高放大 3.2 倍，视觉区分度足够。rgb()的透明度让重叠样点也能看清。如果样点重叠严重，可以用ordilabel()只标记目标样点名，或者使用orditorp()自动避免文字重叠。这张图导出后，图注里写清楚等值线是 GAM 拟合的 A1 梯度，气泡大小代表 Shannon 指数，解释起来非常直观。

另外，如果你有分组信息（比如处理 vs 对照），可以在图边用ordiellipse()加置信椭圆，和 α 多样性叠加起来看，组间差异和多样性关系就很清楚了。

做完整套流程，我最后还有一个习惯：所有脚本从一开始就用 RMarkdown 或纯 R 脚本记录，随机种子、数据预处理细节和 k 值选择理由都写在注释里。这样三个月后回来补分析，或者审稿人要求改参数重跑，都是几分钟的事。生态学数据分析里，可复现性比炫技重要得多。希望这套 NMDS 加 GAM 映射的流程能帮你在自己的数据上少走弯路。

本文还有配套的精品资源，点击获取