从贝叶斯视角看多层模型

从贝叶斯视角看多层模型 | 统计之都

多层模型常被用于处理嵌套数据 (即，具有层次结构的数据)，如，从不同的学校中抽取学生样本 (第一层为学生，第二层为学校，学生样本嵌套于学校中)，或从不同社区中抽取居民样本等。追踪数据也可被视为多层数据中的一种，其中第一层为不同的测量时间点，第二层为个体。为阐述多层模型的原理，以一个仅包含一个变量的两层模型为例，模型对应公式如下：

Level1:yij=θj+eijLevel2:θj=μθ+rj

其中，yij 代表第 j 个群组中的第 i 个被试在变量 y 上的得分，θj 代表第 j 个群组在变量 y 上的均值，μθ 代表 θj 的均值。eij 代表个体水平的误差，服从于 N(0,σ2e) 的分布，rj 代表群体水平的误差，服从于 N(0,σ2θ) 的分布。

对于这样的数据，如果研究者认为 θ1=θ2=⋯=θj=μθ，则可以忽略分组变量进行建模，其对应模型如下： yi=μ+ei 而如果研究者认为各组之间无法相互提供信息，可以对每组分别建立 Level 1 的模型。

多层模型则介于两者之间，认为 θj 来源于 N(μθ,σ2θ) 的正态分布，进而将组间的变异 σ2θ (组间差异带来的误差方差) 与组内的变异 σ2e (个体间差异带来的误差方差) 分离开。

贝叶斯思想

多层模型本质上就蕴含着贝叶斯思想，它将参数 θj 视为随机变量而非定值 μθ。而贝叶斯方法和传统频率学派方法本质的区别是：频率学派将未知参数看作常数，根据样本参数估计总体参数；而贝叶斯方法则将未知参数视为随机变量，分析的目的就是得到参数的分布 (王孟成，邓倩文，毕向阳，2017)。

而在贝叶斯多层模型中，不仅将 θj 参数视为随机变量，所有参数 (如，μθ 参数) 都会被视为随机变量。研究者对于这些随机变量可以提供一定的先验信息，结合先验信息和数据似然函数，就可以得到参数的后验分布，而通过马尔科夫链蒙特卡洛算法，可以从后验分布中迭代地抽取大量的样本近似地反映后验分布，在该算法达到收敛的前提下，可以利用抽取的样本进行参数估计。

下文将通过实例 1 展示贝叶斯多层模型结合先验信息的灵活性，并通过实例 2 演示如何在 Mplus 软件中进行贝叶斯多层建模。

贝叶斯多层模型

Jackman(2009)采用了 18 名棒球运动员在 1970 年赛季的前 45 次击球数据预测每个球员在整个赛季中的平均击中率 θj。分析采用贝叶斯多层模型，在建模时，研究者对 θj 的均值 (所有球员的平均击中率 μθ ) 和方差 σ2θ 提供了一定的先验分布，研究者指出 θj 应当大于 0，且基本小于 0.4，因为上一次有球员的平均击中率大于 0.4 已经可以追溯到 1941 年赛季。研究者进一步猜测平均击中率 μθ 主要分布于 (0.15, 0.3) 之间，基于此，研究者对于 μθ 提供了 N(0.225,(0.15/4)2) 的先验分布，而对于 σ2θ 则提供了 Inv−Gamma(14,0.005) 的先验分布以控制 θj 主要分布于 (0, 0.4) 之间，如下图所示。

贝叶斯估计结果如下图所示 (数据和 R 代码)：

图中第一行 MLE 的结果指采用极大似然估计对于每个个体分别估计其平均击中率，第二行为贝叶斯估计结果，第三行是球员在整个赛季中真实的平均击中率结果。

相比于分组建模，多层建模将 MLE 的结果向总体均值 μθ 拉近，尽管对于每个个体的击中率产生了有偏估计，但是个体间差异更小，预测误差也更小。

实例 2 则以 Heck 和 Thomas(2015) 第三章的数据和模型为例，演示如何在 Mplus 软件中进行贝叶斯多层建模，并灵活地在各个层次间探究变量间的关系。模型如下所示，在组内水平 (Within) 上，研究者希望探究在控制了性别和种族后，员工对于薪水的满意程度 (satpay) 对他们的士气 (morale) 的影响 ( β1 )，而在组间水平(Between) 上，研究者希望探究员工所处部门的整体工资水平 (pctbelow) 对工作士气的影响，及部门的整体工资水平对个体水平上薪水满意度 (satpay) 对员工士气 (morale) 影响的调节作用 (即 β1 是否会随着部门整体工资水平的变化而变化)。

Heck 和 Thomas(2015) 对该模型提供的是极大似然估计，本例中将采用贝叶斯估计，并对参数提供无信息先验分布。Mplus 软件默认会对参数提供无信息先验分布，因此代码中无需额外设定。

首先需要将结果变量士气 (morale) 的方差分解到组内和组间水平(代码和结果)，结果显示，个体间方差为 33.312，群体间方差为 5.506，计算组内相关系数 ICC(Intraclass Correlation) 为 0.142。ICC 指组间方差占总方差的比例，公式如下：

ICC=σ2between/(σ2between+σ2within)

接下来建立假设模型，Mplus 代码如下所示，在模型部分，%Between% 设定群体水平的模型，%Within% 部分设定个体水平的模型。在个体水平部分，用性别、人种和 satpay 预测 morale，并将 satpay 对 morale 的效应赋值给 S。在群体水平部分，用 pctbelow 预测 morale 和 S(即检验 pctbelow 对个体水平上薪水满意度 (satpay) 对员工士气 (morale) 影响 ( β1 ) 的调节作用)。

模型对应公式如下，其中 ϵij,μ0j 和 μ1j 代表误差项:

Level1:yij=β0j+β1jsatpayij+β2femaleij+β3whiteij+ϵijLevel2:β0j=γ00+γ01pctbelowj+μ0jβ1j=γ10+γ11pctbelowj+μ1j

TITLE:	 Model 4: Explaining variation in the level-2 intercept and slope;
DATA:  FILE IS: C:\ Mplus1\ch3new.dat;
	     Format is 5f8.0, 3f8.2;
VARIABLE:	
        Names are deptid morale satpay female white pctbelow lev1wt lev2wt;
        Usevariables are deptid morale satpay female white pctbelow;
        Cluster is deptid;  !分组变量为deptid
        Between = pctbelow; !定义组间变量
        Within = satpay female white; 
Define:	Center satpay female white pctbelow (grandmean);  !中心化处理
ANALYSIS:	
        Type= Twolevel random;   !Twolevel定义两层模型，random定义随机斜率(β1)模型
        Estimator = Bayes;  !定义估计方法
        Biterations = (20000);  !定义模型最小迭代次数为20000次
        Point = Median;  ! 定义点估计报告后验分布的中位数结果
Model:	%Between%
    	morale S on pctbelow;  !用pctbelow预测morale和S
    	S with morale;  !估计morale和S之间的相关
    	%Within%
    	morale on female white;  !用性别和人种预测morale 
    	S | Morale on satpay;  !用satpay预测morale并用S表示satpay对morale的效应，即β1
OUTPUT:	TECH1;
PLOT: TYPE = PLOT2;	!输出各参数的后验分布直方图等

Mplus 软件会在算法迭代收敛后输出估计结果，运行过程如下所示

其中，ITERATION 代表算法迭代次数，每迭代 100 次会输出 PSR(Potential Scale Reduction) 值，当迭代次数大于研究者设定的最小迭代次数且 PSR 值满足收敛标准后，Mplus 软件会停止迭代并输出结果。结果如下所示：

结果显示，员工对于薪水的满意程度将显著的正向预测他们的工作士气 (β1=1.196,p<0.001)，控制变量性别的效应不显著 (𝛽=0.005,p=0.473)、种族的效应显著 (𝛽=0.911,p<0.001)。所处部门的整体工资水平也能显著预测他们的工作士气 (𝛽=−0.026,p<.001)，但未发现整体工资水平对薪水满意度和士气的关系间存在调节作用 (𝛽=0.001,p=0.066)。上述无信息先验下贝叶斯方法的估计结果和书中极大似然法的估计结果类似。此时，组内方差被解释了 47.6%，残差为 17.459，组间方差被解释了 67.5%，残差为 1.791。

采用 R 软件中的 rjags 软件包得到了类似的结果 (具体代码和结果)。

多层模型在社会科学等领域的应用非常广泛，尤其是近年来，越来越多的心理学研究开始采用密集追踪研究探究变量间的关系，而多层模型也是处理密集追踪数据的一种常用分析方法。希望本文的介绍可以帮助大家更好地理解多层模型，如果研究者想要深入了解多层模型及其在 Mplus 软件中的实现方法，也可以参阅 Heck 和 Thomas(2015) 一书。

此外，在进行多层数据分析时，实证数据中容易出现第二层样本量不足 (如，学校、社区的数目较少)。在这种组数较少的情况下，依赖于大样本渐近理论的频率学方法(如，极大似然估计) 容易出现对组间方差的有偏估计，此时采用贝叶斯方法进行模型估计相对更为稳健。而贝叶斯多层模型也可以灵活地结合先验信息 (如，以往研究结果或理论) 进行参数估计，如果没有准确的先验信息，也可以提供无信息先验分布或模糊信息先验分布。

本例中仅展示了在 Mplus 软件中提供无信息先验分布进行贝叶斯多层建模的方法，如果研究者需要提供信息先验分布，可以参阅 Mplus 用户手册 (Muthén & Muthén, 1998-2017)。Mplus 8.3 也可以被用于分析贝叶斯三层模型，如不同时间点嵌套于个体，个体同时也嵌套于群组的数据。尽管 Mplus 软件对于一些特殊、复杂的模型难以实现，且能够提供的先验分布的类型也不如 Stan 和 JAGS 丰富 (如，无法提供双指数先验分布)，但是它易于入门，也基本能够满足大部分研究者的建模需求。

此外，在 R 软件中也可以采用 brms 软件包或 rjags 软件包进行贝叶斯多层建模，它们分别依赖于 Stan 和 JAGS 这两个常用的贝叶斯建模软件，实例 1、2 均给出了采用 rjags 软件包建模的代码。

王孟成, 邓倩文, 毕向阳. (2017). 潜变量建模的贝叶斯方法. 心理科学进展, 25(10), 1682–1695.

Jackman, S. (2009). Bayesian Analysis for the Social Science. John Wiley & Sons, Ltd.

Heck, R. H., & Thomas, S. L. (2015). An Introduction to Multilevel Modeling Techniques: MLM and SEM Approaches Using Mplus. New Jersey: Lawrence Erlbaum Associate Inc.

Muthén, L. K., & Muthén, B. O. (1998-2017). Mplus User’s Guide. Eighth Edition. Los Angeles, CA: Muthén & Muthén.

中山大学心理学系研究生，师从潘俊豪副教授，研究方向为贝叶斯结构方程建模。个人博客 https://www.lijinzhang.xyz/blog.html 张沥今