Compute partial $\eta^2$ for all fixed effects anovas from a lme4 model

Question

Disclamer: I wasn't sure where to post this question: CV or SO, but eventually decided to try here first

I've been asked by one of the reviewers to add effects sizes (preferably $\eta^2_p$ which is standard in my field) for all the $F$ and $t$ tests reported in auxiliary analyses of my paper.

To be specific I'm mentioning a significant 4-way interaction (A:B:C:response) between fixed effects of a anova(lmer.model) call and explaining it with simple interactions and estimated marginal means using emmeans package and pairs() function (reproducible exemplary code is attached below).

I'm aware that there at least a debate on whether or not computing effect sizes for mixed effects models makes sense, but I'd like to satisfy the reviewer without any further discussion.

I've come across some sources mentioning effect sizes in the context of LMEM but I don't feel strong enough with math to understand it.

My question here is twofold:

1. If there is any citable way to produce $\eta^2$, $\eta^2_p$ or $\omega^2$ for anovas $F$ tests how can I do it? An R package / function / script would be great help. Eventually hints for retrieving crucial values from a lmer() object in order to get eff.sizes by hand would be helpful as well

   Eg. I know that $\eta^2 = {SS_{Effect}}/{SS_{Total}}$ and anova(lmer_obj) gives $SS_{Effect}$ but no $SS_{Total}$ and I have no clue how to get it with R

2. If there is no way computing it - what paper (not blog post) can I cite when answering the reviewer why I insisted on skipping effect sizes.

Note that I'm not interested in random effects per se - I defined a maximum justified by design structure of random effects only to control more error variance and my only interest are the fixed effects - namely Fs from the anova table and corresponding marginal means.

Some exemplary results (corresponding to a real data structure by with random values) are as follow:

The A:B:C:response interaction is significant at $F(1; 12082,1)=4,60; p=.032, (\eta^2_p=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

Follow-up comparisons show that the interaction of A:B:C is only valid for response=0: $F(1; 809,32)=4,68; p=.031, (\eta^2_p=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

In my real dataset there are more predicted effects but the idea is still the same. I'm reporting $F$s or $t$s following APA 6 ed. rules and i'm being asked to add some sort of effect sizes to them.

Any help with this matter will be greatly appreciated.

Full reproducible example

###############################################################
#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"))

Answer 1

First, I think @Henrik's answer here is sound and should (at least in excerpts) be stated on this site:

The fact that calculating a global measure of model fit (such as R2) is already riddled with complications and that no simple single number can be found, should be a hint that doing so for a subset of the model parameters (i.e., main-effects or interactions) is even more difficult. Given this, I would not recommend to try finding a measure of standardized effect sizes for mixed models.

This is what he suggests to answer to a reviewer like yours (although the type-I-error argument mainly applies to designs with crossed (random) grouping factors):

Unfortunately, due to the way that variance is partitioned in linear mixed models (e.g., Rights & Sterba, 2019), there does not exist an agreed upon way to calculate standard effect sizes for individual model terms such as main effects or interactions. We nevertheless decided to primarily employ mixed models in our analysis, because mixed models are vastly superior in controlling for Type I errors than alternative approaches and consequently results from mixed models are more likely to generalize to new observations (e.g., Barr, Levy, Scheepers, & Tily, 2013; Judd, Westfall, & Kenny, 2012). Whenever possible, we report unstandardized effect sizes which is in line with general recommendation of how to report effect sizes (e.g., Pek & Flora, 2018).

That being said, personally, I like the $R^2$ measures of Johnson (2014) which are extensions to the corresponding versions of Nakagawa and Schielzeth (2013):

The marginal R squared ($R^2_m$) can be interpreted as the variance explained by all fixed effects in the model, the conditional R squared ($R^2_c$) estimates the variance explained by all fixed effects and all random effects in the model taken together.

Below I show how we can obtain these $R^2$ measures from lmerMod and lmerModLmerTest objects via the MuMIn::r.squaredGLMM() function.

library("MuMIn")

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

Additionally, the semi-partial (marginal) R squared describes the variance explained by each fixed effect adjusted for the other predictors in the model. Jaeger, Edwards, Das, and Sen (2017) show how one can calculate this measure as a Wald statistic of the desired subset of fixed effects.

We can calculate the semi-partial R squared for each fixed effect via the r2glmm::r2beta() function. Note that you should get the r2glmm package from GitHub as the CRAN version is not up to date (accessed 2019-12-08) and misses important functionality.

# 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