rlme4interactionsjplotlmertest

Plotting an interaction with confidence intervals from an lme4 or LmerTest model in R

Using dat (found here), I run the following model:

library(lmerTest)

model <- lmerTest::lmer(eval ~ post + ess + post*ess + (1|ID), data = dat)

The output of summary(model) indicates the interaction term is significant:

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: eval ~ post + ess + ess * post + (1 | ID)
   Data: dat

REML criterion at convergence: 163.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.95714 -0.48596  0.00623  0.49208  1.82729 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 0.33344  0.5774  
 Residual             0.02944  0.1716  
Number of obs: 170, groups:  ID, 85

Fixed effects:
            Estimate Std. Error       df t value             Pr(>|t|)    
(Intercept)  1.50194    0.09082 90.00645  16.538 < 0.0000000000000002 ***
post        -0.24537    0.03658 83.00000  -6.707        0.00000000226 ***
ess          0.15444    0.13076 90.00645   1.181              0.24067    
post:ess     0.15620    0.05267 83.00000   2.965              0.00395 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
         (Intr) post   ess   
post     -0.201              
ess      -0.695  0.140       
post:ess  0.140 -0.695 -0.201

But when I try to plot the interaction using sjPlot using 95% confidence intervals, the resulting intervals don't make the interaction look significant at all...

library(sjPlot)
library(TMB)

plot_model(model, type="int", ci.lvl=0.95)

enter image description here

My two questions:

  1. Why do the estimates and the plotted results appear to tell different stories?
  2. How can I extract confidence intervals for coefficients from the model to create my own graph instead of using plot_model()? I would like to make a bar plot to illustrate the interaction because the variables ess and post are binary.

Note: I'm happy to use lme4 - should get to the same result, it is just not as obvious what coeficients are significant when lme4 objects are summarized and I wanted the question to be very clear.

Solution

  • I'm going to answer your questions in reverse order:

    1. The plot_model() function calls functions from the ggeffects package. Specifically, ggpredict() does a lot of the work. If you go to the following URL, you will find lots of information on how to alter effects plots and getting all sorts of information from the fitted model.

    https://cran.r-project.org/web/packages/ggeffects/vignettes/ggeffects.html

    1. I don't really agree that the interaction does not look significant. The confidence interval doesn't overlap with the other categories mean for the majority of the plot. That might be beside the point, however, as you are currently plotting categorical data in a manner that makes it look continuous. This does not change how the model is being fit under the hood, but it changes the default behaviour of sjPlot. I have fit the model you specified with factors and plotted it how I think it should be plotted below. I don't think this changes much about the plot, but it might change your interpretation. The difference shown here in this plot is consistent with the model summary output and can be stated like this: the difference between the two levels of post is not the same when ess is a 0 vs. a 1. Also, see how the CIs of the two ess categories overlap with the mean when post is 0 but they are significantly different when post is 1.

    Let me know if any of this requires further clarification.

    enter image description here