计算来自 lme4 模型的所有固定效应 anova 的部分η2η2

机器算法验证 r 混合模式 重复测量 lme4-nlme 多层次分析
2022-03-12 09:14:26

Disclamer:我不确定在哪里发布这个问题:CV 或 SO,但最终决定先在这里尝试

一位审稿人要求我为我论文的辅助分析中报告,这是我所在领域的标准) 。ηp2Ft

具体来说,我提到了调用的固定A:B:C:response效果之间的重要 4 向交互(anova(lmer.model)emmeanspairs()

我知道至少存在关于计算混合效果模型的效果大小是否有意义的争论,但我想在没有任何进一步讨论的情况下让审稿人满意。

我遇到了一些在 LMEM 的上下文中提到效果大小的消息来源,但我觉得数学不够强大,无法理解它。

我的问题是双重的:

  1. 如果有任何可引用的方式来为 anovas测试我该怎么做一个包/函数/脚本会很有帮助。最终提示从对象中检索关键值以手动获取 eff.sizes 也会有所帮助η2ηp2ω2FRlmer()

    例如。我知道给出但没有我不知道如何获得它η2=SSEffect/SSTotalanova(lmer_obj)SSEffectSSTotalR

  2. 如果没有办法计算它 - 在回答审稿人为什么我坚持跳过效果大小时,我可以引用什么论文(不是博客文章)。

请注意,我对随机效应本身不感兴趣——我定义了一个由随机效应的设计结构证明的最大值,只是为了控制更多的误差方差,我唯一感兴趣的是固定效应——即来自 anova 表的 Fs 和相应的边际均值。

一些示例性结果(对应于具有随机值的真实数据结构)如下:

A:B:C:response交互作用在 F(1; 12082,1)=4,60 时显 )F(1;12082,1)=4,60;p=.032,(ηp2=xxx

> 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 )F(1;809,32)=4,68;p=.031,(ηp2=xxx

> 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

在我的真实数据集中,有更多的预测效果,但想法仍然相同。我在APA 6 ed s 或 s 。规则,我被要求为它们添加某种效果大小。Ft

任何有关此事的帮助将不胜感激。

完全可重现的例子

###############################################################
#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"))
1个回答

首先,我认为@Henrik 的回答合理的,应该(至少在摘录中)在这个网站上说明:

计算模型拟合的全局度量(例如 R2)已经充满了复杂性,并且无法找到简单的单个数字,这一事实应该暗示对模型参数的子集(即主效应)这样做或交互)更加困难。鉴于此,我不建议尝试为混合模型找到标准化效应大小的度量。

这是他建议对像您这样的审稿人的回答(尽管 I 型错误论点主要适用于具有交叉(随机)分组因素的设计):

不幸的是,由于方差在线性混合模型中的划分方式(例如,Rights & Sterba,2019 年),不存在为单个模型项(例如主效应或交互作用)计算标准效应大小的商定方法。尽管如此,我们还是决定在我们的分析中主要使用混合模型,因为混合模型在控制 I 类错误方面比其他方法要好得多,因此混合模型的结果更有可能推广到新的观察结果(例如,Barr、Levy、Scheepers, & Tily,2013 年贾德、韦斯特法尔和肯尼,2012 年)。只要有可能,我们就会报告非标准化的效应量,这与如何报告效应量的一般建议一致(例如, 佩克和弗洛拉,2018 年)。

话虽如此,我个人喜欢Johnson (2014)度量,它是Nakagawa 和 Schielzeth (2013)相应版本的扩展R2

边际 R 平方( ) 可以解释为模型中所有固定效应解释的方差,条件R 平方( ) 估计模型中所有固定效应和所有随机效应共同解释的方差. Rm2Rc2

下面我展示了我们如何通过函数从对象度量。R2lmerModlmerModLmerTestMuMIn::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