2
$\begingroup$

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!

$\endgroup$
1
  • $\begingroup$ 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. ? $\endgroup$ May 3, 2020 at 11:12

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ @Lenth, if I use lm(), will the EMMs always be the same as raw treatment means? $\endgroup$
    – Cmagelssen
    Sep 30, 2021 at 10:19
  • 1
    $\begingroup$ 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 $\endgroup$
    – Russ Lenth
    Sep 30, 2021 at 13:40

Your Answer

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

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