library(tidyverse)library(emmeans)library(lme4)set.seed(999)n <-1000dat <-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()