I have the following mixed model

subject <- factor(rep(c(1,2,3,4,5,6),each=12))
dep <- c(0.3763244126,0.2185001692,0.4191841742,0.9812978664,0.7429273683,0.6254715486,0.6200958213,0.4693300191,0.1779032899,0.0035873980,0.8821949826,0.4818012617,0.0008013437,0.6280700732,0.7126500814,0.4984349359,0.2457996449,0.3085733312,0.5903398243,0.3704800352,0.8215325437,0.0445236221,0.1849731791,0.3670945817,0.0022268933,0.1630332691,0.9734050406,0.2638539758,0.8550054496,0.9413964085,0.4548943471,0.0440815873,0.5222098769,0.6553600784,0.6853486744,0.0571945074,0.0923124240,0.6976544929,0.9257440316,0.5658967043,0.0636543627,0.1038574059,0.0662497468,0.9165439918,0.5200087291,0.9528015053,0.5347318368,0.1848373057,0.9948602219,0.9633110918,0.1482162909,0.9000614029,0.0898618386,0.7975253051,0.8334557347,0.8629821099,0.0001795699,0.2488384889,0.6382902598,0.1103540359,0.2199716354,0.2737281912,0.5694398067,0.7940423761,0.4906677457,0.5191186895,0.4770589883,0.2823238128,0.2458788699,0.6363522802,0.0306954833,0.6979198116)

f1 <- factor(rep(c("Female","Male","Female","Male","Male","Female"), each=12))
f2 <- c(0.098788608,0.934606288,0.145045152,0.841969882,0.498234471,0.562897249,0.359740488,0.082046687,0.183987342,0.082418820,0.173424633,0.799291329,0.041450568,0.686708743,0.352092230,0.823550310,0.650857094,0.331705763,0.659111451,0.745187314,0.066165065,0.870759966,0.154977488,0.031703774,0.065251788,0.707452073,0.564604314,0.224798417,0.656363138,0.047954841,0.500513114,0.923316812,0.706629266,0.561530974,0.670860932,0.414969178,0.709973062,0.452946384,0.187624344,0.278656351,0.562138433,0.655193272,0.014868182,0.518697012,0.414113229,0.273464316,0.844080831,0.962636550,0.952739605,0.728627219,0.761122951,0.309977150,0.755239042,0.208627128,0.481429897,0.376021223,0.753871400,0.164842337,0.921081061,0.859677311,0.600462073,0.119193708,0.276722102,0.854752641,0.962710853,0.956277061,0.228313179,0.920405764,0.001594131,0.104930433,0.241548888,0.549643015)
f3 <- factor(rep(c("day1","day2","day3","day4"),each=3, times=6))

data <- data.frame(sub=subject, dep=dep, f1=f1, f2=f2, f3=f3)

m <- lmer(dep ~ f1*f2*f3 + (1|sub), data=data)

I use the emmeans package for post-hoc tests and ggplot2 to plot the results. I'm finding some differences between the means calculated by ggplot and the means from emmeans

ggplot(aes(x=f3,y=dep,colour=f1),data=data) + stat_summary(fun.y=mean, geom="point")
emmeans(m, c("f1","f3"))

For example the mean for male in day1 is 0.573, but the emmean is 0.558. Which value best represents my model then?

The data are balanced so I don't understand why the emmeans are different from the ordinary mean. Can someone help me?

  • 3
    $\begingroup$ "Estimated marginal means are based on a model – not directly on data." See here for more information: cran.r-project.org/web/packages/emmeans/vignettes/basics.html $\endgroup$ – Stefan Oct 22 '18 at 19:44
  • $\begingroup$ Thanks @Stefan for your comment. Yes, you are right, I'd forgotten about that. But what exactly does that mean. Is it inaccurate to plot data like in my example with ggplot2 for mixed models? How can I visualize my mixed model effects by plotting the means of day1:male, day1:female etc without needing to calculate estimated marginal means? $\endgroup$ – locus Oct 22 '18 at 19:58
  • $\begingroup$ I would run the model first, then calculate the estimated marginal means via the emmeans() function (incl. standard errors) and take those values for the ggplot call. Alternatively you can also plot emmGrid objects via the emmip() and plot() functions. See also here: cran.r-project.org/web/packages/emmeans/vignettes/… $\endgroup$ – Stefan Oct 22 '18 at 20:06
  • $\begingroup$ Thanks @Stefan, I read those vignettes but I'm still confused why plotting with emmeans is a better representation of my model than simply plot the data. Do you have an opinion about that? I've read many posts in CV and most answers use functions like fixef, coef and predict to plot the data of a mixed model, although I'm not completely sure I understand how to do it. $\endgroup$ – locus Oct 22 '18 at 20:18
  • 2
    $\begingroup$ Since you fitted a model to your data, you want to have means (or predictions) based on the model fit and not your raw data. fixef and coef will extract model-specific values. You can also use the predict function (and a bit more) to calculate emmeans yourself. Have a look at the answer here and the comments: stats.stackexchange.com/questions/253756/… $\endgroup$ – Stefan Oct 22 '18 at 20:50

The data are balanced, but one of your “factors,” namely f2, is actually a covariate, and your data are not balanced with respect to the values of that variable. If you were to take f2 out of the model, I think the EMMs will match the ordinary marginal means, unless I’m overlooking something. But with f2 in the model, the EMMs are like adjusted means in the analysis of covariance, based on predictions when f2 is set at its mean value.

Another way to make them equal is to put cov.adjust = list(f2~f1*f3) in the emmeans() call. This sets it up so that predictions are made with f2 set at its cell means for those two factors. Such EMMs are appropriate when f2 is in the causal path, i.e. it is a mediating variable. That situation is discussed in the vignettes too, I think the one on messy data (could look it up in the vignette-topics index).


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.