研究生态、环境、土壤或微生物数据的同学,大概率都遇到过这样一个困境:变量之间的关系明明是一条因果链,环境因子影响土壤资源,土壤资源影响植物多样性,植物多样性最终影响生态系统功能。这种问题用传统的回归分析一次只能处理一个“因变量”,做多条回归又没法检验整条路径是否合理。等你想到结构方程模型(SEM),打开lavaan一跑,又发现数据不满足多元正态分布、样本量不够大、样方之间还不独立。模型要么运行报错,要么结果不收敛。
正文将围绕这个痛点展开。分段结构方程模型(piecewiseSEM)是解决上述问题的一套成熟方案,它在R语言中有专门实现piecewiseSEM包。与传统SEM不同,分段SEM把全局模型拆成多个局部回归模型,允许你混合使用lm、glm、lme、lmer等不同类型模型,对非正态数据、嵌套结构、空间自相关都有更好的兼容性。读完这篇文章,你会理解分段SEM的基本原理,掌握R语言实现步骤,并能独立完成结果解读。
1. 为什么生态分析需要分段结构方程模型
结构方程模型本质上是一个多方程框架,传统SEM通过最大似然估计同时拟合整个变量之间的方差-协方差矩阵,一次得到所有路径系数和全局拟合优度。这种方法在社会科学、心理学等领域非常成熟,因为那些数据通常来自受控实验或问卷,变量满足或近似满足多元正态分布,样本量也比较大。
但生态学、环境科学、农学等领域的观测数据,经常有以下特点:
第一,许多响应变量并不服从正态分布。比如物种丰富度是计数数据,土壤重金属浓度往往是偏态分布,土壤酶活性可能是正偏态数据。传统SEM对分布假设非常敏感,强行用线性模型拟合会得到错误的显著性检验。
第二,数据存在固有的嵌套或空间结构。比如多个样方属于同一个样地,同一样地内的样方不独立;或者采样点空间距离很近,存在空间自相关。传统SEM通常假设样本之间独立,一旦这个假设被打破,路径系数的标准误会被低估,p值也会偏小,容易得出“显著”但不可靠的结果。
第三,全局估计对样本量要求较高。传统SEM要求每个自由参数有足够多的样本支撑,在样本量较小或变量较多时,模型容易不收敛。而分段SEM逐个拟合局部模型,在样本量不理想的情况下也能工作。
分段结构方程模型的关键思路是:不直接拟合整个因果网络,而是把它拆成若干个局部模型,每个内生变量单独建立一个回归方程。这些局部模型可以是线性回归、广义线性模型、混合效应模型、空间模型等。之后再对整体网络的拟合程度做统一检验。
这意味着它在保住了SEM“多路径、中介效应、因果网络”能力的同时,解决了传统SEM面对生态数据时的“水土不服”。从R实现上看,piecewiseSEM包支持的模型对象包括lm、glm、gls、lme、lmer、glmer、negbin等,覆盖面非常广。如果你正在做生态数据、环境因子对物种组成的影响、生物多样性与生态系统功能关系这类课题,分段SEM几乎是绕不开的分析工具。
2. 分段SEM的核心原理:从全局似然到局部估计
分段SEM不是传统SEM的简单替代,而是一套独立的建模逻辑。理解它的底层原理,比记住几个函数更重要,因为后续所有参数解释都建立在这个原理之上。
2.1 传统SEM与分段SEM的差异
传统SEM用一句话概括:把多个线性方程联立起来,通过极大似然估计同时求解所有未知参数,目标是让模型隐含的协方差矩阵尽量接近观测协方差矩阵。模型整体拟合好不好,通过卡方检验判断。即零假设是“模型隐含协方差矩阵与观测协方差矩阵相同”,当卡方检验的p值大于0.05时,说明模型没有被拒绝。
分段SEM则采用局部估计策略。先根据先验因果假设画出有向无环图(DAG),图中每个箭头代表一个因果关系。然后对每个内生变量分别建立回归模型,拟合方式完全自由:可以是普通最小二乘、广义最小二乘、混合效应模型、广义线性模型等。最后把这些局部模型组合成一个整体结构,再用专门方法检验整体拟合度。
2.2 d分离检验与Fisher's C统计量
分段SEM最核心的检验方法是有向分离检验(d-separation test)。这个名字来自图论中的“有向分离”概念。
假设你有一个因果网络,比如“环境因子A -> 土壤资源B -> 植物多样性C -> 生产力D”。在这个网络中,如果给定中间变量B,A和C应该是条件独立的;如果给定B和C,A和D应该是条件独立的。这些没有直接连线的变量对,称为“缺失路径”。
d分离检验做的事情是:找出图中所有缺失路径,对每一条缺失路径在给定其父节点或祖先节点的条件下做独立性检验。如果某条缺失路径的p值很小,说明变量之间仍然存在显著相关性,最大可能的原因是模型遗漏了一条直接路径。反过来,如果所有缺失路径检验都不显著,说明现有网络结构没有明显遗漏。
把这些独立性检验的p值代入Fisher's C统计量公式:
其中 k 是独立性检验个数,p_i 是每个检验的p值。可以证明,在模型正确且数据满足假设时,C服从自由度为2k的卡方分布。如果整体p值大于0.05,说明模型整体与数据吻合,没有被拒绝。
需要注意的是,这里的“整体p > 0.05”只能说明模型没有明显遗漏路径,不能说明模型是“真实”或“最优”的。因为它只检验了缺失路径,并没有验证路径方向是否正确。方向问题依靠先验知识和多模型比较。
2.3 两种思路的对比
| 维度 | 传统SEM(如lavaan) | 分段SEM(piecewiseSEM) |
|---|---|---|
| 估计方式 | 全局方差-协方差矩阵拟合 | 局部模型分别估计 |
| 分布假设 | 多元正态分布为主 | 支持多种分布,通过glm/glmer扩展 |
| 独立性假设 | 要求样本独立 | 可使用混合效应模型处理嵌套/相关 |
| 样本量要求 | 较高,参数越多越明显 | 相对宽容,局部模型更稳健 |
| 空间自相关 | 难以直接处理 | 可在各模型中添加相关结构 |
| 整体拟合指标 | 卡方、CFI、TLI、RMSEA | Fisher's C、AIC |
| 灵活性 | 模型形式相对固定 | 可组合多种模型类型 |
这张表能直观看出分段SEM的优势集中在生态数据的“不规整”上。但也要注意,传统SEM具备潜变量建模能力,可以处理测量误差,这是分段SEM目前的短板。后续如果变量测量误差较大,仍需要考虑传统SEM或贝叶斯SEM方案。
3. R环境准备与piecewiseSEM安装
本文演示使用的环境是R 4.x + RStudio,操作系统Windows、macOS或Linux均可,代码没有平台依赖。
安装piecewiseSEM包推荐直接从CRAN安装:
install.packages("piecewiseSEM")因为分段SEM会调用大量模型拟合函数,建议同时安装常用依赖包:
install.packages(c("nlme", "lme4", "lmerTest", "MuMIn", "DiagrammeR"))加载包:
library(piecewiseSEM) library(nlme) library(lme4)如果在安装时遇到编译问题,常见原因是R版本较旧或缺少系统编译工具。Windows用户建议安装与当前R版本匹配的Rtools;macOS用户可以检查Xcode Command Line Tools是否完整。对于大多数用户,直接安装预编译二进制包即可,通常不会遇到额外问题。
版本方面,本文演示基于当前CRAN正式版。piecewiseSEM不同版本在输出格式上有细微差异,比如旧版用fisherC()函数,新版部分函数被整合进summary(),但核心API没有变化。只要代码能运行,输出内容解读逻辑一致。
4. 模拟数据与数据结构准备
为了完整演示建模流程,我们构造一份模拟数据。这种做法的好处是:路径系数已知,可以验证piecewiseSEM能否正确还原数据生成过程。真实项目数据通常只存储一份,不方便反复尝试,模拟数据则完全可控。
模拟情景设定如下:
- 样方内环境异质性(env)影响土壤资源可用性(soil)。
- 土壤资源可用性影响植物多样性(plant_div)。
- 植物多样性直接影响生产力(productivity)。
- 环境异质性也能直接影响生产力,构成一条直接路径。
- 数据来源于30个样地,每个样地4个样方,样方之间在样地内存在随机效应。
生成代码如下:
set.seed(123) # 固定随机种子保证结果可复现 n_site <- 30 n_plot <- 4 n <- n_site * n_plot site <- factor(rep(1:n_site, each = n_plot)) # 环境异质性:正态分布观测变量 env <- rnorm(n, mean = 50, sd = 10) # 土壤资源:受环境影响,增加随机误差 soil <- 0.6 * env + rnorm(n, mean = 0, sd = 5) # 植物多样性:受土壤资源影响 plant_div <- 0.7 * soil + rnorm(n, mean = 0, sd = 3) # 生产力:受植物多样性影响 + 环境直接效应 productivity <- 0.5 * plant_div + 0.3 * env + rnorm(n, mean = 0, sd = 4) # 生成数据框 dat <- data.frame(site = site, env = env, soil = soil, plant_div = plant_div, productivity = productivity) # 查看数据结构 head(dat) str(dat)这份数据中,样地(site)是分组变量,样方与样方之间在同一样地内共享环境背景,这里的随机效应没有在生成公式中体现,我们后面示范如何加入随机截距模型。
实际项目里,这一步对应的是数据清洗流程:检查缺失值、确认变量类型、处理异常值、检验共线性。如果原始数据存在大量缺失,建议先用mice或missForest做缺失值插补,再进入模型构建阶段。不要直接拿原始数据建模,否则分段SEM的每个局部模型样本量可能不一致,导致结果可解释性下降。
5. piecewiseSEM完整建模流程示例
下面逐步演示如何用piecewiseSEM构建、拟合、检验分段SEM。
5.1 构建基础分段SEM
先构建没有随机效应的基础模型,把所有变量当作独立样本处理:
model_lm <- psem( lm(soil ~ env, data = dat), lm(plant_div ~ soil, data = dat), lm(productivity ~ plant_div + env, data = dat) ) summary(model_lm)这里psem()是piecewiseSEM的核心函数,它接受一组模型对象,然后把它组合成分段SEM。每个lm()对应一个内生变量,公式中的解释变量就是该内生变量的所有直接父节点。
运行summary()后,屏幕输出主要包含三部分:
- 每个局部模型的路径系数、标准误、p值和R方。
- 缺失路径的d分离检验结果。
- 整体模型拟合指标:Fisher's C、自由度、p值、AIC。
从示例代码看,数据生成时土壤对植物多样性的真实系数为0.7,植物多样性对生产力的真实系数为0.5,环境对生产力的真实直接效应为0.3。模型的估计值应该接近这些真值,这可以验证建模过程是否正确。
5.2 引入随机效应处理嵌套结构
由于原始数据结构中同一site包含4个样方,样方并不完全独立。在真实数据中,如果不考虑这种嵌套结构,路径系数的标准误会偏小。分段SEM的很大价值就在于能直接加入混合效应模型。
将三个局部方程分别改为lme模型,并且加入样地随机截距:
model_mixed <- psem( lme(soil ~ env, random = ~ 1 | site, data = dat), lme(plant_div ~ soil, random = ~ 1 | site, data = dat), lme(productivity ~ plant_div + env, random = ~ 1 | site, data = dat) ) summary(model_mixed)这里使用nlme包中的lme()函数。随机截距的含义是:假设每个样地有自己的背景水平,样地内部样方在截距上有相关性。
如果你更习惯lme4语法,也可以写成:
model_mer <- psem( lmer(soil ~ env + (1 | site), data = dat), lmer(plant_div ~ soil + (1 | site), data = dat), lmer(productivity ~ plant_div + env + (1 | site), data = dat) ) summary(model_mer)两种方式结果基本一致,具体选择取决于你更熟悉哪个函数。如果需要同时处理多个随机效应或更复杂的随机斜率,lme4语法更灵活。
5.3 提取系数与标准化系数
summary()输出已经能给出路径系数,但如果你需要把系数整理成表格,可以使用coefs()和stdCoefs():
# 提取原始尺度系数 raw_coef <- coefs(model_mixed) print(raw_coef) # 提取标准化系数,方便比较不同路径影响大小 std_coef <- stdCoefs(model_mixed) print(std_coef)标准化系数消除了自变量量纲影响,在生态学论文中更常用来比较“哪个路径影响更大”。例如土壤资源对植物多样性的标准化系数是0.7,环境异质性对生产力的标准化系数是0.25,那么可以说前者效应更强。
5.4 模型比较与路径图绘制
如果你同时构建了多个候选模型,可以通过AIC比较。summary()输出中会包含AIC值,也可以用AIC()或MuMIn::AICc()获取:
AIC(model_lm, model_mixed)AIC越小代表模型在拟合与简洁性之间的平衡越好。当候选模型数量较多时,推荐使用AICc(小样本校正版本),避免样本量较小时AIC偏保守。
piecewiseSEM还提供了plot()方法,可以在RStudio Viewer中查看路径图:
plot(model_mixed)这个函数依赖于DiagrammeR包,输出的是有向图。图中每个变量一个节点,路径箭头旁标注了路径系数。如果要在论文中引用,建议导出为PNG或SVG后再调整。
6. 运行结果解读与拟合检验
下面重点解释summary()输出中的每一块内容,因为许多读者第一次接触分段SEM,面对大量输出会感到无从下手。
6.1 局部模型部分
输出首先列出每个内生变量的回归结果。比如:
Response: plant_div Predictor Estimate Std.Error DF Crit.Value P.Value Std.Estimate soil 0.689 0.051 118 13.51 <0.001 0.775这里“predictor”表示该内生变量的预测变量,“estimate”是未标准化系数,“std.estimate”是标准化系数,“crit.value”是t值或z值。DF是自由度,混合效应模型的DF计算方式与普通回归不同,无需过度关注细节。
每个局部模型的底部会给出R方,表示该内生变量被其解释变量解释的比例。R方不是判断模型好坏的唯一标准,但要注意R方太低说明该路径的预测力弱,可能存在重大遗漏变量。
6.2 d分离检验部分
接着输出的是dSep检验,它有点像模型修正指数。例如:
Independ.Claim Test.Type DF Crit.Value P.Value env ~ plant_div + soil ... 0.231 soil ~ productivity + plant_div + env ... 0.544这些p值全部大于0.05,说明这些缺失路径都不显著。如果某个缺失路径p值小于0.05,说明数据中存在一条你尚未纳入模型的显著关系,这时需要回到假设阶段考虑是否加入该路径。
6.3 整体拟合指标部分
输出的最后一段是关键:
Fisher's C = 2.57, df = 4, P-value = 0.632 AIC = 54.31Fisher's C统计量的p值大于0.05,说明模型整体没有被拒绝。这里的“P-value”越大越好,意味着缺失路径不显著,现有因果结构可以解释数据。
还要关注df。它等于2乘以缺失路径数量。如果df为0,说明模型是饱和模型,Fisher's C无法计算或没有意义,这种情况通常发生在所有变量之间都有直接路径时。
6.4 拟合检验失败的处理
如果Fisher's C检验p值小于0.05,模型整体被拒绝,说明图中某条缺失路径实际上存在显著关系。处理顺序建议:
先看dSep检验输出中哪条路径检验的p值最小,这个路径最可能被遗漏。然后基于生态学或机理知识判断这条路径是否应该加入,不要机械地添加。添加路径后重新拟合模型,再看Fisher's C是否改善。
如果添加路径后仍然被拒绝,可能是某个局部模型本身分布假设错误,比如对计数数据用了正态分布,此时需要改成广义线性模型。
7. 常见问题与排查方法
分段SEM在R中运行并不复杂,但实际使用中会遇到几类高频问题,下面用表格总结。
| 问题现象 | 可能原因 | 排查方式 | 解决方案 |
|---|---|---|---|
| 安装piecewiseSEM失败 | 依赖包未安装,或R版本较旧 | 查看报错信息,确认R版本 | 先安装nlme、lme4等依赖包;更新R到最新版;Windows用户安装Rtools |
psem()报错,模型对象不被支持 | 传入了bamlss、brms等不在支持列表中的模型对象 | 查看psem函数文档 | 改用支持的lm、glm、lme、lmer、glmer等 |
| 整体Fisher's C检验p值小于0.05 | 模型遗漏重要路径或局部模型分布错误 | 查看dSep检验结果,对比缺失路径p值 | 根据显著性结果和先验知识谨慎添加路径 |
| 某个缺失路径检验p值很小 | 图中缺少一条直接因果关系 | 检查该路径的生态学依据 | 添加路径并重新拟合,再评估整体拟合变化 |
| 路径系数标准误为NaN或极大 | 解释变量之间存在严重共线性,或样本量太小 | 计算VIF,检查样本量 | 删除高度相关的变量,或合并变量 |
| 标准误偏小,p值过于显著 | 数据存在嵌套或空间结构但模型没有处理 | 检查数据结构,确认是否存在固定分组 | 换成lme/lmer模型,加入随机截距 |
| 非正态计数数据拟合效果差 | 仍在使用lm模型 | 检查响应变量分布 | 使用glm的family=poisson/negative.binomial,或glmer |
| 样点间存在空间自相关 | 模型没有设置空间相关结构 | 做残差空间自相关检验(如Moran's I) | 在gls/lme中加入correlation参数(如corExp、corSpher) |
这些问题的排查逻辑是:先怀疑模型结构,再怀疑分布假设,最后检查数据结构。不要一上来就删变量。
如果在lme()中加入空间相关参数,示例写法如下:
library(nlme) model_sp <- psem( gls(soil ~ env, correlation = corExp(form = ~ x + y), data = dat), lme(plant_div ~ soil, random = ~ 1 | site, correlation = corExp(form = ~ x + y), data = dat), lme(productivity ~ plant_div + env, random = ~ 1 | site, correlation = corExp(form = ~ x + y), data = dat) )这里要求数据框中有坐标变量x和y,并且建模时传入完整数据。空间相关结构对大数据集计算较慢,但如果采样点确实存在空间自相关,忽略它会严重低估标准误。
8. 最佳实践与论文报告建议
分段SEM在论文中越来越多见,但很多使用者只关注最终p值是否显著,忽略了建模过程中的关键细节。下面几条建议可以大幅提升分析的可靠性。
第一,先画因果假设图,再写代码。分段SEM的本质是拟合你预设的因果网络,而不是自动寻找最佳网络。作图的过程会迫使你思考每个箭头的含义、方向、是否存在反向因果。用DiagrammeR先画图,再根据图写psem()中的公式,能避免漏写路径。
第二,不要为了得到“好结果”反复加路径。分段SEM的d分离检验给出的是“缺失路径是否显著”的判断,修改模型应该基于机理假设,而不是机械地追求p值。每添加一条路径,都应该能解释为什么这两个变量存在直接因果关系。否则模型只是过拟合了当前数据,在新数据集上很容易崩溃。
第三,关注局部模型的残差诊断。分段SEM由多个局部回归组成,应用回归诊断的通用标准检查每个模型的残差:残差是否正态、是否异方差、是否有异常值、是否空间相关。你可以分别对每个局部模型做plot()和residuals()检查,不要只看整体拟合指标。
第四,报告完整模型信息。在论文方法部分,至少应报告:观测数据条数、变量定义、每个局部模型的类型及其分布家族、随机效应结构、Fisher's C值、自由度、p值、AIC值。如果进行了模型比较,应列出各候选模型AIC。
第五,用标准化系数汇报效应大小。不同变量的量纲差异很大,未标准化系数不能直接比较效应强弱。stdCoefs()输出的标准化系数更适合作为论文正文中的路径效应量。
第六,做好可重复性管理。在开头设定set.seed(),把数据清洗、模型拟合、结果汇总写成独立R脚本,用RMarkdown或R脚本记录每次分析版本。分段SEM中的随机效应模型在非独立数据中容易受到随机数种子影响,固定种子能保证结果可复现。
第七,关于样本量。分段SEM对样本量的要求比传统SEM宽容,但这不意味着小样本可以随便跑。一般建议每个路径参数至少有10到20个样本支持。如果总样本量只有几十个,即使分段SEM能跑出结果,稳定性也值得怀疑。此时更适合使用简化模型或贝叶斯SEM,配合正则化先验。
9. 总结与后续学习方向
本文围绕分段结构方程模型,重点讲了几个关键点:传统SEM在生态数据场景下为什么容易失效,分段SEM通过局部估计和d分离检验解决了什么问题,以及如何使用R语言piecewiseSEM包完成从数据准备、模型构建、结果解读到模型比较的完整流程。
理解分段SEM最核心的收获,不是会调用几个函数,而是理解“因果网络建模”和“模型诊断”的思路。一个分段SEM模型的价值,取决于你的因果假设质量、变量测量质量和对缺失路径的合理判定。工具本身只是把假设转化成可检验的形式。
这篇内容适合先收藏再动手练习。建议用模拟数据跑通全文流程,再把自己的数据代入,逐项检查输出结果。如果模拟数据的结果与预期不符,优先检查路径方向和数据分布。
下一步可以继续学习三个方向:
一是piecewiseSEM的更多高级用法,比如处理交互项、分段模型中加入二次项、多个内生变量的随机斜率、借用sem.fit()做更复杂的模型比较。
二是如果你需要处理潜变量和测量误差,传统SEM中的lavaan仍然是不可替代的工具,可以把两种方法结合使用:先做测量模型评估潜变量,再通过显变量或因子得分进入分段SEM。
三是贝叶斯SEM,使用brms或stan可以灵活指定先验、处理非正态数据、处理缺失数据,并且能获得参数的后验分布。当数据复杂程度超出分段SEM范围时,贝叶斯路线是一个自然延伸。
在实际项目中,建议把分段SEM当作分析工具箱里的一个选项,而不是所有问题的默认答案。你需要先用领域知识回答“这些变量之间是否存在因果关系”“是否遗漏了关键变量”,再决定用哪种统计框架。数据方法和因果理解是两条腿,一起走,分析结果才真正有意义。