Ytukyg

I've been trying to calculate marginal means for my lmer & glmer in R. I found the emmeans function and I've been trying to understand it and apply it to my model. I found that it's hard to get the means for an interaction, so I'm starting with just additive predictors, but the function doesn't work the way it's presented in examples (e.g. here https://cran.r-project.org/web/packages/emmeans/vignettes/sophisticated.html)

emmeans(Oats.lmer, "nitro") 

nitro    emmean       SE   df  lower.CL  upper.CL

0.0  78.89207 7.294379 7.78  61.98930  95.79484

 0.2  97.03425 7.136271 7.19  80.25029 113.81822

0.4 114.19816 7.136186 7.19  97.41454 130.98179

0.6 124.06857 7.070235 6.95 107.32795 140.80919

what I'm getting is:

emmeans(model2, "VariableA")

VariableA  emmean       SE    df lower.CL upper.CL

0.4657459 2649.742 120.8955 19.07 2396.768 2902.715

Only one line and the variable is averaged instead of split into 0 and 1 (which are the values in the dataset, and maybe the problem is that it's categorical?)
The model I'm running is :

model2 = lmer (rt ~ variableA + variableB + (1 |participant) + (1 |sequence/item), data=memoryData, REML=FALSE)

EDIT: The data file is quite big and I wasn't sure how to extract useful information from it, but here is the structure:

> str(memoryData) 

'data.frame':   3168 obs. of  123 variables:

 $ participant    : int  10 10 10 10 10 10 10 10 10 10 ...

$ variableA      : int  1 1 1 1 1 1 1 1 1 1 ...

 $ variableB      : int  1 1 1 1 1 1 1 1 1 1 ...

$ sequence: int  1 1 1 1 1 1 1 1 1 1 ...

$ item     : int  25 26 27 28 29 30 31 32 33 34 ...

$ accuracy       : int  1 1 1 1 1 1 0 1 1 1 ...

 $ rt             : num  1720 1628 1728 2247 1247 ...

Why is the function not working for me?
And as a further question, is there a way to get these means when I include interaction between variables A and B?

EDIT 2: ok, it did work when I changed it to factor, I guess my method of doing it was incorrect. But I'm still not sure how to calculate it when there is an interaction? Because with this method, R says "NOTE: Results may be misleading due to involvement in interactions"

To see marginal means of interactions, add all variables of the interaction term to emmeans(), and you need to use the at-argument if you want to see the marginal means at different levels of the interaction terms.

Here are some examples, for the average effect of the interaction, and for marginal effects at different levels of the interaction term. The latter has the advantage in terms of visualization.

library(ggeffects)

library(lme4)

library(emmeans)



data("sleepstudy")

sleepstudy$inter <- sample(1:5, size = nrow(sleepstudy), replace = T)

m <- lmer(Reaction ~ Days * inter + (1 + Days | Subject), data = sleepstudy)



# average marginal effect of interaction

emmeans(m, c("Days", "inter"))

#>  Days    inter   emmean      SE    df lower.CL upper.CL

#>   4.5 2.994444 298.3427 8.84715 16.98 279.6752 317.0101

#> 

#> Degrees-of-freedom method: kenward-roger 

#> Confidence level used: 0.95



# marginal effects at different levels of interactions -

# useful for plotting

ggpredict(m, c("Days [3,5,7]", "inter"))

#> 

#> # Predicted values of Reaction 

#> # x = Days 

#> 

#> # inter = 1

#>  x predicted std.error conf.low conf.high

#>  3   279.349     8.108  263.458   295.240

#>  5   304.839     9.818  285.597   324.082

#>  7   330.330    12.358  306.109   354.551

#> 

#> # inter = 2

#>  x predicted std.error conf.low conf.high

#>  3   280.970     7.624  266.028   295.912

#>  5   304.216     9.492  285.613   322.819

#>  7   327.462    11.899  304.140   350.784

#> 

#> # inter = 3

#>  x predicted std.error conf.low conf.high

#>  3   282.591     7.446  267.997   297.185

#>  5   303.593     9.384  285.200   321.985

#>  7   324.594    11.751  301.562   347.626

#> 

#> # inter = 4

#>  x predicted std.error conf.low conf.high

#>  3   284.212     7.596  269.325   299.100

#>  5   302.969     9.502  284.345   321.594

#>  7   321.726    11.925  298.353   345.099

#> 

#> # inter = 5

#>  x predicted std.error conf.low conf.high

#>  3   285.834     8.055  270.046   301.621

#>  5   302.346     9.839  283.062   321.630

#>  7   318.858    12.408  294.540   343.177

#> 

#> Adjusted for:

#> * Subject = 308



emmeans(m, c("Days", "inter"), at = list(Days = c(3, 5, 7), inter = 1:5))

#>  Days inter   emmean        SE    df lower.CL upper.CL

#>     3     1 279.3488  8.132335 23.60 262.5493 296.1483

#>     5     1 304.8394  9.824196 20.31 284.3662 325.3125

#>     7     1 330.3300 12.366296 20.69 304.5895 356.0704

#>     3     2 280.9700  7.630745 18.60 264.9754 296.9646

#>     5     2 304.2160  9.493225 17.77 284.2529 324.1791

#>     7     2 327.4621 11.901431 17.84 302.4420 352.4822

#>     3     3 282.5912  7.445982 16.96 266.8786 298.3038

#>     5     3 303.5927  9.383978 16.98 283.7927 323.3927

#>     7     3 324.5942 11.751239 16.98 299.7988 349.3896

#>     3     4 284.2124  7.601185 18.34 268.2639 300.1609

#>     5     4 302.9694  9.504102 17.85 282.9900 322.9487

#>     7     4 321.7263 11.927612 17.99 296.6666 346.7860

#>     3     5 285.8336  8.076779 23.02 269.1264 302.5409

#>     5     5 302.3460  9.845207 20.48 281.8399 322.8521

#>     7     5 318.8584 12.416642 21.02 293.0380 344.6788

#> 

#> Degrees-of-freedom method: kenward-roger 

#> Confidence level used: 0.95

And a plotting example:

ggpredict(m, c("Days", "inter [1,3,5]")) %>% plot()

You say that "changing the vari[a]ble to factor doesn't help", but I would think this would (as documented in the emmeans FAQ):

 md <- transform(memoryData,

         variableA=factor(variableA),

         variableB=factor(variableB))



 model2 = lmer (rt ~ variableA + variableB + 

    (1 |participant) + (1 |sequence/item), data=md, REML=FALSE)

 emmeans(model2, ~variableA)

 emmeans(model2, ~variableB)

 emmeans(model2, ~variableA + variableB)

If this really doesn't work, then we need a reproducible example ...