I've seen some similar (ex ex2) questions, but hopefully this is not a duplicate. As it is mentioned in one of them, I'm using eemmeans to do pairwise comparisons after my linear mixed effects model. Question: WHAT is the most appropriate effect size estimate for these comparisons? eta ? partial eta?*

  • My model + my post-hoc:
mod1 <- lmer(MY_CONT ~  YEAR * GROUP_2 + (1|ID), data = data, REML = FALSE)

group <- emmeans(mod1,~MY_GROUP|YEAR) %>% pairs(adjust="Tukey")
year <- emmeans(mod1,~YEAR|MY_GROUP) %>% pairs(adjust="Tukey")
  • I've also seen eff_size (such as here) option from the same package, but I couldn't understand from its documentation which estimate it is actually doing. I need some help to comprehend which would be the best estimate for me and how to perform that in R. Thanks in advance!

  • tips on how to report these results would be very appreciated too :)

  • is eff_size equivalent to cohen's d?

  • EDIT for bounty:

Russel kindly answered my issue, but I still have some remaining questions:

  • A) What is the Cohen' D type I'm getting see

  • B) Am I estimating it right?

### the model is:

mod1 <- lmer(CONT_Y ~  MY_GROUP * YEAR + (1|ID), data = dfModels)

### estimate eemmeans:

group <- emmeans(mod1,~ MY_GROUP|YEAR)
year <- emmeans(mod1,~ YEAR|MY_GROUP)

### pairwise comparisons:

group_p <- emmeans(mod1,~ MY_GROUP|YEAR) %>% pairs(adjust="Tukey")
year_p <- emmeans(mod1,~ YEAR|MY_GROUP) %>% pairs(adjust="Tukey")

### correctiong sigma and edf

sigmaValues <- VarCorr(mod1)

sigma <- sqrt((0.25743)^2 + (0.15054)^2)

### calculate Cohen's d: 

# eff1 <- eff_size(emm1, sigma = sigma(mod1), edf = df.residual(mod1)) ### before: 

group_p ### check the lowest df
eff1 <- eff_size(group, sigma = sigma, edf = 60)   ## and for group

year_p ### check the lowest df

eff2 <- eff_size(year, sigma = sigma, edf = 60)   ## and for year
  • C) how can I adapt this for lmers ?
simp <- lm(CONT_Y ~  MY_GROUP * YEAR), data = dfModels)
eff_size(pairs(emm), sigma(fiber.lm), df.residual(fiber.lm), method = "identity")
  • NOTE: I guess it goes without saying, but I DON'T have a background on math, so please bear with me :)

  • data:

data <- structure(list(PARTICIPANTS = c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
                                        3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 
                                        7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 
                                        10L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 
                                        14L, 14L, 14L, 14L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 17L, 
                                        17L, 17L, 17L, 18L, 18L, 18L, 18L, 19L, 19L, 19L, 19L, 20L, 20L, 
                                        20L, 20L, 21L, 21L, 21L, 21L), CONT_Y = c(19.44, 20.07, 19.21, 
                                                                                  16.35, 11.37, 12.82, 19.42, 18.94, 19.59, 20.01, 19.7, 17.92, 
                                                                                  18.78, 19.21, 19.27, 18.46, 19.52, 20.02, 16.19, 19.97, 13.83, 
                                                                                  15.93, 14.79, 21.55, 18.8, 19.42, 19.27, 19.37, 17.14, 14.45, 
                                                                                  17.63, 20.01, 20.28, 17.93, 19.36, 20.15, 16.06, 17.04, 19.16, 
                                                                                  20.1, 16.44, 18.39, 18.01, 19.05, 18.04, 19.69, 19.61, 16.88, 
                                                                                  19.02, 20.42, 18.27, 18.43, 18.08, 17.1, 19.98, 19.43, 19.71, 
                                                                                  19.93, 20.11, 18.41, 20.31, 20.1, 20.38, 20.29, 13.6, 18.92, 
                                                                                  19.05, 19.13, 17.75, 19.15, 20.19, 18.3, 19.43, 19.8, 19.83, 
                                                                                  19.53, 16.14, 21.14, 17.37, 18.73, 16.51, 17.51, 17.06, 19.42
                                        ), CATEGORIES = structure(c(1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 
                                                                    1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 
                                                                    1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 
                                                                    1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 
                                                                    1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 
                                                                    1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L), .Label = c("A", 
                                                                                                                            "B"), class = "factor"), MY_GROUP = structure(c(1L, 2L, 1L, 2L, 
                                                                                                                                                                            1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
                                                                                                                                                                            1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
                                                                                                                                                                            1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
                                                                                                                                                                            1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
                                                                                                                                                                            1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L
                                                                                                                            ), .Label = c("G1", "G2"), class = "factor")), row.names = c(NA, 
                                                                                                                                                                                         -84L), class = c("tbl_df", "tbl", "data.frame"))

### rename column:

data <- data %>%  rename(., YEAR = CATEGORIES)

1 Answer 1


You've gotten a lot of good advice in the comments here and in your discussion with Russ Lenth elsewhere. Here's some advice for a simpler way to deal with the effect-size issue and for what might be a more useful model of your data.

WHAT is the most appropriate effect size estimate ...?

I've found "effect size" to be something that sounds like it should be important but ends up being disappointing or even misleading. I find actual differences in magnitudes easier to think about. That might represent my background in biology and biochemistry rather than in social science, where I gather that "effect size" estimates tend to be expected.

In comments, you've found a good deal of skepticism that there is an "appropriate effect size estimate," particularly in a model with random effects. Cohen's $d$ is the ratio between a difference in some type of outcome estimate and a standard deviation estimate. With random effects it's not at all clear what should be included in that standard deviation estimate, as there are both among-participant and residual variances to evaluate.

Nevertheless, you've been required to report some type of "effect sizes" for your study. Part of research training should involve learning how to deal with demands from supervisors or reviewers in a way that is honest to the data and respectful to the demands (however unrealistic or outdated those demands might be). I think that the Meltzer et al paper that you cited gives you a simple way to meet that demand. They say:

We compute Cohen’s $d$ as: $d = \frac{M_E-M_C}{s}$. For effect sizes of training effects (post – pre) within single groups ..., $M_E$ is the average score of the group post-training, $M_C$ is the average score pre-training, and $s$ is the standard deviation of the pre-training scores pooled across all ... groups. For effect size estimates of differences between groups, $M_E$ is the average difference score (post minus pre) of the experimental group... $M_C$ is the average difference score within the Control group, and $s$ is the standard deviation of difference scores within the Control group.

In your case, A presumably represents pre-training (or its equivalent), B represents post-training, and you have groups G1 and G2. I'd suggest reworking your data into a wide form, both for this task and for later suggestions about your model.

dataWide <- tidyr::pivot_wider(data,names_from=YEAR,values_from=CONT_Y)
## YEAR_B - YEAR_A differences for each PARTICIPANT/GROUP combination
dataWide[,"BAdiff"] <- dataWide$B-dataWide$A
## [1] 2.272576
aggregate(BAdiff ~ MY_GROUP, data = dataWide, FUN = function(x) mean(x)/sd(dataWide$A))
##  MY_GROUP    BAdiff
## 1       G1 0.5033097
## 2       G2 0.2382444

What's shown as BAdiff in the output above are Cohen's $d$ values, as defined by Meltzer et al., for the post-pre (B-A)differences divided by the pooled standard deviation (2.27) of pre (A) values. You can report those Cohen's $d$ values and cite Meltzer et al.

For the other comparison, I'm not sure which would be considered the "control group" in your situation (or if that's even a consideration in your study). You would have to make a choice. The post-pre differences and their standard deviations for your groups are:

aggregate(BAdiff ~ MY_GROUP, data = dataWide, FUN = mean)
##   MY_GROUP    BAdiff
## 1       G1 1.1438095
## 2       G2 0.5414286
aggregate(BAdiff ~ MY_GROUP, data = dataWide, FUN = sd)
##   MY_GROUP   BAdiff
## 1       G1 2.317972
## 2       G2 2.803794

Reporting the "effect sizes" this way allows you to met the demands placed on you in a simple way, supported in the literature, that obviates your other questions (at least for this data set and its completely balanced design): "What is the Cohen' D type I'm getting?" "Am I estimating it right?" "How can I adapt this for lmers?" A reader like me might choose to ignore those Cohen's $d$ values and focus instead on the formal statistical analysis.

Suggestion for your model

Your mod1 seems to violate the assumptions about the distribution of residuals in a way that's troubling. Look at plot(mod1) and qqmath(mod1). The first suggests a rise in residuals with higher predicted values, and the second shows some substantial deviations from normality.

The raw data hint at what might be going on. Try this plot:

ggplot(data = data, mapping = aes(x = YEAR, y = CONT_Y,
    group = MY_GROUP, color = MY_GROUP)) + geom_point() + 
    facet_wrap(facets = vars(PARTICIPANTS)) + geom_line()

It looks like there are typically very big changes in scores between the years when the YEAR_A score is very low (about 15 or lower), but not so much when the YEAR_A score is higher. That suggests you need to take the YEAR_A score into account in some way.

You thus might consider a different way to handle repeated measures over time. It can make sense to use the initial (pre-training, A) values as predictors in a model of later (post-training, B) values. That's particularly helpful when there are multiple later time points, but I think it can help you here too. I got a warning when I tried to do that with lmer, but for your type of study a generalized least squares (GLS) model can account for inter-individual correlations appropriately. See Chapter 7 of Frank Harrell's notes on Regression Modeling Strategies. You have to specify a correlation form; corCompSymm is equivalent to the assumption in repeated-measures ANOVA.

gls1 <- nlme::gls(B ~ A * MY_GROUP,       
         correlation = nlme::corCompSymm(form=~1|PARTICIPANTS), data = dataWide)

GLS is just an extension of linear regression (e.g., lm(B ~ A * MY_GROUP)) that takes the within-PARTICIPANT correlations into account. "GLS is equivalent to applying ordinary least squares to a linearly transformed version of the data." Wikipedia

The residuals seem much better behaved than in your mod1 (see plot(gls1) and qqnorm(gls1,abline=c(0,1))). With a continuous predictor A, explore the emtrends() function in emmeans.

emtrends(gls1, pairwise ~ MY_GROUP, var = "A", mode = "df.error", infer = TRUE)$contrasts
#  contrast estimate    SE df lower.CL upper.CL t.ratio p.value
#  G1 - G2     0.478 0.188 36   0.0969    0.859   2.544  0.0154
# Degrees-of-freedom method: df.error 
# Confidence level used: 0.95 

In GROUP_G1 the YEAR_B values (called B in this model) are more positively associated with the values in YEAR_A (called A in this model) than is the case for GROUP_G2. I think that illustrates the difference between your two groups in a more helpful way than your mod1. You still can report the 4 comparisons of original interest, but I don't think that alone does justice to your data.

  • $\begingroup$ Hi, @EdM , thank you very much. I still have some questions, would you mind answered them? The first is: Why shouldn't I calculate the effect sizes based on the emmeans output if concerning that the effect sizes that I need are for the derived pairwise differences of the model? Concerning the variables, $\endgroup$ Feb 12 at 11:28
  • $\begingroup$ Participants were testes twice (Year A and Year B) in two different tests each year (G1 and G2). Hence, each student has 4 scores: Score 1: Year A - Test G1 Score 2: Year A - Test G2 Score 3: Year B - Test G1 Year B - Test G2 (in which G1 and G2 represent different languages of testing, basically) and I need to know 4 things: A) Is G1 and G2 different in Year A? + Is G2 and G1 different in Year B? + Is G1 and G1 different between Year A and Year B? and finally Is G2 and G2 different between Year A and Year B? (basically, I an effect of group and of year) $\endgroup$ Feb 12 at 11:30
  • $\begingroup$ concerning the model suggestion, I really appreciate it, I was also kindly suggested that before (stats.stackexchange.com/questions/588224/…) but the thing is: I think don't know how to explain this estimation rather than OLS, so I didn't go with that for my knowledge limitation at the moment $\endgroup$ Feb 12 at 11:35
  • $\begingroup$ back to the effect's estimation, when I did: eff1 <- eff_size(group, sigma = sigma, edf = 60) being it :sigma = sqrt((error SD)^2 + (subject SD)^2)) , basically my equation was (using group as ex): diff in group G1 and G2 means of Year A (or B, for the 2nd output)/ total sd, right? being it total sd = sum(model sd + subject sd) Question: but why do I have to take the squared root and ^2? since I didn't input the variance but the SD already? Why couldn't I just sum up the results of VarCorr(mod1)'s Std.Dev for both the model (residuals) and ID ? $\endgroup$ Feb 12 at 11:57
  • 1
    $\begingroup$ @LarissaCury Meltzer et al. (and I) see "effect sizes" as being something that can (and should, in my opinion) be considered separately from formal statistical analysis (via emmeans in your case). They say: "In addition to statistical significance testing, we estimated effect sizes for all outcome measures using a simple traditional approach..." (Emphasis added.) There is no harm in that separation of "effect size" reporting from statistical testing; there is no inference being made on effect sizes. Your choice. I added more about my suggested model, which is just a type of least squares. $\endgroup$
    – EdM
    Feb 12 at 17:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.