简介一份面向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 与 trymaxmetaMDS()的两个参数我见过大量使用者没搞清楚就默认跑完了。第一个是autotransform默认值为 TRUE。它会在计算距离前自动对丰度数据做 Wisconsin 双标准化先把每个物种的丰度除以该物种在所有样点的总量得到相对优势度再把每个样点的所有值按样点总量标准化。这个变换对常见生态数据效果不错但它是「静默地」修改了你的距离矩阵。如果你提前用 Hellinger 转换或其他方法预处理过数据就必须把autotransform设为 FALSE否则等于做了两遍变换结果解释起来会很别扭。第二个是trymax。默认值只有 20意味着只尝试 20 次随机起始。数据量大、群落结构复杂时20 次常常找不到低 stress 的解运行结束会报「结果未收敛」的警告。我一般会在最终分析时把trymax提到 100 甚至 200虽然耗时增加但换来的是更低的 stress 和更稳定的排序。set.seed()也一定要在运行前固定否则每次结果都有细微旋转差异后续的 GAM 映射图和文字描述就对不上了。3. 一般加性模型凭什么能映射环境梯度从样条到等值线3.1 一个公式理解 GAM 映射NMDS 排序图本身只有样点和物种的位置环境变量如何在这个空间中变化需要额外建模。一般加性模型GAMgeneralized 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 是连续数值变量可以直接进 GAMMoisture 是有序因子适合用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$stressdecostand()做 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) testprotest()里的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 映射的流程能帮你在自己的数据上少走弯路。本文还有配套的精品资源点击获取 NMDS与GAM映射：生态学排序图制作与解释实战指南
简介一份面向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 与 trymaxmetaMDS()的两个参数我见过大量使用者没搞清楚就默认跑完了。第一个是autotransform默认值为 TRUE。它会在计算距离前自动对丰度数据做 Wisconsin 双标准化先把每个物种的丰度除以该物种在所有样点的总量得到相对优势度再把每个样点的所有值按样点总量标准化。这个变换对常见生态数据效果不错但它是「静默地」修改了你的距离矩阵。如果你提前用 Hellinger 转换或其他方法预处理过数据就必须把autotransform设为 FALSE否则等于做了两遍变换结果解释起来会很别扭。第二个是trymax。默认值只有 20意味着只尝试 20 次随机起始。数据量大、群落结构复杂时20 次常常找不到低 stress 的解运行结束会报「结果未收敛」的警告。我一般会在最终分析时把trymax提到 100 甚至 200虽然耗时增加但换来的是更低的 stress 和更稳定的排序。set.seed()也一定要在运行前固定否则每次结果都有细微旋转差异后续的 GAM 映射图和文字描述就对不上了。3. 一般加性模型凭什么能映射环境梯度从样条到等值线3.1 一个公式理解 GAM 映射NMDS 排序图本身只有样点和物种的位置环境变量如何在这个空间中变化需要额外建模。一般加性模型GAMgeneralized 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 是连续数值变量可以直接进 GAMMoisture 是有序因子适合用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$stressdecostand()做 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) testprotest()里的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 映射的流程能帮你在自己的数据上少走弯路。本文还有配套的精品资源点击获取