8
$\begingroup$

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"))
$\endgroup$
  • $\begingroup$ checking out the links from googling "site:stat.ethz.ch/pipermail/r-sig-mixed-models "effect size" eta" now ... $\endgroup$ – Ben Bolker Jul 27 '18 at 21:39
  • 3
    $\begingroup$ I suggest to not report standardized effect for mixed models and instead report unstandardized effect sizes or alternatively explain that it is simply not possible. See my longer response with pointers to literature here: afex.singmann.science/forums/topic/… $\endgroup$ – Henrik Jul 28 '18 at 15:10
  • $\begingroup$ @BenBolker I googled that but couldn't find anything useful $\endgroup$ – blazej Jul 30 '18 at 7:03
  • 1
    $\begingroup$ I think @Henrik's answer is definitive. I was going to point you to his afex package, since there are some functions in there that appear to do a bit with partial eta-squared, but since he has responded himself ... $\endgroup$ – Ben Bolker Jul 30 '18 at 12:19
  • $\begingroup$ What is SO you refer to in the disclaimer? $\endgroup$ – Clark Jul 31 '18 at 23:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.