計算部分𝜂2這2eta^2對於來自 lme4 模型的所有固定效應方差分析

Disclamer：我不確定在哪裡發布這個問題：CV 或 SO，但最終決定先在這裡嘗試

一位審閱者要求我添加效果大小（最好這在我的領域是標準的）對於所有和在我的論文的輔助分析中報告的測試。

具體來說，我提到了調用的固定A:B:C:response效果之間的重要 4 向交互（anova(lmer.model)``emmeans``pairs()

我知道至少存在關於計算混合效果模型的效果大小是否有意義的爭論，但我想在沒有任何進一步討論的情況下讓審稿人滿意。

我遇到了一些在 LMEM 的上下文中提到效果大小的消息來源，但我覺得數學不夠強大，無法理解它。

我的問題是雙重的：

1. 如果有任何可引用的生產方式,或者 用於方差分析測試我該怎麼做？一個R包/函數/腳本會很有幫助。最終提示從對lmer()像中檢索關鍵值以手動獲取 eff.sizes 也會有所幫助

例如。我知道並 anova(lmer_obj)給出但不是我不知道如何得到它R 2. 如果沒有辦法計算它 - 在回答審稿人為什麼我堅持跳過效果大小時，我可以引用什麼論文（不是博客文章）。

請注意，我對隨機效應本身不感興趣——我定義了一個由隨機效應的設計結構證明的最大值，只是為了控制更多的誤差方差，我唯一感興趣的是固定效應——即來自 anova 表的 Fs 和相應的邊際均值。

一些示例性結果（對應於具有隨機值的真實數據結構）如下：

交互作用A:B:C:response顯著)

> anova(m0)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F) 
A 1.6717 1.6717 1 134.7 0.4210 0.51755 
B 0.1375 0.1375 1 11860.0 0.0346 0.85238 
C 7.1708 7.1708 1 133.5 1.8058 0.18129 
response 3.7775 3.7775 1 12070.3 0.9513 0.32940 
A:B 1.1291 1.1291 1 11872.8 0.2844 0.59387 
A:C 3.2427 3.2427 1 121.1 0.8166 0.36796 
B:C 4.1048 4.1048 1 121.8 1.0337 0.31130 
A:response 0.0000 0.0000 1 12080.8 0.0000 0.99828 
B:response 4.2350 4.2350 1 12078.9 1.0665 0.30176 
C:response 1.6567 1.6567 1 12082.3 0.4172 0.51834 
A:B:C 8.5249 8.5249 1 131.4 2.1469 0.14525 
A:B:response 0.8765 0.8765 1 132.0 0.2207 0.63926 
A:C:response 1.5150 1.5150 1 119.5 0.3815 0.53797 
B:C:response 0.5921 0.5921 1 122.5 0.1491 0.70005 
A:B:C:response 18.2803 18.2803 1 12082.1 4.6036 0.03192 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

後續比較表明 的交互作用A:B:C僅對 有效response=0：)

> joint_tests(m0, by="response")
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, set emm_options(pbkrtest.limit = 12292) or larger,
but be warned that this may result in large computation time and memory use.
response = no:
model term df1      df2 F.ratio p.value
A            1   858.01   0.161  0.6884
B            1 12177.54   0.530  0.4664
C            1   837.27   0.212  0.6453
A:B          1   785.54   0.365  0.5457
A:C          1   507.23   0.872  0.3509
B:C          1   476.41   0.145  0.7035
A:B:C        1   809.32   4.680  0.0308

response = yes:
model term df1      df2 F.ratio p.value
A            1   184.95   0.361  0.5489
B            1 11853.63   0.601  0.4381
C            1   182.47   3.237  0.0736
A:B          1   177.43   0.000  0.9874
A:C          1   130.20   0.072  0.7882
B:C          1   123.71   1.465  0.2285
A:B:C        1   178.96   0.235  0.6281

在我的真實數據集中，有更多的預測效果，但想法仍然相同。我在報告或s 遵循*APA 6 版。*規則，我被要求為它們添加某種效果大小。

任何有關此事的幫助將不勝感激。

完全可重現的例子

###############################################################
#Simulate data replicating the structure of a real life example
###############################################################

library(AlgDesign) #for generating all levels a factorial design)

set.seed(2)
df <-gen.factorial(c(16,2,2,2,100), factors = "all", 
                  varNames = c("rep", "A", "B", "C", "Subject"))
df$rep <- as.numeric(df$rep)
df$Subject <- as.numeric(df$Subject)
logRT <- rnorm(n=12800, m=7, sd=2)
trialno<- rep(1:128, times = 100)
response <- factor(sample(0:1, 12800, prob = c(0.3, 0.7), replace = T), 
                  labels= c("no", "yes"))

#Simulate some values to drop (set as missings) due to extremly low latency
missingRTs<- as.logical(sample(0:1, 12800, prob = c(0.96, 0.04), replace = T))
logRT[missingRTs==1] <- NA

df <- cbind(df, logRT, trialno, response)
df <- df[complete.cases(df),]

##########################
#Setup model with afex
########################## 

library(afex)
m0 <- mixed(logRT ~ A*B*C*response + (A*B*C*response||Subject), 
           data = df, return = 'merMod', method = "S", expand_re = TRUE)

##########################
#Get results for paper
########################## 

anova(m0)

emm_options(lmerTest.limit = 12292)
em0<-emmeans(m0, ~A*B*C*response)

joint_tests(m0, by="response")
joint_tests(m0, by=c("response", "B"))

首先，我認為@Henrik 的回答是合理的，應該（至少在摘錄中）在這個網站上說明：

計算模型擬合的全局度量（例如 R2）已經充滿了複雜性，並且找不到簡單的單個數字，這一事實應該暗示對模型參數的子集（即主效應）這樣做或交互）更加困難。鑑於此，我不建議嘗試為混合模型找到標準化效應大小的度量。

這是他建議對像您這樣的審稿人的回答（儘管 I 型錯誤論點主要適用於具有交叉（隨機）分組因素的設計）：

不幸的是，由於方差在線性混合模型中的劃分方式（例如，Rights & Sterba，2019 年），不存在為單個模型項（例如主效應或交互作用）計算標準效應大小的商定方法。儘管如此，我們還是決定在我們的分析中主要使用混合模型，因為混合模型在控制 I 類錯誤方面比其他方法要好得多，因此混合模型的結果更有可能推廣到新的觀察結果（例如，Barr、Levy、Scheepers， & Tily，2013 年；賈德、韋斯特法爾和肯尼，2012 年）。只要有可能，我們就會報告非標準化的效應量，這與如何報告效應量的一般建議一致（例如， 佩克和弗洛拉，2018 年）。

話雖如此，就個人而言，我喜歡 $ R^2 $ Johnson (2014)的度量，它們是Nakagawa 和 Schielzeth (2013)相應版本的擴展：

邊際R平方( $ R^2_m $ ) 可以解釋為模型中所有固定效應解釋的方差，條件 R 平方( $ R^2_c $ ) 估計模型中所有固定效應和所有隨機效應共同解釋的方差。

下面我將展示我們如何獲得這些 $ R^2 $ 通過該功能測量lmerMod和測量lmerModLmerTest對象。MuMIn::r.squaredGLMM()

library("MuMIn")

r.squaredGLMM(m0)
# R2m R2c
# [1,] 0.001241324 0.01762294

此外，半部分（邊際）R 平方描述了由針對模型中的其他預測變量調整的每個固定效應所解釋的方差。Jaeger、Edwards、Das 和 Sen (2017)展示瞭如何將此度量計算為所需固定效應子集的 Wald 統計量。

我們可以通過函數計算每個固定效應的半偏 R 平方r2glmm::r2beta()。請注意，您應該r2glmm從 GitHub 獲取包，因為 CRAN 版本不是最新的（訪問時間為 2019-12-08）並且缺少重要功能。

# devtools::install_github('bcjaeger/r2glmm')
library("r2glmm")

r2beta(m0, method = "nsj")
# Effect Rsq upper.CL lower.CL
# 1 Model 0.001 0.004 0.001
# 16 A:B:C:response 0.000 0.001 0.000
# 12 A:B:C 0.000 0.001 0.000
# 4 C 0.000 0.001 0.000
# 8 B:C 0.000 0.001 0.000
# 10 B:response 0.000 0.001 0.000
# 7 A:C 0.000 0.001 0.000
# 5 response 0.000 0.001 0.000
# 2 A 0.000 0.001 0.000
# 14 A:C:response 0.000 0.001 0.000
# 11 C:response 0.000 0.001 0.000
# 6 A:B 0.000 0.001 0.000
# 13 A:B:response 0.000 0.001 0.000
# 15 B:C:response 0.000 0.000 0.000
# 3 B 0.000 0.000 0.000
# 9 A:response 0.000 0.000 0.000

引用自：https://stats.stackexchange.com/questions/358927