Plotting an interaction

Author

Matti Vuorre

library(tidyverse)
library(emmeans)
library(lme4)

set.seed(999)
n <- 1000
dat <- tibble(
  id = sample(letters, n, TRUE),
  x = rnorm(n),
  y = rnorm(n),
  z = rnorm(n),
  d = rnorm(n, x + y + z + x*y)
)
fit <- lmer(d ~ x + x*y*z + (1 | id), data = dat)
boundary (singular) fit: see help('isSingular')
summary(fit)
Linear mixed model fit by REML ['lmerMod']
Formula: d ~ x + x * y * z + (1 | id)
   Data: dat

REML criterion at convergence: 2912

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1910 -0.6753 -0.0383  0.6705  3.3373 

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 0.000    0.000   
 Residual             1.043    1.021   
Number of obs: 1000, groups:  id, 26

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.01612    0.03257   0.495
x            1.01769    0.03227  31.536
y            1.04441    0.03181  32.835
z            0.95667    0.03243  29.502
x:y          0.99425    0.03136  31.707
x:z         -0.03213    0.03158  -1.017
y:z         -0.01824    0.03149  -0.579
x:y:z       -0.03755    0.03275  -1.147

Correlation of Fixed Effects:
      (Intr) x      y      z      x:y    x:z    y:z   
x      0.046                                          
y      0.011  0.066                                   
z     -0.013 -0.042 -0.098                            
x:y    0.065 -0.021  0.042  0.023                     
x:z   -0.044 -0.053  0.017  0.072 -0.010              
y:z   -0.095  0.016  0.015  0.007 -0.029  0.130       
x:y:z  0.010 -0.013 -0.034  0.135  0.010  0.037  0.157
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
emm <- emtrends(
  fit, 
  var = "x", 
  by = c("y", "z"), 
  at = list(y = seq(-2, 2, by = .1))
)
Cannot use mode = "kenward-roger" because *pbkrtest* package is not installed
x <- summary(emm)

x %>% 
  ggplot(aes(y, x.trend)) +
  geom_ribbon(aes(ymin = x.trend - SE, ymax = x.trend + SE)) +
  labs(x = "Value of y", y = "Effect of x on d") +
  geom_line()