I'm using a multiply imputed dataset (mids) to run a linear mixed model in R, and then I tried the emmeans package to obtain my emmeans with p-values. X2 is a factor which I'm interested in, that was multiply imputed as a factor, 0,1.
library(foreign)
library(mice)
library(broom.mixed)
library(lme4)
library(emmeans)
Some of my df1 data (dput):
structure(list(imp = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5), Subj = c(1, 1, 1, 2, 2, 2,
3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9,
10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 14, 14, 14, 15,
15, 15, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20,
20, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7,
7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13,
13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 16, 17, 17, 17, 18, 18,
18, 19, 19, 19, 20, 20, 20, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4,
4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11,
11, 11, 12, 12, 12, 13, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16,
16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 1, 1, 1,
2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8,
9, 9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 14,
14, 14, 15, 15, 15, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19,
19, 20, 20, 20, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5,
6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 12,
12, 12, 13, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 16, 17, 17,
17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 1, 1, 1, 2, 2, 2, 3,
3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10,
10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 14, 14, 14, 15, 15,
15, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20
), X1 = c(0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1,
3, 0, 1, 3, 0, 1, 3), X2 = c(1, 1, 1, 1, 1, 1, NA, NA, NA, 1,
1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, NA, NA, NA, 1, 1, 1, NA, NA, NA, 0, 0, 0, 0, 0, 0, 1, 1,
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1,
1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0,
0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1,
1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0,
0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0), X3 = c(28, 28,
28, 37, 37, 37, 36, 36, 36, 57, 57, 57, 34, 34, 34, 53, 53, 53,
42, 42, 42, 25, 25, 25, 29, 29, 29, 21, 21, 21, 52, 52, 52, 23,
23, 23, 39, 39, 39, 49, 49, 49, 39, 39, 39, 30, 30, 30, 65, 65,
65, 35, 35, 35, 34, 34, 34, 42, 42, 42, 28, 28, 28, 37, 37, 37,
36, 36, 36, 57, 57, 57, 34, 34, 34, 53, 53, 53, 42, 42, 42, 25,
25, 25, 29, 29, 29, 21, 21, 21, 52, 52, 52, 23, 23, 23, 39, 39,
39, 49, 49, 49, 39, 39, 39, 30, 30, 30, 65, 65, 65, 35, 35, 35,
34, 34, 34, 42, 42, 42, 28, 28, 28, 37, 37, 37, 36, 36, 36, 57,
57, 57, 34, 34, 34, 53, 53, 53, 42, 42, 42, 25, 25, 25, 29, 29,
29, 21, 21, 21, 52, 52, 52, 23, 23, 23, 39, 39, 39, 49, 49, 49,
39, 39, 39, 30, 30, 30, 65, 65, 65, 35, 35, 35, 34, 34, 34, 42,
42, 42, 28, 28, 28, 37, 37, 37, 36, 36, 36, 57, 57, 57, 34, 34,
34, 53, 53, 53, 42, 42, 42, 25, 25, 25, 29, 29, 29, 21, 21, 21,
52, 52, 52, 23, 23, 23, 39, 39, 39, 49, 49, 49, 39, 39, 39, 30,
30, 30, 65, 65, 65, 35, 35, 35, 34, 34, 34, 42, 42, 42, 28, 28,
28, 37, 37, 37, 36, 36, 36, 57, 57, 57, 34, 34, 34, 53, 53, 53,
42, 42, 42, 25, 25, 25, 29, 29, 29, 21, 21, 21, 52, 52, 52, 23,
23, 23, 39, 39, 39, 49, 49, 49, 39, 39, 39, 30, 30, 30, 65, 65,
65, 35, 35, 35, 34, 34, 34, 42, 42, 42, 28, 28, 28, 37, 37, 37,
36, 36, 36, 57, 57, 57, 34, 34, 34, 53, 53, 53, 42, 42, 42, 25,
25, 25, 29, 29, 29, 21, 21, 21, 52, 52, 52, 23, 23, 23, 39, 39,
39, 49, 49, 49, 39, 39, 39, 30, 30, 30, 65, 65, 65, 35, 35, 35,
34, 34, 34, 42, 42, 42), X4 = c(2, 2, 2, 1, 1, 1, 1, 1, 1, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1,
1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2,
2, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 2,
2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2,
2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1,
1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1,
1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1,
1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 2,
2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1,
1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2,
2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1,
1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 2,
2, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1,
1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1), X5 = c(2, 2, 2, 1,
1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1,
1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2,
2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1,
1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1,
1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1,
1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
Dependant = c(32, 11, 25, 30, 10, 8, 34, NA, NA, 16, 20,
20.5, 13, 14, 15, 16, 17, 12, 21, NA, 20, 41, 25, 26, 24,
28, 33, 16, 6, 3, 7, 0, NA, 17, 3, NA, 40, 25, 21, 35, NA,
NA, 17, 5, 4, 27, 8, 4, 22, 20, NA, 18, 23, 10, 13, 14, 4,
15, 12, 14, 32, 11, 25, 30, 10, 8, 34, 22.74679565, 9.078754233,
16, 20, 20.5, 13, 14, 15, 16, 17, 12, 21, 23.89331785, 20,
41, 25, 26, 24, 28, 33, 16, 6, 3, 7, 0, 9.563034676, 17,
3, 8.171110057, 40, 25, 21, 35, 13.53468633, 21.4990267,
17, 5, 4, 27, 8, 4, 22, 20, 10.75448127, 18, 23, 10, 13,
14, 4, 15, 12, 14, 32, 11, 25, 30, 10, 8, 34, 23.678318,
6.38024107, 16, 20, 20.5, 13, 14, 15, 16, 17, 12, 21, 15.33803817,
20, 41, 25, 26, 24, 28, 33, 16, 6, 3, 7, 0, 26.37095017,
17, 3, 5.8830194, 40, 25, 21, 35, 17.5314952, 13.24750912,
17, 5, 4, 27, 8, 4, 22, 20, 16.64639694, 18, 23, 10, 13,
14, 4, 15, 12, 14, 32, 11, 25, 30, 10, 8, 34, 24.23722065,
14.13450552, 16, 20, 20.5, 13, 14, 15, 16, 17, 12, 21, 21.8956743,
20, 41, 25, 26, 24, 28, 33, 16, 6, 3, 7, 0, 13.31806771,
17, 3, 12.46910086, 40, 25, 21, 35, 5.346067919, 14.71408964,
17, 5, 4, 27, 8, 4, 22, 20, 9.644988209, 18, 23, 10, 13,
14, 4, 15, 12, 14, 32, 11, 25, 30, 10, 8, 34, 16.25464311,
14.55203722, 16, 20, 20.5, 13, 14, 15, 16, 17, 12, 21, 20.15197466,
20, 41, 25, 26, 24, 28, 33, 16, 6, 3, 7, 0, 7.169741223,
17, 3, 15.89606111, 40, 25, 21, 35, 8.178602917, -2.653330003,
17, 5, 4, 27, 8, 4, 22, 20, 17.11276677, 18, 23, 10, 13,
14, 4, 15, 12, 14, 32, 11, 25, 30, 10, 8, 34, 28.30217443,
7.983293096, 16, 20, 20.5, 13, 14, 15, 16, 17, 12, 21, -2.445506679,
20, 41, 25, 26, 24, 28, 33, 16, 6, 3, 7, 0, 0.897162825,
17, 3, 12.38497675, 40, 25, 21, 35, 13.76105932, 23.24251938,
17, 5, 4, 27, 8, 4, 22, 20, 9.744084982, 18, 23, 10, 13,
14, 4, 15, 12, 14)), class = "data.frame", row.names = c(NA,
-360L), variable.labels = c(imp = "", Subj = "", X1 = "Temps",
X2 = "PsychoT1", X3 = "T0_age", X4 = "Group", X5 = "Sexe", Dependant = "Dependant"
), codepage = 65001L)
df.mids1<-as.mids(long=df1, .imp="imp")
yPsy1 <- with(data = df.mids1, exp = lme4::lmer(Dependant ~ factor(X1) +
factor(X2) +
X3 +
factor(X4) +
factor(X5) +
factor(X1) * factor(X2) +
factor(X2) * factor(X4) +
factor(X1) * factor(X4) +
factor(X1) * factor(X2) * factor(X4) +
(1| Subj)), reml=TRUE)
mice.res <- summary(pool(yPsy1))
mice.res
#> term estimate std.error statistic
#> 1 (Intercept) 25.9340792 7.8864719 3.28842599
#> 2 factor(X1)1 -11.0168803 6.1641328 -1.78725552
#> 3 factor(X1)3 -12.8614164 6.5151402 -1.97408128
#> 4 factor(X2)1 7.6964991 7.2210293 1.06584516
#> 5 X3 -0.1129837 0.1455105 -0.77646422
#> 6 factor(X4)2 -6.6860358 7.3265611 -0.91257490
#> 7 factor(X5)2 6.7636729 4.4682356 1.51372343
#> 8 factor(X1)1:factor(X2)1 -2.5189278 7.1551292 -0.35204505
#> 9 factor(X1)3:factor(X2)1 -1.8445553 7.8818619 -0.23402533
#> 10 factor(X2)1:factor(X4)2 -4.2784766 8.7986350 -0.48626595
#> 11 factor(X1)1:factor(X4)2 8.7468803 7.6738824 1.13982465
#> 12 factor(X1)3:factor(X4)2 6.1926579 7.7724684 0.79674276
#> 13 factor(X1)1:factor(X2)1:factor(X4)2 -0.6267509 9.4682741 -0.06619485
#> 14 factor(X1)3:factor(X2)1:factor(X4)2 6.0615966 9.6874709 0.62571507
#> df p.value
#> 1 36.93801 0.00221774
#> 2 42.12344 0.08109302
#> 3 34.16544 0.05650184
#> 4 31.15342 0.29468469
#> 5 34.29064 0.44279816
#> 6 36.25076 0.36749560
#> 7 41.83535 0.13761496
#> 8 40.67191 0.72662209
#> 9 25.35173 0.81684491
#> 10 36.72275 0.62966797
#> 11 39.75094 0.26118111
#> 12 37.82968 0.43057262
#> 13 34.49832 0.94760483
#> 14 30.61215 0.53613943
emmeans(yPsy1, pairwise ~ X1|X2|X4)
#> $emmeans
#> X2 = 0.0, X4 = 1:
#> X1 emmean SE df lower.CL upper.CL
#> 0 25.0 6.58 18.0 11.152 38.8
#> 1 13.9 6.56 18.1 0.172 27.7
#> 3 12.1 6.18 21.7 -0.732 24.9
#>
#> X2 = 0.6, X4 = 1:
#> X1 emmean SE df lower.CL upper.CL
#> 0 25.0 6.58 18.0 11.152 38.8
#> 1 13.9 6.56 18.1 0.172 27.7
#> 3 12.1 6.18 21.7 -0.732 24.9
#>
#> X2 = 0.8, X4 = 1:
#> X1 emmean SE df lower.CL upper.CL
#> 0 25.0 6.58 18.0 11.152 38.8
#> 1 13.9 6.56 18.1 0.172 27.7
#> 3 12.1 6.18 21.7 -0.732 24.9
#>
#> X2 = 1.0, X4 = 1:
#> X1 emmean SE df lower.CL upper.CL
#> 0 32.7 3.65 21.4 25.071 40.3
#> 1 19.1 3.85 18.1 11.037 27.2
#> 3 18.0 3.97 16.1 9.540 26.4
#>
#> X2 = 0.0, X4 = 2:
#> X1 emmean SE df lower.CL upper.CL
#> 0 18.3 4.68 20.9 8.540 28.0
#> 1 16.0 5.01 16.9 5.437 26.6
#> 3 11.6 4.68 20.8 1.869 21.4
#>
#> X2 = 0.6, X4 = 2:
#> X1 emmean SE df lower.CL upper.CL
#> 0 18.3 4.68 20.9 8.540 28.0
#> 1 16.0 5.01 16.9 5.437 26.6
#> 3 11.6 4.68 20.8 1.869 21.4
#>
#> X2 = 0.8, X4 = 2:
#> X1 emmean SE df lower.CL upper.CL
#> 0 18.3 4.68 20.9 8.540 28.0
#> 1 16.0 5.01 16.9 5.437 26.6
#> 3 11.6 4.68 20.8 1.869 21.4
#>
#> X2 = 1.0, X4 = 2:
#> X1 emmean SE df lower.CL upper.CL
#> 0 21.7 3.42 23.1 14.630 28.8
#> 1 16.3 3.72 17.3 8.434 24.1
#> 3 19.2 3.50 21.7 11.981 26.5
#>
#> Results are averaged over the levels of: X5
#> Confidence level used: 0.95
#>
#> $contrasts
#> X2 = 0.0, X4 = 1:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 11.02 6.16 30.2 1.787 0.1909
#> X10 - X13 12.86 6.52 25.1 1.974 0.1394
#> X11 - X13 1.84 6.50 25.4 0.284 0.9567
#>
#> X2 = 0.6, X4 = 1:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 11.02 6.16 30.2 1.787 0.1909
#> X10 - X13 12.86 6.52 25.1 1.974 0.1394
#> X11 - X13 1.84 6.50 25.4 0.284 0.9567
#>
#> X2 = 0.8, X4 = 1:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 11.02 6.16 30.2 1.787 0.1909
#> X10 - X13 12.86 6.52 25.1 1.974 0.1394
#> X11 - X13 1.84 6.50 25.4 0.284 0.9567
#>
#> X2 = 1.0, X4 = 1:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 13.54 3.63 25.1 3.725 0.0028
#> X10 - X13 14.71 3.92 17.7 3.754 0.0040
#> X11 - X13 1.17 3.86 18.9 0.303 0.9508
#>
#> X2 = 0.0, X4 = 2:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 2.27 4.57 24.7 0.497 0.8733
#> X10 - X13 6.67 4.34 29.6 1.537 0.2885
#> X11 - X13 4.40 4.76 20.4 0.924 0.6320
#>
#> X2 = 0.6, X4 = 2:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 2.27 4.57 24.7 0.497 0.8733
#> X10 - X13 6.67 4.34 29.6 1.537 0.2885
#> X11 - X13 4.40 4.76 20.4 0.924 0.6320
#>
#> X2 = 0.8, X4 = 2:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 2.27 4.57 24.7 0.497 0.8733
#> X10 - X13 6.67 4.34 29.6 1.537 0.2885
#> X11 - X13 4.40 4.76 20.4 0.924 0.6320
#>
#> X2 = 1.0, X4 = 2:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 5.42 3.83 21.4 1.412 0.3525
#> X10 - X13 2.45 3.56 28.9 0.689 0.7720
#> X11 - X13 -2.96 3.90 19.8 -0.760 0.7311
#>
#> Results are averaged over the levels of: X5
#> P value adjustment: tukey method for comparing a family of 3 estimates
I've had no problem obtaining the model coefficient + associated stats from the mice package using summary(pool(yPsy1)). Everything looks good.
I ran this with my Large Mira object (obtained with the previous model) to obtain my emmeans:
emmeans(yPsy1, pairwise ~ X1|X2|X4)
THE PROBLEM: Although X2 is a factor (0, 1), emmeans gives me additional levels (0, 0.6, 0.8, 1). There are even more levels when I considered my 25 imputes datasets.
This variable was imputed, but I inspected my 25 imputed dataset (df.mids) and can't see these values. I also checked the individual 25 frames in my Mira model element, and all I see is 0 and 1. Also, the mixed linear model outputs give me only 2 levels, which is normal.
Could you tell me where I went wrong and how to fix this?
I was redirected from Cross Validated where someone suggested that emmeans may be considering X2 as a numeric value. So, rather then specifying as.factor in my lmer(), I should try modifying my variable in my imputed dataframe. So I did :
Same initial data set
df.mids2<-as.mids(long=df2, .imp="imp")
df.mids2$X2<-as.factor(df.mids2$X2)
yPsy2 <- with(data = df.mids2, exp = lme4::lmer(Dependant ~ factor(X1) + factor(X2) + X3 +factor(X4) +factor(X5) + factor(X1) * factor(X2) +
factor(X2) * factor(X4) + factor(X1) * factor(X4) +
factor(X1) * factor(X2) * factor(X4) + (1| Subj)), reml=TRUE)
mice.res <- summary(pool(yPsy2))
mice.res
#> term estimate std.error statistic
#> 1 (Intercept) 25.9340792 7.8864719 3.28842599
#> 2 factor(X1)1 -11.0168803 6.1641328 -1.78725552
#> 3 factor(X1)3 -12.8614164 6.5151402 -1.97408128
#> 4 factor(X2)1 7.6964991 7.2210293 1.06584516
#> 5 X3 -0.1129837 0.1455105 -0.77646422
#> 6 factor(X4)2 -6.6860358 7.3265611 -0.91257490
#> 7 factor(X5)2 6.7636729 4.4682356 1.51372343
#> 8 factor(X1)1:factor(X2)1 -2.5189278 7.1551292 -0.35204505
#> 9 factor(X1)3:factor(X2)1 -1.8445553 7.8818619 -0.23402533
#> 10 factor(X2)1:factor(X4)2 -4.2784766 8.7986350 -0.48626595
#> 11 factor(X1)1:factor(X4)2 8.7468803 7.6738824 1.13982465
#> 12 factor(X1)3:factor(X4)2 6.1926579 7.7724684 0.79674276
#> 13 factor(X1)1:factor(X2)1:factor(X4)2 -0.6267509 9.4682741 -0.06619485
#> 14 factor(X1)3:factor(X2)1:factor(X4)2 6.0615966 9.6874709 0.62571507
#> df p.value
#> 1 36.93801 0.00221774
#> 2 42.12344 0.08109302
#> 3 34.16544 0.05650184
#> 4 31.15342 0.29468469
#> 5 34.29064 0.44279816
#> 6 36.25076 0.36749560
#> 7 41.83535 0.13761496
#> 8 40.67191 0.72662209
#> 9 25.35173 0.81684491
#> 10 36.72275 0.62966797
#> 11 39.75094 0.26118111
#> 12 37.82968 0.43057262
#> 13 34.49832 0.94760483
#> 14 30.61215 0.53613943
emmeans(yPsy2, pairwise ~ X1|X2|X4)
#> $emmeans
#> X2 = 0.0, X4 = 1:
#> X1 emmean SE df lower.CL upper.CL
#> 0 25.0 6.58 18.0 11.152 38.8
#> 1 13.9 6.56 18.1 0.172 27.7
#> 3 12.1 6.18 21.7 -0.732 24.9
#>
#> X2 = 0.6, X4 = 1:
#> X1 emmean SE df lower.CL upper.CL
#> 0 25.0 6.58 18.0 11.152 38.8
#> 1 13.9 6.56 18.1 0.172 27.7
#> 3 12.1 6.18 21.7 -0.732 24.9
#>
#> X2 = 0.8, X4 = 1:
#> X1 emmean SE df lower.CL upper.CL
#> 0 25.0 6.58 18.0 11.152 38.8
#> 1 13.9 6.56 18.1 0.172 27.7
#> 3 12.1 6.18 21.7 -0.732 24.9
#>
#> X2 = 1.0, X4 = 1:
#> X1 emmean SE df lower.CL upper.CL
#> 0 32.7 3.65 21.4 25.071 40.3
#> 1 19.1 3.85 18.1 11.037 27.2
#> 3 18.0 3.97 16.1 9.540 26.4
#>
#> X2 = 0.0, X4 = 2:
#> X1 emmean SE df lower.CL upper.CL
#> 0 18.3 4.68 20.9 8.540 28.0
#> 1 16.0 5.01 16.9 5.437 26.6
#> 3 11.6 4.68 20.8 1.869 21.4
#>
#> X2 = 0.6, X4 = 2:
#> X1 emmean SE df lower.CL upper.CL
#> 0 18.3 4.68 20.9 8.540 28.0
#> 1 16.0 5.01 16.9 5.437 26.6
#> 3 11.6 4.68 20.8 1.869 21.4
#>
#> X2 = 0.8, X4 = 2:
#> X1 emmean SE df lower.CL upper.CL
#> 0 18.3 4.68 20.9 8.540 28.0
#> 1 16.0 5.01 16.9 5.437 26.6
#> 3 11.6 4.68 20.8 1.869 21.4
#>
#> X2 = 1.0, X4 = 2:
#> X1 emmean SE df lower.CL upper.CL
#> 0 21.7 3.42 23.1 14.630 28.8
#> 1 16.3 3.72 17.3 8.434 24.1
#> 3 19.2 3.50 21.7 11.981 26.5
#>
#> Results are averaged over the levels of: X5
#> Confidence level used: 0.95
#>
#> $contrasts
#> X2 = 0.0, X4 = 1:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 11.02 6.16 30.2 1.787 0.1909
#> X10 - X13 12.86 6.52 25.1 1.974 0.1394
#> X11 - X13 1.84 6.50 25.4 0.284 0.9567
#>
#> X2 = 0.6, X4 = 1:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 11.02 6.16 30.2 1.787 0.1909
#> X10 - X13 12.86 6.52 25.1 1.974 0.1394
#> X11 - X13 1.84 6.50 25.4 0.284 0.9567
#>
#> X2 = 0.8, X4 = 1:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 11.02 6.16 30.2 1.787 0.1909
#> X10 - X13 12.86 6.52 25.1 1.974 0.1394
#> X11 - X13 1.84 6.50 25.4 0.284 0.9567
#>
#> X2 = 1.0, X4 = 1:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 13.54 3.63 25.1 3.725 0.0028
#> X10 - X13 14.71 3.92 17.7 3.754 0.0040
#> X11 - X13 1.17 3.86 18.9 0.303 0.9508
#>
#> X2 = 0.0, X4 = 2:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 2.27 4.57 24.7 0.497 0.8733
#> X10 - X13 6.67 4.34 29.6 1.537 0.2885
#> X11 - X13 4.40 4.76 20.4 0.924 0.6320
#>
#> X2 = 0.6, X4 = 2:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 2.27 4.57 24.7 0.497 0.8733
#> X10 - X13 6.67 4.34 29.6 1.537 0.2885
#> X11 - X13 4.40 4.76 20.4 0.924 0.6320
#>
#> X2 = 0.8, X4 = 2:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 2.27 4.57 24.7 0.497 0.8733
#> X10 - X13 6.67 4.34 29.6 1.537 0.2885
#> X11 - X13 4.40 4.76 20.4 0.924 0.6320
#>
#> X2 = 1.0, X4 = 2:
#> contrast estimate SE df t.ratio p.value
#> X10 - X11 5.42 3.83 21.4 1.412 0.3525
#> X10 - X13 2.45 3.56 28.9 0.689 0.7720
#> X11 - X13 -2.96 3.90 19.8 -0.760 0.7311
#>
#> Results are averaged over the levels of: X5
#> P value adjustment: tukey method for comparing a family of 3 estimates
I also tried emmeans(yPsy1, pairwise ~ factor(X1)|factor(X2)|factor(X4)) and the results did not change.
Simply convert the variables to factors before fitting the model:
df1 <- datawizard::to_factor(df1, c("X1", "X2", "X4", "X5"))
# impute missing values
df.mids1 <- as.mids(long = df1, .imp = "imp")
yPsy2 <- with(data = df.mids1, exp = lme4::lmer(Dependant ~ X1 + X2 + X3 + X4 + X5 +
X1 * X2 + X2 * X4 + X1 * X4 + X1 * X2 * X4 + (1 | Subj)))
mice.res <- summary(pool(yPsy2))
mice.res
emmeans::emmeans(yPsy2, pairwise ~ X1 | X2 | X4)
A "longer" solution, which requires a hack: the ggeffects packages automatically detects the variable type used in the model and (should) treat it as intended:
library(ggeffects)
library(easystats)
library(lme4)
library(mice)
# ... data set
# impute missing values
df.mids1 <- as.mids(long = df1, .imp = "imp")
# fit models
models <- lapply(1:5, function(i) {
lme4::lmer(Dependant ~ factor(X1) + factor(X2) +
X3 +
factor(X4) +
factor(X5) +
factor(X1) * factor(X2) +
factor(X2) * factor(X4) +
factor(X1) * factor(X4) +
factor(X1) * factor(X2) * factor(X4) +
(1 | Subj), data = mice::complete(df.mids1, action = i))
})
# pooled estimates
pool_parameters(models)
#> # Fixed Effects
#>
#> Parameter | Coefficient | SE | 95% CI | Statistic | df | p
#> -----------------------------------------------------------------------------------------------------
#> (Intercept) | 25.93 | 7.89 | [ 9.95, 41.91] | 3.29 | 36.94 | 0.002
#> X1 [1] | -11.02 | 6.16 | [-23.46, 1.42] | -1.79 | 42.12 | 0.081
#> X1 [3] | -12.86 | 6.52 | [-26.10, 0.38] | -1.97 | 34.17 | 0.057
#> X2 [1] | 7.70 | 7.22 | [ -7.03, 22.42] | 1.07 | 31.15 | 0.295
#> X3 | -0.11 | 0.15 | [ -0.41, 0.18] | -0.78 | 34.29 | 0.443
#> X4 [2] | -6.69 | 7.33 | [-21.54, 8.17] | -0.91 | 36.25 | 0.367
#> X5 [2] | 6.76 | 4.47 | [ -2.25, 15.78] | 1.51 | 41.84 | 0.138
#> X1 [1] × factor(X2)1 | -2.52 | 7.16 | [-16.97, 11.93] | -0.35 | 40.67 | 0.727
#> X1 [3] × factor(X2)1 | -1.84 | 7.88 | [-18.07, 14.38] | -0.23 | 25.35 | 0.817
#> X2 [1] × factor(X4)2 | -4.28 | 8.80 | [-22.11, 13.55] | -0.49 | 36.72 | 0.630
#> X1 [1] × factor(X4)2 | 8.75 | 7.67 | [ -6.77, 24.26] | 1.14 | 39.75 | 0.261
#> X1 [3] × factor(X4)2 | 6.19 | 7.77 | [ -9.54, 21.93] | 0.80 | 37.83 | 0.431
#> X1 [1] × factor(X2)1 × factor(X4)2 | -0.63 | 9.47 | [-19.86, 18.60] | -0.07 | 34.50 | 0.948
#> X1 [3] × factor(X2)1 × factor(X4)2 | 6.06 | 9.69 | [-13.71, 25.83] | 0.63 | 30.61 | 0.536
#>
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#> using a Wald distribution approximation.
# do the same for predictions...
predictions <- lapply(1:5, function(i) {
m <- lme4::lmer(Dependant ~ factor(X1) + factor(X2) +
X3 +
factor(X4) +
factor(X5) +
factor(X1) * factor(X2) +
factor(X2) * factor(X4) +
factor(X1) * factor(X4) +
factor(X1) * factor(X2) * factor(X4) +
(1 | Subj), data = mice::complete(df.mids1, action = i))
predict_response(m, c("X1", "X2", "X4"), margin = "marginalmeans")
})
# pooled estimated marginal means
pooled_pr <- pool_predictions(predictions)
pooled_pr
#> # Predicted values of Dependant
#>
#> X2: 0
#> X4: 1
#>
#> X1 | Predicted | 95% CI
#> -----------------------------
#> 0 | 24.97 | 12.08, 37.85
#> 1 | 13.95 | 1.09, 26.81
#> 3 | 12.10 | -0.02, 24.23
#>
#> X2: 0
#> X4: 2
#>
#> X1 | Predicted | 95% CI
#> -----------------------------
#> 0 | 18.28 | 9.10, 27.46
#> 1 | 16.01 | 6.19, 25.83
#> 3 | 11.61 | 2.43, 20.79
#>
#> X2: 1
#> X4: 1
#>
#> X1 | Predicted | 95% CI
#> -----------------------------
#> 0 | 32.66 | 25.50, 39.82
#> 1 | 19.13 | 11.58, 26.68
#> 3 | 17.96 | 10.17, 25.74
#>
#> X2: 1
#> X4: 2
#>
#> X1 | Predicted | 95% CI
#> -----------------------------
#> 0 | 21.70 | 15.00, 28.40
#> 1 | 16.28 | 8.98, 23.58
#> 3 | 19.25 | 12.39, 26.11
#>
#> Adjusted for:
#> * X3 = 38.50
# This is a hack! The `test_predictions` function needs information from the
# original model object (like, linear model, mixed model etc.), which is not
# available in the pooled object. Thus, we simply fit a single model again
# to have "m" in the environment - "m" is the name we used in the loop above.
m <- lme4::lmer(Dependant ~ factor(X1) + factor(X2) +
X3 +
factor(X4) +
factor(X5) +
factor(X1) * factor(X2) +
factor(X2) * factor(X4) +
factor(X1) * factor(X4) +
factor(X1) * factor(X2) * factor(X4) +
(1 | Subj), data = mice::complete(df.mids1, action = 1))
# test predictions (pairwise comparisons), grouped by X2 and X4
test_predictions(pooled_pr, by = c("X2", "X4"), engine = "ggeffects")
#> # Pairwise comparisons
#>
#> X2/X4 = 0.1
#>
#> X1 | Contrast | 95% CI | p
#> --------------------------------------
#> 0-1 | 11.02 | -7.70, 29.74 | 0.242
#> 0-3 | 12.86 | -5.33, 31.05 | 0.161
#> 1-3 | 1.84 | -16.33, 20.02 | 0.839
#>
#> X2/X4 = 1.1
#>
#> X1 | Contrast | 95% CI | p
#> --------------------------------------
#> 0-1 | 13.54 | 2.84, 24.24 | 0.014
#> 0-3 | 14.71 | 3.83, 25.58 | 0.009
#> 1-3 | 1.17 | -9.98, 12.32 | 0.833
#>
#> X2/X4 = 0.2
#>
#> X1 | Contrast | 95% CI | p
#> --------------------------------------
#> 0-1 | 2.27 | -11.55, 16.09 | 0.742
#> 0-3 | 6.67 | -6.68, 20.01 | 0.319
#> 1-3 | 4.40 | -9.42, 18.22 | 0.525
#>
#> X2/X4 = 1.2
#>
#> X1 | Contrast | 95% CI | p
#> --------------------------------------
#> 0-1 | 5.42 | -4.77, 15.60 | 0.290
#> 0-3 | 2.45 | -7.41, 12.31 | 0.619
#> 1-3 | -2.96 | -13.27, 7.34 | 0.565
Created on 2024-09-02 with reprex v2.1.1