I have a dataset with some growth measurements pre and post vaccination and 3 different groups (3 different types of vaccine) so there are both a within (pre/post) and between factor (group).
- Problem 1: The number of observations varies between people (some people can have 1 observation and other up to 5)
- Problem 2: Some participants have pre intervention but no post intervention data.
- Problem 3: Time between vaccine and growth measurement varies (ie. some people will have data collected 10 days after vaccination, some other 15 days after etc..)
With all those problems I am a bit confused about how should I answer my questions: Is there a change in BMI z-score between groups after vaccination?
- zBMI: BMI zscore
- tipo_vacuna: vaccine 1, vaccine 2, vaccine 3.
- tiempo_peso_vacuna1: time between first vaccine and zBMI measurement (negative if before -vaccination and positive if after vaccination)
- Index1: number of observation per participants
- study_id: participants
The idea was to apply a linear mixed model such as:
df1_lmm =lmer(zBMI ~ tiempo_peso_vacuna1 * tipo_vacuna + (1|study_id), data=df1)
summary(df1_lmm)
anova (df1_lmm)
However I have a lot of doubts:
- Is the model appropriate to answer the research question?
- Does it control for zBMI pre intervention?
- Does the random intercept accounts for the number of observations and their closeness in time to solve the autocorrelation problem?
Any insight would be much appreciated, thank you.
Reproducible sample:
df1<-structure(list(study_id = structure(c(3, 3, 3, 3, 3, 3, 4, 4,
4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7,
7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10,
11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 13, 13, 13, 13,
13, 13, 14, 14, 14, 14, 14, 14, 16, 16, 16, 16, 16, 16, 18, 18,
18, 18, 18, 18, 19, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 20,
21, 21, 21, 21), format.spss = "F4.0", display_width = 11L),
ageyear_vacuna1 = structure(c(0.780287474332649, 0.780287474332649,
0.780287474332649, 0.780287474332649, 0.780287474332649,
0.780287474332649, 0.815879534565366, 0.815879534565366,
0.815879534565366, 0.815879534565366, 0.815879534565366,
0.815879534565366, 0.432580424366872, 0.432580424366872,
0.432580424366872, 0.432580424366872, 0.432580424366872,
0.432580424366872, 1.54414784394251, 1.54414784394251, 1.54414784394251,
1.54414784394251, 1.54414784394251, 1.54414784394251, 0.240930869267625,
0.240930869267625, 0.240930869267625, 0.240930869267625,
0.240930869267625, 0.240930869267625, 0.684462696783025,
0.684462696783025, 0.684462696783025, 0.684462696783025,
0.684462696783025, 0.684462696783025, 0.0711841204654346,
0.0711841204654346, 0.0711841204654346, 0.0711841204654346,
0.0711841204654346, 0.0711841204654346, 0.355920602327173,
0.355920602327173, 0.355920602327173, 0.355920602327173,
0.355920602327173, 0.355920602327173, 0.407939767282683,
0.407939767282683, 0.407939767282683, 0.407939767282683,
0.407939767282683, 0.407939767282683, 1.79603011635866, 1.79603011635866,
1.79603011635866, 1.79603011635866, 1.79603011635866, 1.79603011635866,
0.213552361396304, 0.213552361396304, 0.213552361396304,
0.213552361396304, 0.213552361396304, 0.213552361396304,
0.550308008213552, 0.550308008213552, 0.550308008213552,
0.550308008213552, 0.550308008213552, 0.550308008213552,
0.292950034223135, 0.292950034223135, 0.292950034223135,
0.292950034223135, 0.292950034223135, 0.292950034223135,
0.785763175906913, 0.785763175906913, 0.785763175906913,
0.785763175906913, 0.785763175906913, 0.785763175906913,
0.380561259411362, 0.380561259411362, 0.380561259411362,
0.380561259411362, 0.380561259411362, 0.380561259411362,
NA, NA, NA, NA, NA, NA, 0.303901437371663, 0.303901437371663,
0.303901437371663, 0.303901437371663), format.spss = "F8.2", display_width = 17L),
height_m = structure(c(157, 157, 157, 157, 157, 157, 165,
165, 165, 165, 165, 165, 150, 150, 150, 150, 150, 150, 175,
175, 175, 175, 175, 175, 171, 171, 171, 171, 171, 171, 148,
148, 148, 148, 148, 148, 165, 165, 165, 165, 165, 165, 164,
164, 164, 164, 164, 164, 169, 169, 169, 169, 169, 169, 158,
158, 158, 158, 158, 158, 164, 164, 164, 164, 164, 164, 171,
171, 171, 171, 171, 171, 160, 160, 160, 160, 160, 160, 166,
166, 166, 166, 166, 166, 165, 165, 165, 165, 165, 165, 170,
170, 170, 170, 170, 170, 165, 165, 165, 165), label = "height mothers", format.spss = "F6.2", display_width = 14L),
weight_m = structure(c(47.3, 47.3, 47.3, 47.3, 47.3, 47.3,
96, 96, 96, 96, 96, 96, 50, 50, 50, 50, 50, 50, 66, 66, 66,
66, 66, 66, 60, 60, 60, 60, 60, 60, 48, 48, 48, 48, 48, 48,
59, 59, 59, 59, 59, 59, 76.6, 76.6, 76.6, 76.6, 76.6, 76.6,
79, 79, 79, 79, 79, 79, 61, 61, 61, 61, 61, 61, 62, 62, 62,
62, 62, 62, 64, 64, 64, 64, 64, 64, 69, 69, 69, 69, 69, 69,
56, 56, 56, 56, 56, 56, 55, 55, 55, 55, 55, 55, 60, 60, 60,
60, 60, 60, 62, 62, 62, 62), label = "weight mothers", format.spss = "F6.2", display_width = 16L),
weight_m_2 = structure(c(48, 48, 48, 48, 48, 48, 96, 96,
96, 96, 96, 96, 50, 50, 50, 50, 50, 50, 64, 64, 64, 64, 64,
64, 65, 65, 65, 65, 65, 65, 46, 46, 46, 46, 46, 46, 75, 75,
75, 75, 75, 75, 69, 69, 69, 69, 69, 69, 84, 84, 84, 84, 84,
84, 64, 64, 64, 64, 64, 64, 64.5, 64.5, 64.5, 64.5, 64.5,
64.5, 64, 64, 64, 64, 64, 64, 63, 63, 63, 63, 63, 63, 52,
52, 52, 52, 52, 52, 60, 60, 60, 60, 60, 60, 62, 62, 62, 62,
62, 62, 61, 61, 61, 61), format.spss = "F6.2", display_width = 13L),
gwg_m = structure(c(12, 12, 12, 12, 12, 12, 4, 4, 4, 4, 4,
4, 6, 6, 6, 6, 6, 6, 10, 10, 10, 10, 10, 10, 12, 12, 12,
12, 12, 12, 8, 8, 8, 8, 8, 8, 20, 20, 20, 20, 20, 20, 8,
8, 8, 8, 8, 8, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
15, 12, 12, 12, 12, 12, 12, 18, 18, 18, 18, 18, 18, 7, 7,
7, 7, 7, 7, 17, 17, 17, 17, 17, 17, 20, 20, 20, 20, 20, 20,
10, 10, 10, 10, 10, 10, 9, 9, 9, 9), format.spss = "F4.2", display_width = 15L),
birth_date = structure(c("2020/04/10", "2020/04/10", "2020/04/10",
"2020/04/10", "2020/04/10", "2020/04/10", "2020/04/03", "2020/04/03",
"2020/04/03", "2020/04/03", "2020/04/03", "2020/04/03", "2020/10/19",
"2020/10/19", "2020/10/19", "2020/10/19", "2020/10/19", "2020/10/19",
"2019/06/30", "2019/06/30", "2019/06/30", "2019/06/30", "2019/06/30",
"2019/06/30", "2020/12/06", "2020/12/06", "2020/12/06", "2020/12/06",
"2020/12/06", "2020/12/06", "2020/05/26", "2020/05/26", "2020/05/26",
"2020/05/26", "2020/05/26", "2020/05/26", "2021/01/30", "2021/01/30",
"2021/01/30", "2021/01/30", "2021/01/30", "2021/01/30", "2020/10/20",
"2020/10/20", "2020/10/20", "2020/10/20", "2020/10/20", "2020/10/20",
"2020/11/02", "2020/11/02", "2020/11/02", "2020/11/02", "2020/11/02",
"2020/11/02", "2019/05/11", "2019/05/11", "2019/05/11", "2019/05/11",
"2019/05/11", "2019/05/11", "2020/12/21", "2020/12/21", "2020/12/21",
"2020/12/21", "2020/12/21", "2020/12/21", "2020/09/06", "2020/09/06",
"2020/09/06", "2020/09/06", "2020/09/06", "2020/09/06", "2020/11/21",
"2020/11/21", "2020/11/21", "2020/11/21", "2020/11/21", "2020/11/21",
"2020/05/14", "2020/05/14", "2020/05/14", "2020/05/14", "2020/05/14",
"2020/05/14", "2020/08/29", "2020/08/29", "2020/08/29", "2020/08/29",
"2020/08/29", "2020/08/29", "2020/11/16", "2020/11/16", "2020/11/16",
"2020/11/16", "2020/11/16", "2020/11/16", "2020/11/06", "2020/11/06",
"2020/11/06", "2020/11/06"), format.spss = "A10", display_width = 10L),
sex_infant = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 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, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L), .Label = c("0", "1"), class = "factor"), delivery_method = structure(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, 1, 1, 1, 1, 1, 1, 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, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2,
0, 0, 0, 0), format.spss = "F1.0", display_width = 11L),
gestational_age = structure(c(41.4, 41.4, 41.4, 41.4, 41.4,
41.4, 41, 41, 41, 41, 41, 41, 40.86, 40.86, 40.86, 40.86,
40.86, 40.86, 38, 38, 38, 38, 38, 38, 41, 41, 41, 41, 41,
41, 41.29, 41.29, 41.29, 41.29, 41.29, 41.29, 41.14, 41.14,
41.14, 41.14, 41.14, 41.14, 41, 41, 41, 41, 41, 41, 40, 40,
40, 40, 40, 40, 40.14, 40.14, 40.14, 40.14, 40.14, 40.14,
40.71, 40.71, 40.71, 40.71, 40.71, 40.71, 40.71, 40.71, 40.71,
40.71, 40.71, 40.71, 36, 36, 36, 36, 36, 36, 40, 40, 40,
40, 40, 40, 40.71, 40.71, 40.71, 40.71, 40.71, 40.71, 40,
40, 40, 40, 40, 40, 37, 37, 37, 37), format.spss = "F5.2", display_width = 16L),
tipo_vacuna = structure(c(3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L,
2L, 4L, 4L, 4L, 4L, 4L, 4L, 2L, 2L, 2L, 2L, 2L, 2L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 3L, 3L,
3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L), .Label = c("NA", "1", "2", "3"), class = "factor"),
fecha_vacuna_1 = structure(c("2021/01/20", "2021/01/20",
"2021/01/20", "2021/01/20", "2021/01/20", "2021/01/20", "2021/01/26",
"2021/01/26", "2021/01/26", "2021/01/26", "2021/01/26", "2021/01/26",
"2021/03/26", "2021/03/26", "2021/03/26", "2021/03/26", "2021/03/26",
"2021/03/26", "2021/01/14", "2021/01/14", "2021/01/14", "2021/01/14",
"2021/01/14", "2021/01/14", "2021/03/04", "2021/03/04", "2021/03/04",
"2021/03/04", "2021/03/04", "2021/03/04", "2021/01/31", "2021/01/31",
"2021/01/31", "2021/01/31", "2021/01/31", "2021/01/31", "2021/02/25",
"2021/02/25", "2021/02/25", "2021/02/25", "2021/02/25", "2021/02/25",
"2021/02/27", "2021/02/27", "2021/02/27", "2021/02/27", "2021/02/27",
"2021/02/27", "2021/03/31", "2021/03/31", "2021/03/31", "2021/03/31",
"2021/03/31", "2021/03/31", "2021/02/25", "2021/02/25", "2021/02/25",
"2021/02/25", "2021/02/25", "2021/02/25", "2021/03/09", "2021/03/09",
"2021/03/09", "2021/03/09", "2021/03/09", "2021/03/09", "2021/03/26",
"2021/03/26", "2021/03/26", "2021/03/26", "2021/03/26", "2021/03/26",
"2021/03/08", "2021/03/08", "2021/03/08", "2021/03/08", "2021/03/08",
"2021/03/08", "2021/02/25", "2021/02/25", "2021/02/25", "2021/02/25",
"2021/02/25", "2021/02/25", "2021/01/15", "2021/01/15", "2021/01/15",
"2021/01/15", "2021/01/15", "2021/01/15", "", "", "", "",
"", "", "2021/02/25", "2021/02/25", "2021/02/25", "2021/02/25"
), format.spss = "A10", display_width = 10L), fecha_vacuna_2 = structure(c("2021/02/22",
"2021/02/22", "2021/02/22", "2021/02/22", "2021/02/22", "2021/02/22",
"2021/02/18", "2021/02/18", "2021/02/18", "2021/02/18", "2021/02/18",
"2021/02/18", "2021/04/26", "2021/04/26", "2021/04/26", "2021/04/26",
"2021/04/26", "2021/04/26", "2021/02/04", "2021/02/04", "2021/02/04",
"2021/02/04", "2021/02/04", "2021/02/04", "2021/06/02", "2021/06/02",
"2021/06/02", "2021/06/02", "2021/06/02", "2021/06/02", "2021/02/22",
"2021/02/22", "2021/02/22", "2021/02/22", "2021/02/22", "2021/02/22",
"2021/05/29", "2021/05/29", "2021/05/29", "2021/05/29", "2021/05/29",
"2021/05/29", "2021/06/02", "2021/06/02", "2021/06/02", "2021/06/02",
"2021/06/02", "2021/06/02", "2021/07/06", "2021/07/06", "2021/07/06",
"2021/07/06", "2021/07/06", "2021/07/06", "2021/05/30", "2021/05/30",
"2021/05/30", "2021/05/30", "2021/05/30", "2021/05/30", "2021/06/05",
"2021/06/05", "2021/06/05", "2021/06/05", "2021/06/05", "2021/06/05",
"2021/04/23", "2021/04/23", "2021/04/23", "2021/04/23", "2021/04/23",
"2021/04/23", "2021/05/31", "2021/05/31", "2021/05/31", "2021/05/31",
"2021/05/31", "2021/05/31", "2021/05/31", "2021/05/31", "2021/05/31",
"2021/05/31", "2021/05/31", "2021/05/31", "2021/02/09", "2021/02/09",
"2021/02/09", "2021/02/09", "2021/02/09", "2021/02/09", "",
"", "", "", "", "", "2021/03/18", "2021/03/18", "2021/03/18",
"2021/03/18"), format.spss = "A10", display_width = 18L),
Index1 = structure(c(1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6,
1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1,
2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2,
3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3,
4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4,
5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4), format.spss = "F4.0"),
date = structure(c("", "", "", "", "", "", "2020/04/03",
"2020/11/11", "2021/03/03", "2021/04/07", "", "", "2020/10/25",
"2021/02/19", "", "", "", "", "2019/07/05", "2020/07/01",
"2021/01/15", "", "", "", "2020/12/06", "2021/02/11", "2021/03/11",
"", "", "", "2020/06/05", "2021/01/26", "", "", "", "", "2021/01/30",
"2021/02/17", "2021/03/01", "2021/04/07", "", "", "2020/10/20",
"2021/02/22", "2021/03/25", "", "", "", "2020/11/16", "2021/03/09",
"2021/04/08", "", "", "", "2019/05/20", "", "", "", "", "",
"2020/12/21", "2021/02/22", "", "", "", "", "2020/09/06",
"2021/03/04", "", "", "", "", "2020/11/27", "2021/01/25",
"2021/02/23", "2021/03/22", "", "", "2020/05/14", "2020/11/16",
"", "", "", "", "2020/09/08", "2020/12/30", "2021/02/04",
"2021/03/03", "", "", "", "", "", "", "", "", "2020/11/19",
"2021/01/11", "2021/03/08", ""), format.spss = "A10", display_width = 10L),
weight = structure(c(NA, NA, NA, NA, NA, NA, 3.6, 9.7, 10,
11.6, NA, NA, 4.4, 8.2, NA, NA, NA, NA, 2.84, 10.87, 12.5,
NA, NA, NA, 3.56, 6.3, 7.11, NA, NA, NA, 2.84, 8.02, NA,
NA, NA, NA, 3.95, 3.9, 4.97, 6.18, NA, NA, 3.4, 7.7, 8.4,
NA, NA, NA, 4.13, 8.88, 9.2, NA, NA, NA, 3.5, NA, NA, NA,
NA, NA, 3.22, 4.79, NA, NA, NA, NA, 4.6, 8.36, NA, NA, NA,
NA, 1.98, 3.38, 4.32, 5.18, NA, NA, 3.42, 6.12, NA, NA, NA,
NA, 3.28, 6.74, 7.36, 7.68, NA, NA, NA, NA, NA, NA, NA, NA,
3.52, 5.12, 7.2, NA), format.spss = "F5.2", display_width = 12L),
length = structure(c(NA, NA, NA, NA, NA, NA, 49, 68, 72,
73, NA, NA, 52.5, 64, NA, NA, NA, NA, 48.5, 77, 87.5, NA,
NA, NA, 53, 61.5, 65, NA, NA, NA, 48, 68.5, NA, NA, NA, NA,
54, 57, 59.5, 63, NA, NA, 52, 64, 68, NA, NA, NA, 54, 67,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 50, 56.5, NA, NA,
NA, NA, 53, 71.5, NA, NA, NA, NA, NA, 51, NA, 56, NA, NA,
49, 65.5, NA, NA, NA, NA, 51, 63, 64.5, 68, NA, NA, NA, NA,
NA, NA, NA, NA, 51.5, 57, 62, NA), format.spss = "F5.2", display_width = 15L),
headcirc = structure(c(NA, NA, NA, NA, NA, NA, 34.5, 45,
49, 49, NA, NA, 36.5, 42, NA, NA, NA, NA, 34.5, 47.5, 49,
NA, NA, NA, 35, 41, 42, NA, NA, NA, 34, 43, NA, NA, NA, NA,
36, 37.5, 38.5, 41, NA, NA, 35, 42.5, 44, NA, NA, NA, 37,
43.5, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 33, 37.7, NA,
NA, NA, NA, 37.5, 44, NA, NA, NA, NA, NA, 37.5, NA, 40.5,
NA, NA, 33, 41, NA, NA, NA, NA, 35.5, 41.5, NA, 42, NA, NA,
NA, NA, NA, NA, NA, NA, 35, 38, 41, NA), format.spss = "F5.2", display_width = 15L),
zhaz = structure(c(NA, NA, NA, NA, NA, NA, -0.47, -0.73,
-1.08, -1.22, NA, NA, 1.24, 0.84, NA, NA, NA, NA, -1.19,
0.5, 1.72, NA, NA, NA, 1.65, 1.21, 1.58, NA, NA, NA, -1.91,
-0.98, NA, NA, NA, NA, 2.17, 2.11, 2.49, 1.95, NA, NA, 1.12,
-0.06, 0.88, NA, NA, NA, 0.87, 1.31, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, 0.46, -0.38, NA, NA, NA, NA, 1.65, 1.9,
NA, NA, NA, NA, NA, -3.16, NA, -2.79, NA, NA, -0.08, -0.18,
NA, NA, NA, NA, 0.09, 0.38, 0.03, 0.92, NA, NA, NA, NA, NA,
NA, NA, NA, 0.08, -0.28, -0.05, NA), label = "Length/height-for-age z-score", format.spss = "F8.2", display_width = 10L),
BMI = structure(c(NA, NA, NA, NA, NA, NA, 14.993752603082,
20.977508650519, 19.2901234567901, 21.7676862450741, NA,
NA, 15.9637188208617, 20.01953125, NA, NA, NA, NA, 12.073546604315,
18.3336144375105, 16.3265306122449, NA, NA, NA, 12.6735493058028,
16.656751933373, 16.8284023668639, NA, NA, NA, 12.3263888888889,
17.0920134263946, NA, NA, NA, NA, 13.5459533607682, 12.0036934441367,
14.0385565991102, 15.5706727135299, NA, NA, 12.5739644970414,
18.798828125, 18.1660899653979, NA, NA, NA, 14.1632373113855,
19.7816885720651, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
12.88, 15.0050904534419, NA, NA, NA, NA, 16.375934496262,
16.3528778913394, NA, NA, NA, NA, NA, 12.9950019223376, NA,
16.5178571428571, NA, NA, 14.2440649729279, 14.264902977682,
NA, NA, NA, NA, 12.6105344098424, 16.9816074577979, 17.6912445165555,
16.6089965397924, NA, NA, NA, NA, NA, NA, NA, NA, 13.2717504006033,
15.7586949830717, 18.7304890738814, NA), format.spss = "F8.2", display_width = 10L),
zBMI = structure(c(NA, NA, NA, NA, NA, NA, 1.15, 2.29, 1.59,
3.09, NA, NA, 2.02, 1.99, NA, NA, NA, NA, -1.02, 1.08, 0.19,
NA, NA, NA, -0.59, 0.12, -0.08, NA, NA, NA, -0.9, -0.12,
NA, NA, NA, NA, 0.11, -1.58, -0.66, -0.66, NA, NA, -0.68,
1.07, 0.59, NA, NA, NA, 0.4, 1.67, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, -0.38, -0.56, NA, NA, NA, NA, 2.03, -0.71,
NA, NA, NA, NA, NA, -2.11, NA, -0.1, NA, NA, 0.71, -1.91,
NA, NA, NA, NA, -0.54, 0.2, 0.53, -0.2, NA, NA, NA, NA, NA,
NA, NA, NA, -0.11, -0.09, 1.27, NA), label = "BMI-for-age z-score", format.spss = "F8.2", display_width = 10L),
zwaz = structure(c(NA, NA, NA, NA, NA, NA, 0.51, 1.34, 0.56,
1.67, NA, NA, 2.02, 1.94, NA, NA, NA, NA, -1.25, 1.08, 1.1,
NA, NA, NA, 0.43, 0.75, 0.83, NA, NA, NA, -1.59, -0.67, NA,
NA, NA, NA, 1.17, -0.05, 0.83, 0.59, NA, NA, 0.11, 0.77,
0.95, NA, NA, NA, 0.69, 2.02, 1.77, NA, NA, NA, -0.11, NA,
NA, NA, NA, NA, -0.03, -0.61, NA, NA, NA, NA, 2.3, 0.53,
NA, NA, NA, NA, -3.26, -3.29, -2.46, -1.72, NA, NA, 0.4,
-1.5, NA, NA, NA, NA, -0.31, 0.36, 0.42, 0.37, NA, NA, NA,
NA, NA, NA, NA, NA, -0.03, -0.21, 0.91, NA), label = "Weight-for-age z-score", format.spss = "F8.2", display_width = 10L),
zwhz = structure(c(NA, NA, NA, NA, NA, NA, 1.52, 2.31, 1.42,
2.83, NA, NA, 1.27, 1.9, NA, NA, NA, NA, -0.76, 1.11, 0.38,
NA, NA, NA, -1.38, -0.18, -0.27, NA, NA, NA, -0.42, -0.09,
NA, NA, NA, NA, -0.91, -3.36, -2.01, -1.14, NA, NA, -1.17,
1.1, 0.64, NA, NA, NA, -0.39, 1.64, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, -0.45, -0.36, NA, NA, NA, NA, 1.55, -0.57,
NA, NA, NA, NA, NA, -0.59, NA, 0.79, NA, NA, 0.87, -1.83,
NA, NA, NA, NA, -0.93, 0.21, 0.6, -0.09, NA, NA, NA, NA,
NA, NA, NA, NA, -0.48, 0.08, 1.3, NA), label = "Weight-for-length/height z-score", format.spss = "F8.2", display_width = 10L),
zhcz = structure(c(NA, NA, NA, NA, NA, NA, 0.03, 0.69, 2.55,
2.25, NA, NA, 1.75, 1.09, NA, NA, NA, NA, -0.37, 1.1, 1.15,
NA, NA, NA, 0.42, 1.33, 1.13, NA, NA, NA, -1.22, -1.25, NA,
NA, NA, NA, 1.21, 1.07, 1.08, 1.33, NA, NA, 0.42, 0.64, 1.11,
NA, NA, NA, 0.97, 1.42, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, -0.74, -0.54, NA, NA, NA, NA, 2.39, 0.62, NA, NA, NA,
NA, NA, -0.78, NA, -0.05, NA, NA, -0.74, -0.98, NA, NA, NA,
NA, 0.57, 0.69, NA, -0.21, NA, NA, NA, NA, NA, NA, NA, NA,
-0.11, -0.4, 0.32, NA), label = "hc/age SD", format.spss = "F8.2", display_width = 10L),
tiempo_peso_vacuna1 = c(NA, NA, NA, NA, NA, NA, -298, -76,
36, 71, NA, NA, -152, -35, NA, NA, NA, NA, -559, -197, 1,
NA, NA, NA, -88, -21, 7, NA, NA, NA, -240, -5, NA, NA, NA,
NA, -26, -8, 4, 41, NA, NA, -130, -5, 26, NA, NA, NA, -135,
-22, 8, NA, NA, NA, -647, NA, NA, NA, NA, NA, -78, -15, NA,
NA, NA, NA, -201, -22, NA, NA, NA, NA, -101, -42, -13, 14,
NA, NA, -287, -101, NA, NA, NA, NA, -129, -16, 20, 47, NA,
NA, NA, NA, NA, NA, NA, NA, -98, -45, 11, NA), tiempo_peso_vacuna2 = structure(c(NA,
NA, NA, NA, NA, NA, -321, -99, 13, 48, NA, NA, -183, -66,
NA, NA, NA, NA, -580, -218, -20, NA, NA, NA, -178, -111,
-83, NA, NA, NA, -262, -27, NA, NA, NA, NA, -119, -101, -89,
-52, NA, NA, -225, -100, -69, NA, NA, NA, -232, -119, -89,
NA, NA, NA, -741, NA, NA, NA, NA, NA, -166, -103, NA, NA,
NA, NA, -229, -50, NA, NA, NA, NA, -185, -126, -97, -70,
NA, NA, -382, -196, NA, NA, NA, NA, -154, -41, -5, 22, NA,
NA, NA, NA, NA, NA, NA, NA, -119, -66, -10, NA), format.spss = "F8.2", display_width = 21L),
agemos_vacuna1 = structure(c(9.36344969199179, 9.36344969199179,
9.36344969199179, 9.36344969199179, 9.36344969199179, 9.36344969199179,
9.79055441478439, 9.79055441478439, 9.79055441478439, 9.79055441478439,
9.79055441478439, 9.79055441478439, 5.19096509240246, 5.19096509240246,
5.19096509240246, 5.19096509240246, 5.19096509240246, 5.19096509240246,
18.5297741273101, 18.5297741273101, 18.5297741273101, 18.5297741273101,
18.5297741273101, 18.5297741273101, 2.8911704312115, 2.8911704312115,
2.8911704312115, 2.8911704312115, 2.8911704312115, 2.8911704312115,
8.2135523613963, 8.2135523613963, 8.2135523613963, 8.2135523613963,
8.2135523613963, 8.2135523613963, 0.854209445585216, 0.854209445585216,
0.854209445585216, 0.854209445585216, 0.854209445585216,
0.854209445585216, 4.27104722792608, 4.27104722792608, 4.27104722792608,
4.27104722792608, 4.27104722792608, 4.27104722792608, 4.8952772073922,
4.8952772073922, 4.8952772073922, 4.8952772073922, 4.8952772073922,
4.8952772073922, 21.5523613963039, 21.5523613963039, 21.5523613963039,
21.5523613963039, 21.5523613963039, 21.5523613963039, 2.56262833675565,
2.56262833675565, 2.56262833675565, 2.56262833675565, 2.56262833675565,
2.56262833675565, 6.60369609856263, 6.60369609856263, 6.60369609856263,
6.60369609856263, 6.60369609856263, 6.60369609856263, 3.51540041067762,
3.51540041067762, 3.51540041067762, 3.51540041067762, 3.51540041067762,
3.51540041067762, 9.42915811088296, 9.42915811088296, 9.42915811088296,
9.42915811088296, 9.42915811088296, 9.42915811088296, 4.56673511293635,
4.56673511293635, 4.56673511293635, 4.56673511293635, 4.56673511293635,
4.56673511293635, NA, NA, NA, NA, NA, NA, 3.64681724845996,
3.64681724845996, 3.64681724845996, 3.64681724845996), format.spss = "F8.2", display_width = 16L),
agemos_vacuna2 = structure(c(10.4476386036961, 10.4476386036961,
10.4476386036961, 10.4476386036961, 10.4476386036961, 10.4476386036961,
10.5462012320329, 10.5462012320329, 10.5462012320329, 10.5462012320329,
10.5462012320329, 10.5462012320329, 6.20944558521561, 6.20944558521561,
6.20944558521561, 6.20944558521561, 6.20944558521561, 6.20944558521561,
19.2197125256674, 19.2197125256674, 19.2197125256674, 19.2197125256674,
19.2197125256674, 19.2197125256674, 5.84804928131417, 5.84804928131417,
5.84804928131417, 5.84804928131417, 5.84804928131417, 5.84804928131417,
8.93634496919918, 8.93634496919918, 8.93634496919918, 8.93634496919918,
8.93634496919918, 8.93634496919918, 3.90965092402464, 3.90965092402464,
3.90965092402464, 3.90965092402464, 3.90965092402464, 3.90965092402464,
7.39219712525667, 7.39219712525667, 7.39219712525667, 7.39219712525667,
7.39219712525667, 7.39219712525667, 8.08213552361396, 8.08213552361396,
8.08213552361396, 8.08213552361396, 8.08213552361396, 8.08213552361396,
24.6406570841889, 24.6406570841889, 24.6406570841889, 24.6406570841889,
24.6406570841889, 24.6406570841889, 5.45379876796715, 5.45379876796715,
5.45379876796715, 5.45379876796715, 5.45379876796715, 5.45379876796715,
7.52361396303901, 7.52361396303901, 7.52361396303901, 7.52361396303901,
7.52361396303901, 7.52361396303901, 6.27515400410678, 6.27515400410678,
6.27515400410678, 6.27515400410678, 6.27515400410678, 6.27515400410678,
12.5503080082136, 12.5503080082136, 12.5503080082136, 12.5503080082136,
12.5503080082136, 12.5503080082136, 5.38809034907598, 5.38809034907598,
5.38809034907598, 5.38809034907598, 5.38809034907598, 5.38809034907598,
NA, NA, NA, NA, NA, NA, 4.33675564681725, 4.33675564681725,
4.33675564681725, 4.33675564681725), format.spss = "F8.2", display_width = 16L),
agew = structure(c(NA, NA, NA, NA, NA, NA, 0, 222, 334, 369,
NA, NA, 6, 123, NA, NA, NA, NA, 5, 367, 565, NA, NA, NA,
0, 67, 95, NA, NA, NA, 10, 245, NA, NA, NA, NA, 0, 18, 30,
67, NA, NA, 0, 125, 156, NA, NA, NA, 14, 127, 157, NA, NA,
NA, 9, NA, NA, NA, NA, NA, 0, 63, NA, NA, NA, NA, 0, 179,
NA, NA, NA, NA, 6, 65, 94, 121, NA, NA, 0, 186, NA, NA, NA,
NA, 10, 123, 159, 186, NA, NA, NA, NA, NA, NA, NA, NA, 13,
66, 122, NA), format.spss = "F8.2", display_width = 10L)), row.names = c(NA,
-100L), class = c("tbl_df", "tbl", "data.frame"))
Best Answer
Your data is complex and it may take several iterations of analysis and validation to come up with a good model. We cannot know whether a model fits a dataset from looking at its formula.
You don't mention this explicitly but the sample data suggests that the "participants" are babies since birth date ranges from May 2019 to Jan 2021. Is BMI a meaningful concept for babies?
In a way. The model assumes the relationship between zBMI (BMI z-score) and time is linear, with a different slope and intercept for each vaccine type. The model also includes an offset (a random intercept) for each baby, so the vaccine-specific line can be adjusted up or down. To an extent, the random intercepts account for the initial zBMI measurements.
On a related note, it's not clear what the BMI z-score is; I guess it's a standardization of BMI (by baby?), which is itself a function of length and weight, ie. growth. Since the outcome zBMI has a complex definition, the coefficients in the model might be difficult to interpret. Why not model BMI directly?
Since all measurements of a baby share the same random intercept, they are not independent: every two measurements are equally correlated, no matter how much time has elapsed between measurements. To let the correlation decay with time, you have to specify a different covariance structure for the within-baby measurements. You can do this with a mixed effects model of course. You can also use Generalized Least Squares (GLS); see
nlme::glsand in particular thecorrelationargument.You haven't included any covariates in the model; this is probably a mistake. Even if you are only interested in the relationship between vaccine type and BMI, you should include all predictors associated with a baby's growth: age, gestational age, length (I gather that with babies we measure length rather than height), weight, sex, etc. If these covariates explain variability in BMI, the bigger model will have smaller residual error and it will be possible to detect (smaller) differences due to vaccine type, if there are any. It might also help to log transform covariates, length and weight in particular.
Some babies received a second shot. You can't ignore this as you do now. A simple way to deal with this extra complexity is to exclude data points collected after a second vaccine shot.
However, the biggest challenge might be whether & how to incorporate post-treatment effects for the three vaccine types. The proposed model allows each vaccine to have a different linear relationship between zBMI and time; however, the slope and intercept of these lines are the same before and after the vaccine is administered. This is a strong assumption, particularly if there are more measurements before than after vaccination. I fit your model on the sample data to demonstrate why.
The two panels show observed and predicted BMI z-scores for 14 babies. On a high level, the model fails to capture a lot of the variability between the babies given the same vaccine (random slopes might help). More subtly, all babies contribute to the slope and intercept estimates of their vaccine type, whether they have been measured after vaccination or not. I suspect that, in addition to covariates, the model should include "change after vaccination" effects for each vaccine type. How to represent such an effect might be the most difficult question.