# Calculating confidence intervals of marginal means in linear mixed models

I'm using different R packages (effects, ggeffects, emmeans, lmer) to calculate confidence intervals of marginal means in a linear mixed model. My problem is that the effects package produces smaller CIs compared to other methods. Here is an example:

library(effects)
library(ggeffects)
library(emmeans)
library(lmerTest)
library(ggplot2)

ham$$Age_scale <- scale(ham$$Age)

# fit model without the intercept
fit <- lmer(Informed.liking ~ -1 + Product + Age_scale +
(1|Consumer/Product), data=ham)


Model summary:

Random effects:
Groups           Name        Variance Std.Dev.
Product:Consumer (Intercept) 3.1464   1.7738
Consumer         (Intercept) 0.3908   0.6252
Residual                     1.6991   1.3035
Number of obs: 648, groups:  Product:Consumer, 324; Consumer, 81

Fixed effects:
Estimate Std. Error        df t value Pr(>|t|)
Product1  5.809e+00  2.327e-01 3.110e+02  24.960   <2e-16 ***
Product2  5.105e+00  2.327e-01 3.110e+02  21.936   <2e-16 ***
Product3  6.093e+00  2.327e-01 3.110e+02  26.180   <2e-16 ***
Product4  5.926e+00  2.327e-01 3.110e+02  25.464   <2e-16 ***
Age_scale 9.621e-03  1.311e-01 7.900e+01   0.073    0.942


Then, I used different methods to calculate the CIs for the marginal means for Product:

term = 'Product'

## ggpredict
c0 <- as.data.frame(ggpredict(fit, terms = term))
c0 <- c0[,4:5]

## confint.merMod
c1 <- confint(fit, method='profile')
c1 <- c1[4:7,]

## confint.merMod
c2 <- confint(fit,method='Wald')
c2 <- c2[4:7,]

## confint.merMod
c3 <- confint(fit,method='boot')
c3 <- c3[4:7,]

## effect
c4 <- with(effect(term,fit),cbind(lower,upper))

## emmeans,'kenward-roger'
c5 <- with(summary(emmeans(fit,spec=term)),cbind(lower.CL,upper.CL))


I put all the CIs together:

tmpf <- function(method,val) {
data.frame(method=method,
v=LETTERS[1:4],
setNames(as.data.frame(tail(val,4)),
c("lwr","upr")))
}

allCI <- rbind(tmpf("ggpredict",c0),
tmpf('profile',c1),
tmpf('wald',c2),
tmpf('boot',c3),
tmpf("effects",c4),
tmpf("emmeans",c5))

ggplot(allCI,aes(v,ymin=lwr,ymax=upr,colour=method))+
geom_linerange(position=position_dodge(width=0.8),size=1) + theme_bw()


The confidence intervals produced by the effects package are substantially smaller than others. I noticed someone asked a similar question before, but its answers showed that CIs produced byeffects are very similar to others. I wonder what's going on here. Did I miss something?

UPDATE: @Daniel ran the exact same code and found no deviation for the effects package. ggeffect() and ggemmeans() from version 0.8.0 of the ggeffects package also produced very similar CIs.

However, I still got smaller CIs for effects (version 4.0.1) and ggeffects (version 0.7.0), at .95 confidence level.

## effect
eff = effect(term,fit)
c4 <- with(eff,cbind(lower,upper))

## emmeans,'kenward-roger'
c5 <- with(summary(emmeans(fit,spec=term)),cbind(lower.CL,upper.CL))

## ggeffect
c6 <- ggeffect(fit, terms = term)
c6 <- c6[,4:5]

packageVersion('ggeffects')
‘0.7.0’
packageVersion('effects')
‘4.0.1’
eff\$confidence.level
0.95


packageVersion('ggeffects')
‘0.8.0’
packageVersion('effects')
‘4.1.0’

## effect
eff = effect(term,fit)
c4 <- with(eff,cbind(lower,upper))

## emmeans,'kenward-roger'
c5 <- with(summary(emmeans(fit,spec=term)),cbind(lower.CL,upper.CL))

## ggeffect
c6 <- ggeffect(fit, terms = term)
c6 <- c6[,4:5]

## ggemmeans (only the 0.8.0 version supports this)
c7 <- ggemmeans(fit, terms = term)
c7 <- c7[,4:5]


• Good to know where this error comes from. The smaller CI for ggeffects 0.7.0 in your example seem plausible (and would also have occured with ggeffects 0.8.0), because the ggeffect() function simply wraps effect() - and if the source of this "bug" was the older version 4.0.1 of the effect-package, ggeffect() will also have this "bug". – Daniel Jan 10 '19 at 20:24

I have just checked your example (using the just released version 0.8.0 of ggeffects), however, instead of using effect() or emmeans() directly, I use the functions from ggeffects, which actually wrap around these functions (ggeffect() and ggemmeans()).

## effect
c4 <- ggeffect(fit, terms = term)
c4 <- c4[,4:5]

## emmeans,'kenward-roger'
c5 <- ggemmeans(fit, terms = term)
c5 <- c5[,4:5]


Here, the CI's are similar.

Also, after exactly running your examples, the plot is identical, with no deviation for effects.

## effect
c4 <- with(effect(term,fit),cbind(lower,upper))

## emmeans,'kenward-roger'
c5 <- with(summary(emmeans(fit,spec=term)),cbind(lower.CL,upper.CL))


Have you set some options for the effects-package that might change the ci-level?

It seems to me that all these methods produce confidence intervals, not prediction intervals. If you use type = "re" in ggpredict(), the intervals are much larger (except for this model, because the variance in the random effects is almost zero).

• Thanks for your answer. I found that updating the effects package from 4.0.1 to 4.1.0 resolved the problem. – Han Zhang Jan 10 '19 at 19:42