贝叶斯混合效应模型：基于brms 的应用教程

潘晚坷 温秀娟 金海洋

（1 南京师范大学心理学院，南京 210097）

（2 广州医科大学附属脑科医院，广州 510370）

（3 Department of Psychology,Division of Science,New York University Abu Dhabi,PO Box 129188,Saadiyat Island,Abu Dhabi,United Arab Emirates）

1 引言

近年来，心理学研究所使用的统计方法日渐丰富。随着新的统计方法的推广与普及，心理学研究者逐渐意识到了一些传统统计方法的局限性，实验心理学中最常用的重复测量方差分析就是其中之一。首先，它处理的是汇总后（aggregated）的数据，比如每个被试在每个条件下的平均反应时。汇总后的数据忽略了每个被试个体内部的变异，这可能会导致无法准确估计自变量的效应（Haaf &Rouder,2017;Veenman et al.,2022）。其次，它通常只能考虑一个随机因素（Judd et al.,2017），而无法同时考虑多个随机因素。例如，它一般把被试作为随机因素，借此通过被试样本的数据推断被试总体的效应。类似地，重复测量方差分析可以把刺激材料作为随机因素，通过实验中使用的刺激材料推断刺激材料总体的效应，但它无法同时考虑被试和刺激材料这两个随机因素。最后，当某个被试在某一条件下的汇总数据缺失时（例如在纵向研究中某些被试未能参加某个时间点的测试），通常只能删除存在缺失值的被试的所有数据，这会造成数据进一步流失，可能导致有偏差的估计（Magezi,2015;Matuschek et al.,2017）。

相比之下，线性混合效应模型（Linear Mixedeffects Models，LMM），亦称多水平模型（Multilevel Models）或层级模型（Hierarchical Models），可以更好地处理这些问题。混合效应模型可以考虑被试在各个条件内的变化性及数据中可能存在的相互依存关系和层级结构，并且能够更好地处理数据缺失等问题（Singmann &Kellen,2019;Tuerlinckx et al.,2006）。例如，在以学生为被试的实验中，学生来自不同的班级，不同的班级属于不同的学校；同一个班级或同一个学校学生的表现一般更相似。混合效应模型可以更好地考虑学生、班级和学校之间的层级关系，对实验效应的估计更准确。统计工具（如 R 语言）的成熟使得在心理学研究中使用混合效应模型分析数据变得更加容易，越来越多的研究也开始使用该模型（Bono et al.,2021）。

与（频率论的，Frequentist）线性混合效应模型相比，贝叶斯线性混合效应模型（Bayesian Linear Mixed-effects Models，BLMM）因 其在基于贝叶斯统计推断（即依照贝叶斯定理，通过结合先验信息（prior）与模型似然性（likelihood）对参数的后验（posterior）进行统计推断）以及处理复杂模型等方面的优势而受到了越来越多的关注（Bürkner,2017;Sorensen &Vasishth,2016）。首先，与频率论统计相比，贝叶斯分析的一大优势是可以量化参数估计的不确定性，即它的结果是一个分布（例如后验概率密度分布），其宽窄描述了对参数分析的不确定性。其次，贝叶斯分析的可信区间（credible interval）可以理解为一个拥有一定概率包含真值（true value）的区间，比如95%的可信区间表示有95%的可能性包含真值（Morey et al.,2016）。但该理解通常被误认为是频率论置信区间（confidence interval）的属性（Morey et al.,2016）（95%置信区间表示在通过无数样本所获得的无数个区间中，大约有95%的区间包含了真值（Lambert,2018））。其次，基于频率论的统计方法一般不能为支持虚无假设提供证据（Cohen,1994），而贝叶斯分析中的贝叶斯因子（Bayes Factor，BF）或实际等价区域（Region of Practical Equivalence，ROPE）可以（Makowski et al.,2019）。最后，如果混合效应模型试图包含复杂的随机效应结构，基于频率论的模型通常难以收敛，从而无法给出可靠的估计（Matuschek et al.,2017），而贝叶斯模型能够更好地处理该结构（Sorensen &Vasishth,2016）。此外，贝叶斯分析还可以通过整合研究者的先验知识提供更合理的参数估计。

贝叶斯混合效应模型兼具（频率论）混合效应模型和贝叶斯分析的优点，是心理学研究中极具潜力和值得推广的统计方法。本文致力于向读者提供一个相对容易理解和上手的教程。

2 模型介绍

2.1 混合效应模型

混合效应模型（Baayen et al.,2008）是一般线性模型的扩展形式，得名于它同时包含了固定效应（fixed effects）和随机效应（random effects）。固定效应通常是研究者所感兴趣的效应，即自变量在总体水平的效应。与之相对的，随机效应是自变量在随机因素（例如被试或刺激材料）各个水平的特有效应。例如，在探究图片类型对情绪唤醒的影响的实验中，每个被试受到图片类型的影响可能并不完全相同。其中，图片类型对所有被试的平均影响为固定效应，而每个被试受到图片类型的影响与平均影响的差值为随机效应。

为更好地解释固定效应与随机效应，我们首先回顾最基本的一般线性模型：

其中，Yi,j是因变量，比如被试的脑电波幅（详见3 应用示例），表示第j名被试在第i试次的观测数据。Xi,j为自变量，表示研究者感兴趣的实验操作，可以是组间或者组内变量。参数β0是截距，表示所有被试的均值，比如所有被试的平均脑电波幅；参数β1为斜率，表示自变量对因变量的影响，比如不同图片类型对被试脑电波幅的影响（本文中涉及截距和斜率的解释仅适用于应用逐次差分对比编码（Successive Differences Contrast Coding，即本文所使用的编码）且只有两个水平的自变量（Venables &Ripley,2002）。当使用不同的编码时，截距和斜率的解释也会有所不同，详见 Schad et al.,2020）。；∈i,j为残差，表示模型中不能被解释的部分。

需要注意的是，使用一般线性模型需要满足残差独立性和正态性的前提预设，即不同观测数据的残差应该相互独立并且服从正态分布。然而，在被试内实验设计中，每个被试完成不同图片类型的多个试次，因此他们在不同试次上的表现并不是完全独立的。此时使用一般线性模型无法处理观测数据间的依存性。此外，公式（1）的模型没有考虑不同被试间的个体差异，比如，不同被试的平均脑电波幅（即每个被试的截距）可能存在差异，而忽略这些差异可能会掩盖图片类型的真实效应。为更好地考虑这些方面，可以在公式1的基础上增加随机截距：

与公式（1）相比，公式（2）主要区分了固定截距（β0）和随机截距（S0,j）。其中，固定截距可以理解为所有被试截距的平均值，即所有被试的平均脑电波幅；而随机截距可以理解为每个被试的截距与固定截距的差值，即每个被试各个条件的平均脑电波幅与所有被试平均脑电波幅的差值。虽然加入随机截距可以考虑到各被试平均脑电波幅的差异，但仍然没有考虑到自变量对因变量的影响在不同被试上的差异。比如，不同被试受到图片类型的影响并不一定相同。为更好地解释自变量的影响在不同被试间的差异，可以进一步在混合模型中加入随机斜率：

与公式（2）相比，公式（3）主要区分了固定斜率（β1）和随机斜率（S1,j）。其中，固定斜率是所有被试斜率的平均值，即所有被试受到图片类型的平均影响；随机斜率是每个被试斜率与平均斜率的差值，即每个被试受到图片类型的特有影响。需要注意的是，由于随机截距与随机斜率都受到同一个随机因素（被试）的影响，因此两者可能存在一定的相关性，比如对于平均脑电波幅（随机截距）更高的被试，其在不同图片类型下的脑电波幅的差异（随机斜率）可能更大。此时通过预设两种随机效应来自同一个多元正态分布可以更好地考虑他们之间的相关性。

总的来说，混合效应模型是一般线性模型的扩展，它能够借助随机截距和随机斜率考虑随机因素（如被试和实验材料）对实验效应的影响，可以更好地处理数据中的依存和层级关系，弥补重复测量方差分析的不足，提高对实验效应估计的准确性。

2.2 贝叶斯混合效应模型

在构建混合效应模型的过程中，为了在控制合理的一类错误率的同时保持良好的统计效力，研究者通常需要考虑最大化的模型（maximal model）（Barr et al.,2013），即把所有可能的随机截距和随机斜率以及它们之间的相关系数加入模型中。然而，由于模型参数增多，可能导致模型难以收敛，从而给出不可靠的估计（Bates et al.,2015）。相比之下，贝叶斯混合效应模型可以更好地处理这种复杂的模型。更重要的是，相较于频率论混合效应模型，贝叶斯混合效应模型具有考虑先验信息和提供更直观的参数后验分布等优势。这些优点使得贝叶斯混合效应模型更受欢迎。

贝叶斯定理（Bayes’ theorem）在贝叶斯数据分析中起着重要作用，它借助数据更新先验信息，从而得到后验信息：

其中，p（θ|Y）为参数的后验分布（posterior distribution），即连续型参数对应的概率密度函数或离散型参数对应的概率质量函数，表示在观测到数据Y的基础上参数θ的概率分布。对于贝叶斯混合效应模型来说，θ包括所有的固定效应β、随机效应S、随机效应间的相关性和残差∈（见公式3）所对应的分布的所有参数。p（θ|Y）可以由公式右侧部分计算得来，其中p（θ）为参数的先验分布（prior distribution），表示研究者在观测到数据前对参数θ的先验知识。p（Y|θ）为似然，描述了当参数θ为不同的值时，观测到当前数据的似然值。值得注意的是，p（Y）并不包含θ，这是因为它是一个边际似然（marginal likelihood），是通过对所有可能的θ的似然值与相应的先验概率乘积求和或积分而获得的（当θ为离散变量时：p（Y）=∑θp（Y|θ）*p（θ）；当θ为连续变量时，p（Y）=∫θp（Y|θ）*p（θ）dθ，（Kruschke &Liddell,2018）。p（Y）是一个“归一化”因子（normalizing factor），能够确保后验分布是一个有效的概率分布（使后验分布的面积为1）。

基于贝叶斯定理，在对贝叶斯混合效应模型进行参数估计时，研究者首先需要为各参数设置合理的先验分布，之后通过数据对先验进行更新，从而获得参数的后验分布。需要注意的是，随着模型变得复杂，贝叶斯定理中的边际似然（公式4右侧分母部分）将变得难以计算。此时通常通过马尔科夫链蒙特卡洛（Markov Chain Monte Carlo,MCMC,Lambert,2018;Morey et al.,2011）算法推测参数的后验分布（Bürkner,2017），其具体原理可参考Plummer（2023）。与MCMC相关的重要参数会在下文中结合示例予以介绍。

3 应用示例

我们将借助一个假想的心理学实验展示如何使用brms工具包进行贝叶斯混合效应模型的数据分析。具体来说，本节将从数据模拟和预处理、模型构建、先验设定、模型的拟合与检验和统计推断五个方面介绍贝叶斯分析最基本的工作流程（Gelman et al.,2020）。具体代码见https://osf.io/cgxwn/。

3.1 数据模拟与预处理

该假想研究的目的是探究抑郁症人群加工不同情绪刺激的神经基础。实验中，40名抑郁症患者和40名健康对照组被试观看30张积极图片和30张中性图片，期间我们采集了他们的脑电数据。因变量是晚期正电位（late positive potentials,LPP）的波幅。该实验为2（组别group：抑郁症组depression、对照组 control）×2（图片类型type：积极positive、中性neutral）混合实验设计，其中组别为被试间因素，图片类型为被试内因素。

实验数据使用R的faux工具包生成（DeBruine,2021），具体代码见线上补充材料数据模拟部分（https://osf.io/3g79b）。简单来说，设定（图1和表1）图片类型的主效应为6.5（即积极图片诱发的LPP比中性图片高6.5），组别的主效应为0.1（即对照组LPP比抑郁组高0.1）和组别与图片类型的交互作用为0.1。

表1 模拟数据表述统计结果

图1 模拟数据柱状图

首先需要确保数据集的结构符合工具包的要求。重复测量方差分析通常需要“宽”格式数据，即每一行是一个被试的数据。而混合效应模型的分析通常需要 “长”格式数据，即每一行是一个试次的数据。例如，R语言中的brms工具包需要使用长数据格式进行混合效应分析。与大部分实验程序（如E-prime）生成的数据类似，我们通过faux工具包生成的模拟数据也是以长格式保存，符合brms工具包输入数据的要求。该模拟数据集的变量名为df_simu，我们可以通过head(df_simu,10)查看前10行数据（图2）。

此外，我们还需要对数据集df_simu中的称名（自）变量进行编码。这里我们使用逐次差分对比编码对组别（group）和图片类型（type）这两个自变量进行编码（图3），即使用R中的MASS::contr.sdif()函数。通过此设定，贝叶斯模型的输出结果刚好对应大部分研究者熟悉的主效应和交互作用（见公式1处解释）。

3.2 构建模型

完成数据的预处理后，研究者需要构建模型。模型构建主要有两种方法。一种是从包含了研究者最感兴趣效应的简单模型开始，逐渐添加更多的效应以寻找能够最优解释数据的模型。另一种方法是根据实验设计构建最大化的模型（Barr et al.,2013）。鉴于第二种方法相对简单，且更多地被认知和心理学研究者所采用（Schad et al.,2021），本文将主要介绍如何根据实验设计构建最大化模型（对第一种途径感兴趣的读者可参考Schad et al.,2021）。

在实际分析中建立贝叶斯线性混合效应模型时，可以直接根据实验设计设定最大化的模型。需要注意的是，本文为了更清晰地介绍如何理解和设定固定效应和随机效应，将以三个模型为例，分别解释如何结合实验设计使用brms设定固定效应、随机截距和随机斜率。

首先建立一个仅包含固定效应的模型（model_1，见图4）。具体来说，我们将组别和图片类型两个自变量以及他们的交互作用作为固定效应。该模型假定所有被试的平均波幅相同，且图片类型对每个被试的影响也相同。

图4 固定效应模型

在图4的代码中，LPP为因变量，即被试观看不同图片时LPP的波幅，group*type表示组别和图片类型两个自变量的主效应以及他们的交互作用。df_simu为数据集。其他参数是模型拟合时所基于的MCMC采样算法的相关参数。为了保证后验分布的有效性，MCMC通常会同时运行多个独立的采样过程，每个独立的采样过程称为一条链（chain）。只有当各个独立链的采样“融合”在一起时（见图6），MCMC获得的参数后验分布才更可能是有效的（McElreath,2020）。在model_1中，我们设置了4条链（chains=4，默认为4）。与之相对应的是cores参数，它可以设定运行MCMC时使用的计算机中央处理器核心（Central Processing Unit cores）的个数。设置cores=4（默认为1）表示使用四个中央处理器核心同时分别各运行一条马尔可夫链（共四条链），借此提高数据分析的速度。iter参数（iteration的简写）指的是每条链的迭代次数，即对每个参数进行多少次的采样。warmup参数（或称为burn in）设定了在每条链采样的开始阶段用于“校准”采样器（sampler）以提升采样效率的迭代次数（详见https://mc-stan.org/docs/reference-manual/hmc-algorithm-parameters.html），该阶段的采样不会被纳入最终的参数概率分布。thin参数为“稀释比例（thinning rate，默认为1）”，表示每隔thin个采样保留一个样本，研究者可以通过设定较大的thin来减少电脑内存的需求（详见https：//mc-stan.org/docs/referencemanual/effective-sample-size.html#thinningsamples）。总的来说，每条链最后保留的采样数量与iter、warmup和thin相关，最后保留的采样数等于（iter-warmup）/thin。例如，在model_1中每一条链会保留每个参数的4000个样本，即（5000-1000）/1，最终四条链一共会保留16000个样本。

通过summary()函数可以展示模型拟合的相关参数和结果，详见图5。

输出结果的第一部分展示了模型设置的信息。例如，“Family：gaussian，Links：mu=identity；sigma=identity”表示模型预设自变量和因变量间的关系是线性的，且残差是一个正态分布。对于其他非线性关系的模型，研究者可以设定相应的残差分布和链接函数（详见第4部分）。

Population-Level Effects部分为固定效应的结果，通常是研究者最为关注的部分。其中，Intercept是模型的固定截距，表示所有被试的平均LPP波幅。group2M1表示在平均图片类型的效应之后，健康对照组的LPP波幅减去抑郁症组LPP波幅的差值，即组别的主效应。与此类似，type2M1是在平均组别的效应之后，积极图片诱发的LPP波幅减去中性图片诱发的LPP波幅的差值，即图片类型的主效应。group2M1：type2M1是组别和图片类型交互作用，即（LPP健康对照组，积极图片-LPP健康对照组，中性图片）-（LPP抑郁症组，积极图片-LPP抑郁症组，中性图片）。

在输出的结果中，Estimate表示参数的后验分布的均值，Est.Error表示后验分布的标准差。l-95%CI和u-95%CI表示参数的95%可信区间的两个边界，例如group2M1的95%的可信区间是［0.11,0.46］。Rhat表示不同链条采样的融合程度。Rhat越接近1，表示不同链条采样融合得越好，即模型的收敛情况越好。如果Rhat大于1.1（或更严格的标准下，如大于1.01），则说明模型的收敛情况并不理想（Gelman &Rubin,1992;Vehtari et al.,2021）。此时，研究者不应该使用获得的后验分布进行进一步分析，而应该检查模型设定和数据是否存在问题，并进行相应调整。ESS为有效样本量，是对MCMC算法采样到的有效且独立的样本数量的估计。ESS越高表明能用来有效估计后验分布的样本数量越大。通常建议Bulk-ESS应是MCMC链数量的100倍及以上。例如，当使用四条MCMC链时，Bulk-ESS应至少为400（Vehtari et al.,2021）。为了确保足够的Bulk-ESS，研究者可以适当增加iter，例如本文使用的iter为5000。Tail-ESS表示后验分布两端的有效样本量（effective sample size）。通常情况下，过小的Bulk-ESS和Tail-ESS会伴随着过大的Rhat。设定合理的模型和先验分布一般可以解决大部分Rhat过大的情况。最后，研究者还可以通过可视化的方式检查各条链的融合情况（图6）。

图6 model_1 的MCMC 各条链采样轨迹图（又称“毛毛虫”图）

图6描述了model_1中MCMC不同链抽样的轨迹图，反映了不同链的融合情况。不同链参数的抽样很好地 “融合”在一起，说明其对应的参数后验分布相对可靠。由于融合的链看起来像一条毛毛虫，因此轨迹图又称“毛毛虫”图。如果一条链的抽样与其他链的抽样存在一定差异（通常伴随着较大的Rhat），那这些后验抽样也很可能无法有效表征参数的后验分布（详见线上补充材料第三部分MCMC链不收敛示例，https：//osf.io/3g79b），此时，研究者需要重新考虑模型的合理性（Roy,2020）。

如前所述，仅包含固定效应的模型（model_1）假定所有被试的平均脑电波幅都是相同的，即所有被试共享同一个截距值。但这种情况在现实中不太可能发生。更可能的情况是，每个被试的平均脑电波幅各不相同。因此，我们需要在模型中增加随机截距来考虑被试在平均脑电波幅上的个体差异。具体来说，我们需要在model_1的基础上增加基于被试（by-subject）的随机截距（对应公式2中的S0,i），即(1|subj)。其中subj为被试编号，1表示随机截距，详见图7。

Population-Level Effects为model_2固定效应的结果，其相关参数的含义与model_1中的相同（详见3.2 部分）。Group-Level Effects为随机效应结构的结果，其中sd(Intercept)为随机截距的标准差，反映了LPP平均波幅在不同被试间的离散程度。虽然model_2考虑了不同被试平均LPP波幅的个体差异，但它仍然假定所有被试受到实验处理（即图片类型）的影响是相同的。然而在实际情况中，每个被试受到图片类型的影响更可能是不同的。因此，我们需要设置随机斜率以考虑图片类型对不同被试的不同影响。具体来说，我们可以在model_2模型的基础上加入随机斜率（对应公式3中的S1,i，即 (type|subj)，其中subj为被试编号，type为图片类型的随机斜率。与此同时，我们也在模型中加入了随机截距与随机斜率间的相关，详见图8。

图8 模型3 代码及结果

与model_2结果相比，model_3中的Group-Level Effects部分额外展示了基于被试的随机斜率和它与随机截距相关性的后验分布情况。其中，sd(type2M1)为基于被试的图片类型随机斜率的标准差，反映了不同被试LPP波幅受到图片类型不同影响的离散程度。cor(Intercept,type2M1)为随机截距和随机斜率的相关，该相关结果表明两者之间存在一定的正相关，即平均LPP波幅更大的被试，其对积极图片的LPP波幅也比中性图片更大。

至此，我们完成了对该模拟数据的最大化模型的构建，即根据实验设计设定了所有合理的固定效应和随机效应结构。在该示例中，我们从仅包含固定效应的线性模型出发，在不同的模型中依次展示了如何添加随机截距和随机斜率。需要强调的是，在实际的贝叶斯数据分析中，研究者可以根据实验设计直接设定最大化的模型，而无需在不同的模型中依次添加不同的随机效应。

如前所述，贝叶斯混合效应模型还可以同时考虑多个随机因素，以及不同随机因素间的层级关系。例如，在分析包含“学生—班级—学校”层级关系的数据时，研究者可以通过设定“（1 ｜学校班级学生）”的随机效应结构来考虑同一个班级和同一个学校学生间的相似性等（详见Bürkner,2018）。

3.3 先验的设置

完成构建模型后，研究者需要根据当前的先验知识自定义模型的先验，以在数据分析中整合先验知识。由于先前没有自定义先验分布（例如model_3），因此模型使用的均是brms默认的。这些先验分布通常是弱信息先验（weakly informative prior），即相对实验数据来讲，基本不能为参数估计提供有效信息的先验（Seaman et al.,2012）。使用这些默认的弱信息先验非但不能发挥贝叶斯数据分析的最大作用，反而可能会对数据分析的结果产生负面影响（Seaman et al.,2012）。例如，贝叶斯因子很容易受到先验分布的影响，不合理的先验会导致有偏差的贝叶斯因子。因此，为贝叶斯分析设定合理的先验至关重要。

设定先验分布的第一步是了解可以对一个模型的哪些参数设置先验。这里可以借助get_prior()函数得到相应的信息，详见图9。

图9 先验设置

如结果所示，我们可以为model_3模型的以下参数设定先验：固定截距（class为Intercept）；固定斜率（class为b），包含组别主效应、图片类型主效应和他们的交互作用（分别对应coef为group2M1，type2M1和group2M1:type2M1）；基于被试（group为subj）的随机效应的标准差（class为sd），包含随机截距（coef为Intercept）和随机斜率（coef为type2M1），基于被试的随机效应间的相关（class为cor），以及残差（class为sigma）。由于sd和sigma不能为负数，因此这些参数的最小值为0（即列lb，为low bound的缩写）。

在实际数据分析中，研究者一般可以根据以往研究的结果为当前研究的模型设定先验分布。此时，研究者可以通过先验预测检验（prior predictive check）来测试这些先验的合理性，并进行相应调整，进而为模型确定合理的先验。在进行先验预测检验时，需要考虑领域知识，即研究者对研究领域的专业知识和经验，它可以帮助研究者确定先验分布的合理范围，并检查先验分布是否与领域知识一致。如果先验分布与领域知识相矛盾，则需要重新考虑先验分布的选择。

如果可以从先前类似的研究获得参数的后验概率分布，研究者可以将这些概率分布作为当前研究模型的先验分布。即使无法获得参数准确的概率分布，研究者仍然可以根据先前研究结果以及已有知识为模型设定先验分布。如果无法从先前研究中获取参数的概率分布作为模型的先验，我们可以首先尝试设置弱信息的先验分布，详见图10。

图10 定义先验分布

在图10的代码中，我们把固定截距和固定斜率分别设定为一个均值为0、标准差为10的正态分布（在代码中表示为normal(0,10)）。对于prior(normal(0,10),class=Intercept)，我们可以简单地把它理解为，我们设定所有被试的平均值脑电波幅大致有95%的概率落在-20和20之间（即两个标准差之内）。对于平均脑电波幅来讲，这是一个相对较大的区间。因此，它对应的先验分布是弱信息先验分布。类似地，所有固定斜率（包含组别和图片类型的主效应以及他们的交互作用）也大致落在相同的区间。此外需要注意的是，尽管使用正态分布作为残差sigma和随机效应的标准差sd设置先验，但由于这些参数不能小于0，因此他们实际对应的先验分布是截断正态分布（truncated normal distribution）。最后，随机效应间的相关性设置为弱信息先验prior(lkj(2),class=cor)，其中lkj为LKJ相关性矩阵（LKJ-Correlation prior matrix，Lewandowski et al.，2009）。

为检验先验分布的合理性，我们可以使用先验预测检验考查仅基于先验信息生成的因变量（即LPP的波幅）的情况。进行先验预测检验需要对先验分布进行采样，即在模型中设定先验分布prior=prior_01的基础上加入sample_prior="only"，通过该设置可以仅从先验分布中进行采样（忽略观测数据的似然值信息）。其余参数的设定和我们在分析数据时设定相同（例如model_3），具体代码如图11所示。

图11 prior_01 先验预测检验代码

通过上述代码，我们可以获得16,000组（（5000-1000）*4）来自先验分布的抽样。其中，每一组都包含模型的所有参数，即固定截距，固定斜率和随机效应等。每一组先验参数都可以生成一套新的与数据集（df_simu）结构相同的模拟观测数据（即因变量）。通过考查这些模拟观测数据的合理性（比如检查每一套模拟数据的极值）可以判断先验的设定是否合适。例如在一个以反应时为因变量的stroop实验中，如果那些由先验分布生成的模拟数据集中大部分数据集的最大值都超过预期（比如10秒），就很可能说明模型某些参数先验的设定不太合理（比如截距或标准差过大）。在当前模拟实验中，尽管我们无法根据先前研究获取精确的先验，但是可以利用现有的知识，即试次水平的脑电事件相关电位的波幅通常在一定的范围之内，来确定先验。因此，我们计算16,000套由先验分布生成的模拟观测数据各自的最小值、均值和最大值，并通过观察它们的范围是否符合我们的先验知识来判断先验分布的合理性。

上述分析可以通过brms提供的pp_check()函数实现。图12中的代码展示了prior=prior_01的先验预测的结果。

图12 prior_01 先验预测检验结果代码

图13 基于先验prior_01 生成的16,000 套预测数据（y(pred)）的最小值、最大值和均值的分布

结果显示，基于prior_01生成的每一套模拟观测数据的最小值大部分落在［-75,0］，均值大部分落在［-25,25］，而最大值大部分落在［0,80］。根据先验知识，试次水平的脑电波幅通常很少会低于-50（μV）或高于50（μV）（伪迹除外）。因此，我们可以适当地调整参数的先验分布以符合我们的先验知识，即减小参数先验分布的不确定性。

在接下来的尝试中，我们将固定截距和固定斜率的先验设置成均值为0，标准差为5的正态分布，并进行与之前相同的先验预测检验（见图14）。

图14 prior_02 先验预测检验

经过同样的步骤，我们可以获得基于先验prior_02生成的16,000套模拟观测数据（y(pred)）的最小值、均值和最大值的分布（图15）。结果显示模拟观测数据的最小值大部分落在［-60,0］，均值大部分落在了［-10,10］，最大值大部分落在［0,60］。与prior_01的先验预测结果（图 13）相比，prior_02预测的结果更多落在了［-50,50］内，能够相对更合理地描述我们的先验知识。需要注意的是，确定合理的先验分布通常需要重复多次类似的过程，且并不是先验分布的信息越强就越好，它需要合理反映研究者的先验知识。

图15 基于先验prior_02 生成的每一套预测数据（y(pred)）的最小值、最大值和均值的分布

3.4 模型拟合

在确定先验后就可以开始拟合模型，该模型包括所有可能的随机截距、随机斜率和它们之间的相关，并通过prior=prior_02为模型设定先前确定的先验。此外，我们设定sample_prior="yes"，以此保留仅依靠先验信息获得的采样（即先验分布的采样）和同时依靠先验和数据似然值信息获得的采样（即后验分布的采样），这两种采样将会被用于计算贝叶斯因子（详见3.5部分）。具体代码见图16。

图16 模型4 代码及结果

完成模型拟合后，我们可以借助后验预测检验（posterior predictive check）检查模型对实际观测数据的“描述充分性”（descriptive adequacy,Gelman et al.,2013）。与先验预测检验（3.3小节）相似，后验预测检验也可以根据参数值生成对应的模拟数据，区别在于后验预测检验模拟过程使用的参数来自参数后验分布。我们同样可以使用brms的pp_check（）函数来绘制后验预测检验的结果（图17）。

图17 后验预测检验

图18展示了观测数据的均值（黑色线条为中性图片和积极图片类型LPP的均值）和依据模型生成的预测数据（即后验预测分布）的均值的分布（浅灰色为模型预测的数据均值）。如图所示，大部分后验预测数据的均值与观测数据的均值接近，说明构建的模型较好地描述了观测数据。如果观测数据的统计量和模拟生成的预测数据的统计量间存在较大的差异，研究者则需要慎重考虑构建的模型是否忽略了某些重要因素，如未使用合适的残差分布或链接函数等（Schad et al.，2021）。

图18 后验预测分布图（黑色线条T（y）为观测数据的均值；灰色部分T（yrep）为后验预测数据均值的分布

3.5 统计推断：贝叶斯因子

在实际应用中，研究者通常感兴趣的是实验条件之间是否存在差异，例如本研究中的组别（组间因素）和图片类型（被试内因素）是否会对LPP产生影响。在贝叶斯分析框架下，我们可以基于后验分布进行贝叶斯推断和假设检验。进行假设检验的常见方式之一是计算贝叶斯因子。贝叶斯因子可以量化基于不同的假设或模型（例如零假设和备择假设）观测到当前数据的似然比值（Rouder et al.,2018），研究者可以根据其数值大小推断观测到的数据更可能支持哪种假设或模型。一种流行的计算贝叶斯因子的方法是 Savage-Dickey density ratio，该方法通过计算在参数取值为0时两个模型的似然值的比值来估算贝叶斯因子，详见Wagenmakers等（2010）。通过brms提供的hypothesis()函数可以实现该贝叶斯因子的计算。

首先，利用model_4检验图片类型是否会对LPP产生影响，代码见图19。

图19 组内比较

其中，"type2M1=0"表示图片类型的零假设，即LPP在积极图片与中性图片上的差异为0。hypothesis $ Evid.Ratio描述了贝叶斯因子BF01，表示基于零假设（H0）观测到当前数据的似然与基于备择假设（H1）观测到当前数据的似然的比值。为了依据BF01进行统计推断，研究者通常需要预先设定贝叶斯因子的判断标准（胡传鹏等,2018）。根据贝叶斯因子判定标准(Lee&Wagenmakers,2014)，可以将BF01支持零假设的证据强调分为十类（表2）。上述结果显示，BF01=1.63*10-15，远远小于1/100，说明存在极强的证据支持图片类型的效应存在的备择假设，因此我们推断LPP的波幅在不同图片类型之间存在差异。

表2 贝叶斯因子决策标准

其次，我们检验LPP在不同组别间的差异是否为0，代码见图20。

图20 组间比较

"group2M1=0"表示组别（group）的零假设，即抑郁组与对照组被试在LPP上的表现差异为0。结果显示，BF01=8.90 >3，这表明存在中等强度的证据支持组别的零假设，从而我们推断组别之间的LPP波幅不存在差异。

同时，我们还可以检验交互作用是否为0，代码见图21。

图21 交互项比较

"group2M1:type2M1=0"表示组别和图片类型交互作用的零假设，即抑郁组中积极图片和中性图片LPP波幅差异与对照组中积极图片和中性图片LPP波幅差异的差异（即差的差）为0。结果显示，BF01=11.75>10，这表明存在较强的证据支持不存在组别和图片类型的交互作用。

最后，我们可以根据model_4参数分析的结果（固定效应）和贝叶斯因子分析的结果报告：“贝叶斯混合效应模型的结果显示，存在极强的证据表明，积极图片诱发的LPP波幅高于中性图片，β图片类型=6.48，BF01=1.63*10-15；同时，存在较强的证据表明，两组的LPP波幅不存在显著差异，β组别=0.27，BF01=8.90；且存在中等的证据表明图片类型和组别不存在交互作用，β图片类型，组别=0，BF01=11.75。”

4 小结

本文首先介绍了贝叶斯混合效应模型在处理数据层级结构和提供更直观的统计结果等方面的独特优势，并概述了其基本原理。更重要的是结合模拟数据和代码示例，介绍了固定效应和随机效应，如何构建贝叶斯混合效应模型，如何利用先验预测检验设置合理的先验分布，以及如何借助贝叶斯因子进行统计推断。

本文主要关注了贝叶斯线性混合效应模型，它是其他扩展模型的基础。一类常用的扩展模型是贝叶斯广义线性混合效应模型（Bayesian Generalized Linear Mixed-effects Models），它能够更好地处理心理学实验中其他常见的因变量，比如二元选择、反应时、阅读时间、李克特评分等（Dixon,2008;Rousselet &Wilcox,2020）。具体来说，研究者可以在贝叶斯广义线性混合效应模型中为因变量设定更合理的残差分布（residual distribution）和链接函数（link function），进而获得更准确的估计结果（Vuorre,2017）。例如，对于二项选择（如正确标记为1，错误为0）的因变量，brm()函数可以通过参数family=bernoulli("logit")设置模型的残差分布为伯努利分布且使用logit链接函数（Bürkner,2017）。再比如对于反应时，研究者可以通过family=“lognormal” 为其设置log转换（Wagenmakers &Brown,2007）。此外，贝叶斯分层模型还可以应用于复杂的认知模型，比如信号检测论、漂移扩散模型、强化学习模型和多项式加工树模型等（Ahn et al.,2017；Schad et al.,2020;Tso et al.,2021）。这些分析可以通过R语言的工具包（如hBayesDM）、Stan和JAGS等其他工具来完成（Ahn et al.,2017;Kruschke,2014）。总的来说，贝叶斯混合效应模型凭借其强大的灵活性，能够很好地适用于多样的心理学实验研究中。

贝叶斯混合效应模型的诸多优点使其逐渐流行于心理学研究的数据分析之中（Bono et al.,2021）。但研究者在实践中使用这种方法时仍会面临一些困难。比如，在初学时可能不清楚如何设置合理的先验，因而会采用工具包提供的默认先验。需要注意的是，这些默认先验并不一定适用于所有的研究，而使用不合适的先验可能造成模型难以收敛等潜在问题（Veenman et al.,2022）。如3.3小节所介绍，研究者可以根据以往的研究和自己的知识自定义先验，并且可以借助先验预测检验来审查这些先验的合理性。其次，当使用贝叶斯因子进行假设检验时，研究者一般无法通过单个模型获得所有感兴趣的效应。例如model_4结果只包含了图片类型和组别的主效应，以及他们的交互作用。相应地，我们也只能计算这三个效应的贝叶斯因子。如果研究者同时对其他效应（如简单效应）感兴趣，则需要构建新的模型并设置不同的对比编码来检验不同的假设（Schad et al.,2020）。最后，拟合贝叶斯混合效应模型相对更加耗时。随着数据量的变大、模型复杂程度的提升、以及采样迭代数量的增长，贝叶斯混合效应模型的拟合时间也会相应地呈指数性增长。对此，研究者可以通过设置合理的先验和使用高性能的计算机提升模型的拟合速度。

本文借助心理学研究示例，通过具体代码向读者展示了如何使用贝叶斯混合效应模型进行数据分析。随着开放科学的兴起（胡传鹏等,2016;Nosek,2019;Open Science Collaboration,2015;Jin et al.,2023），越来越多的文章公开其数据和代码，这无疑为研究者提供了更多结合具体实例的贝叶斯混合效应模型的应用教程。此外，我们也鼓励希望对该方法有更加深入和全面了解的读者阅读贝叶斯分析相关的优秀书籍（Gelman et al.,2013；Kruschke，2014；Lambert，2018；McElreath，2020）。希望研究者可以对贝叶斯混合效应模型有更加深入的了解，能够广泛地将其应用于多样的心理学研究之中。