vegan包envfit单列与整表分析差异及单独分析合理性问询

问题背景 #

我在R语言中使用vegan包的envfit函数，尝试将环境变量绘制至苔藓群落的NMDS图时遇到问题：NMDS图中群落按「2015年多年冻土（浅蓝）、2015年无多年冻土（浅橙）、2021年多年冻土（深蓝）、2021年无多年冻土（深橙）」分组。

单变量分析过程与结果 #

此前从两个环境数据表（combo.env.15.21和veg.TOT2.15.21）中选取7个目标变量，单独对每个变量执行envfit分析，代码如下：

elevationNMDS<-envfit(live.moss.mds~Elev_start, data=combo.env.15.21, perm=999, choices=c(1,2), display="sites", na.rm=TRUE, arrow.mul = 1, add = TRUE)
plot(elevationNMDS, col = "green", cex = 1.0, label = "elevation")

latitudeNMDS<-envfit(live.moss.mds~Lat_start, data=combo.env.15.21, perm=999, choices=c(1,2), display="sites", na.rm=TRUE, arrow.mul = 1, add = TRUE)
plot(latitudeNMDS, col = "magenta", cex = 1.0, label = "latitude")

orgNMDS<-envfit(live.moss.mds~avg_organic_depth, data=veg.TOT2.15.21, perm=999, choices=c(1,2), display = "sites", na.rm = TRUE, arrow.mul = 1, add = TRUE)
plot(orgNMDS, col="saddlebrown", lwd=10, cex = 1.0, labels = "organic")

moistNMDS<-envfit(live.moss.mds~moisture_class_num, data=veg.TOT2.15.21, perm=999, choices=c(1,2), display = "sites", na.rm = TRUE, arrow.mul = 1, add = TRUE)
plot(moistNMDS, col="purple",cex = 1.0, labels = "moisture")

actNMDS<-envfit(live.moss.mds~active_layer_depth, data=veg.TOT2.15.21, perm=999, choices=c(1,2), display = "sites", na.rm = TRUE, arrow.mul = 1, add = TRUE)
plot(actNMDS, col="turquoise",cex = 1.5, labels = " ")

permNMDS<-envfit(live.moss.mds~permafrost_num, data=combo.env.15.21, perm=0, choices=c(1,2), display = "sites", na.rm = TRUE)
plot(permNMDS, col="royalblue",cex = 1.5, labels = "permafrost")

yearNMDS<-envfit(live.moss.mds~year, data=combo.env.15.21, perm=0, choices=c(1,2), display = "sites", na.rm = TRUE)
plot(yearNMDS, col="black",cex = 1.5, labels = "year")

此时year箭头呈水平方向、permafrost箭头呈垂直方向，符合群落分组的逻辑预期。

整表分析过程与结果 #

之后将两个数据表通过唯一ID合并，筛选出目标7个变量生成envdataTEST，代码如下：

envdataTEST <- merge(combo.env.15.21,veg.TOT2.15.21,by.x ="X", by.y = "unique_id")
rownames(envdataTEST)<-combo.env.15.21[,2]
envdataTEST<-envdataTEST[c(3,5,11,56,58,116,123)]

再对整表执行envfit分析，代码如下：

allenvNMDS <- envfit(live.moss.mds, envdataTEST, perm=999, choices=c(1,2), display = "sites", na.rm = TRUE)
plot(allenvNMDS, col = "red")

结果箭头完全不符合预期，year和permafrost箭头方向错误。

现咨询：vegan包的envfit函数在单独分析单列变量与整表分析时，是否执行不同操作？若单独分析的结果更符合数据逻辑，是否可以采用该分析方式？

问题解答 #

1. 单变量与整表分析的操作差异 #

• 单变量分析：每次仅拟合单个环境变量与NMDS排序轴的关系，本质是基于该变量的单变量线性回归（或相关分析，依变量类型而定），结果仅反映该变量自身与排序轴的独立关联，不受其他环境变量的干扰。
• 整表分析：同时拟合所有输入环境变量与NMDS排序轴的关系，此时函数执行的是多元线性回归，每个变量的系数是控制其他变量后的偏相关效应。如果变量间存在共线性（比如year或permafrost_num与其他环境变量存在关联），会显著改变单个变量的回归系数，导致箭头方向、长度与单变量分析结果出现偏差，这也是你遇到的核心问题。
• 补充：你单变量分析中对year和permafrost_num设置了perm=0（不做置换检验），而整表分析用了perm=999，但这不是箭头方向差异的主要原因，核心是多元拟合的偏效应。

2. 是否可以采用单变量分析的结果 #

• 如果单变量分析的结果更符合你的生态学逻辑，且你的研究目标是展示每个环境变量独立对群落结构的解释能力，而非变量间的交互或控制效应，完全可以采用单变量分析的结果。
• 但需要注意：单变量分析无法排除变量间共线性的干扰，若多个变量存在较强相关，单独拟合的结果可能存在重复解释。如果要呈现更严谨的多元效应，建议先做以下处理：
  • 用cor()函数计算变量间的相关系数，或用VIF（方差膨胀因子）检验共线性；
  • 移除高度相关的变量后再进行整表分析；
  • 或选择生态学意义最关键的变量单独展示，避免结果冗余。

内容的提问来源于stack exchange，提问作者jane.f