2
$\begingroup$

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"))
Improve this question
$\endgroup$

1 Answer 1

Reset to default
0
$\begingroup$

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?

Does the model control for pre-treatment zBMI?

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?

Does the random intercept account for the number of observations and their closeness in time to solve the autocorrelation problem?

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::gls and in particular the correlation argument.

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.

enter image description here

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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks for your detailed answer. $\endgroup$
    – Charlotte Juton
    Jul 13, 2022 at 15:31
  • $\begingroup$ zBMI is a standardised measure of babies' or child's BMI regardless of his age or sex, the tables can be found on the WHO's website. You said that every two measurements are equally correlated so why should I let the correlation decay with time? I agree with you regarding the covariates and the fact that we can't ignore the second vaccine and should exclude growth measurement taken. With this in mind, I was thinking about creating a dummy variable 'on/off treatment' to add to the model to account for the change after first vaccination alongside the time and group. Would that make any sense? $\endgroup$
    – Charlotte Juton
    Jul 13, 2022 at 15:31
  • 1
    $\begingroup$ Your question mentions autocorrelation, so I point out that with the random intercept the correlation structure of within subject measurements is simpler. (Of course, more complex is not always better.) A dummy variable for each treatment after vaccination may be a reasonable approach. These dummy variables would correspond to an offset in the intercept? Then there is the question of whether before vaccination all babies should have the same E{zBMI | covariates}. $\endgroup$
    – dipetkov
    Jul 13, 2022 at 20:12
  • $\begingroup$ The idea with the additional dummy variable that would be to create a triple interaction: zBMI ~ groupon/off treatmenttime, is that possible and accurate? does it represent the offset that you're talking about? thank you $\endgroup$
    – Charlotte Juton
    Jul 15, 2022 at 9:02
  • $\begingroup$ Three-way interactions can be hard to interpret. Making the interaction even more complex doesn't get around the fact that the models implies there are differences between babies on different vaccines before they are vaccinated. I find this assumption strange but I also guess your data is observational; otherwise you'll have zBMI measured at time 0. $\endgroup$
    – dipetkov
    Jul 15, 2022 at 9:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.