★ FEATURED ARTICLE 简介这份PDF教程面向生态学、环境科学及社会科学领域的数据分析人员聚焦R语言中非度量多维标尺排序NMDS与一般加性模型GAM映射的完整实现路径。内容围绕vegan包的metaMDS()函数展开涵盖物种组成、环境因子与样点分类三类数据集的准备方式并演示如何借助ggplot2绘制带凸包边界的排序图再通过mgcv包的gam()函数将环境变量以非线性光滑项映射到NMDS轴上配合gam.check()评估模型合理性。资源包为单个PDF文件大小约45.43MB图文并茂地呈现代码与结果解读便于对照练习。目前已有1206人学习下载适合希望掌握降维排序与非线性环境效应分析、提升生态统计可视化能力的中级R语言用户参考。1. 从一张排序图说起NMDS 为什么在生态数据里比 PCA 更抗造做过群落数据的人大概都有过这种经历辛辛苦苦跑完 PCA第一轴解释率只有 12%图上一团糊审稿人还追问“你这排序能反映梯度吗”。问题往往不在数据而在方法选错了。物种组成矩阵里全是零、丰度分布严重右偏、样点之间还夹着缺失值这种数据天然不满足 PCA 对欧氏距离和线性的假设。非度量多维标尺排序NMDS走的是另一条路——它不要求距离矩阵有欧氏性质只保留样点之间的排位顺序把“谁和谁更像”这件事压缩到低维空间里。换句话说它关心的是秩序不是绝对距离。这份 R 语言教程围绕的就是这条链路用vegan跑 NMDS用mgcv的一般加性模型GAM把环境因子以等值线形式映射回排序空间再用ggplot2把分组凸包、样点、环境梯度叠在一张图上。它适合做群落生态、微生物组、环境监测的从业者尤其是手里有“物种矩阵 环境因子 分组信息”三张表、想在一张排序图上同时讲清分组差异和环境驱动的人。下面按数据准备、NMDS 执行、GAM 映射、避坑、进阶五段拆开讲每一步都给可抄的代码和参数解释。2. 三张表的准备与 metaMDS 执行从原始矩阵到排序坐标2.1 物种矩阵、环境因子、分组表的格式约定这套流程对输入格式有硬性要求格式错了后面全崩。第一张是物种组成数据每一列一个物种每一行一个样点单元格是丰度或相对丰度行名必须是样点 ID。第二张是环境因子同样行是样点、列是环境变量行名要和物种矩阵严格对齐。第三张是样点分类表两列即可样点 ID 和分组标签。三张表共用一个样点 ID 体系这是后面能把 GAM 预测面映射回排序图的前提。常见做法是把三张表都读成 data.frame然后强制统一行名。我一般会先做一次all(rownames(species) rownames(env))的断言不通过就先停下来对齐别急着往下跑。# 读取三张表row.names 1 表示第一列是样点ID species - read.csv(species.csv, row.names 1, check.names FALSE) env - read.csv(env.csv, row.names 1, check.names FALSE) group - read.csv(group.csv, row.names 1, check.names FALSE) # 关键断言三张表的样点顺序必须一致 stopifnot(all(rownames(species) rownames(env))) stopifnot(all(rownames(species) rownames(group))) # 物种矩阵里不能有全零样点否则 Bray-Curtis 会出 NaN species - species[rowSums(species) 0, ]check.names FALSE是为了防止 R 把物种名里的特殊字符自动改写物种名一旦被改后面匹配就乱了。rowSums(species) 0这步过滤掉空样点因为 Bray-Curtis 距离在空样点之间没有定义metaMDS会直接报错或产出 NaN 坐标。2.2 metaMDS 的参数选择距离、迭代与维度metaMDS()是vegan包的核心函数它内部做了两件容易被忽略的事一是用多次随机起始配置反复跑取 stress 最小的解二是对数据做平方根或其他转换来降低高丰度物种的权重。默认距离是 Bray-Curtis这对群落数据基本够用但如果你的数据是存在-缺失矩阵改成jaccard更合适。library(vegan) # 设置随机种子保证每次跑出来的起始配置一致结果可复现 set.seed(123) nmds_result - metaMDS( species, distance bray, # 群落数据默认用Bray-Curtis k 2, # 降到2维方便画图 trymax 100, # 最多尝试100次随机起始直到收敛 autotransform TRUE,# 自动对数据进行转换降低高丰度权重 trace FALSE # 不打印每次迭代过程 ) # 查看stress值和收敛情况 nmds_result$stress stressplot(nmds_result)k 2是为了出二维排序图如果只为看 stress 可以试k 3但三维图不好读。trymax 100是保险值正常数据几十次内就收敛。stress值是判断排序可信度的核心指标小于 0.1 算好0.1 到 0.2 可接受超过 0.2 这张图基本不能用来下结论。stressplot()画的是排序距离和原始距离的散点点越贴近那条非线性拟合线说明排序保留的秩序越多。2.3 提取坐标并合并分组信息metaMDS跑完后坐标要用scores()提取默认返回的是样点在 NMDS1、NMDS2 上的得分。这一步要把坐标和分组表拼起来拼的时候靠行名匹配不能靠行顺序否则分组会错位。# 提取样点坐标 NMDS.point - as.data.frame(scores(nmds_result, display sites)) NMDS.point$sample_id - rownames(NMDS.point) # 用样点ID合并分组而不是按行顺序cbind NMDS.point - merge(NMDS.point, group, by.x sample_id, by.y row.names) colnames(NMDS.point)[ncol(NMDS.point)] - grp # 检查合并后行数是否和原始样点数一致 nrow(NMDS.point) nrow(species)用merge而不是cbind是血泪经验只要中途过滤过样点行顺序就可能和分组表对不上cbind会静默地贴错标签图上分组全乱但代码不报错。合并后做一次行数断言能挡住大部分低级错误。3. 凸包绘制与 GAM 环境因子映射把梯度叠到排序图上3.1 用 chull 提取分组边界并绘制凸包分组凸包的作用是让读者一眼看出各组的空间范围。chull()是 R 基础函数输入一组二维坐标返回构成最小凸多边形的点索引。对每个分组单独取凸包点再合并成一个数据框就能用geom_polygon一次性画出所有分组的边界。library(ggplot2) # 对每个分组提取凸包边界点 grp.a - NMDS.point[NMDS.point$grp A, ][chull(NMDS.point[NMDS.point$grp A, c(NMDS1, NMDS2)]), ] grp.b - NMDS.point[NMDS.point$grp B, ][chull(NMDS.point[NMDS.point$grp B, c(NMDS1, NMDS2)]), ] grp.c - NMDS.point[NMDS.point$grp C, ][chull(NMDS.point[NMDS.point$grp C, c(NMDS1, NMDS2)]), ] hull.data - rbind(grp.a, grp.b, grp.c) # 基础排序图凸包 样点 p_base - ggplot() geom_polygon(data hull.data, aes(x NMDS1, y NMDS2, fill grp), alpha 0.7) geom_point(data NMDS.point, aes(x NMDS1, y NMDS2, shape grp), size 1.5) scale_fill_manual(values c(#fce8b7, #D2E2FE, #b6e7da), name 采样区域) scale_shape_manual(values c(15, 16, 17), name 采样区域) theme_bw()chull返回的是点的行索引所以NMDS.point[...][chull(...), ]这种写法是先按分组筛选再在筛选结果里取凸包点。alpha 0.7让凸包半透明避免盖住下面的样点。scale_fill_manual和scale_shape_manual分别控制填充色和点形状两套映射用同一个分组变量图例会自动合并。3.2 GAM 拟合环境因子与排序轴的非线性关系环境因子不是直接画在排序图上的而是先对每个环境变量和 NMDS 轴建 GAM再用模型预测值生成网格上的等值面。mgcv包的gam()里s()是光滑项用来拟合非线性关系k参数控制光滑度的上限太小会欠拟合太大会过拟合。library(mgcv) # 以NMDS1为响应环境变量为光滑项建GAM gam1 - gam(NMDS1 ~ s(env_var1, k 5) s(env_var2, k 5), data cbind(NMDS.point, env)) summary(gam1) gam.check(gam1)k 5是常见起点数据点少时用 3 到 4点多可以到 6 到 8。summary()看每个光滑项的edf有效自由度和p值edf接近 1 说明关系接近线性越大越弯。gam.check()会输出残差诊断图重点看 QQ 图和残差是否随拟合值有系统趋势有明显喇叭形或弯曲就说明模型设定有问题。3.3 生成预测网格并映射等值线要把 GAM 的结果画成等值线需要先在 NMDS1、NMDS2 平面上生成一张规则网格对网格上每个点用 GAM 预测环境值再把预测结果交给stat_contour或geom_text_contour。这一步的坐标范围要覆盖所有样点否则等值线会在边缘断掉。# 生成覆盖样点范围的网格 x_seq - seq(min(NMDS.point$NMDS1), max(NMDS.point$NMDS1), length.out 100) y_seq - seq(min(NMDS.point$NMDS2), max(NMDS.point$NMDS2), length.out 100) grid - expand.grid(NMDS1 x_seq, NMDS2 y_seq) # 用GAM预测网格上每个点的环境值 grid$z - predict(gam1, newdata grid) # 把预测面叠到排序图上 p_map - p_base stat_contour(data grid, aes(x NMDS1, y NMDS2, z z), color grey40) geom_text_contour(data grid, aes(x NMDS1, y NMDS2, z z), skip 1, breaks seq(80, 120, by 1), check_overlap TRUE, size 2.5)length.out 100决定网格密度太稀等值线会锯齿太密计算慢但图更顺。predict(gam1, newdata grid)用的是 GAM 的预测接口注意newdata的列名必须和建模时一致。stat_contour画线geom_text_contour加数值标签skip 1表示每隔一条线标一个数breaks控制标注的数值范围check_overlap TRUE防止标签叠在一起。geom_text_contour来自metR包需要单独library(metR)。4. 避坑与排查NMDS 和 GAM 联用最容易翻车的五处4.1 stress 值偏高却强行解读现象nmds_result$stress跑到 0.25 以上图还是照画结论照下。原因通常是k设得太小、数据里噪声样点太多或者距离度量选错。解决先把k提到 3 看 stress 是否明显下降如果降幅大说明二维压得太狠再检查是否有极端样点用goodness()看每个样点的拟合优度把长期偏离的样点单独核查。stress 超过 0.2 的图我一般只用来做探索不进正文。4.2 分组凸包把不存在的区域也包进去现象某个分组的凸包覆盖了另一组样点所在的区域视觉上两组“重叠”了。原因chull取的是最小凸多边形如果组内样点分布呈凹形凸包会把凹口填平。解决凸包只适合组内分布大致凸的情况凹形分布改用geom_encircle或直接不画凸包只画点加置信椭圆。别为了图好看硬套凸包。4.3 GAM 预测面外推到无数据区域现象等值线延伸到排序图角落那里根本没有样点。原因网格范围按坐标轴极值生成但样点只占其中一部分GAM 在无数据区做的是外推预测值不可信。解决把网格范围收缩到样点坐标的分位数区间比如 5% 到 95%或者用geom_point的分布手动裁掉空白区。外推的等值线看着漂亮但审稿人一问就露馅。4.4 环境因子和排序轴方向对不上现象GAM 拟合出来的等值线梯度方向和预期相反。原因NMDS 轴的符号是任意的metaMDS每次跑出来的 NMDS1 可能整体翻转。解决跑之前set.seed()固定随机数跑完后如果方向不对手动把 NMDS1 或 NMDS2 乘以 -1 再合并环境数据。别在 GAM 里改改坐标本身更干净。4.5 物种矩阵未做转换导致高丰度物种主导现象排序图被几个高丰度物种拽着走稀有物种的信息全被淹没。原因metaMDS默认autotransform TRUE会做 Wisconsin 双重转换但如果手动设了autotransform FALSE又没自己做转换高丰度物种的权重会过大。解决保持autotransform TRUE或者手动做sqrt()、log1p()转换后再跑。转换方式要在方法里写清楚不同转换出来的排序图不能直接比。5. 进阶用 envfit 交叉验证 GAM 的梯度方向GAM 映射出的等值线是模型预测不是直接观测所以需要一条独立的证据链来验证梯度方向是否合理。vegan的envfit()做的是另一件事把环境因子向量拟合到排序坐标上返回每个环境因子的方向和显著性。把envfit的箭头和 GAM 的等值线叠在同一张图上如果箭头指向等值线升高的方向说明两条路径结论一致如果相反大概率是坐标轴翻转或 GAM 过拟合。# envfit 拟合环境因子向量 ef - envfit(nmds_result, env, permutations 999) ef # 提取显著环境因子的箭头坐标 ef_scores - as.data.frame(scores(ef, display vectors)) ef_scores$var - rownames(ef_scores) # 叠到排序图上 p_final - p_map geom_segment(data ef_scores, aes(x 0, y 0, xend NMDS1, yend NMDS2), arrow arrow(length unit(0.2, cm)), color red) geom_text(data ef_scores, aes(x NMDS1 * 1.1, y NMDS2 * 1.1, label var), color red, size 3)permutations 999是置换检验次数999 够用要更稳可以到 9999。scores(ef, display vectors)返回的是箭头终点坐标起点固定在原点。箭头长度代表环境因子与排序轴的相关系数强度方向代表梯度方向。把envfit箭头和 GAM 等值线放一起看是我每次出图前的固定动作——两条独立方法指向同一个方向这张图才敢往文章里放。还有一个常被忽略的验证点是 GAM 的concurvity()。当两个环境因子高度相关时各自的s()项会互相抢解释力edf和p值都不可信。跑一下concurvity(gam1, full TRUE)如果某个项的 concurvity 超过 0.8就考虑删掉一个或改用张量积te()。这个检查不花时间但能挡住“两个因子都显著但其实是一个在起作用”的误判。从那以后我每次跑完 NMDS 都强制走一遍先看 stress再看goodness()排异常样点然后envfit和 GAM 交叉验证梯度方向最后concurvity查共曲。这套流程跑下来图能不能用、结论站不站得住心里基本有数。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
