I am trying to execute linear-mixed models with 3 factors as groups. However running the syntax, I just get 2 terms: time:grup_int2 or time:grup_int3 as you can see in the picture. There should be another term for time:grup_int1. My reference factor group is the level 1, however I don't know how to relevel in the loop, and I intended to post a example with a database provided by the package. The problem is with the database provided by the package works... and I get 3 terms
My code with genes_wide dataframe (plenty of NAs)
#pivot wider to make horizontal data
genes_wide <- genes1 %>% pivot_wider(names_from = gen, values_from = dCt)
genes_wide$time [genes_wide$time == "1"] <- 1
genes_wide$time [genes_wide$time == "3"] <- 3
genes_wide$time<-as.numeric(genes_wide$time)
# if they need to be relevel
#genes_wide$grup_int <- relevel(genes_wide$grup_int, "2")
# names of variables
genes.names <- colnames(genes_wide)[13:67]
no.genes <- length(genes.names)
# create a named list to hold the fitted models
gene.list <- as.list(1:no.genes)
names(gene.list) <- genes.names
# loop over gene names
for(i in genes.names){
# create temporary data matrix and model formula
#1st: adjusting per several variables (edad0 = numeric, grup_int = factor, sexo = factor, time = numeric, id = numeric (need to change?)
#2nd: adjusting just per group
tmp <- genes_wide[, c(i, "grup_int","time", "id", "edad0", "sexo")]
fml <- as.formula( paste( i, "~", paste(c("time", "time:grup_int", "edad0", "sexo"), collapse="+")))
gene.list[[i]] <- lme(fml, random= ~ time|id , method="REML", data=tmp, control = lmeControl(opt = "optim"), na.action = na.omit)
}
# List of data.frames/tibble
result <- lapply(gene.list, function(x) tidy(x))
# Solution 2
result <- dplyr::bind_rows(result, .id = "var")
result2 <- result %>% dplyr::select( var, term, p.value)
**
Example database
** However running the exact syntax with an example database (+ 5 variables to iterate) from the nlme package, the result is different
#adding variables to df problem
BodyWeight$var1 <- rnorm(176, 47 ,3)
BodyWeight$var2 <- rnorm(176, 50 ,10)
BodyWeight$var3 <- rnorm(176, 68, 7)
BodyWeight$var4 <- rnorm(176, 150 ,10)
BodyWeight$var5 <- rnorm(176, 140, 7)
# names of variables
var.names <- colnames(BodyWeight)[5:9]
no.var <- length(var.names)
# create a named list to hold the fitted models
var.list <- as.list(1:no.var)
names(var.list) <- var.names
# loop over gene names
for(i in var.names){
#2nd: adjusting just per group
tmp <- BodyWeight[, c(i, "Diet","Time", "Rat", "weight")]
fml <- as.formula( paste( i, "~", paste(c("Time:Diet", "weight"), collapse="+")))
var.list[[i]] <- lme(fml, random= ~ Time|Rat , method="REML", data=tmp, control = lmeControl(opt = "optim"), na.action = na.omit)
}
# List of data.frames/tibble
result <- lapply(var.list, function(x) tidy(x))
# Solution
result <- dplyr::bind_rows(result, .id = "var")
result2 <- result %>% dplyr::select( var, term, p.value)
As you can see I get 3 terms:
- Time:Diet1
- Time:Diet2
- Time:Diet3
My original database is made up by more variables, but time is numeric, as well as grouping variables are factor (grup_int and Diet respectively). The original database contains NAs I was wondering if it is possible to relevel to get the interaction remaining
UPDATE I don't know how to put in words, what is happening when I get this:
What does it mean? time:grup_int1 is significant but time:grup_int2 not? Is time:grup_int1 the meaning for statistical differences between group 1 ant the others groups (2 and 3) across time ? It doesn't make any sense to be significant for one but not for the others. How can I get an overall p-value comparison between groups?
The df used for the loop
structure(list(weight = c(240, 250, 255, 260, 262, 258, 266,
266, 265, 272, 278, 225, 230, 230, 232, 240, 240, 243, 244, 238,
247, 245, 245, 250, 250, 255, 262, 265, 267, 267, 264, 268, 269,
260, 255, 255, 265, 265, 268, 270, 272, 274, 273, 275, 255, 260,
255, 270, 270, 273, 274, 273, 276, 278, 280, 260, 265, 270, 275,
275, 277, 278, 278, 284, 279, 281, 275, 275, 260, 270, 273, 274,
276, 271, 282, 281, 284, 245, 255, 260, 268, 270, 265, 265, 267,
273, 274, 278, 410, 415, 425, 428, 438, 443, 442, 446, 456, 468,
478, 405, 420, 430, 440, 448, 460, 458, 464, 475, 484, 496, 445,
445, 450, 452, 455, 455, 451, 450, 462, 466, 472, 555, 560, 565,
580, 590, 597, 595, 595, 612, 618, 628, 470, 465, 475, 485, 487,
493, 493, 504, 507, 518, 525, 535, 525, 530, 533, 535, 540, 525,
530, 543, 544, 559, 520, 525, 530, 540, 543, 546, 538, 544, 553,
555, 548, 510, 510, 520, 515, 530, 538, 535, 542, 550, 553, 569
), Time = c(1, 8, 15, 22, 29, 36, 43, 44, 50, 57, 64, 1, 8, 15,
22, 29, 36, 43, 44, 50, 57, 64, 1, 8, 15, 22, 29, 36, 43, 44,
50, 57, 64, 1, 8, 15, 22, 29, 36, 43, 44, 50, 57, 64, 1, 8, 15,
22, 29, 36, 43, 44, 50, 57, 64, 1, 8, 15, 22, 29, 36, 43, 44,
50, 57, 64, 1, 8, 15, 22, 29, 36, 43, 44, 50, 57, 64, 1, 8, 15,
22, 29, 36, 43, 44, 50, 57, 64, 1, 8, 15, 22, 29, 36, 43, 44,
50, 57, 64, 1, 8, 15, 22, 29, 36, 43, 44, 50, 57, 64, 1, 8, 15,
22, 29, 36, 43, 44, 50, 57, 64, 1, 8, 15, 22, 29, 36, 43, 44,
50, 57, 64, 1, 8, 15, 22, 29, 36, 43, 44, 50, 57, 64, 1, 8, 15,
22, 29, 36, 43, 44, 50, 57, 64, 1, 8, 15, 22, 29, 36, 43, 44,
50, 57, 64, 1, 8, 15, 22, 29, 36, 43, 44, 50, 57, 64), Rat = structure(c(4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L,
11L, 11L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L,
13L, 13L, 13L, 13L, 13L, 13L, 13L, 15L, 15L, 15L, 15L, 15L, 15L,
15L, 15L, 15L, 15L, 15L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L,
14L, 14L, 14L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L,
16L), levels = c("2", "3", "4", "1", "8", "5", "6", "7", "11",
"9", "10", "12", "13", "15", "14", "16"), class = c("ordered",
"factor")), Diet = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), levels = c("1", "2", "3"), class = "factor"),
var1 = c(46.804, 45.792, 43.572, 50.587, 52.268, 45.402,
45.996, 48.566, 47.933, 48.336, 43.964, 52.441, 45.077, 45.699,
51.348, 47.888, 50.245, 48.111, 42.373, 44.231, 46.461, 45.639,
44.993, 50.391, 47.997, 48.368, 44.939, 47.805, 49.195, 51.419,
50.242, 45.856, 48.115, 48.053, 44.318, 50.493, 40.48, 47.418,
49.441, 49.414, 49.776, 49.043, 41.456, 48.52, 46.079, 43.058,
43.649, 43.967, 44.799, 48.081, 50.728, 51.207, 50.01, 47.346,
43.816, 52.575, 47.962, 49.061, 45.438, 43.822, 43.022, 48.864,
41.669, 45.58, 42.206, 42.788, 43.988, 46.636, 45.651, 42.502,
42.93, 49.196, 48.314, 46.412, 49.436, 47.193, 48.806, 45.142,
48.4, 48.419, 47.937, 42.605, 46.027, 45.325, 43.811, 45.922,
47.262, 46.566, 47.016, 46.111, 50.573, 47.436, 47.859, 47.839,
48.239, 46.762, 46.581, 55.9, 42.816, 49.636, 47.514, 39.489,
46.513, 49.323, 51.504, 45.558, 42.402, 45.153, 48.632, 47.908,
46.033, 44.963, 49.312, 45.733, 51.149, 45.755, 47.247, 46.674,
49.121, 43.411, 43.735, 44.777, 44.942, 48.207, 50.962, 47.1,
44.167, 50.007, 48.531, 46.749, 47.69, 47.396, 45.074, 49.203,
42.869, 47.713, 46.186, 52.191, 46.809, 47.578, 45.365, 53.228,
45.896, 45.428, 42.702, 47.328, 46.108, 49.857, 45.128, 43.772,
46.479, 51.573, 51.255, 45.33, 45.651, 48.521, 49.006, 51.906,
48.924, 53.693, 48.555, 48.367, 49.659, 47.839, 47.829, 48.212,
45.867, 51.965, 48.275, 49.091, 44.54, 44.978, 50.431, 48.915,
51.569, 46.008), var2 = c(53.273, 35.7, 44.406, 62.644, 42.971,
49.393, 52.585, 30.572, 52.156, 58.566, 44.631, 50.967, 48.583,
60.878, 50.712, 48.741, 47.584, 54.774, 62.288, 49.952, 33.225,
39.251, 53.158, 46.74, 46.796, 40.343, 58.491, 48.12, 41.854,
55.315, 44.553, 55.674, 47.404, 56.085, 53.443, 56.382, 74.641,
44.822, 73.191, 46.409, 66.453, 60.929, 51.394, 56.812, 53.244,
48.993, 54.763, 35.624, 63.346, 38.389, 48.288, 61.422, 45.707,
50.213, 42.629, 64.381, 39.218, 58.923, 61.569, 39.279, 54.748,
63.149, 54.639, 56.494, 66.907, 36.303, 52.297, 46.288, 62.852,
53.063, 43.096, 70.983, 49.175, 43.663, 44.456, 55.17, 51.383,
56.945, 51.471, 49.053, 63.062, 46.386, 48.094, 43.282, 35.509,
48.079, 48.223, 34.605, 42.669, 40.146, 53.089, 40.911, 39.396,
39.507, 41.481, 41.729, 28.367, 62.781, 59.154, 50.788, 84.282,
58.214, 53.404, 46.312, 36.983, 61.712, 53.109, 37.285, 50.834,
32.306, 47.969, 36.232, 46.744, 37.876, 50.369, 59.684, 40.427,
44.2, 46.955, 56.388, 46.667, 46.373, 49.932, 60.3, 27.732,
56.772, 54.569, 45.164, 37.747, 57.827, 49.998, 23.109, 53.339,
51.024, 58.875, 55.315, 40.637, 37.781, 47.209, 57.51, 46.956,
66.715, 51.034, 52.54, 60.342, 34.439, 58.08, 64.298, 56.374,
64.829, 49.392, 48.939, 38.746, 40.048, 46.856, 61.697, 51.253,
54.024, 61.166, 45.834, 43.578, 39.212, 67.632, 35.208, 23.676,
50.577, 46.832, 42.866, 57.808, 53.292, 42.474, 62.779, 34.747,
55.367, 27.928, 57.481), var3 = c(67.483, 67.977, 70.749,
76.851, 77.351, 59.821, 67.203, 72.635, 67.767, 67.627, 64.729,
70.388, 71.889, 61.95, 54.182, 71.049, 56.33, 68.982, 70.615,
68.416, 60.661, 73.205, 68.086, 72.611, 62.082, 74.026, 69.073,
67.488, 79.656, 79.612, 71.016, 61.284, 76.297, 62.713, 79.335,
68.155, 67.184, 73.889, 61.59, 68.273, 57.056, 66.95, 60.504,
77.121, 73.056, 60.739, 75.891, 74.363, 82.829, 86.902, 59.38,
77.63, 67.613, 72.666, 68.903, 71.631, 68.86, 57.653, 65.119,
60.825, 63.301, 64.945, 65.89, 80.092, 79.951, 65.475, 66.491,
72.355, 71.909, 70.339, 83.508, 59.686, 62.102, 67.291, 52.908,
70.508, 63.2, 66.62, 75.138, 72.311, 66.394, 64.812, 74.063,
51.087, 55.095, 56.906, 55.88, 71.916, 81.805, 67.339, 66.459,
70.005, 82.891, 65.456, 86.611, 65.401, 72.647, 59.198, 73.558,
66.536, 69.83, 77.101, 65.494, 62.129, 71.86, 85.406, 70.619,
79.225, 71.806, 62.86, 79.615, 65.049, 63.472, 63.46, 80.215,
74.033, 58.949, 69.519, 66.39, 68.751, 69.356, 75.201, 54.815,
50.592, 76.019, 62.177, 66.677, 69.495, 82.428, 70.743, 59.458,
75.326, 69.615, 67.157, 69.793, 66.205, 75.328, 76.311, 64.36,
62.207, 74.16, 51.799, 66.378, 67.964, 66.937, 59.779, 59.735,
64.992, 66.817, 60.595, 65.388, 50.404, 61.742, 65.789, 77.047,
56.86, 69.756, 66.338, 57.149, 63.407, 73.601, 66.042, 67.421,
51.208, 47.265, 69.465, 63.471, 75.845, 74.938, 69.009, 81.807,
81.195, 64.288, 77.085, 64.698, 65.747), var4 = c(143.048,
139.381, 152.849, 163.726, 155.712, 146.957, 140.231, 148.321,
141.917, 156.873, 165.845, 158.824, 157.305, 144.521, 145.253,
134.892, 153.597, 124.482, 153.096, 134.954, 171.238, 147.43,
150.983, 129.945, 150.652, 142.259, 156.631, 152.215, 152.196,
153.54, 139.712, 168.495, 152.043, 148.538, 164.067, 172.445,
144.175, 144.015, 149.093, 135.921, 156.787, 148.292, 144.958,
137.332, 148.185, 151.678, 160.774, 168.479, 175.839, 148.211,
163.747, 145.416, 131.502, 154.997, 152.35, 125.955, 167.116,
141.734, 135.958, 162.752, 139.482, 149.315, 144.425, 144.654,
160.05, 153.767, 148.532, 143.102, 135.553, 143.305, 153.021,
127.982, 165.047, 144.142, 136.568, 145.793, 143.913, 158.002,
154.043, 167.63, 165.472, 143.681, 166.417, 152.012, 145.843,
144.448, 138.782, 141.475, 150.717, 140.765, 164.657, 133.655,
151.908, 160.084, 138.25, 155.93, 142.756, 170.935, 128.346,
149.986, 154.099, 141.688, 150.616, 149.999, 139.195, 160.921,
148.261, 144.611, 154.183, 157.2, 130.512, 155.004, 125.064,
124.813, 150.074, 145.615, 152.238, 144.241, 149.079, 159.354,
153.398, 146.609, 155.032, 145.522, 143.296, 139.502, 148.944,
143.446, 155.094, 148.208, 154.417, 153.249, 154.544, 140.746,
156.798, 146.069, 134.905, 145.318, 137.198, 142.087, 137.117,
151.782, 147.245, 148.583, 152.535, 151.265, 152.272, 162.014,
150.055, 154.343, 139.717, 148.752, 140.933, 146.761, 145.311,
146.361, 155.264, 154.984, 168.623, 129.31, 147.492, 139.828,
161.22, 145.89, 139.146, 138.104, 148.543, 155.836, 143.594,
159.32, 163.648, 152.058, 160.412, 133, 155.305, 144.094),
var5 = c(135.794, 136.92, 131.548, 132.091, 143.587, 130.393,
133.678, 143.005, 136.036, 128.161, 134.036, 144.237, 125.767,
133.32, 137.047, 140.3, 137.92, 144.883, 140.901, 140.632,
136.315, 138.138, 140.855, 124.721, 143.153, 144.514, 142.092,
144.765, 139.205, 134.905, 141.815, 122.992, 138.249, 125.918,
125.236, 136.503, 137.462, 130.757, 135.141, 135.413, 134.991,
137.143, 153.329, 148.779, 138.249, 133.978, 132.454, 139.955,
122.84, 135.984, 140.941, 134.203, 146.091, 146.766, 136.259,
151.474, 139.683, 132.913, 142.39, 134.625, 124.717, 145.74,
135.647, 135.119, 144.182, 133.961, 130.387, 141.314, 129.26,
127.154, 154.74, 142.327, 142.23, 139.433, 148.245, 132.123,
127.374, 134.518, 133.422, 143.66, 129.306, 138.31, 147.787,
137.415, 138.686, 143.896, 140.973, 128.592, 139.53, 146.29,
145, 129.926, 132.338, 155.193, 132.286, 140.86, 135.343,
137.037, 140.147, 131.872, 137.858, 142.137, 131.587, 121.455,
142.264, 138.533, 138.769, 136.296, 141.055, 133.475, 148.136,
136.175, 143.947, 152.026, 137.649, 142.466, 144.819, 153.745,
138.95, 138.146, 133.404, 147.928, 145.928, 146.421, 145.156,
134.289, 127.955, 148.223, 129.987, 139.295, 141.035, 150.867,
136.412, 139.705, 142.013, 129.377, 143.521, 145.595, 142.936,
137.703, 133.344, 144.46, 143.001, 137.017, 149.985, 144.242,
158.029, 151.835, 133.262, 123.893, 136.338, 143.951, 136.71,
136.007, 147.102, 134.472, 128.058, 145.521, 142.379, 135.869,
144.135, 141.297, 149.133, 134.008, 138.566, 117.773, 151.701,
140.288, 144.17, 140.886, 151.585, 136.438, 150.399, 140.535,
142.601, 136.449)), row.names = c(NA, 176L), class = c("nfnGroupedData",
"nfGroupedData", "groupedData", "data.frame"), outer = ~Diet, formula = weight ~
Time | Rat, labels = list(x = "Time", y = "Body weight"), units = list(
x = "(days)", y = "(g)"), FUN = structure(function (x)
max(x, na.rm = TRUE), source = "function(x) max(x, na.rm = TRUE)"), order.groups = TRUE)