What's the appropriate effect size estimate (and power analysis) for post-hoc regression?

Question

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)
sigmaValues

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)

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
sd(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:

library(ggplot2)
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.