企业官网建设流程全解析

简介：这份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查共曲。这套流程跑下来，图能不能用、结论站不站得住，心里基本有数。希望帮到你。

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