后台经常有人问我Stata到底能不能做HLM当然能而且从命令成熟度和输出友好度来看Stata可以说是做多层线性模型最顺手的工具之一。HLM多层线性模型在Stata中对应的一套语句核心就是mixed命令。这篇文章我就把从数据准备、空模型、随机截距、随机斜率到结果解读的完整语句流程给你捋一遍顺便把我踩过的坑也一并交代清楚。这套内容适合谁正在处理嵌套结构数据的同学比如教育研究里学生嵌在班级和学校、组织研究里员工嵌在企业、医疗研究里病人嵌在医院或者追踪研究里重复测量嵌在个体身上。你只要发现自己的数据天然带“层”的概念并且对Stata的基本操作有一定熟悉度这篇文章可以直接当你的操作手册。1. 多层线性模型到底在解决什么问题1.1 不是所有数据都是“独立”的很多人在入门统计时学的经典回归默认一个前提每一个观测值之间相互独立误差项互不相关。这个前提在嵌套数据结构下基本不成立。我给你举个最直观的例子。研究学生成绩的影响因素你从10所学校抽了500个学生。这500个学生并不是真正独立的个体——同一个学校的学生共享同一个校园环境、同样的教师风格、相似的生源背景。他们的成绩天然存在“抱团”现象来自优质学校的学生整体偏高来自薄弱学校的学生整体偏低。如果直接用普通最小二乘回归等于硬生生把这种组内相似性忽略掉。同样的逻辑出现在很多场景。追踪研究里同一个人在不同时间点的测量存在自相关组织研究里同一个团队的员工彼此影响临床研究里同一个医院的患者接受相同的诊疗流程。这些数据的共同点就是存在“个体”和“层”两个层面的信息。1.2 普通回归在多层级数据上会出什么乱子用普通回归硬解嵌套数据最直接的问题有两个。第一个是标准误被低估。因为你把500个学生的数据当成500个独立样本去算标准误但真实的有效信息量可能没有那么大。同一个学校的学生之间有相关性等价于你的样本打了折扣。标准误偏小p值就偏小本来不显著的结果可能变得显著这就是I类错误膨胀。第二个问题是生态谬误。普通回归把组间效应和组内效应混在一起估计。假设学生家庭收入与成绩的关系在个体层面是促进的但学校层面可能因为“名校聚集效应”呈现负相关。普通回归给出的系数是这两个效应的混合物解释起来非常尴尬。多层线性模型的核心逻辑就是显式承认数据存在层级结构把方差拆成“层内”和“层间”两部分。量化这个问题有一个指标叫ICC组内相关系数计算公式是组间方差除以总方差。ICC越接近0说明组间差异越小你用普通回归问题不大ICC大于0.05甚至0.1就必须认真考虑层级结构了。注意ICC的计算是HLM入门第一道坎很多人跑完mixed命令直接跳过estat icc这一步数据到底适不适合做多层模型都没验证后面的结论就很难站稳脚跟。2. Stata做HLM的语句体系从xtmixed到mixed2.1 命令为什么有两个版本Stata里做多层线性模型最早的主命令是xtmixed很多教材和网上资料还在用它。但Stata 13之后官方就推荐改用mixed命令了。两者核心算法是一样的只是mixed的后续命令更完善比如margins、pwcompare这些后估计工具都能和mixed顺畅配合xtmixed则相对受限。我的习惯是直接用mixed。如果老代码里看到xtmixed直接换成mixed也能跑输出结果基本一致但能享受更完整的后估计功能。2.2 mixed命令的30秒语法速记先记住最核心的模板mixed 因变量 自变量 || 组变量:, 选项竖线后面跟的是层级结构。嵌套多层的写法是连续用空格分隔多个括号mixed y x || school: || class:这表示个体嵌套在班级、班级嵌套在学校三水平的模型。每一层括号后面可以指定随机系数。如果只写组变量名表示这一层只有随机截距。如果想允许某个自变量在不同组间有随机斜率这样写mixed y x1 x2 || school: x1意思是截距和x1的斜率都随学校变化。默认情况下随机截距和随机斜率之间是允许相关的对应方差协方差矩阵的不限定结构。估计方法上mixed默认用REML限制最大似然估计老版本资料里说的“reml选项”现在已经是默认行为。想改用最大似然估计就加mle选项。3. 手把手实操从空模型到完整模型的语句开展3.1 第一步不是跑模型而是先处理数据结构和变量很多人拿到数据直接敲mixed命令结果跑完发现组变量里的某些组只有一两个观测随机效应方差瞎估计再回头看数据清理浪费时间。我建议建模前先把数据形态捋清楚。先检查每组样本量bysort school: gen n_per_group _N summarize n_per_group这步能快速发现哪些组的样本量过少。一般来说一组只有一两个观测的组对随机效应的贡献很大但非常不稳定需要在建模前评估是否合并或处理。接下来是关键的一步中心化。这是HLM实操里新手最容易犯糊涂的地方。level-1变量个体层建议做组均值中心化group-mean centering也就是用变量值减去该个体所在组的均值。这样做的好处有两个第一组均值中心化后的变量只保留组内变异估计出的系数是纯粹的组内效应第二把中心化后的个体变量和组均值作为两个变量放进去就能同时估计组内效应和组间效应把混杂问题拆开。Stata写法bysort school: egen x1_group_mean mean(x1) gen x1_c x1 - x1_group_meanlevel-2变量学校层建议做总均值中心化grand-mean centering因为它本身只有组间变异不用再按组分。直接减去全样本均值egen w1_grand_mean mean(w1) gen w1_c w1 - w1_grand_mean提示为什么要中心化核心目的是让截距有实际解释意义。如果学生家庭收入x1不中心化截距代表x10时的预测成绩但收入0根本不现实。中心化之后截距变成“组均值处”的预测成绩解释起来顺得多。3.2 空模型一切HLM的起点空模型也叫零模型模型里不加入任何自变量只分解方差。这是判断“要不要用HLM”的关键步骤。mixed y || school: estat icc记一下输出结果var(_cons)学校层面的方差分量也就是组间方差var(Residual)个体层面的残差方差也就是组内方差estat icc直接给出ICC和标准误如果ICC很小比如小于0.01说明组间几乎没差异用普通回归就行。如果ICC介于0.05到0.2之间这是很典型的嵌套数据结构继续做HLM完全合理。ICC超过0.2说明组间差异非常大研究组间影响因素就很有价值。空模型的另一个用途是为后续模型的方差分解提供基准。后续每加入一批变量看组间方差缩小了多少就说明这批变量“解释”了多大比例的组间差异。3.3 随机截距模型把个体层变量放进来空模型跑通后开始加level-1变量。mixed y x1_c x2_c || school: , reml这里我用的是x1_c和x2_c也就是已经做过组均值中心化的变量。模型允许不同学校的基线水平不同但x1和x2对y的效应在所有学校是一样的——这就是随机截距模型也叫随机效应ANCOVA模型。随机截距模型跑出来的结果里固定效应部分看x1_c和x2_c的系数。因为已经做了组均值中心化这些系数可以直接解释为“个体层面每增加一个单位y在组内平均变化多少”。随机效应部分关注两点第一加入x变量之后组间方差是否明显缩小第二剩余组间方差是否仍然显著。如果加入level-1变量后组间方差微弱缩小说明这个学校的差异主要由个体构成不同来解释而不是学校本身的机制差异。3.4 完整模型把学校层变量放进来level-1变量整理完之后再加入level-2变量比如学校资源投入来解释组间差异。mixed y x1_c x2_c w1_c || school: , reml这里需要注意w1_c是学校层的变量放进去是解释随机截距的变异来源相当于用学校层面的特征去预测学校的基线水平。它的系数解释是学校层面的w1每增加一个单位该校学生的平均预测成绩变化多少。完整模型跑完后可以对比一下随机效应的变化。原来的空模型中学校方差是某个值加入个体层和学校层变量后学校方差如果大幅下降说明这些变量确实在解释学校之间的差异。有一个便捷操作是先把模型存起来再做似然比检验estimates store m0 mixed y x1_c x2_c w1_c || school: , reml estimates store m1 lrtest m0 m1如果p值显著说明加入的变量整体上显著改善了模型拟合。不过要注意lrtest要求两个模型在同一个样本上估计所以数据里有缺失值的话务必保证两个模型用完全相同的样本。3.5 随机斜率模型允许效应因组而异随机截距模型默认“x对y的效应在各组一样”。这个假设有时候站不住脚。比如研究教学方式x对学生成绩y的影响有些学校可能因为执行力度好教学方式的效果更强有些学校则效果平平。这时候x的斜率就不该是固定值。Stata写法mixed y x1_c x2_c w1_c || school: x1_c, covariance(unstructured) remlcovariance(unstructured)是默认选项表示随机截距和随机斜率之间可以存在相关性。也就是说基线水平高的学校x1的效果可能更强或更弱。随机斜率模型跑完后一定要看两个东西第一随机斜率方差是否显著。可以用estat recovariance查看随机效应的方差协方差矩阵也可以和随机截距模型做似然比检验estimates store m_randint mixed y x1_c x2_c w1_c || school: x1_c, covariance(unstructured) reml estimates store m_randslope lrtest m_randint m_randslope似然比检验显著说明x1的斜率确实在各学校间有差异保留随机斜率是合理的。第二注意模型是否收敛。加了随机斜率之后待估参数变多模型更容易出现收敛问题。如果提示不收敛后面在常见问题里我会详细讲调试思路。实操心得随机斜率模型不是越多越好。我见过有人一口气把五六个自变量全部设成随机斜率结果模型不收敛收敛了运行时间也感人。一般只对自己核心关注的变量设随机斜率其他控制变量固定即可。4. Stata输出怎么读别只看星号4.1 固定效应部分mixed命令的固定效应输出和普通回归长得差不多有Coef.、Std. Err.、z、P|z|、95%置信区间。注意这里是z值而不是t值因为REML框架下用的是近似正态分布。当组数很少的时候这个近似可能不够好有研究者会建议用更复杂的小样本校正但Stata自带的mixed输出不自动做这个校正。固定效应系数的解读跟普通回归类似但加上“控制了组间差异”这个前提。比如x1_c的系数是2.3含义是在同一个学校内部x1_c每增加一个单位y平均增加2.3个单位。4.2 随机效应部分与ICC随机效应部分有var(_cons)和var(Residual)再往下还可能有随机斜率的方差和协方差。这一块的收入经常被新手忽略但随机效应的诊断价值非常高。var(_cons)是组间方差var(Residual)是组内方差。计算ICC的公式就是ICC var(_cons) / (var(_cons) var(Residual))Stata的estat icc会直接帮你算出来同时报告标准误。在完整模型里看ICC的变化也很有意义如果完整模型的ICC比空模型小了很多说明加入的变量解释了相当比例的组间差异。4.3 模型比较的指标HLM报告里经常会看到三种指标似然比检验、AIC、BIC。似然比检验适用于嵌套模型比较用lrtest命令即可。要比较固定效应结构不同但随机效应结构相同的模型时需要用mle重新估计因为REML的似然值不能直接用于比较固定效应不同的嵌套模型。记住这个规则随机效应结构不同 → 保持相同的固定效应用reml固定效应结构不同 → 保持相同的随机效应用mleAIC和BIC的规则是数值越小越好但两者惩罚项不同BIC对模型复杂度惩罚更重。实际选型时这两个指标要结合学术领域的惯例用不能盲目追求最低。5. 常见问题与排查技巧实录5.1 模型不收敛怎么办这是HLM实操里最普遍的坑。模型不收敛的表现是Stata提示not concave或convergence not achieved大概率是下面几个原因。第一个原因是迭代次数不够。解决办法是调大迭代次数mixed y x1_c || school: x1_c, iterate(5000)第二个原因是随机效应方差接近0或者模型过于复杂待估参数过多导致无法收敛。这种时候反过来做减法比如去掉随机斜率、去掉相关性结构改成covariance(independent)或者直接用方差分量结构都能缓解收敛压力。第三个原因是变量量纲问题。x变量如果数值特别大比如收入几万几十万会让优化过程困难。可以把变量缩放一下再放入模型比如除以1000输出结果的解释再相应调整。第四个是用快速可靠的初始值。Stata的mixed命令有startvalues()选项可以指定初始值不过我实测下来多数情况调迭代次数和简化模型就够用了。5.2 随机效应方差估计为0有时候跑完mixed发现var(_cons)几乎是0或恰好是0LL值也有点怪。这说明数据里组间变异极小或者模型过度抽取了组间变异随机部分已经无方差可分。类比理解就好像你给每个人发一套校服结果发现所有人身材几乎一样那“身高这个组间因素”的方差自然趋近于0。这时候要做的是重新审视变量选择和数据分组而不是硬塞随机效应。组数太少也是方差估计为0的一个原因。一般的经验法则是第二层组数至少要有20到30个低于10个时随机效应的方差估计非常不可靠。如果确实现实条件受限可以考虑直接把组作为固定效应放进模型或者用带稳健标准误的普通回归作为替代方案。注意组数太少时REML估计通常比ML更稳定一些但这并不意味着REML能创造奇迹。真实组数太少时随机效应求稳不如把组别当作固定效应做出结果再和HLM结果互相印证。5.3 缺失值处理mixed命令对缺失值的默认处理方式是无缺失样本分析也就是只要模型中任何一个变量缺失这个样本就会被剔除。如果数据里的缺失值比较多这样做会损失不少样本。我建议建模前自己先做缺失值诊断misstable summarize y x1 x2 w1如果缺失模式有规律比如某个学校的变量全缺失直接在模型里剔除这个学校会更干净。如果缺失比例不高无缺失样本分析问题不大。如果缺失较多就慎重考虑要不要用多重插补插补后的数据直接放进mixed也是支持的。5.4 中心化到底该怎么选中心化是HLM里面最容易混乱的操作我单独拿出来说一次。level-1变量的中心化有两种主流选择组均值中心化和总均值中心化。组均值中心化是把个体观测减去所在组的均值优点是能干净分离组内效应和组间效应且与随机斜率模型配合良好。缺点是如果你还想看组间效应就必须把组均值作为level-2变量放回模型。总均值中心化是减去全样本均值计算简单、解释方便但它不改变原理上组内和组间效应的混淆因此对研究问题需要更加谨慎。实际操作我的建议是如果做随机斜率和交叉层级交互优先组均值中心化如果只是想减少共线性、让结果更容易解释总均值中心化足够省事。5.5 随机斜率变化太大但方差又不显著有一种情况随机斜率的方差看起来数值不小但检验不显著。原因可能是随机斜率估计方差的标准误偏大或者样本量不够。这时候不要硬下结论说斜率不随机再看一下学校数量。如果学校数量不到20随机斜率的统计检验功效很低。有条件就继续收集学校数据没条件就如实报告“样本不足以检验斜率随机性”。6. 最后分享一点我的实操习惯做HLM这几年我自己的一套工作流已经固定下来先算ICC确认数据确实存在层级结构再跑空模型中间任何一步都不急着跳级每加入一组变量就保存一次模型估计最后统一做似然比检验汇报结果时固定效应、随机效应方差、ICC、AIC/BIC四个表格一个都不少这样无论审稿人还是领导追问模型细节都有据可查。另外一个小技巧是善用estat icc和estat recovariance这两个后估计命令前者帮你快速回答“数据有没有组间差异”后者帮你检查随机斜率模型的方差协方差结构是否合理。这两个命令的输出直接放进论文的表格里也够规范。还有一点提醒Stata里HLM的后估计能力虽然不错但如果后续要画复杂的交互效应图margins和marginsplot一定要配合用。随机斜率模型的交互效应可视化这两个命令能省掉你一半手动计算的功夫。希望这篇东西能帮你把HLM在Stata里的语句流程走通。数据建模这种事最大的障碍不是命令记不住而是不知道每一步到底在做什么。命令只有那几十个字母背后的模型逻辑想通了一切就顺了。