用stata写多层线性模型这事说难真不难说简单也藏着不少坑。我最早接触HLM多层线性模型也叫混合模型、随机效应模型的时候被各种命令绕得晕头转向上网搜了一圈要么是理论推导堆满屏要么是只给一段命令不解释为什么换了个数据就不会用了。后来跑了几十组实际项目数据踩了不少坑才算把stata里做HLM这条路的语句逻辑理顺了。这篇博文就是把我在stata里跑多层线性模型的完整语句套路、参数选择的逻辑、以及那些文档里不常写的实操教训一次性讲清楚。不管你是教育、心理、管理还是医学背景的数据分析者只要数据结构里有学生嵌套在班级里员工嵌套在公司里患者嵌套在医院里这种层级关系下文这套操作都能直接用。1. 建模思路与命令选型先想清楚再敲键盘1.1 什么时候必须用HLM平时用普通回归regress分析的时候默认所有观测是独立的。但实际数据里它根本不独立——同一个班的学生共享一个班主任同一个医院的患者面对同样的医疗资源这种组内相关性如果不管标准误就偏低p值容易被低估本来不显著的变量可能就假显著了。多层线性模型的核心价值就是把这部分组间差异拆出来同时刻画个体层和组层变量各自的影响。判断一个项目要不要用HLM我一般看两个硬指标。第一数据是不是有明确的层级结构至少两层个体层嵌套在组层里第二组内相关系数ICC有没有意义。ICC说的是总变异里有多少比例来自组间差异。用stata跑一个空模型看estat icc的结果如果ICC只有0.01那基本可以当普通回归处理如果ICC在0.1以上尤其还在0.2、0.3这种水平说明组间差异很大必须做多水平建模。我在教育类数据里见过ICC高达0.4的情况班级层面的效应大到不能忽视这时候用普通OLS结论基本是偏的。1.2 stata里做主流的几条技术路线stata里做HLM市面上能看到的方案大概四类mixed、gllamm、menl以及runmlwin。这几条路线各有各的使用场景但主流推荐还是mixed命令。mixed新版叫mixed老版本叫xtmixed是stata官方命令跑线性多层模型的首选。它支持连续结局变量固定效应、随机截距、随机斜率都能处理估计方法有最大似然mle和限制性最大似然reml输出结果规范配合后验预测、estat系列命令基本能满足绝大多数实证需求。gllamm是第三方命令需要安装功能很强大可以处理非连续结局比如二分类、计数、有序分类变量的多水平模型模型灵活度高但语法比较晦涩运行速度也不算快新手不友好。menl是非线性混合模型的官方命令处理的是非线性函数形式的模型比如组间参数本身就是非线性关系的场景用的情况相对少。runmlwin是在stata里调MLwiN的接口适合跑多水平多变量模型、复杂交叉分类模型前提是电脑装了MLwiN它对跨软件操作不太熟悉的人来说有点门槛。我的建议很直接常规的连续结局变量直接用mixed其他几条路等你碰上了再说不要在选工具上耗太多时间。命令只是一个入口理解和建模思路才是重点。1.3 为什么我推荐从mixed命令入门mixed命令的语法骨架长这样mixed 因变量 固定效应变量 || 层级变量: 随机效应变量, 估计方法选项这个语法的设计逻辑很符合写模型的直觉。竖线前面的部分和普通回归的写法一样指定因变量和自变量竖线后面的部分专门描述层级嵌套和随机效应。它把固定部分和随机部分分开层次感很清楚。而且mixed命令是基于mata写的高效算法处理几千个组、几万条记录的常见数据完全没问题。我处理过一组近10万条患者数据、嵌套在500多家医院的例子mixed跑起来也就几十秒。另一个重要原因是mixed命令配套的后验分析功能很完善。estat icc能直接给组内相关系数estat lrtest可以做似然比检验estat recovr能显示随机效应的方差协方差矩阵这些在做模型诊断和结果汇报时都是必用的。换一条技术路线这些配套功能不一定齐还得自己去算很麻烦。2. 完整语句流程从数据准备到结果输出2.1 数据结构检查与变量中心化处理在写任何mixed语句之前数据整理这步不可跳过。我见过很多人在模型结果异常时回去查数据结构发现变量格式不对白白浪费半天。要做的事无非三件。第一确保层级变量存在并且确认它的类型是数值型。ID变量必须是干净的数值。层级变量可以用一个简单tab命令检查分组信息tab school_id看分组数量和每组样本量。如果group数量太少小于20的情况就要小心后面随机效应方差的估计会不稳定。第二检查变量缺省情况。多层模型对缺失值的处理比较敏感必要时先看一下每个变量的缺失比例misstable summarize如果缺失比例高的变量又很重要建议考虑多重插补不要直接listwise删除否则样本量骤减随机效应方差估计会很不稳。第三中心化。这是多水平模型里特别关键的一步。如果自变量是连续变量且包含有意义的0点比如年龄、收入可不中心化但如果0点没有实际含义比如SES得分最好做组均值中心化或总均值中心化。组均值中心化能把个体层的效应解释为组内相对水平的效应在跨层交互模型里尤其常用。用stata实现很简单* 总均值中心化 center ses, prefix(c_) * 组均值中心化 bysort school_id: egen mean_ses mean(ses) gen gm_ses ses - mean_ses我自己的经验是跑交互项的时候连续变量一定提前中心化否则交互项和主效应之间的共线性会让你怀疑人生。2.2 空模型零模型语句与组内相关系数计算空模型是HLM的起点也叫截距模型里面只有一个因变量和一个层级变量不放任何解释变量。它的作用是把总变异拆成组间变异和个体变异两块。以教育数据为例学生math成绩嵌套在学校school_id里空模型语句如下mixed math || school_id:, mle跑完看两个地方。一个是随机效应的方差结果——school_id的方差是组间方差var(_cons)残差的方差是组内方差。另一个是直接用estat icc得到组内相关系数estat iccICC算式就是组间方差除以组间方差组内方差。假设组间方差是2.0组内方差是8.0ICC就是2.0 / (2.0 8.0) 0.2意思是math成绩的变异有20%来自学校之间的差异。我拿到ICC后怎么决定下一步如果是小样本探索性研究ICC超过0.05就可以考虑多水平模型如果是严谨的实证项目一般要到0.1以上才说明组层机制不可忽略。空模型还承担另一个职责给后续模型提供基准对数似然值。后面每加一组变量谁优谁劣靠的就是对数似然值和信息准则的比较。2.3 逐步加入固定效应空模型确认了层级结构存在接下来就是把个体层变量放进去。一般顺序是先加个体层变量确认哪些在组内有显著效应再加组层变量再考虑随机斜率和跨层交互。我拿学生数学成绩举例。第一层放学生个体变量比如性别female0/1、家庭SES得分sesmixed math female c_ses || school_id:, mle注意我用的c_ses是中心化后的变量。female作为二分类变量可以直接放如果是多分类变量建议用ib3.grade这种因子变量写法明确指定参照组。模型跑完后看固定效应系数的解读方式female的系数表示在同一个学校内部女生相对男生的平均成绩差异c_ses的系数表示在同一所学校内SES每高一个单位学生成绩预期变化多少。随机截距部分输出的var(_cons)代表控制了这些个体变量后学校之间仍然存在的平均成绩差异。然后加学校层变量比如学校资源指数fundingmixed math female c_ses funding || school_id:, mle这个时候funding的系数解释是控制学生个体特征后学校资源每增加一个单位学校平均成绩水平的变化。这就是跨层固定效应的直接体现。再加随机斜率的时候要小心。随机斜率说的是某个个体层变量的效应在不同组之间是否不同。比如我怀疑家庭SES对成绩的影响在不同学校不一样就可以让c_ses的斜率随机化mixed math female c_ses funding || school_id: c_ses, mle在竖线后的school_id:后面写上c_ses就表示c_ses的系数允许在不同学校之间变化。跑出来会多一个随机效应项比如var(c_ses)和var(_cons)之间的协方差。这里马上会遇到一个常见问题——随机斜率模型经常不收敛。我在后面第4节会专门讲排查方法。2.4 随机斜率与跨层交互上面这步已经把随机斜率加进来了但真正的跨层交互还没做。跨层交互的目的是检验某个个体层变量的效应是否取决于某个组层变量比如家庭SES对成绩的影响是否因学校资源水平而变。在stata里生成交互的方式是用井字号运算符mixed math female c_ses funding c_ses#c.funding || school_id: c_ses, mle这里c_ses#c.funding表示连续变量SES和连续变量funding的交互项。跑出来之后交互项系数是核心关注的对象。如果为正说明家庭SES对成绩的正面影响在资源更好的学校里更强烈如果为负则说明资源好的学校反而在缩小SES带来的成绩差距。解释交互项时必须注意这时的主效应系数含义已经变了。c_ses的系数代表funding等于0如果中心化了就是funding平均水平时SES的效应funding的系数代表c_ses等于0时资源的影响。所以连续变量中心化在这种模型里几乎必不可少否则主效应系数会无法解释。随机斜率和交互项同时放进模型模型方程会变得复杂。除了看固定效应系数我还习惯把随机效应的方差协方差矩阵导出来看一眼estat recovr这个命令会展示各个随机效应之间的方差和协方差。如果随机斜率方差很小或者协方差部分出现了奇怪的极端值就要考虑这个随机斜率是不是有必要保留。用似然比检验可以给这个判断提供统计依据把带随机斜率的模型和不带随机斜率的模型用一个命令对比estat lrtest这个检验的p值如果大于0.05说明加随机斜率没有显著改善模型那就不如把随机斜率去掉保持模型简洁。2.5 模型比较与结果输出常用语句模型比较这件事在stata里最常用的是信息准则和似然比检验。信息准则在mixed运行结果中会直接给出AIC和BIC。模型的AIC/BIC越小越好注意BIC对参数个数的惩罚更重所以两个准则结论不一致时我会多看BIC防止过度拟合。对于嵌套模型一个模型是另一个模型的特殊情形用似然比检验最规范。做法是先保存一个模型estimates store model1再跑一个更复杂的模型mixed math female c_ses funding c_ses#c.funding || school_id: c_ses, mle estimates store model2 lrtest model1 model2lrtest输出一个p值如果p值小于0.05说明模型2显著优于模型1新加的效应或随机参数有保留价值。这里必须强调一个前提lrtest只对mle估计方法有效如果你用了reml两个模型的固定效应部分不同时不能直接用这个检验比较。结果输出方面我强烈建议用esttab把多个模型放到一张表里方便写论文或做汇报esttab model1 model2 using hlm_results.rtf, /// b(3) se(3) star(* 0.05 ** 0.01 *** 0.001) /// title(HLM Model Comparison) replace这样能一次性输出系数、标准误和显著性标记比自己一个个抄到Word里省事得多。变量标签如果没设定建议先把标签设置好label variable c_ses SES (centered)3. 语句背后的统计细节为什么这样写坑在哪3.1 随机效应方差结构选择mixed命令里有一个covariance()选项默认是identity也就是unstructured意思是允许随机截距和随机斜率之间存在协方差。实际操作里到底用哪种方差结构学问不少。如果模型只有随机截距不需要额外设置var(_cons)就是组间方差。如果同时有随机截距和随机斜率默认的unstructured会估计var(_cons)、var(斜率)和它们之间的协方差这个协方差是有实际含义的——它告诉我们组的平均水平和个体层变量的斜率之间是否存在关联。比如截距高、斜率低的学校协方差会是负值说明高平均水平的学校个体差异反而更小。这个信息对教育政策的解读很有价值。但问题是当随机效应比较多时unstructured参数数量增长很快模型容易不收敛。我处理这种情况的办法是先跑unstructured如果迭代十几轮都收敛不了就换成independentcovariance(independent)强制截距与斜率不相关减少协方差参数再不收敛就干脆把随机斜率项去掉只保留随机截距。不要在一棵树上吊死建模的目的是解释数据不是追求复杂。3.2 极大似然与限制性极大似然如何选mle和reml的区别很多人记不住。一句话说mle估计的是固定效应和方差成分的时候把固定效应当作已知的reml在估计方差成分的时候先调整掉固定效应占用的自由度对方差成分的估计更无偏。所以经验法则是需要比较固定效应嵌套模型时用mle更关心方差成分的准确估计时用reml。但这里有个关键限制reml下不能比较固定效应部分不同的模型因为reml的似然值里包含了对固定效应的调整两组不同固定效应的模型之间不再有可比性。因此我的常规操作是模型探索阶段用mle方便各种比较模型定下来了用reml重新估计一遍作为最终结果。样本量大的时候mle和reml的估计结果差异通常很小但在二三十个组的小样本场景下差异能大到影响结论必须认真对待。stata命令切换很简单mixed math female c_ses funding || school_id:, reml3.3 自由度、p值与显著性mixed命令输出的固定效应表里有z值和p值。这个z检验是大样本近似当组数和每组样本量不够大的时候p值可能偏小。很多人不知道的是stata里可以通过dof选项来调整小样本自由度。比如用以下语句做组数减一或者Satterthwaite调整mixed math female c_ses funding || school_id:, mle dfmethod(satterthwaite)还有一个选项是dfmethod(kroger)。我在做医学项目时经常会看到审稿人要求报告自由度调整所以如果你的组数偏少比如只有30个医院建议用dfmethod选项而不是直接默认z检验。这个细节在初级教程里很少被提到但现实中很重要。3.4 迭代收敛与参数估计没有收敛的模型结果不能直接采用这是基本的职业底线。mixed默认最大迭代次数是1600如果迭代次数超了说明模型很可能设定有问题。常见的原因包括随机效应的初始值太差、变量量纲差异过大、模型过于复杂超过数据能支撑的参数数量。针对量纲问题最直接的办法是把连续变量标准化或者做尺度变换。我在stata里常用egen z_ses std(ses)用标准化后的变量重跑模型收敛概率会大幅提高。还有一种情况是数值精度到了但被判定为没收敛可以尝试指定迭代技巧mixed ..., technique(bhhh)或者放宽收敛容限mixed ..., tolerance(1e-7)如果这些手段都用上了还是不收敛那大概率是模型设定与数据结构不匹配老老实实简化模型才是正解。4. 常见问题与排查技巧实录4.1 模型不收敛的几种典型处理法我在实操中遇到的不收敛大致分三类。第一类是随机斜率模型不收敛解决办法是先将斜率变量的方差初始值调整或者改用independent协方差结构或者去掉随机斜率。第二类是复杂模型在某个变量上出现完美共线这时候stata往往会提示perfect prediction查一下因子变量设置、分类变量的参照水平编码是否对。第三类是数据里层级ID编码不连续导致矩阵运算异常把ID变量重新编码egen new_id group(school_id)运行新ID变量代替旧ID。这类问题容易忽略只要层级ID是字符串或者有缺失后面数据合并时会出现各种诡异结果。4.2 固定效应与聚类稳健标准误另一个容易混淆的问题是HLM里要不要用聚类稳健标准误。答案是混合模型本质上已经通过随机效应处理了组内相关性所以一般不需要额外标注vce(cluster)。但如果你担心随机效应模型对异方差敏感可以在mixed里加mixed math female c_ses || school_id:, mle vce(robust)沙盒试验里它通常能给出更保守的标准误。对比一下随机效应标准差稳健后的结果如果差异很大说明模型可能没有正确捕捉方差结构这时候值得回头审查随机效应的设定而不是简单选用稳健标准误了事。我自己做报告的时候会把普通标准误和稳健标准误的结果都放上去让读者判断稳健性。4.3 三分层与嵌套结构处理现实场景经常不只有两层。比如学生嵌套在班级里班级又嵌套在学校里。stata的mixed命令支持多层嵌套竖线表示层级嵌套关系mixed math female c_ses || class_id: || school_id:, mle这里有两个竖线表示class_id嵌套在school_id之内。输出结果会多出class_id层和school_id层的方差成分。各层的方差分别代表班级层面和学校层面的组间变异。需要提醒的是层数越多数据量和组数量要求越高。三层模型至少要保证每层有足够的组数否则最高层级的方差会估计不准甚至出现负方差。4.4 输出表格合并与论文排版跑模型不是终点把模型报告清楚才是。上面提到用esttab输出结果我再补一个实用技巧估计完模型后用以下命令把模型之间的关系存好方便随时调出来estimates table model1 model2, star stats(N ll aic bic)这句命令会在stata结果窗口里横向对比多个模型的系数、星星和拟合指标用来快速判断哪个模型更合适。除此之外写论文时经常要报告ICC和随机效应方差可以手动用estat icc里的数值或者结合返回值return list这里有回归后保存的全部标量比如方差值、对数似然值等可以进一步自动化汇总到表格里。我写过一个小循环批量处理多个因变量并直接生成Word表格效率很高。如果你经常做这类多水平分析建议花点时间研究esttab的选项这比手动复制粘贴省出大量时间。4.5 常见问题速查表问题现象可能原因解决措施模型不收敛随机效应结构太复杂简化模型、用independent协方差、去掉随机斜率变量被omitted共线性或分类变量参照组设置不当检查因子变量语法考虑中心化或删去高相关变量ICC几乎为0组间差异太小可考虑放弃HLM用普通回归即可p值过于显著组数少导致z检验失真使用dfmethod调整自由度随机效应方差为负模型过度拟合或数据与模型不匹配精简随机部分检查数据结构esttab表格变量名不友好变量标签未设置提前用label variable设置中文或好读的标签reml下lrtest提示无效固定效应不同的模型不能比较改用mle估计后再做似然比检验分组ID为字符串矩阵运算无法执行用egen new_id group(id)转换数值型这张表是我自己的排查清单建议收藏起来碰到异常先对照一遍。4.6 我踩过的典型坑一个完整排查案例最后分享一个我最近真实处理的案例。某管理类数据600多个员工嵌套在45个部门研究员工绩效和工作满意度的关系。我一开始直接上了随机斜率模型mixed perf satis || dept: satis, mle结果迭代2000次都不收敛。我当时第一反应是变量量纲问题——绩效是0到100分满意度是1到5分。把满意度标准化之后再跑egen z_satis std(satis) mixed perf z_satis || dept: z_satis, mle这个问题就消失了。后来我又发现部门数量45个其实有点少随机斜率方差估计出来较大但不显著结合estat lrtest发现随机斜率项没有显著改善模型最终模型去掉了随机斜率只保留随机截距mixed perf z_satis || dept:, mle这个案例说明很多模型问题不是出在数据质量而是模型复杂度和数据结构不匹配学会做减法也是建模的一部分。另外部门数量偏少时我还建议用bootstrap或鲁棒标准误做敏感性分析看结论是否稳健这在管理类实证中越来越常见。说点实在的stata里的HLM语句本身并不长难得是理解每条语句在模型里做了什么、随机效应怎么设置、模型估计结果是否可信。我建议你拿到自己的数据后先老老实实跑空模型算ICC再逐步加变量控制固定效应最后再决定要不要随机斜率、要不要交互。每一步都用lrtest或信息准则检验一下比一次性堆一个超复杂模型靠谱得多。如果中间发现结果不稳定不要急着换工具先检查和简化模型设定往往问题就出在最不起眼的那一步。