结构方程模型做过的人都知道,最耗时的往往不是统计本身,而是反复改模型。去年我帮某课题组审核一份投稿时,对方把量表总分直接当成显变量丢进模型,又说不清测量模型和结构模型的边界,根源就在于图形界面软件改起来太痛苦。我后来专门用R里的lavaan包重整理了整套分析流程——从潜变量、复合变量,到分组、嵌套、时间序列,再到非递归模型,一份脚本跑到底,再也不用来回画线。
这篇文章适合两类读者:一类是正被图形界面软件的各种弹窗和画线折磨的人,另一类是刚接触lavaan,想把常用SEM场景一口气串起来的人。我尽量不堆理论公式,重点讲实际操作中怎么选、怎么写语法、会碰到哪些坑,以及每类模型背后的识别逻辑和取舍思路。
1. 为什么我最后选了lavaan:从图形界面到代码流的迁移
1.1 图形界面软件的三个"卡脖子"时刻
先说三个让我彻底放弃图形界面软件的真实瞬间。第一个瞬间是我要给一个模型加两条路径,结果发现所有约束都要从头逐一点一遍,点完还要检查有没有漏掉,模型复杂到我自己都记不清哪些画过哪些没画过。第二个瞬间是分组分析,同一个模型要在界面里复制好几份,换一组数据又得重新连一遍,做测量不变性检验时简直像在做手工劳动。第三个瞬间是非递归模型,界面里画了一条从A到B、又一条从B到A的反馈线,软件直接报错,报错信息还看不出来到底是识别问题还是收敛问题。
不是说图形界面一无是处,它对初学者友好,画出来的路径图也确实直观。但一旦模型开始涉及多组比较、时间序列、潜增长、反馈回路这些稍复杂的设计,界面操作的成本会指数级上升。更重要的是,图形界面很容易让人忽略模型背后的识别条件和参数约束——你只是在"画"模型,而不是在"写"一个可复现的分析过程。
1.2 lavaan的语法骨架,其实只有三个符号
lavaan和那些软件的底层逻辑没有本质区别,都是对协方差矩阵做拟合,只是它把模型变成了一段可以被反复修改、批量执行的文本。初次接触会很别扭,但熬过第一周之后,你会发现自己再也不想回去。
它的核心运算符其实就几个:
=~表示潜变量的测量关系,读作"潜变量 =~ 它的指标",指向外行话叫"因子载荷"。~表示回归路径,左边是被解释变量,右边是解释变量。~~表示方差或协方差,比如潜变量之间的相关、残差之间的相关。~1表示截距,多组分析时做截距不变性会用。
一个最朴素的工作流长这样:
library(lavaan) # 测量模型:员工满意度由5个题项反映 model_cfa <- ' satisfaction =~ S1 + S2 + S3 + S4 + S5 ' fit_cfa <- cfa(model_cfa, data = df_sim) summary(fit_cfa, fit.measures = TRUE, standardized = TRUE)就是这么简单,模型写进字符串,cfa()跑起来,summary()看结果。后面的所有复杂模型,包括潜增长、交叉滞后、非递归,都是在这个字符串里做文章。代码就是模型,模型就是代码,这才是lavaan最大的价值。
2. 潜变量建模:先把测量模型做扎实,再谈路径
2.1 先跑CFA:验证性因子分析的基本盘
我见过的绝大多数SEM翻车案例,都不是死在结构路径上,而是死在测量模型上。测量模型没做好,后面一切拟合指标都只是数字游戏。
所以无论最终要检验多复杂的结构模型,我永远先单独跑一遍CFA。下面是一份我常用的完整测量模型,对应着模拟项目X的员工调研数据:
model_cfa <- ' # 四个潜变量,每个由若干题项反映 satisfaction =~ S1 + S2 + S3 + S4 + S5 stress =~ P1 + P2 + P3 + P4 commitment =~ C1 + C2 + C3 + C4 turnover =~ T1 + T2 + T3 ' fit_cfa <- cfa(model_cfa, data = df_sim) summary(fit_cfa, fit.measures = TRUE, standardized = TRUE)跑完之后我会先看两个东西。第一是标准化因子载荷,通常要求每个题项在所属潜变量上的载荷大于0.5,理想情况大于0.7。如果某个题项载荷低到0.3几,我会考虑它是不是该被剔除,或者是不是存在严重的交叉载荷。第二是整体拟合指标,尤其是CFI、TLI、RMSEA、SRMR这四个。
| 拟合指标 | 参考阈值 | 我的实际判断标准 |
|---|---|---|
| CFI | > 0.90 可接受,> 0.95 较好 | 低于0.90基本不看结构模型 |
| TLI | 同CFI | TLI会惩罚复杂模型,比CFI略保守 |
| RMSEA | < 0.08 可接受,< 0.05 较好 | 如果置信区间上限超过0.08,要警惕 |
| SRMR | < 0.08 | 对错误模型很敏感,和小样本更相关 |
这里有个很容易被忽略的细节:拟合指标达标的模型不一定就是好模型。有一次我跑出来的CFA拟合指标非常漂亮,但仔细看修改指数,发现好几个人气极高的题项残差被建议两两相关,而且这些题项在语义上都脱不开干系——说白了就是同一个意思换了种问法。这种局部依赖如果靠修改指数硬加残差相关,拟合指标当然会上去,但模型的概念结构已经变形了。
2.2 从CFA到完整SEM:结构路径怎么写
CFA确认测量模型没问题之后,才轮到结构模型。做法就是在原模型语法后面加上回归路径,把潜变量之间的关系显式写出来:
model_sem <- ' # 测量模型部分 satisfaction =~ S1 + S2 + S3 + S4 + S5 stress =~ P1 + P2 + P3 + P4 commitment =~ C1 + C2 + C3 + C4 turnover =~ T1 + T2 + T3 # 结构模型部分 commitment ~ satisfaction + stress turnover ~ commitment + stress ' fit_sem <- sem(model_sem, data = df_sim) summary(fit_sem, fit.measures = TRUE, standardized = TRUE)在写结构模型之前,先在心里过一遍识别问题。SEM的识别可以简化成三句话:样本协方差矩阵提供的信息要大于等于你要求模型估计的参数个数;每个潜变量至少要有两个或三个指标,否则载荷没法稳定估计;每个潜变量的尺度必须固定,要么固定某个载荷为1,要么固定潜变量方差为1。cfa()和sem()默认会固定每个潜变量的第一个载荷为1,这就是在解决尺度问题。
2.3 拟合指标不能只看P值,还要看这些
很多初学者拿到summary()输出,先看卡方P值,一看P值小于0.05就慌,觉得模型被拒绝了。实际上卡方检验对样本量极其敏感,样本量一大,好模型也照样显著。所以我基本不看它的P值,而是看模型比较时会用到的卡方差异检验——比如做嵌套模型比较时,那个检验才有实际意义。
结构模型跑完之后,比拟合指标更值得关注的是路径系数的方向和量级。我习惯同时保存标准化系数和置信区间,如果某条路径方向跟理论假设相反、效应量又小到几乎可以忽略,哪怕它P值小于0.01,我也会怀疑是不是存在抑制效应或多重共线性。这种情况在加入多个相关潜变量时很常见。
另一个建议是,路径系数的解读要带着"研究设计"一起看。截面数据里做出来的"效应"本质上还是相关关系,lavaan不会替你做因果推断。
3. 复合变量:合成分数和潜变量是两回事
3.1 什么时候该用复合变量
复合变量是SEM里最容易引发争执的概念。潜变量理论上是"背后那个共同原因",反映型指标的共同方差是它留下的痕迹;而复合变量更像是一个"组合分数",多个指标按一定权重合成一个新变量,指标之间的关系是构成型的,也就是指标合在一起定义了这个概念。打个比方,潜变量像温度计里的汞柱——温度是背后原因,汞柱高度是结果;复合变量像体质指数——身高和体重组合到一起才构成这个概念,缺一个都不完整。
我在实际项目里用复合变量,一般就这三种情况:一是指标数量很多而样本量有限,潜变量模型估不动,用复合变量降低参数负担;二是有明确的理论或实务依据,说明几个维度应该等权重或按特定权重相加,而不是让数据自己估计载荷;三是把交互项或二次项构造好之后作为显变量放进模型,这类高阶项用复合变量处理比用潜变量交互简单太多。
3.2 lavaan里定义复合变量的两种写法
如果你用的是较新版本的lavaan,可以直接用<~符号定义构成型复合变量:
model_comp <- ' # 工作压力复合指数:由四个维度构成 stress_comp <~ P1 + P2 + P3 + P4 # 结构路径 turnover =~ T1 + T2 + T3 turnover ~ stress_comp + satisfaction '注意一个关键差异:用=~定义的是反映型潜变量,载荷是数据估计出来的;用<~定义的是构成型复合变量,你默认它是由指标"组合"出来的。如果你的lavaan版本不支持<~,也可以绕行——把所有载荷固定成同一个值,比如全部固定为1,构造一个"等权重复合变量":
model_comp_fix <- ' stress_comp =~ 1*P1 + 1*P2 + 1*P3 + 1*P4 turnover =~ T1 + T2 + T3 turnover ~ stress_comp '但我要提醒一句:这个"绕行方案"里的stress_comp依然会被当成一个潜变量,它会试图解释指标之间的共同方差,而不是完整保留每个指标的独特信息。如果四个指标内部一致性本来就不高,这种写法会得到糟糕的拟合。它跟直接计算rowMeans(df[, c("P1","P2","P3","P4")])得到的均分并不是一回事。
3.3 复合变量的识别风险与规避
复合变量最大的坑是识别不足。一个构成型复合变量本身没有测量尺度,如果它既没有指向任何内生变量,自身又没有外生锚,模型根本就估计不出来。实际建模时,至少要让复合变量指向一个结果变量,或者被某个外生变量影响,这样它才有了可以锚定的结构关系。在我上面的例子里,stress_comp指向turnover,这个锚就足够了。
如果你把复合变量当成控制变量丢进模型,理论依据又不充分,我建议还是用潜变量或者直接用显变量原始题项。复合变量不是万金油,它只在你确实需要"组合"这个操作时才该出现。
4. 分组与嵌套数据:测量不变性和多水平边界
4.1 多组SEM:三步完成测量不变性检验
拿模拟项目X的数据为例,如果我想比较两个部门之间模型路径是否有差异,首先得确认一个问题:这个量表在两组人眼里是不是在测量同一个东西?这就是测量不变性检验。lavaan做这件事非常顺手。
第一步跑形态不变性模型,也就是所有参数完全自由:
fit_configural <- sem(model_sem, data = df_sim, group = "dept")第二步在形态不变性基础上,加上载荷等值约束:
fit_metric <- sem(model_sem, data = df_sim, group = "dept", group.equal = "loadings")第三步再去掉截距的自由估计,强制截距也等值:
fit_scalar <- sem(model_sem, data = df_sim, group = "dept", group.equal = c("loadings", "intercepts"))然后做嵌套模型比较:
lavTestLRT(fit_configural, fit_metric, fit_scalar)如果fit_metric和fit_configural差异不显著,说明两组人对题项的敏感度差不多,弱等值成立;如果fit_scalar也不显著,说明两组的截距也一致,可以放心比较潜均值。注意,很多文章只做到载荷等值就去比较路径系数,这其实还不够严。截距不等值意味着系统性的响应风格差异,比较出来的组间路径差异可能是假的。
4.2 嵌套数据:lavaan的cluster到底做了什么
嵌套数据,典型场景是员工嵌套在门店里,学生嵌套在学校里。同一个门店里的人会共享某些环境因素,观测不再独立,标准误就可能低估,检验就会变得虚高。lavaan对这类数据的支持有一个明确的边界:
fit_nested <- sem(model_sem, data = df_sim, cluster = "store_id") summary(fit_nested, fit.measures = TRUE)cluster参数做的事情是稳健标准误修正,它会把各门店内部的聚类依赖考虑进去,让标准误、P值和置信区间变得更可靠。但它不改变系数估计本身,也不会估计门店层的随机截距或随机斜率。换句话说,lavaan调整的是"推断的稳健性",不是"结构的层级性"。
如果你的数据聚类性很强,比如组内相关系数ICC超过0.2甚至更高,又或者在理论层面你需要回答"门店层面的氛围如何影响员工个人满意度和离职意愿",那么lavaan的这个cluster就不够用了。真正的多水平SEM需要在模型里同时写层1和层2的表达式,这类功能通常要在专门的商业软件或者贝叶斯建模框架里才能跑舒服。
4.3 一条退而求其次的处理路径
如果你手头只有lavaan,又不得不面对明显嵌套的数据,我会给出一个降级方案:先用icc()或类似方法算一下关键变量的组内相关,确认嵌套到底有多严重;如果是小区域分组且组数不多,可以在模型里加入一组组虚拟变量作为控制,这至少能吸收组间均值差异,比完全忽略嵌套要强得多。还有一个更简单的处置是组均值中心化——把每个人的原始值减去所在组的均值,剥离掉组间共同成分,但这个方法会改变变量的解释,最好在方法部分写清楚。
这种做法不算完美,但它在"工具受限"时比假装数据完全独立更诚实。
5. 时间数据:潜增长曲线和交叉滞后模型
5.1 潜增长曲线模型:三次测量能回答什么
纵向数据在lavaan里最经典的用法是潜增长曲线模型。它把个体随时间的变化拆成两个潜变量:一个是截距因子i,代表个体的初始水平;一个是斜率因子s,代表个体的变化速率。编码方式直接用测量的时间点。
假设项目X在三个时间点测量了满意度,变量名分别是S_t1、S_t2、S_t3:
model_lgm <- ' i =~ 1*S_t1 + 1*S_t2 + 1*S_t3 s =~ 0*S_t1 + 1*S_t2 + 2*S_t3 ' fit_lgm <- growth(model_lgm, data = df_sim) summary(fit_lgm, fit.measures = TRUE, standardized = TRUE)i的载荷全部固定为1,意思是不管哪个时间点,初始水平对观测的贡献都一样;s的载荷从0、1到2,这是把三次测量的时间间隔编码成等距的线性趋势。这里有个实操细节:时间编码的零点可以平移,只影响截距的解释,不影响模型拟合。但时间间隔必须和你的设计一致,如果第二次和第三次测量间隔是第一次的两倍,就用0、2、4而不是0、1、2。
潜增长模型里我最关心的是两个潜变量的方差。截距方差显著,说明个体起点确实有差异;斜率方差显著,说明个体的变化速率不一样——这就为下一步寻找"哪些因素解释变化速率差异"提供了动机。
5.2 交叉滞后模型:时序关系的证据边界
纵向数据的另一个高频需求是检验两个变量之间谁先影响谁,交叉滞后模型就是为此设计的。最简单的两波次版本是这样的:
model_clpm <- ' S_t2 ~ S_t1 + T_t1 T_t2 ~ T_t1 + S_t1 S_t1 ~~ T_t1 S_t2 ~~ T_t2 S_t1 ~~ S_t1 T_t1 ~~ T_t1 S_t2 ~~ S_t2 T_t2 ~~ T_t2 '这个模型里,T_t1到S_t2的路径衡量的是"前测离职意向对后测满意度的预测",S_t1到T_t2则反过来。两波次只能看到这样的单步预测,要谈更稳定的交叉滞后关系,至少需要三波次数据。另外,传统CLPM有个毛病:它把个体间差异和个体内变化混在一起估计。更严格的替代是随机截距交叉滞后模型,在普通CLPM基础上加上两个ri潜变量吸收稳定的个体差异,但这对数据和收敛性的要求也更高。
我得强调一句:交叉滞后路径显著,只能说"时序上的先后关系与模型一致",离因果还有距离。漏掉第三个共同原因,或者测量时间间隔设置不合理,都可能导致完全相反的结论。
5.3 纵向数据的缺失值处理
纵向调查必然伴随样本流失,这是跑时间模型时绕不过去的坎。lavaan默认遇到缺失值会做成列删除,这对纵向数据非常浪费,也可能引入严重偏倚。正确做法是在拟合函数里加上缺失处理参数:
fit_lgm <- growth(model_lgm, data = df_sim, missing = "fiml")missing = "fiml"用的是全信息最大似然,它利用所有可用信息来估计参数,前提是数据缺失机制总体上不是完全不可忽略的。只要缺失和"已观测到的变量"有关,FIML通常比粗暴删除要好得多。
6. 非递归模型:反馈回路的识别与收敛实战
6.1 什么时候模型会出现"回路"
前面讲的标准SEM都是递归模型,也就是变量之间只有单向箭头,可以按因果关系排序。但现实里很多问题天生带回路。比如工作压力和工作满意度,常说压力越大满意度越低,但也可能是满意度低的人对压力的感受更强烈;再比如组织承诺和离职意向,承诺低会让人想走,可一旦有离职念头,对组织的承诺也会进一步崩坏。这种双向箭头构成的反馈回路,就是非递归模型。
在lavaan里,写双向路径很简单,难的是让模型可识别、能收敛、有稳定的解。
6.2 阶条件与秩条件:识别靠的是排除限制
非递归模型识别的核心一句话:每个方程都必须有一个"自己的工具变量"。所谓工具变量,在SEM语境下就是那些只在其中一个方程里出现、在另一个方程里被排除掉的外生变量。用模拟项目X的数据举例,假设我要同时估压力和离职意向的互相影响,我写:
model_nonrec <- ' stress =~ P1 + P2 + P3 + P4 turnover =~ T1 + T2 + T3 turnover ~ stress + support stress ~ turnover + workload '这里我用了两个额外变量:support(支持感)放在了离职方程里,workload(工作量)放在了压力方程里。识别条件要求,第一个方程里有一个不出现在第二个方程的外生变量,第二个方程里也有一个不出现在第一个方程的外生变量。这就是所谓的排除限制——你把某些路径显式排除在某个方程之外,才可能同时分辨出两条相反路径的效应。
判断识别有两种层面的检查。先看阶条件,即每个方程被排除的变量数要达到内生变量总数减去1;再看秩条件,即这些被排除的变量不能和该方程的因变量毫无关系,否则就是弱工具。很多时候模型写出来阶条件满足,跑起来却收敛得很勉强,多半是秩条件出了问题——那个"独有变量"其实不解释目标变异的任何部分。
6.3 lavaan跑非递归模型的三个实操细节
第一个细节是收敛失败。反馈回路的模型非常容易不收敛,因为迭代过程中两条路径互相牵制。我的办法是给关键路径一个合理的初始值,比如:
turnover ~ support + 0.3*stress stress ~ workload + (-0.3)*turnover初始值能让优化器少走弯路,但不建议随便改,最好先跑一遍单方向模型,用估计结果作为非递归模型的起点。
第二个细节是结果出现负方差或异常大的标准化系数。这时候先别急着改模型,去检查你的工具变量和因变量之间的相关方向。如果support和turnover几乎不相关,秩条件可能已经在倒塌边缘,这时候任何估计值都不可信。
第三个细节是,双向路径估计出来如果一条显著一条不显著,不要马上把所有不显著路径删掉再重估。正确的做法是把双向模型简化成递归模型,或者把双向路径换成~~相关,做一次模型比较,看拟合有没有显著恶化,再下结论。
7. 空间数据在lavaan里的边界:能做什么,不该硬做什么
7.1 把空间当成"分组嵌套"来处理
空间数据在SEM里是个比较尴尬的话题,因为lavaan本身并没有空间自回归参数。如果你手头的数据来自不同城市、不同区域或不同门店,最方便的做法是把空间当成一种分组嵌套。和第四章的嵌套逻辑一样,用cluster = "region_id"去调整标准误,避免同一区域内的样本依赖造成推断偏误。
这种方案适用的场景是:你关心的是个体层面的关系,空间只是干扰因素。比如员工满意度模型里,同一家门店的员工会受同样的店长和排班制度影响,你要的是把这种"店效应"从推断误差里剔除掉。这时候cluster修正完全够用,也足够透明。
7.2 空间滞后变量:手动构造后纳入模型
另一种处理路径是把"空间影响"显式构造为变量,再放进SEM里。典型做法是先根据经纬度或区域邻接关系计算空间权重矩阵,然后对某个关键变量构造空间滞后值,也就是"邻居平均值",最后把这个滞后变量作为解释变量加入模型。
这部分严格来说不再是lavaan的空间功能,而是你在建模之前先完成了空间数据处理。用这种方式,模型可以回答"邻居区域的某项指标是否会影响本区域的个体结果",但完全无法处理"邻近区域之间的相互反馈"这类内生空间效应,因为那已经超出了SEM的范围,需要空间计量模型来处理。
7.3 真正的空间SEM该用别的工具
我必须把边界划清楚:lavaan不是空间计量工具。如果你需要估计空间自回归系数,也就是让模型显式地表达出"区域之间的依赖强度",lavaan做不了。这种情况应该转向专门的空间计量包或者贝叶斯建模框架。空间SEM是一个比自己动手构造滞后变量复杂得多的体系,涉及空间权重矩阵的设定、空间误差和空间滞后的选择、以及空间溢出效应的解释,硬塞进lavaan里只会得到一堆无法自圆其说的数字。
反过来,如果你只是想在SEM里照顾到空间结构,避免标准误失真,cluster配合手动构造的空间变量,已经是lavaan能在不引入其他工具的前提下做到的最合理方案。至于研究方法部分怎么写,我会诚实地写"标准误按区域聚类进行调整"和"模型纳入空间滞后控制变量",而不是写"做了空间结构方程模型"。
我自己跑这些模型最大的体会是:lavaan强在"模型即文本",所以千万不要一上来就追求跑通完整模型。我会先写测量模型,跑CFA,保存结果,再逐步加结构路径,再加组间约束,每一步都留下一个R脚本文件。另一个很实用的经验是,凡是遇到非递归模型或复合变量,先拿张白纸把模型画一遍,数清楚每个方程有多少个外生变量、哪些路径被排除在外,确认识别之后再写语法,这样能省掉大量调试时间。
最后分享一个小技巧:建模之前先用inspect(fit_lgm, "sampstat")看一眼样本协方差和相关矩阵,确认变量的方向和量级跟你的理论预期一致。很多时候模型跑不出来,不是语法的错,而是数据里变量的编码方向反了,或者你想当然以为高分的含义却恰好相反。这一步十秒钟就能做完,却能帮你少走最远的弯路。