如何使用 Stata 进行多层次回归分析？
原创
每个人都有意或无意地生活在特定社会结构之下，人类行为的异质性决定了社会结构的多样性和复杂性。当数据具有多层次结构时，常规的在同一层次上所开展的多元回归分析便不再适用，需要更为适配的多层次建模方法（Multilevel models）。使用多层次回归分析可以帮助我们进行正确的推断、探索群体或组效应、估计组效应的同时估计组层面自变量的影响，以及推断组的总体。
1. 何谓多层次结构数据？
作为社会动物，人类行为的本质在于彼此间无处不在的互动，而人类彼此间的互动又往往是在一个个局部群落环境下发生的，由此便有了“人以群分”之说。对于“人以群分”，简单的理解通常为特定时点下根据某一分类因素而对个体进行的区分，个体汇集在分类因素所包含的各个水平或类别之下。例如，社区与居民，不同社区内部有各自的居民；再比如班级与学生，通常一所学校内同一年级有许多班级，每一个班级内部是学生。因此，我们总能够将微观层面彼此不同的个体按照某一宏观特征因素进行划分，这时便构成了一个具有嵌套关系的多层次结构（a Multilevel Structure）。
下图展示了一个典型的“学生（Student）—学校（School）” 的双层嵌套结构。 分析这样的数据需要进行多层级建模，多层次建模是定量社会科学研究中用于对具有复杂层次结构的数据（Data with complex hierarchical structures）进行建模分析的基础技术之一。
2. 多层次结构数据示例
载入一份包含 65 所学校、 4059 名学生的样本数据集。
该数据集包含了个体层面（层1）学生的总成绩（y）、阅读成绩（x1）以及性别（x3，女性=1），还包含学校层面（层2）的学校类型（x3，男女混合=1、男子学校=2、女子学校=3）。
clear all
use "$raw_data_multilevel\schools", clear
levelsof school
return list
dis r(r) //65
/* Observations: 4,059
Variables: 7 11 Nov 2009 11:36
----------------------------------------------------------------
Variable Storage Display Value
name type format label Variable label
----------------------------------------------------------------
school float %9.0g School id
student float %9.0g Student id
y float %9.0g Score
x1 float %9.0g Reading test
x2 float %9.0g Female=1
x3 float %10.0g x3 Type of school */接下来，可将每所学校内所有学生的成绩以及每所学校的学生平均成绩进行绘制。
由图可知，学生成绩除了在学校内部存在差异，在学校间也存在差异。因此，若仅在个体层面或学校层面开展传统的回归分析，很可能存在较为严重的估计问题。就个体层面而言，忽略学校层次意味着假设学校内部的学生成绩差异在学校之间并无较大差异；但在现实中，学校内部的个体往往具有群聚特征，使得学校间存在差异，这便会产生异方差问题。就学校层面而言，忽略个体层次则可能导致生态学谬误的出现，即加总层面得到的关系无法直接推论至个体。
bysort school: egen y_mean=mean(y)
twoway scatter y school, msize(small) mcolor(gray) msymbol(oh) ///
|| connected y_mean school, lwidth(medthick) lcolor(black) msymbol(none) ///
||, ytitle("学生得分", size(medium)) xtitle("学校 id", size(medium)) ///
legend(off) plotregion(fcolor(white) lcolor(white)) graphregion(fcolor(white) lcolor(white)) 
3. 不同层 2 单元下的 OLS 回归
基于上面的示例数据，我们希望探讨阅读成绩（x1）对总成绩（y）的影响。在层级数据结构下，若仍使用传统的基于单一层级的回归分析方法，可对每所学校分别进行回归（实际上就是根据学校 id 分类的子样本回归）。由于有 65 所学校，因而可以得到 65 组子样本回归结果，每组回归结果包含了每所学校的截距和斜率系数。
全样本和子样本的样本回归方程表示如下：
下图对这些估计结果进行了图示，显而易见的是，黑色粗线所代表的全样本回归无法代表各所学校的情况。因此，需要使用多层次回归建模方法将组间变异和组内变异同时纳入到分析之中。
*计算每所学校的斜率和截距，保存至数据集
statsby intercept = _b[_cons] slope = _b[x1], by(school) saving("$temp_data_multilevel\OLS", replace): reg y x1
*与原始数据集进行合并
use "$raw_data_multilevel\schools", clear
merge school using "$temp_data_multilevel\OLS"
drop _merge
gen y_fitted = intercept + slope*x1, before(y)
sort school x1
separate y_fitted, by(school) //seperate 命令能够按照 “by()” 中的唯一标识生成新变量
twoway line y_fitted1-y_fitted65 x1 ///
|| lfit y x1, lwidth(thick) lcolor(black) ///
legend(off) ytitle("y") ///
plotregion(fcolor(white) lcolor(white)) graphregion(fcolor(white) lcolor(white)) 
4. 多层次回归
多层次回归分析的建模逻辑并不复杂，其实质是组内方程的部分或全部参数（截距项或斜率项）作为因变量，进而用组间方程加以解释。因此，多层次回归分析也是一种类似于元分析的“回归的回归”。仍以上面的数据为例，考虑更多学生层面解释变量（假设有 k 个，数据集中有 2 个）的组内方程可表示如下：
我们注意到，式（3）除了加入更多自变量外，与式（2）并无差异，两者都是个体层面的组内方程，将位于第 j 所学校（层 2）的第 i 名学生（层1） 的成绩表示为学生个体变量（比如数据集中的阅读成绩和性别）和随机扰动项的函数。多层次回归模型的关键则在于探讨层1的回归系数（如 \beta_{kj}）如何在层 2 单位间发生变化，因而可进一步表示为层 2 解释变量和随机扰动项的函数：
假设层 2 有 p 个解释变量（数据集中有 x_3），k = 0, 1, 2, ..., m；显然，当 k=0 时，式（4）为截距项的层 2 回归方程，当 k≠0 时，式（4）为斜率项的层 2 回归方程。基于此，可以将式（4）带入式（3），进而得到完整的回归模型。
更为重要的是如何对 \gamma_{kp} 进行解释，其含义为某个层 2 变量（比如数据集中的学校类型 x_3）如何对层1中的截距项 \beta_{0j} 或特定变量（x_k）的回归系数（\beta_{kj}）产生影响，体现为层 2 变量如何影响层1 因变量。
4.1 变截距模型（零模型）
零模型是多层次回归分析的起点。
顾名思义，所谓零模型，即层 1 和层 2 均未加入自变量，只有截距项。
mixed y || school: , mle nolog
/*
Mixed-effects ML regression Number of obs = 4,059
Group variable: school Number of groups = 65
Obs per group:
min = 2
avg = 62.4
max = 198
Wald chi2(0) = .
Log likelihood = -14851.502 Prob > chi2 = .
------------------------------------------------------------------------------
y | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_cons | -.1317107 .5362734 -0.25 0.806 -1.182787 .9193659
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects parameters | Estimate Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
school: Identity |
var(_cons) | 16.86388 3.28458 11.51248 24.7028
-----------------------------+------------------------------------------------
var(Residual) | 84.77541 1.897109 81.13751 88.57642
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 498.72 Prob >= chibar2 = 0.0000 */由结果可计算组内相关系数（Intra-class Correlation Coeficient, ICC）：
ICC 越小则组内相关性越低，表明组内个体间差异较大，因而 “组内方差” 也越大，其统计学意义在于组内个体之间具有更高的独立性，这时是否使用多层次模型便成为一个经验判别问题。通常将 0.059 作为阈值（Cohen, 1988）。示例中，组间方差 = 16.86388，组内方差 = 84.77541，由此可得 ICC ≈ 0.17，表明使用多层次回归分析是合适的。
4.2 变截距模型（层 1 加入自变量）
更为一般化的变截距模型可在在层 1 加入更多个体变量。
mixed y x1 || school: , mle nolog
/*
Mixed-effects ML regression Number of obs = 4,059
Group variable: school Number of groups = 65
Obs per group:
min = 2
avg = 62.4
max = 198
Wald chi2(1) = 2042.57
Log likelihood = -14024.799 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
y | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
x1 | .5633697 .0124654 45.19 0.000 .5389381 .5878014
_cons | .0238706 .4002255 0.06 0.952 -.760557 .8082982
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects parameters | Estimate Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
school: Identity |
var(_cons) | 9.212859 1.85304 6.211364 13.66476
-----------------------------+------------------------------------------------
var(Residual) | 56.57267 1.266255 54.14451 59.10973
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 403.27 Prob >= chibar2 = 0.0000 */4.3 变截距、变斜率模型
进一步，考虑层 1 的截距项和解释变量的斜率，但在层 2 不考虑其他解释变量。
此时 ICC = [var(_x1_) + var(_\_cons_)] / var(_Residual_)
mixed y x1 || school: x1, mle nolog covariance(unstructure)
/*
Mixed-effects ML regression Number of obs = 4,059
Group variable: school Number of groups = 65
Obs per group:
min = 2
avg = 62.4
max = 198
Wald chi2(1) = 779.79
Log likelihood = -14004.613 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
y | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
x1 | .556729 .0199368 27.92 0.000 .5176535 .5958044
_cons | -.115085 .3978346 -0.29 0.772 -.8948264 .6646564
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects parameters | Estimate Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
school: Unstructured |
var(x1) | .0145358 .0045773 .0078415 .0269452
var(_cons) | 9.04472 1.831021 6.08243 13.44971
cov(x1,_cons) | .1804042 .0691523 .0448682 .3159402
-----------------------------+------------------------------------------------
var(Residual) | 55.36531 1.249281 52.97012 57.86881
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 443.64 Prob > chi2 = 0.0000
Note: LR test is conservative and provided only for reference. */4.4 变斜率模型
当仅考虑层 1 自变量时，与常规的 OLS 回归结果几乎没有差异。
若纳入层 2 自变量，则意味着出现层际交互效应。
mixed y x1 || _all: R.x1, mle nolog
/*
Mixed-effects ML regression Number of obs = 4,059
Group variable: _all Number of groups = 1
Obs per group:
min = 4,059
avg = 4,059.0
max = 4,059
Wald chi2(1) = 2186.09
Log likelihood = -14226.433 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
y | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
x1 | .5950551 .0127269 46.76 0.000 .5701108 .6199995
_cons | -.011948 .1263914 -0.09 0.925 -.2596706 .2357745
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects parameters | Estimate Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
_all: Identity |
var(R.x1) | 1.70e-10 4.91e-08 4.0e-257 7.3e+236
-----------------------------+------------------------------------------------
var(Residual) | 64.84141 1.439323 62.08088 67.7247
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 0.00 Prob >= chibar2 = 1.0000 */
reg y x1
/*
Source | SS df MS Number of obs = 4,059
-------------+---------------------------------- F(1, 4057) = 2185.01
Model | 141749.015 1 141749.015 Prob > F = 0.0000
Residual | 263191.301 4,057 64.8733797 R-squared = 0.3500
-------------+---------------------------------- Adj R-squared = 0.3499
Total | 404940.316 4,058 99.7881508 Root MSE = 8.0544
------------------------------------------------------------------------------
y | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
x1 | .5950551 .0127301 46.74 0.000 .5700972 .6200131
_cons | -.011948 .1264225 -0.09 0.925 -.2598056 .2359095
------------------------------------------------------------------------------ */--- to be revised ---
初稿日期：2024年09月23日
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
腾讯云开发者
