3
$\begingroup$

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()

Confidence Intervals for Marginal Means

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

enter image description here

After upgrading to the latest versions, however, all methods produced similar CIs.

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]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ 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". $\endgroup$ – Daniel Jan 10 at 20:24
2
$\begingroup$

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.

enter image description here

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))

enter image description here

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).

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

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.