貝葉斯混合效應模型：基于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）。希望研究者可以對貝葉斯混合效應模型有更加深入的了解，能夠廣泛地將其應用于多樣的心理學研究之中。