貝葉斯混合效應(yīng)模型：基于brms 的應(yīng)用教程

潘晚坷 溫秀娟 金海洋

（1 南京師范大學(xué)心理學(xué)院，南京 210097）

（2 廣州醫(yī)科大學(xué)附屬腦科醫(yī)院，廣州 510370）

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

1 引言

近年來，心理學(xué)研究所使用的統(tǒng)計方法日漸豐富。隨著新的統(tǒng)計方法的推廣與普及，心理學(xué)研究者逐漸意識到了一些傳統(tǒng)統(tǒng)計方法的局限性，實驗心理學(xué)中最常用的重復(fù)測量方差分析就是其中之一。首先，它處理的是匯總后（aggregated）的數(shù)據(jù)，比如每個被試在每個條件下的平均反應(yīng)時。匯總后的數(shù)據(jù)忽略了每個被試個體內(nèi)部的變異，這可能會導(dǎo)致無法準(zhǔn)確估計自變量的效應(yīng)（Haaf &Rouder,2017;Veenman et al.,2022）。其次，它通常只能考慮一個隨機因素（Judd et al.,2017），而無法同時考慮多個隨機因素。例如，它一般把被試作為隨機因素，借此通過被試樣本的數(shù)據(jù)推斷被試總體的效應(yīng)。類似地，重復(fù)測量方差分析可以把刺激材料作為隨機因素，通過實驗中使用的刺激材料推斷刺激材料總體的效應(yīng)，但它無法同時考慮被試和刺激材料這兩個隨機因素。最后，當(dāng)某個被試在某一條件下的匯總數(shù)據(jù)缺失時（例如在縱向研究中某些被試未能參加某個時間點的測試），通常只能刪除存在缺失值的被試的所有數(shù)據(jù)，這會造成數(shù)據(jù)進一步流失，可能導(dǎo)致有偏差的估計（Magezi,2015;Matuschek et al.,2017）。

相比之下，線性混合效應(yīng)模型（Linear Mixedeffects Models，LMM），亦稱多水平模型（Multilevel Models）或?qū)蛹壞Ｐ停℉ierarchical Models），可以更好地處理這些問題。混合效應(yīng)模型可以考慮被試在各個條件內(nèi)的變化性及數(shù)據(jù)中可能存在的相互依存關(guān)系和層級結(jié)構(gòu)，并且能夠更好地處理數(shù)據(jù)缺失等問題（Singmann &Kellen,2019;Tuerlinckx et al.,2006）。例如，在以學(xué)生為被試的實驗中，學(xué)生來自不同的班級，不同的班級屬于不同的學(xué)校；同一個班級或同一個學(xué)校學(xué)生的表現(xiàn)一般更相似。混合效應(yīng)模型可以更好地考慮學(xué)生、班級和學(xué)校之間的層級關(guān)系，對實驗效應(yīng)的估計更準(zhǔn)確。統(tǒng)計工具（如 R 語言）的成熟使得在心理學(xué)研究中使用混合效應(yīng)模型分析數(shù)據(jù)變得更加容易，越來越多的研究也開始使用該模型（Bono et al.,2021）。

與（頻率論的，F(xiàn)requentist）線性混合效應(yīng)模型相比，貝葉斯線性混合效應(yīng)模型（Bayesian Linear Mixed-effects Models，BLMM）因 其在基于貝葉斯統(tǒng)計推斷（即依照貝葉斯定理，通過結(jié)合先驗信息（prior）與模型似然性（likelihood）對參數(shù)的后驗（posterior）進行統(tǒng)計推斷）以及處理復(fù)雜模型等方面的優(yōu)勢而受到了越來越多的關(guān)注（Bürkner,2017;Sorensen &Vasishth,2016）。首先，與頻率論統(tǒng)計相比，貝葉斯分析的一大優(yōu)勢是可以量化參數(shù)估計的不確定性，即它的結(jié)果是一個分布（例如后驗概率密度分布），其寬窄描述了對參數(shù)分析的不確定性。其次，貝葉斯分析的可信區(qū)間（credible interval）可以理解為一個擁有一定概率包含真值（true value）的區(qū)間，比如95%的可信區(qū)間表示有95%的可能性包含真值（Morey et al.,2016）。但該理解通常被誤認(rèn)為是頻率論置信區(qū)間（confidence interval）的屬性（Morey et al.,2016）（95%置信區(qū)間表示在通過無數(shù)樣本所獲得的無數(shù)個區(qū)間中，大約有95%的區(qū)間包含了真值（Lambert,2018））。其次，基于頻率論的統(tǒng)計方法一般不能為支持虛無假設(shè)提供證據(jù)（Cohen,1994），而貝葉斯分析中的貝葉斯因子（Bayes Factor，BF）或?qū)嶋H等價區(qū)域（Region of Practical Equivalence，ROPE）可以（Makowski et al.,2019）。最后，如果混合效應(yīng)模型試圖包含復(fù)雜的隨機效應(yīng)結(jié)構(gòu)，基于頻率論的模型通常難以收斂，從而無法給出可靠的估計（Matuschek et al.,2017），而貝葉斯模型能夠更好地處理該結(jié)構(gòu)（Sorensen &Vasishth,2016）。此外，貝葉斯分析還可以通過整合研究者的先驗知識提供更合理的參數(shù)估計。

貝葉斯混合效應(yīng)模型兼具（頻率論）混合效應(yīng)模型和貝葉斯分析的優(yōu)點，是心理學(xué)研究中極具潛力和值得推廣的統(tǒng)計方法。本文致力于向讀者提供一個相對容易理解和上手的教程。

2 模型介紹

2.1 混合效應(yīng)模型

混合效應(yīng)模型（Baayen et al.,2008）是一般線性模型的擴展形式，得名于它同時包含了固定效應(yīng)（fixed effects）和隨機效應(yīng)（random effects）。固定效應(yīng)通常是研究者所感興趣的效應(yīng)，即自變量在總體水平的效應(yīng)。與之相對的，隨機效應(yīng)是自變量在隨機因素（例如被試或刺激材料）各個水平的特有效應(yīng)。例如，在探究圖片類型對情緒喚醒的影響的實驗中，每個被試受到圖片類型的影響可能并不完全相同。其中，圖片類型對所有被試的平均影響為固定效應(yīng)，而每個被試受到圖片類型的影響與平均影響的差值為隨機效應(yīng)。

為更好地解釋固定效應(yīng)與隨機效應(yīng)，我們首先回顧最基本的一般線性模型：

其中，Yi,j是因變量，比如被試的腦電波幅（詳見3 應(yīng)用示例），表示第j名被試在第i試次的觀測數(shù)據(jù)。Xi,j為自變量，表示研究者感興趣的實驗操作，可以是組間或者組內(nèi)變量。參數(shù)β0是截距，表示所有被試的均值，比如所有被試的平均腦電波幅；參數(shù)β1為斜率，表示自變量對因變量的影響，比如不同圖片類型對被試腦電波幅的影響（本文中涉及截距和斜率的解釋僅適用于應(yīng)用逐次差分對比編碼（Successive Differences Contrast Coding，即本文所使用的編碼）且只有兩個水平的自變量（Venables &Ripley,2002）。當(dāng)使用不同的編碼時，截距和斜率的解釋也會有所不同，詳見 Schad et al.,2020）。；∈i,j為殘差，表示模型中不能被解釋的部分。

需要注意的是，使用一般線性模型需要滿足殘差獨立性和正態(tài)性的前提預(yù)設(shè)，即不同觀測數(shù)據(jù)的殘差應(yīng)該相互獨立并且服從正態(tài)分布。然而，在被試內(nèi)實驗設(shè)計中，每個被試完成不同圖片類型的多個試次，因此他們在不同試次上的表現(xiàn)并不是完全獨立的。此時使用一般線性模型無法處理觀測數(shù)據(jù)間的依存性。此外，公式（1）的模型沒有考慮不同被試間的個體差異，比如，不同被試的平均腦電波幅（即每個被試的截距）可能存在差異，而忽略這些差異可能會掩蓋圖片類型的真實效應(yīng)。為更好地考慮這些方面，可以在公式1的基礎(chǔ)上增加隨機截距：

與公式（1）相比，公式（2）主要區(qū)分了固定截距（β0）和隨機截距（S0,j）。其中，固定截距可以理解為所有被試截距的平均值，即所有被試的平均腦電波幅；而隨機截距可以理解為每個被試的截距與固定截距的差值，即每個被試各個條件的平均腦電波幅與所有被試平均腦電波幅的差值。雖然加入隨機截距可以考慮到各被試平均腦電波幅的差異，但仍然沒有考慮到自變量對因變量的影響在不同被試上的差異。比如，不同被試受到圖片類型的影響并不一定相同。為更好地解釋自變量的影響在不同被試間的差異，可以進一步在混合模型中加入隨機斜率：

與公式（2）相比，公式（3）主要區(qū)分了固定斜率（β1）和隨機斜率（S1,j）。其中，固定斜率是所有被試斜率的平均值，即所有被試受到圖片類型的平均影響；隨機斜率是每個被試斜率與平均斜率的差值，即每個被試受到圖片類型的特有影響。需要注意的是，由于隨機截距與隨機斜率都受到同一個隨機因素（被試）的影響，因此兩者可能存在一定的相關(guān)性，比如對于平均腦電波幅（隨機截距）更高的被試，其在不同圖片類型下的腦電波幅的差異（隨機斜率）可能更大。此時通過預(yù)設(shè)兩種隨機效應(yīng)來自同一個多元正態(tài)分布可以更好地考慮他們之間的相關(guān)性。

總的來說，混合效應(yīng)模型是一般線性模型的擴展，它能夠借助隨機截距和隨機斜率考慮隨機因素（如被試和實驗材料）對實驗效應(yīng)的影響，可以更好地處理數(shù)據(jù)中的依存和層級關(guān)系，彌補重復(fù)測量方差分析的不足，提高對實驗效應(yīng)估計的準(zhǔn)確性。

2.2 貝葉斯混合效應(yīng)模型

在構(gòu)建混合效應(yīng)模型的過程中，為了在控制合理的一類錯誤率的同時保持良好的統(tǒng)計效力，研究者通常需要考慮最大化的模型（maximal model）（Barr et al.,2013），即把所有可能的隨機截距和隨機斜率以及它們之間的相關(guān)系數(shù)加入模型中。然而，由于模型參數(shù)增多，可能導(dǎo)致模型難以收斂，從而給出不可靠的估計（Bates et al.,2015）。相比之下，貝葉斯混合效應(yīng)模型可以更好地處理這種復(fù)雜的模型。更重要的是，相較于頻率論混合效應(yīng)模型，貝葉斯混合效應(yīng)模型具有考慮先驗信息和提供更直觀的參數(shù)后驗分布等優(yōu)勢。這些優(yōu)點使得貝葉斯混合效應(yīng)模型更受歡迎。

貝葉斯定理（Bayes’ theorem）在貝葉斯數(shù)據(jù)分析中起著重要作用，它借助數(shù)據(jù)更新先驗信息，從而得到后驗信息：

其中，p（θ|Y）為參數(shù)的后驗分布（posterior distribution），即連續(xù)型參數(shù)對應(yīng)的概率密度函數(shù)或離散型參數(shù)對應(yīng)的概率質(zhì)量函數(shù)，表示在觀測到數(shù)據(jù)Y的基礎(chǔ)上參數(shù)θ的概率分布。對于貝葉斯混合效應(yīng)模型來說，θ包括所有的固定效應(yīng)β、隨機效應(yīng)S、隨機效應(yīng)間的相關(guān)性和殘差∈（見公式3）所對應(yīng)的分布的所有參數(shù)。p（θ|Y）可以由公式右側(cè)部分計算得來，其中p（θ）為參數(shù)的先驗分布（prior distribution），表示研究者在觀測到數(shù)據(jù)前對參數(shù)θ的先驗知識。p（Y|θ）為似然，描述了當(dāng)參數(shù)θ為不同的值時，觀測到當(dāng)前數(shù)據(jù)的似然值。值得注意的是，p（Y）并不包含θ，這是因為它是一個邊際似然（marginal likelihood），是通過對所有可能的θ的似然值與相應(yīng)的先驗概率乘積求和或積分而獲得的（當(dāng)θ為離散變量時：p（Y）=∑θp（Y|θ）*p（θ）；當(dāng)θ為連續(xù)變量時，p（Y）=∫θp（Y|θ）*p（θ）dθ，（Kruschke &Liddell,2018）。p（Y）是一個“歸一化”因子（normalizing factor），能夠確保后驗分布是一個有效的概率分布（使后驗分布的面積為1）。

基于貝葉斯定理，在對貝葉斯混合效應(yīng)模型進行參數(shù)估計時，研究者首先需要為各參數(shù)設(shè)置合理的先驗分布，之后通過數(shù)據(jù)對先驗進行更新，從而獲得參數(shù)的后驗分布。需要注意的是，隨著模型變得復(fù)雜，貝葉斯定理中的邊際似然（公式4右側(cè)分母部分）將變得難以計算。此時通常通過馬爾科夫鏈蒙特卡洛（Markov Chain Monte Carlo,MCMC,Lambert,2018;Morey et al.,2011）算法推測參數(shù)的后驗分布（Bürkner,2017），其具體原理可參考Plummer（2023）。與MCMC相關(guān)的重要參數(shù)會在下文中結(jié)合示例予以介紹。

3 應(yīng)用示例

我們將借助一個假想的心理學(xué)實驗展示如何使用brms工具包進行貝葉斯混合效應(yīng)模型的數(shù)據(jù)分析。具體來說，本節(jié)將從數(shù)據(jù)模擬和預(yù)處理、模型構(gòu)建、先驗設(shè)定、模型的擬合與檢驗和統(tǒng)計推斷五個方面介紹貝葉斯分析最基本的工作流程（Gelman et al.,2020）。具體代碼見https://osf.io/cgxwn/。

3.1 數(shù)據(jù)模擬與預(yù)處理

該假想研究的目的是探究抑郁癥人群加工不同情緒刺激的神經(jīng)基礎(chǔ)。實驗中，40名抑郁癥患者和40名健康對照組被試觀看30張積極圖片和30張中性圖片，期間我們采集了他們的腦電數(shù)據(jù)。因變量是晚期正電位（late positive potentials,LPP）的波幅。該實驗為2（組別group：抑郁癥組depression、對照組 control）×2（圖片類型type：積極positive、中性neutral）混合實驗設(shè)計，其中組別為被試間因素，圖片類型為被試內(nèi)因素。

實驗數(shù)據(jù)使用R的faux工具包生成（DeBruine,2021），具體代碼見線上補充材料數(shù)據(jù)模擬部分（https://osf.io/3g79b）。簡單來說，設(shè)定（圖1和表1）圖片類型的主效應(yīng)為6.5（即積極圖片誘發(fā)的LPP比中性圖片高6.5），組別的主效應(yīng)為0.1（即對照組LPP比抑郁組高0.1）和組別與圖片類型的交互作用為0.1。

表1 模擬數(shù)據(jù)表述統(tǒng)計結(jié)果

圖1 模擬數(shù)據(jù)柱狀圖

首先需要確保數(shù)據(jù)集的結(jié)構(gòu)符合工具包的要求。重復(fù)測量方差分析通常需要“寬”格式數(shù)據(jù)，即每一行是一個被試的數(shù)據(jù)。而混合效應(yīng)模型的分析通常需要 “長”格式數(shù)據(jù)，即每一行是一個試次的數(shù)據(jù)。例如，R語言中的brms工具包需要使用長數(shù)據(jù)格式進行混合效應(yīng)分析。與大部分實驗程序（如E-prime）生成的數(shù)據(jù)類似，我們通過faux工具包生成的模擬數(shù)據(jù)也是以長格式保存，符合brms工具包輸入數(shù)據(jù)的要求。該模擬數(shù)據(jù)集的變量名為df_simu，我們可以通過head(df_simu,10)查看前10行數(shù)據(jù)（圖2）。

此外，我們還需要對數(shù)據(jù)集df_simu中的稱名（自）變量進行編碼。這里我們使用逐次差分對比編碼對組別（group）和圖片類型（type）這兩個自變量進行編碼（圖3），即使用R中的MASS::contr.sdif()函數(shù)。通過此設(shè)定，貝葉斯模型的輸出結(jié)果剛好對應(yīng)大部分研究者熟悉的主效應(yīng)和交互作用（見公式1處解釋）。

3.2 構(gòu)建模型

完成數(shù)據(jù)的預(yù)處理后，研究者需要構(gòu)建模型。模型構(gòu)建主要有兩種方法。一種是從包含了研究者最感興趣效應(yīng)的簡單模型開始，逐漸添加更多的效應(yīng)以尋找能夠最優(yōu)解釋數(shù)據(jù)的模型。另一種方法是根據(jù)實驗設(shè)計構(gòu)建最大化的模型（Barr et al.,2013）。鑒于第二種方法相對簡單，且更多地被認(rèn)知和心理學(xué)研究者所采用（Schad et al.,2021），本文將主要介紹如何根據(jù)實驗設(shè)計構(gòu)建最大化模型（對第一種途徑感興趣的讀者可參考Schad et al.,2021）。

在實際分析中建立貝葉斯線性混合效應(yīng)模型時，可以直接根據(jù)實驗設(shè)計設(shè)定最大化的模型。需要注意的是，本文為了更清晰地介紹如何理解和設(shè)定固定效應(yīng)和隨機效應(yīng)，將以三個模型為例，分別解釋如何結(jié)合實驗設(shè)計使用brms設(shè)定固定效應(yīng)、隨機截距和隨機斜率。

首先建立一個僅包含固定效應(yīng)的模型（model_1，見圖4）。具體來說，我們將組別和圖片類型兩個自變量以及他們的交互作用作為固定效應(yīng)。該模型假定所有被試的平均波幅相同，且圖片類型對每個被試的影響也相同。

圖4 固定效應(yīng)模型

在圖4的代碼中，LPP為因變量，即被試觀看不同圖片時LPP的波幅，group*type表示組別和圖片類型兩個自變量的主效應(yīng)以及他們的交互作用。df_simu為數(shù)據(jù)集。其他參數(shù)是模型擬合時所基于的MCMC采樣算法的相關(guān)參數(shù)。為了保證后驗分布的有效性，MCMC通常會同時運行多個獨立的采樣過程，每個獨立的采樣過程稱為一條鏈（chain）。只有當(dāng)各個獨立鏈的采樣“融合”在一起時（見圖6），MCMC獲得的參數(shù)后驗分布才更可能是有效的（McElreath,2020）。在model_1中，我們設(shè)置了4條鏈（chains=4，默認(rèn)為4）。與之相對應(yīng)的是cores參數(shù)，它可以設(shè)定運行MCMC時使用的計算機中央處理器核心（Central Processing Unit cores）的個數(shù)。設(shè)置cores=4（默認(rèn)為1）表示使用四個中央處理器核心同時分別各運行一條馬爾可夫鏈（共四條鏈），借此提高數(shù)據(jù)分析的速度。iter參數(shù)（iteration的簡寫）指的是每條鏈的迭代次數(shù)，即對每個參數(shù)進行多少次的采樣。warmup參數(shù)（或稱為burn in）設(shè)定了在每條鏈采樣的開始階段用于“校準(zhǔn)”采樣器（sampler）以提升采樣效率的迭代次數(shù)（詳見https://mc-stan.org/docs/reference-manual/hmc-algorithm-parameters.html），該階段的采樣不會被納入最終的參數(shù)概率分布。thin參數(shù)為“稀釋比例（thinning rate，默認(rèn)為1）”，表示每隔thin個采樣保留一個樣本，研究者可以通過設(shè)定較大的thin來減少電腦內(nèi)存的需求（詳見https：//mc-stan.org/docs/referencemanual/effective-sample-size.html#thinningsamples）。總的來說，每條鏈最后保留的采樣數(shù)量與iter、warmup和thin相關(guān)，最后保留的采樣數(shù)等于（iter-warmup）/thin。例如，在model_1中每一條鏈會保留每個參數(shù)的4000個樣本，即（5000-1000）/1，最終四條鏈一共會保留16000個樣本。

通過summary()函數(shù)可以展示模型擬合的相關(guān)參數(shù)和結(jié)果，詳見圖5。

輸出結(jié)果的第一部分展示了模型設(shè)置的信息。例如，“Family：gaussian，Links：mu=identity；sigma=identity”表示模型預(yù)設(shè)自變量和因變量間的關(guān)系是線性的，且殘差是一個正態(tài)分布。對于其他非線性關(guān)系的模型，研究者可以設(shè)定相應(yīng)的殘差分布和鏈接函數(shù)（詳見第4部分）。

Population-Level Effects部分為固定效應(yīng)的結(jié)果，通常是研究者最為關(guān)注的部分。其中，Intercept是模型的固定截距，表示所有被試的平均LPP波幅。group2M1表示在平均圖片類型的效應(yīng)之后，健康對照組的LPP波幅減去抑郁癥組LPP波幅的差值，即組別的主效應(yīng)。與此類似，type2M1是在平均組別的效應(yīng)之后，積極圖片誘發(fā)的LPP波幅減去中性圖片誘發(fā)的LPP波幅的差值，即圖片類型的主效應(yīng)。group2M1：type2M1是組別和圖片類型交互作用，即（LPP健康對照組，積極圖片-LPP健康對照組，中性圖片）-（LPP抑郁癥組，積極圖片-LPP抑郁癥組，中性圖片）。

在輸出的結(jié)果中，Estimate表示參數(shù)的后驗分布的均值，Est.Error表示后驗分布的標(biāo)準(zhǔn)差。l-95%CI和u-95%CI表示參數(shù)的95%可信區(qū)間的兩個邊界，例如group2M1的95%的可信區(qū)間是［0.11,0.46］。Rhat表示不同鏈條采樣的融合程度。Rhat越接近1，表示不同鏈條采樣融合得越好，即模型的收斂情況越好。如果Rhat大于1.1（或更嚴(yán)格的標(biāo)準(zhǔn)下，如大于1.01），則說明模型的收斂情況并不理想（Gelman &Rubin,1992;Vehtari et al.,2021）。此時，研究者不應(yīng)該使用獲得的后驗分布進行進一步分析，而應(yīng)該檢查模型設(shè)定和數(shù)據(jù)是否存在問題，并進行相應(yīng)調(diào)整。ESS為有效樣本量，是對MCMC算法采樣到的有效且獨立的樣本數(shù)量的估計。ESS越高表明能用來有效估計后驗分布的樣本數(shù)量越大。通常建議Bulk-ESS應(yīng)是MCMC鏈數(shù)量的100倍及以上。例如，當(dāng)使用四條MCMC鏈時，Bulk-ESS應(yīng)至少為400（Vehtari et al.,2021）。為了確保足夠的Bulk-ESS，研究者可以適當(dāng)增加iter，例如本文使用的iter為5000。Tail-ESS表示后驗分布兩端的有效樣本量（effective sample size）。通常情況下，過小的Bulk-ESS和Tail-ESS會伴隨著過大的Rhat。設(shè)定合理的模型和先驗分布一般可以解決大部分Rhat過大的情況。最后，研究者還可以通過可視化的方式檢查各條鏈的融合情況（圖6）。

圖6 model_1 的MCMC 各條鏈采樣軌跡圖（又稱“毛毛蟲”圖）

圖6描述了model_1中MCMC不同鏈抽樣的軌跡圖，反映了不同鏈的融合情況。不同鏈參數(shù)的抽樣很好地 “融合”在一起，說明其對應(yīng)的參數(shù)后驗分布相對可靠。由于融合的鏈看起來像一條毛毛蟲，因此軌跡圖又稱“毛毛蟲”圖。如果一條鏈的抽樣與其他鏈的抽樣存在一定差異（通常伴隨著較大的Rhat），那這些后驗抽樣也很可能無法有效表征參數(shù)的后驗分布（詳見線上補充材料第三部分MCMC鏈不收斂示例，https：//osf.io/3g79b），此時，研究者需要重新考慮模型的合理性（Roy,2020）。

如前所述，僅包含固定效應(yīng)的模型（model_1）假定所有被試的平均腦電波幅都是相同的，即所有被試共享同一個截距值。但這種情況在現(xiàn)實中不太可能發(fā)生。更可能的情況是，每個被試的平均腦電波幅各不相同。因此，我們需要在模型中增加隨機截距來考慮被試在平均腦電波幅上的個體差異。具體來說，我們需要在model_1的基礎(chǔ)上增加基于被試（by-subject）的隨機截距（對應(yīng)公式2中的S0,i），即(1|subj)。其中subj為被試編號，1表示隨機截距，詳見圖7。

Population-Level Effects為model_2固定效應(yīng)的結(jié)果，其相關(guān)參數(shù)的含義與model_1中的相同（詳見3.2 部分）。Group-Level Effects為隨機效應(yīng)結(jié)構(gòu)的結(jié)果，其中sd(Intercept)為隨機截距的標(biāo)準(zhǔn)差，反映了LPP平均波幅在不同被試間的離散程度。雖然model_2考慮了不同被試平均LPP波幅的個體差異，但它仍然假定所有被試受到實驗處理（即圖片類型）的影響是相同的。然而在實際情況中，每個被試受到圖片類型的影響更可能是不同的。因此，我們需要設(shè)置隨機斜率以考慮圖片類型對不同被試的不同影響。具體來說，我們可以在model_2模型的基礎(chǔ)上加入隨機斜率（對應(yīng)公式3中的S1,i，即 (type|subj)，其中subj為被試編號，type為圖片類型的隨機斜率。與此同時，我們也在模型中加入了隨機截距與隨機斜率間的相關(guān)，詳見圖8。

圖8 模型3 代碼及結(jié)果

與model_2結(jié)果相比，model_3中的Group-Level Effects部分額外展示了基于被試的隨機斜率和它與隨機截距相關(guān)性的后驗分布情況。其中，sd(type2M1)為基于被試的圖片類型隨機斜率的標(biāo)準(zhǔn)差，反映了不同被試LPP波幅受到圖片類型不同影響的離散程度。cor(Intercept,type2M1)為隨機截距和隨機斜率的相關(guān)，該相關(guān)結(jié)果表明兩者之間存在一定的正相關(guān)，即平均LPP波幅更大的被試，其對積極圖片的LPP波幅也比中性圖片更大。

至此，我們完成了對該模擬數(shù)據(jù)的最大化模型的構(gòu)建，即根據(jù)實驗設(shè)計設(shè)定了所有合理的固定效應(yīng)和隨機效應(yīng)結(jié)構(gòu)。在該示例中，我們從僅包含固定效應(yīng)的線性模型出發(fā)，在不同的模型中依次展示了如何添加隨機截距和隨機斜率。需要強調(diào)的是，在實際的貝葉斯數(shù)據(jù)分析中，研究者可以根據(jù)實驗設(shè)計直接設(shè)定最大化的模型，而無需在不同的模型中依次添加不同的隨機效應(yīng)。

如前所述，貝葉斯混合效應(yīng)模型還可以同時考慮多個隨機因素，以及不同隨機因素間的層級關(guān)系。例如，在分析包含“學(xué)生—班級—學(xué)?！睂蛹夑P(guān)系的數(shù)據(jù)時，研究者可以通過設(shè)定“（1 ｜學(xué)校班級學(xué)生）”的隨機效應(yīng)結(jié)構(gòu)來考慮同一個班級和同一個學(xué)校學(xué)生間的相似性等（詳見Bürkner,2018）。

3.3 先驗的設(shè)置

完成構(gòu)建模型后，研究者需要根據(jù)當(dāng)前的先驗知識自定義模型的先驗，以在數(shù)據(jù)分析中整合先驗知識。由于先前沒有自定義先驗分布（例如model_3），因此模型使用的均是brms默認(rèn)的。這些先驗分布通常是弱信息先驗（weakly informative prior），即相對實驗數(shù)據(jù)來講，基本不能為參數(shù)估計提供有效信息的先驗（Seaman et al.,2012）。使用這些默認(rèn)的弱信息先驗非但不能發(fā)揮貝葉斯數(shù)據(jù)分析的最大作用，反而可能會對數(shù)據(jù)分析的結(jié)果產(chǎn)生負(fù)面影響（Seaman et al.,2012）。例如，貝葉斯因子很容易受到先驗分布的影響，不合理的先驗會導(dǎo)致有偏差的貝葉斯因子。因此，為貝葉斯分析設(shè)定合理的先驗至關(guān)重要。

設(shè)定先驗分布的第一步是了解可以對一個模型的哪些參數(shù)設(shè)置先驗。這里可以借助get_prior()函數(shù)得到相應(yīng)的信息，詳見圖9。

圖9 先驗設(shè)置

如結(jié)果所示，我們可以為model_3模型的以下參數(shù)設(shè)定先驗：固定截距（class為Intercept）；固定斜率（class為b），包含組別主效應(yīng)、圖片類型主效應(yīng)和他們的交互作用（分別對應(yīng)coef為group2M1，type2M1和group2M1:type2M1）；基于被試（group為subj）的隨機效應(yīng)的標(biāo)準(zhǔn)差（class為sd），包含隨機截距（coef為Intercept）和隨機斜率（coef為type2M1），基于被試的隨機效應(yīng)間的相關(guān)（class為cor），以及殘差（class為sigma）。由于sd和sigma不能為負(fù)數(shù)，因此這些參數(shù)的最小值為0（即列l(wèi)b，為low bound的縮寫）。

在實際數(shù)據(jù)分析中，研究者一般可以根據(jù)以往研究的結(jié)果為當(dāng)前研究的模型設(shè)定先驗分布。此時，研究者可以通過先驗預(yù)測檢驗（prior predictive check）來測試這些先驗的合理性，并進行相應(yīng)調(diào)整，進而為模型確定合理的先驗。在進行先驗預(yù)測檢驗時，需要考慮領(lǐng)域知識，即研究者對研究領(lǐng)域的專業(yè)知識和經(jīng)驗，它可以幫助研究者確定先驗分布的合理范圍，并檢查先驗分布是否與領(lǐng)域知識一致。如果先驗分布與領(lǐng)域知識相矛盾，則需要重新考慮先驗分布的選擇。

如果可以從先前類似的研究獲得參數(shù)的后驗概率分布，研究者可以將這些概率分布作為當(dāng)前研究模型的先驗分布。即使無法獲得參數(shù)準(zhǔn)確的概率分布，研究者仍然可以根據(jù)先前研究結(jié)果以及已有知識為模型設(shè)定先驗分布。如果無法從先前研究中獲取參數(shù)的概率分布作為模型的先驗，我們可以首先嘗試設(shè)置弱信息的先驗分布，詳見圖10。

圖10 定義先驗分布

在圖10的代碼中，我們把固定截距和固定斜率分別設(shè)定為一個均值為0、標(biāo)準(zhǔn)差為10的正態(tài)分布（在代碼中表示為normal(0,10)）。對于prior(normal(0,10),class=Intercept)，我們可以簡單地把它理解為，我們設(shè)定所有被試的平均值腦電波幅大致有95%的概率落在-20和20之間（即兩個標(biāo)準(zhǔn)差之內(nèi)）。對于平均腦電波幅來講，這是一個相對較大的區(qū)間。因此，它對應(yīng)的先驗分布是弱信息先驗分布。類似地，所有固定斜率（包含組別和圖片類型的主效應(yīng)以及他們的交互作用）也大致落在相同的區(qū)間。此外需要注意的是，盡管使用正態(tài)分布作為殘差sigma和隨機效應(yīng)的標(biāo)準(zhǔn)差sd設(shè)置先驗，但由于這些參數(shù)不能小于0，因此他們實際對應(yīng)的先驗分布是截斷正態(tài)分布（truncated normal distribution）。最后，隨機效應(yīng)間的相關(guān)性設(shè)置為弱信息先驗prior(lkj(2),class=cor)，其中l(wèi)kj為LKJ相關(guān)性矩陣（LKJ-Correlation prior matrix，Lewandowski et al.，2009）。

為檢驗先驗分布的合理性，我們可以使用先驗預(yù)測檢驗考查僅基于先驗信息生成的因變量（即LPP的波幅）的情況。進行先驗預(yù)測檢驗需要對先驗分布進行采樣，即在模型中設(shè)定先驗分布prior=prior_01的基礎(chǔ)上加入sample_prior="only"，通過該設(shè)置可以僅從先驗分布中進行采樣（忽略觀測數(shù)據(jù)的似然值信息）。其余參數(shù)的設(shè)定和我們在分析數(shù)據(jù)時設(shè)定相同（例如model_3），具體代碼如圖11所示。

圖11 prior_01 先驗預(yù)測檢驗代碼

通過上述代碼，我們可以獲得16,000組（（5000-1000）*4）來自先驗分布的抽樣。其中，每一組都包含模型的所有參數(shù)，即固定截距，固定斜率和隨機效應(yīng)等。每一組先驗參數(shù)都可以生成一套新的與數(shù)據(jù)集（df_simu）結(jié)構(gòu)相同的模擬觀測數(shù)據(jù)（即因變量）。通過考查這些模擬觀測數(shù)據(jù)的合理性（比如檢查每一套模擬數(shù)據(jù)的極值）可以判斷先驗的設(shè)定是否合適。例如在一個以反應(yīng)時為因變量的stroop實驗中，如果那些由先驗分布生成的模擬數(shù)據(jù)集中大部分?jǐn)?shù)據(jù)集的最大值都超過預(yù)期（比如10秒），就很可能說明模型某些參數(shù)先驗的設(shè)定不太合理（比如截距或標(biāo)準(zhǔn)差過大）。在當(dāng)前模擬實驗中，盡管我們無法根據(jù)先前研究獲取精確的先驗，但是可以利用現(xiàn)有的知識，即試次水平的腦電事件相關(guān)電位的波幅通常在一定的范圍之內(nèi)，來確定先驗。因此，我們計算16,000套由先驗分布生成的模擬觀測數(shù)據(jù)各自的最小值、均值和最大值，并通過觀察它們的范圍是否符合我們的先驗知識來判斷先驗分布的合理性。

上述分析可以通過brms提供的pp_check()函數(shù)實現(xiàn)。圖12中的代碼展示了prior=prior_01的先驗預(yù)測的結(jié)果。

圖12 prior_01 先驗預(yù)測檢驗結(jié)果代碼

圖13 基于先驗prior_01 生成的16,000 套預(yù)測數(shù)據(jù)（y(pred)）的最小值、最大值和均值的分布

結(jié)果顯示，基于prior_01生成的每一套模擬觀測數(shù)據(jù)的最小值大部分落在［-75,0］，均值大部分落在［-25,25］，而最大值大部分落在［0,80］。根據(jù)先驗知識，試次水平的腦電波幅通常很少會低于-50（μV）或高于50（μV）（偽跡除外）。因此，我們可以適當(dāng)?shù)卣{(diào)整參數(shù)的先驗分布以符合我們的先驗知識，即減小參數(shù)先驗分布的不確定性。

在接下來的嘗試中，我們將固定截距和固定斜率的先驗設(shè)置成均值為0，標(biāo)準(zhǔn)差為5的正態(tài)分布，并進行與之前相同的先驗預(yù)測檢驗（見圖14）。

圖14 prior_02 先驗預(yù)測檢驗

經(jīng)過同樣的步驟，我們可以獲得基于先驗prior_02生成的16,000套模擬觀測數(shù)據(jù)（y(pred)）的最小值、均值和最大值的分布（圖15）。結(jié)果顯示模擬觀測數(shù)據(jù)的最小值大部分落在［-60,0］，均值大部分落在了［-10,10］，最大值大部分落在［0,60］。與prior_01的先驗預(yù)測結(jié)果（圖 13）相比，prior_02預(yù)測的結(jié)果更多落在了［-50,50］內(nèi)，能夠相對更合理地描述我們的先驗知識。需要注意的是，確定合理的先驗分布通常需要重復(fù)多次類似的過程，且并不是先驗分布的信息越強就越好，它需要合理反映研究者的先驗知識。

圖15 基于先驗prior_02 生成的每一套預(yù)測數(shù)據(jù)（y(pred)）的最小值、最大值和均值的分布

3.4 模型擬合

在確定先驗后就可以開始擬合模型，該模型包括所有可能的隨機截距、隨機斜率和它們之間的相關(guān)，并通過prior=prior_02為模型設(shè)定先前確定的先驗。此外，我們設(shè)定sample_prior="yes"，以此保留僅依靠先驗信息獲得的采樣（即先驗分布的采樣）和同時依靠先驗和數(shù)據(jù)似然值信息獲得的采樣（即后驗分布的采樣），這兩種采樣將會被用于計算貝葉斯因子（詳見3.5部分）。具體代碼見圖16。

圖16 模型4 代碼及結(jié)果

完成模型擬合后，我們可以借助后驗預(yù)測檢驗（posterior predictive check）檢查模型對實際觀測數(shù)據(jù)的“描述充分性”（descriptive adequacy,Gelman et al.,2013）。與先驗預(yù)測檢驗（3.3小節(jié)）相似，后驗預(yù)測檢驗也可以根據(jù)參數(shù)值生成對應(yīng)的模擬數(shù)據(jù)，區(qū)別在于后驗預(yù)測檢驗?zāi)M過程使用的參數(shù)來自參數(shù)后驗分布。我們同樣可以使用brms的pp_check（）函數(shù)來繪制后驗預(yù)測檢驗的結(jié)果（圖17）。

圖17 后驗預(yù)測檢驗

圖18展示了觀測數(shù)據(jù)的均值（黑色線條為中性圖片和積極圖片類型LPP的均值）和依據(jù)模型生成的預(yù)測數(shù)據(jù)（即后驗預(yù)測分布）的均值的分布（淺灰色為模型預(yù)測的數(shù)據(jù)均值）。如圖所示，大部分后驗預(yù)測數(shù)據(jù)的均值與觀測數(shù)據(jù)的均值接近，說明構(gòu)建的模型較好地描述了觀測數(shù)據(jù)。如果觀測數(shù)據(jù)的統(tǒng)計量和模擬生成的預(yù)測數(shù)據(jù)的統(tǒng)計量間存在較大的差異，研究者則需要慎重考慮構(gòu)建的模型是否忽略了某些重要因素，如未使用合適的殘差分布或鏈接函數(shù)等（Schad et al.，2021）。

圖18 后驗預(yù)測分布圖（黑色線條T（y）為觀測數(shù)據(jù)的均值；灰色部分T（yrep）為后驗預(yù)測數(shù)據(jù)均值的分布

3.5 統(tǒng)計推斷：貝葉斯因子

在實際應(yīng)用中，研究者通常感興趣的是實驗條件之間是否存在差異，例如本研究中的組別（組間因素）和圖片類型（被試內(nèi)因素）是否會對LPP產(chǎn)生影響。在貝葉斯分析框架下，我們可以基于后驗分布進行貝葉斯推斷和假設(shè)檢驗。進行假設(shè)檢驗的常見方式之一是計算貝葉斯因子。貝葉斯因子可以量化基于不同的假設(shè)或模型（例如零假設(shè)和備擇假設(shè)）觀測到當(dāng)前數(shù)據(jù)的似然比值（Rouder et al.,2018），研究者可以根據(jù)其數(shù)值大小推斷觀測到的數(shù)據(jù)更可能支持哪種假設(shè)或模型。一種流行的計算貝葉斯因子的方法是 Savage-Dickey density ratio，該方法通過計算在參數(shù)取值為0時兩個模型的似然值的比值來估算貝葉斯因子，詳見Wagenmakers等（2010）。通過brms提供的hypothesis()函數(shù)可以實現(xiàn)該貝葉斯因子的計算。

首先，利用model_4檢驗圖片類型是否會對LPP產(chǎn)生影響，代碼見圖19。

圖19 組內(nèi)比較

其中，"type2M1=0"表示圖片類型的零假設(shè)，即LPP在積極圖片與中性圖片上的差異為0。hypothesis $ Evid.Ratio描述了貝葉斯因子BF01，表示基于零假設(shè)（H0）觀測到當(dāng)前數(shù)據(jù)的似然與基于備擇假設(shè)（H1）觀測到當(dāng)前數(shù)據(jù)的似然的比值。為了依據(jù)BF01進行統(tǒng)計推斷，研究者通常需要預(yù)先設(shè)定貝葉斯因子的判斷標(biāo)準(zhǔn)（胡傳鵬等,2018）。根據(jù)貝葉斯因子判定標(biāo)準(zhǔn)(Lee&Wagenmakers,2014)，可以將BF01支持零假設(shè)的證據(jù)強調(diào)分為十類（表2）。上述結(jié)果顯示，BF01=1.63*10-15，遠(yuǎn)遠(yuǎn)小于1/100，說明存在極強的證據(jù)支持圖片類型的效應(yīng)存在的備擇假設(shè)，因此我們推斷LPP的波幅在不同圖片類型之間存在差異。

表2 貝葉斯因子決策標(biāo)準(zhǔn)

其次，我們檢驗LPP在不同組別間的差異是否為0，代碼見圖20。

圖20 組間比較

"group2M1=0"表示組別（group）的零假設(shè)，即抑郁組與對照組被試在LPP上的表現(xiàn)差異為0。結(jié)果顯示，BF01=8.90 >3，這表明存在中等強度的證據(jù)支持組別的零假設(shè)，從而我們推斷組別之間的LPP波幅不存在差異。

同時，我們還可以檢驗交互作用是否為0，代碼見圖21。

圖21 交互項比較

"group2M1:type2M1=0"表示組別和圖片類型交互作用的零假設(shè)，即抑郁組中積極圖片和中性圖片LPP波幅差異與對照組中積極圖片和中性圖片LPP波幅差異的差異（即差的差）為0。結(jié)果顯示，BF01=11.75>10，這表明存在較強的證據(jù)支持不存在組別和圖片類型的交互作用。

最后，我們可以根據(jù)model_4參數(shù)分析的結(jié)果（固定效應(yīng)）和貝葉斯因子分析的結(jié)果報告：“貝葉斯混合效應(yīng)模型的結(jié)果顯示，存在極強的證據(jù)表明，積極圖片誘發(fā)的LPP波幅高于中性圖片，β圖片類型=6.48，BF01=1.63*10-15；同時，存在較強的證據(jù)表明，兩組的LPP波幅不存在顯著差異，β組別=0.27，BF01=8.90；且存在中等的證據(jù)表明圖片類型和組別不存在交互作用，β圖片類型，組別=0，BF01=11.75。”

4 小結(jié)

本文首先介紹了貝葉斯混合效應(yīng)模型在處理數(shù)據(jù)層級結(jié)構(gòu)和提供更直觀的統(tǒng)計結(jié)果等方面的獨特優(yōu)勢，并概述了其基本原理。更重要的是結(jié)合模擬數(shù)據(jù)和代碼示例，介紹了固定效應(yīng)和隨機效應(yīng)，如何構(gòu)建貝葉斯混合效應(yīng)模型，如何利用先驗預(yù)測檢驗設(shè)置合理的先驗分布，以及如何借助貝葉斯因子進行統(tǒng)計推斷。

本文主要關(guān)注了貝葉斯線性混合效應(yīng)模型，它是其他擴展模型的基礎(chǔ)。一類常用的擴展模型是貝葉斯廣義線性混合效應(yīng)模型（Bayesian Generalized Linear Mixed-effects Models），它能夠更好地處理心理學(xué)實驗中其他常見的因變量，比如二元選擇、反應(yīng)時、閱讀時間、李克特評分等（Dixon,2008;Rousselet &Wilcox,2020）。具體來說，研究者可以在貝葉斯廣義線性混合效應(yīng)模型中為因變量設(shè)定更合理的殘差分布（residual distribution）和鏈接函數(shù)（link function），進而獲得更準(zhǔn)確的估計結(jié)果（Vuorre,2017）。例如，對于二項選擇（如正確標(biāo)記為1，錯誤為0）的因變量，brm()函數(shù)可以通過參數(shù)family=bernoulli("logit")設(shè)置模型的殘差分布為伯努利分布且使用logit鏈接函數(shù)（Bürkner,2017）。再比如對于反應(yīng)時，研究者可以通過family=“l(fā)ognormal” 為其設(shè)置log轉(zhuǎn)換（Wagenmakers &Brown,2007）。此外，貝葉斯分層模型還可以應(yīng)用于復(fù)雜的認(rèn)知模型，比如信號檢測論、漂移擴散模型、強化學(xué)習(xí)模型和多項式加工樹模型等（Ahn et al.,2017；Schad et al.,2020;Tso et al.,2021）。這些分析可以通過R語言的工具包（如hBayesDM）、Stan和JAGS等其他工具來完成（Ahn et al.,2017;Kruschke,2014）?？偟膩碚f，貝葉斯混合效應(yīng)模型憑借其強大的靈活性，能夠很好地適用于多樣的心理學(xué)實驗研究中。

貝葉斯混合效應(yīng)模型的諸多優(yōu)點使其逐漸流行于心理學(xué)研究的數(shù)據(jù)分析之中（Bono et al.,2021）。但研究者在實踐中使用這種方法時仍會面臨一些困難。比如，在初學(xué)時可能不清楚如何設(shè)置合理的先驗，因而會采用工具包提供的默認(rèn)先驗。需要注意的是，這些默認(rèn)先驗并不一定適用于所有的研究，而使用不合適的先驗可能造成模型難以收斂等潛在問題（Veenman et al.,2022）。如3.3小節(jié)所介紹，研究者可以根據(jù)以往的研究和自己的知識自定義先驗，并且可以借助先驗預(yù)測檢驗來審查這些先驗的合理性。其次，當(dāng)使用貝葉斯因子進行假設(shè)檢驗時，研究者一般無法通過單個模型獲得所有感興趣的效應(yīng)。例如model_4結(jié)果只包含了圖片類型和組別的主效應(yīng)，以及他們的交互作用。相應(yīng)地，我們也只能計算這三個效應(yīng)的貝葉斯因子。如果研究者同時對其他效應(yīng)（如簡單效應(yīng)）感興趣，則需要構(gòu)建新的模型并設(shè)置不同的對比編碼來檢驗不同的假設(shè)（Schad et al.,2020）。最后，擬合貝葉斯混合效應(yīng)模型相對更加耗時。隨著數(shù)據(jù)量的變大、模型復(fù)雜程度的提升、以及采樣迭代數(shù)量的增長，貝葉斯混合效應(yīng)模型的擬合時間也會相應(yīng)地呈指數(shù)性增長。對此，研究者可以通過設(shè)置合理的先驗和使用高性能的計算機提升模型的擬合速度。

本文借助心理學(xué)研究示例，通過具體代碼向讀者展示了如何使用貝葉斯混合效應(yīng)模型進行數(shù)據(jù)分析。隨著開放科學(xué)的興起（胡傳鵬等,2016;Nosek,2019;Open Science Collaboration,2015;Jin et al.,2023），越來越多的文章公開其數(shù)據(jù)和代碼，這無疑為研究者提供了更多結(jié)合具體實例的貝葉斯混合效應(yīng)模型的應(yīng)用教程。此外，我們也鼓勵希望對該方法有更加深入和全面了解的讀者閱讀貝葉斯分析相關(guān)的優(yōu)秀書籍（Gelman et al.,2013；Kruschke，2014；Lambert，2018；McElreath，2020）。希望研究者可以對貝葉斯混合效應(yīng)模型有更加深入的了解，能夠廣泛地將其應(yīng)用于多樣的心理學(xué)研究之中。