# Estimated marginal means and arithmetic means are different

I have the following data set where I am testing how fast people can recognise words in different conditions - pre-test (before learning), immediate (straight after learning) and delayed (a long time after learning).

https://pastebin.com/QM3J8XG2


I first inspect my mean RTs. Subjects are faster during the delayed condition and slowest during immediate

df %>%
drop_na %>%
dplyr::group_by(test.time) %>%
dplyr::summarise(meanRT = mean(RT))

test.time meanRT
<fct>      <dbl>
1 delayed    1340.
2 immediate  1484.
3 pretest    1386.


I fit a linear mixed effect models to investigate the effect of test time on RTs and conduct post hoc contrasts using the emmeans package

library(lmerTest)
libary(emmeans)

mod.1 <- lmerTest::lmer(RT ~ test.time + (1 | Ppt.No),
data = df %>%
drop_na)

emm1 = emmeans(mod.1, specs = pairwise ~ test.time)

emm1

$emmeans test.time emmean SE df lower.CL upper.CL delayed 1343 82.9 43.2 1176 1510 immediate 1483 83.0 43.5 1315 1650 pretest 1402 83.0 43.3 1235 1569 Degrees-of-freedom method: kenward-roger Confidence level used: 0.95$contrasts
contrast            estimate   SE   df t.ratio p.value
delayed - immediate   -140.0 36.2 2454 -3.865  0.0003
delayed - pretest      -59.0 36.0 2455 -1.639  0.2295
immediate - pretest     80.9 36.3 2455  2.227  0.0668

Degrees-of-freedom method: kenward-roger
P value adjustment: tukey method for comparing a family of 3 estimates


The contrasts are what I would expect by looking at my arithmetic means - delayed is significantly different from immediate, pretest is significantly different from immediate, and delayed and pretest are not significantly different.

However my estimated marginal means are vastly different from my arithmetic means - why is this? I have read on a few sites that it might be because my design is unbalanced? I investigated this and found that some participants have far fewer observations than others in my sample. Is this a valid reason for my estimated marginal means to be so different? And should I just stick to reporting arithmetic means when describing my post hoc contrasts?

Any help is appreciated, thank you!

• I do not use these "contrasts" but a quick suggestion: maybe you could restrict your analysis to a sample with the same number of observations per person, and see if you get the same issue. ? May 3, 2020 at 11:12

The reason is that estimated marginal means depend on the model, not on the data.

To illustrate, consider this toy example with 5 subjects and two treatments:

foo = data.frame(subj = factor(c(1,1, 2,2,2,   3, 4,4, 5,5,5)),
trt = factor( c(1,1, 1,1,1,   2, 2,2, 2,2,2)),
y = c(3,2, 4,4,5,   7, 9,9, 10,12,15))


First, here are the raw treatment means:

> with(foo, tapply(y, trt, mean))
1        2
3.60000 10.33333


Now, consider fitting a simple one-way fixed-effects model, and obtaining the EMMs:

> foo.lm = lm(y ~ trt, data = foo)
> emmeans(foo.lm, "trt")
trt emmean    SE df lower.CL upper.CL
1      3.6 0.995  9     1.35     5.85
2     10.3 0.908  9     8.28    12.39

Confidence level used: 0.95


These EMMs are in fact the same as the raw treatment means.

Next, do the same but use a mixed model with random subject effects:

> foo.lmer = lmer(y ~ trt + (1|subj), data = foo)
> emmeans(foo.lmer, "trt")
trt emmean   SE   df lower.CL upper.CL
1     3.45 1.58 2.70    -1.90     8.81
2     9.68 1.35 3.06     5.43    13.93

Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95


These results are different from the ones obtained from foo.lm --- again, because EMMs depend on the model.

Please note that these discrepancies are not a peculiarity of the emmeans package. They are discrepancies in the estimated fixed effects of the two models. For the simple model, we have

> coef(foo.lm)
(Intercept)        trt2
3.600000    6.733333


Interpreting these coefficients (in terms of the default contrasts "contr.treatment" used for coding the factor), the treatment predictions are 3.6000 for trt = 1 (the reference level), and 3.6000 + 6.7333 = 10.3333 for trt = 2. For the mixed model,

> fixef(foo.lmer)
(Intercept)        trt2
3.452249    6.227216


implying predictions of 3.4522 for trt = 1 and 3.4522 + 6.2272 = 9.6794 for trt = 2.

If the mixed model is the right model for the data, then the right post hoc means would be the ones based on that model. It would be wrong to report the raw means and use them as the basis of a post hoc analysis.

• @Lenth, if I use lm(), will the EMMs always be the same as raw treatment means? Sep 30, 2021 at 10:19
• Only if there is only one factor, or if the data are balanced with more than one factor. Read the "basics" vignette for an example Sep 30, 2021 at 13:40