分析接上一篇文章：R語言之重複測量方差分析——ezANOVA的使用與解析

寫在前言

開篇咱依舊先吹吹水，話說當前心理學研究面臨著可重複性危機（Reproducibility Crisis），對NHST(Null Hypothesis Significance Test)的濫用或誤用似乎成為眾矢之的。記得在我剛開學那會，各路研究方法學的教授一上來就是對NHST和p-hacking進行批判（通過修改、剔除不顯著的數據來操縱p值達到顯著水平），著實讓我懷疑人生。

p值：我TM做錯什麼了……

實際上，NHST中用於驗證假設的p值並不能解釋假設本身。根據美國統計協會（ASA）在2016年的一份聲明中對p值的解釋：在統計檢驗當中，p值只能解釋數據與假設之間的關係，並不能衡量研究假設為真時的概率(Wasserstein & Lazar, 2016)。此外，翻翻以前的教科書我驚訝地發現，NHST關於p值的定義是：假定 H_{0} 為真時，以完全相同的條件無數次地重複當前實驗/測量/抽樣，得到的結果與H_{0}一致或者極端相反的概率。而現實生活中，「無數次」這個定語往往被忽略，許多研究者（包括我之前）都單純地誤認為p值是一次檢驗中拒絕H_{0}時犯錯誤的概率。

僅上述兩點，就足以讓我們所有人反思：到底我們的研究是在取悅自己還是真的在推動學科的發展？

而在這一背景之下，陸續有研究者推薦使用貝葉斯因子（Bayes Factor）來替代NHST中的p值。

什麼是貝葉斯統計(Bayesian Statistics)？ 這裡推薦閱讀一個大佬的回答：

怎樣用非數學語言講解貝葉斯定理（Bayess theorem）？ - 馬同學的回答 - 知乎

如何理解貝葉斯推斷和beta分佈？

貝葉斯因子是貝葉斯統計中用來進行模型比較和假設檢驗的方法。在貝葉斯統計框架下的假設檢驗中，相當於我們根據當前收集到的數據來檢驗某個理論模型為真的可能性。因此，貝葉斯因子代表的是當前數據對H_{0}與H_{1}支持的強度之間的比率。相比之下，p值在研究中反映的只是樣本均值之間的差別有無統計學意義，並不表示其差別大小。

一般來說，我們用 BF_{10} 來表示數據支持備擇假設（alternative hypothesis，即 H_{1} ）的程度。其中下標的「10」即代表 H_{1}與H_{0}，所以同理，BF_{01} 則代表數據支持原假設（null hypothesis，即 H_{0} ）的程度。

BF_{10} = frac{p(data|H_{1})}{p(data|H_{0})} = frac{1}{BF_{01}}

對於貝葉斯因子大小的意義一般參照Jeffreys (1961)給出的表述進行解釋

當前有不少用於做貝葉斯統計的工具，如 WinBUGs 、JAGS 、Stan 等（後兩者在R-Studio中均有相應的包），本篇以R-Studio中的『BayesFactor』包以及最近幾年非常流行的JASP(JASP - A Fresh Way to Do Statistics)對上一篇文章的例子進行貝葉斯因子的計算。JASP的運算是基於R的代碼，可以看作具有可視化操作的R，對於剛入門想了解貝葉斯因子又不熟悉R語言的朋友們來說是個不錯的嘗試。

下面分別運行R和JASP的結果，最後再對結果進行解釋。

首先以R為例，接著上一篇的內容，對數據支持主效應和交互作用的程度進行分析。

> library(BayesFactor) #載入包
> bf <- anovaBF( Score ~ Familiarity * Density + Subject, data = data, whichRandom = "Subject" )
> summary(bf)
Bayes factor analysis
--------------
[1] Familiarity + Subject : 7214.359 ±0.9%
[2] Density + Subject : 3104.927 ±0.73%
[3] Familiarity + Density + Subject : 29442574316 ±1.15%
[4] Familiarity + Density + Familiarity:Density + Subject : 2.178987e+19 ±2.14%

Against denominator:
Score ~ Subject
---
Bayes factor type: BFlinearModel, JZS

結果輸出了4個模型：

  1. 只考慮Familiarty主效應的模型
  2. 只考慮Density主效應的模型
  3. 同時考慮Familiarity和Density主效應的模型
  4. 以及即考慮了Familiarity和Density的主效應，又考慮兩者交互作用的模型。

右邊即 BF_{10} ，其數值大小表明數據對每個模型的支持程度。這裡可以看到，包含交互作用的模型[4]獲得最多證據，與上一步的重複測量方差分析結果一致。我們可以畫個圖更直觀地呈現上述4個模型的對比。

>plot(bf)

數據明顯更支持包含交互作用的模型

如果想知道每個主效應具體的BF_{10}是多少，我們可以執行如下運算：

> bf_Familiarity <- bf[3]/bf[2]
> bf_Density <- bf[3]/bf[1]
> bf_int <- bf[4]/bf[3]

> bf_Familiarity
Bayes factor analysis
--------------
[1] Familiarity + Density + Subject : 9482534 ±1.36%

Against denominator:
Score ~ Density + Subject
---
Bayes factor type: BFlinearModel, JZS

> bf_Density
Bayes factor analysis
--------------
[1] Familiarity + Density + Subject : 4081107 ±1.46%

Against denominator:
Score ~ Familiarity + Subject
---
Bayes factor type: BFlinearModel, JZS

> bf_int
Bayes factor analysis
--------------
[1] Familiarity + Density + Familiarity:Density + Subject : 740080349 ±2.43%

Against denominator:
Score ~ Familiarity + Density + Subject
---
Bayes factor type: BFlinearModel, JZS

而JASP的操作就更加簡單了，與SPSS十分相似。

依次執行：ANOVA ——Bayesian Repeated Measures ANOVA，在彈出的界面依次設置好自變數以及自變數的各個水平，然後將左側的數據傳送到右邊的Cells，Bayes Factor默認 BF_{10} ，Order默認與零假設（Null model）進行比較，然後單擊OK。

接下來看輸出結果

Model Comparison

左邊的Models代表所有納入分析的模型，與R的內容一致。P(M)代表所有模型都被賦予了相同的先驗概率。P(M|data)代表根據數據結果得到的後驗概率。可以看到所有的後驗離散分佈都集中在交互作用的模型上。BF_{M}代表從先驗到後驗分佈的概率的變化。從數據增加模型概率的層面上看，明顯是交互作用的模型獲得了數據的支持。error%可以簡單理解為貝葉斯因子的變異係數。

ps：這裡如果仔細觀察，我們會發現JASP計算出來的BF_{10}結果跟上面R的不一樣，這是由於兩次運算受到不同的蒙特卡羅抽樣雜訊的影響。貝葉斯因子的計算隨著每次抽樣的變化而變化，所以即使相同的操作在JASP中重新運算一編，出來的數值依然會有輕微的變化，但並不影響其結果。

說了那麼多主要還是看BF_{10}，所有納入考慮的模型與零假設相比都顯著。但包含了交互作用的模型顯然獲得最多證據的支持。與R一樣，如果我們想知道數值具體有多大，我們可以用上述兩個模型的BF_{10}相除:

BF_{10(Familiarity*Density)} = 2.182e+19 / 2.993e+10 = 7.29e+08

同理，如果我們想知道數據支持兩個主效應顯著的證據有多大，則分別用Familiarty+Density，與Familiarty和Density的 BF_{10} 相除即可:

BF_{10(Familiarity)} = 2.993e+10 / 7105.341 = 4.21e+06

BF_{10(Density)} = 2.993e+10 / 3108.819 = 7.11e+03

綜上所述，R和JASP的結果均表明，相比只考慮主效應的模型，我們的數據有7.29e+08倍（約7億）支持考慮了交互作用的模型。所以，在文章中報告交互作用時我們可以這樣子寫（主效應同理）：

F (1.22, 10.97) = 62.88, p< 0.001, eta_{p}^{2}= 0.875, BF_{10} > 1000，根據Jeffreys (1961)提出的分類標準，極強的證據支持我們的研究存在交互作用。

後續的簡單效應檢驗可以用『BayesFactor』中的ttestBF( )對交互作用的各個水平進行兩兩對比。 至於在JASP中如何實現，胡傳鵬等(2018)在文章中也詳細地講解了其原理與操作過程，這裡也推薦大家閱讀。

關於『BayesFactor』包的其它函數及其使用，可參考：

richarddmorey.github.io

參考文獻

Theory of ProbabilityTheory of probability. Oxford: UK Oxford University Press.

Wagenmakers, E. J., Love, J., Marsman, M., Jamil, T., Ly, A., Verhagen, J., ... & Meerhoff, F. (2018). Bayesian inference for psychology. Part II: Example applications with JASP. Psychonomic bulletin & review, 25(1), 58-76. 胡傳鵬, 孔祥禎, & 彭凱平. (2018). 貝葉斯因子及其在 JASP 中的實現.心理科學進展,26(6), 951-965.

作者聲明：

本文原創性歸作者所有，轉載請註明來源和作者。內容只用於交流與學習，禁止用於任何商業用途。

推薦閱讀：

相關文章