Skip to content

Multi Level Modeling Example Script using lmer4()

kirkvanacore edited this page May 24, 2021 · 1 revision
require(lmer4) # for modeling 
require(sjPlot) # for tables


# NULL MODEL 
MAP_Mnull  <-  lmer(SPRING_MAP_RIT ~  (1 | School/class_names),
              data = orland,
              na.action = na.omit)
summary(MAP_Mnull)
tab_model(MAP_Mnull)
tab_model(MAP_Mnull,
          title = c("MLM Predicting Spring MAP Scores"),
          pred.labels = c("Intercept"
                         
          )
          ,
          dv.labels = c("NULL Model"
          ),
          CSS = list( 
            css.table = 'font-family: Open Sans'
          ),
          #p.style = "asterisk",
          wrap.labels = 25,
          digits = 2)


# MODEL 1: add pretest
class(orland$FALL_MAP_RIT)
orland$FALL_MAP_RIT_mc <- orland$FALL_MAP_RIT - mean(orland$FALL_MAP_RIT, na.rm = T)
summary(orland$FALL_MAP_RIT_mc)

MAP_M1  <-  lmer(SPRING_MAP_RIT ~ FALL_MAP_RIT_mc + (1 | School/class_names),
                   data = orland,
                   na.action = na.omit)
summary(MAP_M1)
tab_model(MAP_Mnull, MAP_M1)

tab_model(MAP_Mnull,MAP_M1,
          title = c("MLM Predicting Spring MAP Scores"),
          pred.labels = c("Intercept",
                          "Fall MAP RIT (mean centered)"
                          
          )
          ,
          dv.labels = c("NULL Model",
                        "Model 1"
          ),
          CSS = list( 
            css.table = 'font-family: Open Sans'
          ),
          #p.style = "asterisk",
          wrap.labels = 25,
          digits = 2)

# MODEL 2: add treatment
MAP_M2  <-  lmer(SPRING_MAP_RIT ~ 
                   FALL_MAP_RIT_mc +
                   Treatment +
                   (1 | School/class_names),
                 data = orland,
                 na.action = na.omit)
summary(MAP_M2)
tab_model(MAP_Mnull,MAP_M1,MAP_M2,
          title = c("MLM Predicting Spring MAP Scores"),
          pred.labels = c("Intercept",
                          "Fall MAP RIT (mean centered)",
                          "Treatment: Core5"
                          
          )
          ,
          dv.labels = c("NULL Model",
                        "Model 1",
                        "Model 2"
          ),
          CSS = list( 
            css.table = 'font-family: Open Sans'
          ),
          #p.style = "asterisk",
          wrap.labels = 25,
          digits = 2)


# MODEL 3: add grade
MAP_M3  <-  lmer(SPRING_MAP_RIT ~ 
                   FALL_MAP_RIT_mc +
                   Treatment +
                   Grade +
                   (1 | School/class_names),
                 data = orland,
                 na.action = na.omit)
summary(MAP_M3)
tab_model(MAP_Mnull, MAP_M1, MAP_M2, MAP_M3)
tab_model(MAP_Mnull,MAP_M1,MAP_M2,MAP_M3,
          title = c("MLM Predicting Spring MAP Scores"),
          pred.labels = c("Intercept",
                          "Fall MAP RIT (mean centered)",
                          "Treatment: Core5",
                          "Grade 1",
                          "Grade 2",
                          "Grade 3",
                          "Grade 4",
                          "Grade 5"
          )
          ,
          dv.labels = c("NULL Model",
                        "Model 1",
                        "Model 2",
                        "Model 3"
          ),
          CSS = list( 
            css.table = 'font-family: Open Sans'
          ),
          #p.style = "asterisk",
          wrap.labels = 25,
          digits = 2,
          show.aic = T)

# MODEL 4: add interaction bw treatment and grade
MAP_M4  <-  lmer(SPRING_MAP_RIT ~ 
                   FALL_MAP_RIT_mc +
                   Treatment*Grade +
                   (1 | School/class_names),
                 data = orland,
                 na.action = na.omit)
summary(MAP_M4)
tab_model(MAP_Mnull, MAP_M1, MAP_M2, MAP_M3, MAP_M4)
tab_model(MAP_M4)
tab_model(MAP_Mnull,MAP_M1,MAP_M2,MAP_M3,MAP_M4,
          title = c("MLM Predicting Spring MAP Scores - Total Sample"),
          pred.labels = c("Intercept",
                          "Fall MAP RIT (mean centered)",
                          "Treatment: Core5",
                          "Grade 1",
                          "Grade 2",
                          "Grade 3",
                          "Grade 4",
                          "Grade 5",
                          "Treatment X Grade 1",
                          "Treatment X Grade 2",
                          "Treatment X Grade 3",
                          "Treatment X Grade 4",
                          "Treatment X Grade 5"
          )
          ,
          dv.labels = c("NULL Model",
                        "Model 1",
                        "Model 2",
                        "Model 3",
                        "Model 4"
          ),
          CSS = list( 
            css.table = 'font-family: Open Sans'
          ),
          #p.style = "asterisk",
          wrap.labels = 25,
          digits = 2,
          show.aic = T)


# MODEL 5: add Disability
MAP_M5  <-  lmer(SPRING_MAP_RIT ~ 
                   FALL_MAP_RIT_mc +
                   Treatment +
                   Grade +
                   SLD +
                   SLI +
                   DD +
                   AU +
                   OHI +
                   VI +
                   (1 | School/class_names),
                 data = orland,
                 na.action = na.omit)
summary(MAP_M5)
tab_model(MAP_Mnull, MAP_M1, MAP_M2, MAP_M5)
tab_model(MAP_M5)

# MODEL 6: add interaction bw treatment and Disability
MAP_M6  <-  lmer(SPRING_MAP_RIT ~ 
                   FALL_MAP_RIT_mc +
                   Treatment +
                   Grade +
                   Treatment*SLD +
                   Treatment*SLI +
                   Treatment*DD +
                   Treatment*AU +
                   Treatment*OHI +
                   Treatment*VI +
                   (1 | School/class_names),
                 data = orland,
                 na.action = na.omit)
summary(MAP_M6)
tab_model(MAP_M5, MAP_M6)
tab_model(MAP_M5, MAP_M6,
          title = c("MLM Predicting Spring MAP Scores - Total Sample"),
          pred.labels = c("Intercept",
                          "Fall MAP RIT (mean centered)",
                          "Treatment: Core5",
                          "Grade 1",
                          "Grade 2",
                          "Grade 3",
                          "Grade 4",
                          "Grade 5",
                          "Specific Learning Disability",
                          "Speech and Language Impairment",
                          "Developmental Delay",
                          "Autism",
                          "Other Health Impairment",
                          "Visual Imparement",
                          "Treatment X Specific Learning Disability",
                          "Treatment X Speech and Language Impairment",
                          "Treatment X Developmental Delay",
                          "Treatment X Autism",
                          "Treatment X Other Health Impairment",
                          "Treatment X Visual Imparement"
          )
          ,
          dv.labels = c(
                        "Model 5",
                        "Model 6"
          ),
          CSS = list( 
            css.table = 'font-family: Open Sans'
          ),
          #p.style = "asterisk",
          wrap.labels = 25,
          digits = 2,
          show.aic = T)

rm(list=setdiff(ls(), "orland"))

#### Subgroups ####
# IDEA Groups

# Specific Learning Disability
table(orland$SLD)
table(orland$Treatment, orland$SLD)
table(orland$Pullout_Dummy, orland$SLD)
MAP_M_SLD1  <-  lmer(SPRING_MAP_RIT ~ 
                   FALL_MAP_RIT_mc +
                   Treatment +
                   (1 | School/class_names),
                 data = orland[orland$SLD == 1, ],
                 na.action = na.omit)
summary(MAP_M_SLD1)
tab_model(MAP_M_SLD1)


MAP_M_SLD2  <-  lmer(SPRING_MAP_RIT ~ 
                       FALL_MAP_RIT_mc +
                       Treatment +
                       Grade +
                       (1 | School/class_names),
                     data = orland[orland$SLD == 1, ],
                     na.action = na.omit)
summary(MAP_M_SLD2)
tab_model(MAP_M_SLD2)


MAP_M_SLD3  <-  lmer(SPRING_MAP_RIT ~ 
                       FALL_MAP_RIT_mc +
                       Treatment*Grade +
                       (1 | School/class_names),
                     data = orland[orland$SLD == 1, ],
                     na.action = na.omit)
summary(MAP_M_SLD3)
tab_model(MAP_M_SLD1, MAP_M_SLD2, MAP_M_SLD3)


tab_model(MAP_M_SLD1, MAP_M_SLD2, MAP_M_SLD3,
          title = c("MLM Predicting Spring MAP Scores - Specific Learning Disability"),
          pred.labels = c("Intercept",
                          "Fall MAP RIT (mean centered)",
                          "Treatment: Core5",
                          "Grade 1",
                          "Grade 2",
                          "Grade 3",
                          "Grade 4",
                          "Grade 5",
                          "Treatment X Grade 1",
                          "Treatment X Grade 2",
                          "Treatment X Grade 3",
                          "Treatment X Grade 4"
          )
          ,
          dv.labels = c(
                        "Model 1",
                        "Model 2",
                        "Model 3"
          ),
          CSS = list( 
            css.table = 'font-family: Open Sans'
          ),
          #p.style = "asterisk",
          wrap.labels = 25,
          digits = 2,
          show.aic = T)
#### Total Sample ####

# NULL MODEL 
MAP_Mnull  <-  lmer(SPRING_MAP_RIT ~  (1 | School/class_names),
              data = orland,
              na.action = na.omit)
summary(MAP_Mnull)
tab_model(MAP_Mnull)
tab_model(MAP_Mnull,
          title = c("MLM Predicting Spring MAP Scores"),
          pred.labels = c("Intercept"
                         
          )
          ,
          dv.labels = c("NULL Model"
          ),
          CSS = list( 
            css.table = 'font-family: Open Sans'
          ),
          #p.style = "asterisk",
          wrap.labels = 25,
          digits = 2)


# MODEL 1: add pretest
class(orland$FALL_MAP_RIT)
orland$FALL_MAP_RIT_mc <- orland$FALL_MAP_RIT - mean(orland$FALL_MAP_RIT, na.rm = T)
summary(orland$FALL_MAP_RIT_mc)

MAP_M1  <-  lmer(SPRING_MAP_RIT ~ FALL_MAP_RIT_mc + (1 | School/class_names),
                   data = orland,
                   na.action = na.omit)
summary(MAP_M1)
tab_model(MAP_Mnull, MAP_M1)

tab_model(MAP_Mnull,MAP_M1,
          title = c("MLM Predicting Spring MAP Scores"),
          pred.labels = c("Intercept",
                          "Fall MAP RIT (mean centered)"
                          
          )
          ,
          dv.labels = c("NULL Model",
                        "Model 1"
          ),
          CSS = list( 
            css.table = 'font-family: Open Sans'
          ),
          #p.style = "asterisk",
          wrap.labels = 25,
          digits = 2)

# MODEL 2: add treatment
MAP_M2  <-  lmer(SPRING_MAP_RIT ~ 
                   FALL_MAP_RIT_mc +
                   Treatment +
                   (1 | School/class_names),
                 data = orland,
                 na.action = na.omit)
summary(MAP_M2)
tab_model(MAP_Mnull,MAP_M1,MAP_M2,
          title = c("MLM Predicting Spring MAP Scores"),
          pred.labels = c("Intercept",
                          "Fall MAP RIT (mean centered)",
                          "Treatment: Core5"
                          
          )
          ,
          dv.labels = c("NULL Model",
                        "Model 1",
                        "Model 2"
          ),
          CSS = list( 
            css.table = 'font-family: Open Sans'
          ),
          #p.style = "asterisk",
          wrap.labels = 25,
          digits = 2)


# MODEL 3: add grade
MAP_M3  <-  lmer(SPRING_MAP_RIT ~ 
                   FALL_MAP_RIT_mc +
                   Treatment +
                   Grade +
                   (1 | School/class_names),
                 data = orland,
                 na.action = na.omit)
summary(MAP_M3)
tab_model(MAP_Mnull, MAP_M1, MAP_M2, MAP_M3)
tab_model(MAP_Mnull,MAP_M1,MAP_M2,MAP_M3,
          title = c("MLM Predicting Spring MAP Scores"),
          pred.labels = c("Intercept",
                          "Fall MAP RIT (mean centered)",
                          "Treatment: Core5",
                          "Grade 1",
                          "Grade 2",
                          "Grade 3",
                          "Grade 4",
                          "Grade 5"
          )
          ,
          dv.labels = c("NULL Model",
                        "Model 1",
                        "Model 2",
                        "Model 3"
          ),
          CSS = list( 
            css.table = 'font-family: Open Sans'
          ),
          #p.style = "asterisk",
          wrap.labels = 25,
          digits = 2,
          show.aic = T)

# MODEL 4: add interaction bw treatment and grade
MAP_M4  <-  lmer(SPRING_MAP_RIT ~ 
                   FALL_MAP_RIT_mc +
                   Treatment*Grade +
                   (1 | School/class_names),
                 data = orland,
                 na.action = na.omit)
summary(MAP_M4)
tab_model(MAP_Mnull, MAP_M1, MAP_M2, MAP_M3, MAP_M4)
tab_model(MAP_M4)
tab_model(MAP_Mnull,MAP_M1,MAP_M2,MAP_M3,MAP_M4,
          title = c("MLM Predicting Spring MAP Scores - Total Sample"),
          pred.labels = c("Intercept",
                          "Fall MAP RIT (mean centered)",
                          "Treatment: Core5",
                          "Grade 1",
                          "Grade 2",
                          "Grade 3",
                          "Grade 4",
                          "Grade 5",
                          "Treatment X Grade 1",
                          "Treatment X Grade 2",
                          "Treatment X Grade 3",
                          "Treatment X Grade 4",
                          "Treatment X Grade 5"
          )
          ,
          dv.labels = c("NULL Model",
                        "Model 1",
                        "Model 2",
                        "Model 3",
                        "Model 4"
          ),
          CSS = list( 
            css.table = 'font-family: Open Sans'
          ),
          #p.style = "asterisk",
          wrap.labels = 25,
          digits = 2,
          show.aic = T)


# MODEL 5: add Disability
MAP_M5  <-  lmer(SPRING_MAP_RIT ~ 
                   FALL_MAP_RIT_mc +
                   Treatment +
                   Grade +
                   SLD +
                   SLI +
                   DD +
                   AU +
                   OHI +
                   VI +
                   (1 | School/class_names),
                 data = orland,
                 na.action = na.omit)
summary(MAP_M5)
tab_model(MAP_Mnull, MAP_M1, MAP_M2, MAP_M5)
tab_model(MAP_M5)

# MODEL 6: add interaction bw treatment and Disability
MAP_M6  <-  lmer(SPRING_MAP_RIT ~ 
                   FALL_MAP_RIT_mc +
                   Treatment +
                   Grade +
                   Treatment*SLD +
                   Treatment*SLI +
                   Treatment*DD +
                   Treatment*AU +
                   Treatment*OHI +
                   Treatment*VI +
                   (1 | School/class_names),
                 data = orland,
                 na.action = na.omit)
summary(MAP_M6)
tab_model(MAP_M5, MAP_M6)
tab_model(MAP_M5, MAP_M6,
          title = c("MLM Predicting Spring MAP Scores - Total Sample"),
          pred.labels = c("Intercept",
                          "Fall MAP RIT (mean centered)",
                          "Treatment: Core5",
                          "Grade 1",
                          "Grade 2",
                          "Grade 3",
                          "Grade 4",
                          "Grade 5",
                          "Specific Learning Disability",
                          "Speech and Language Impairment",
                          "Developmental Delay",
                          "Autism",
                          "Other Health Impairment",
                          "Visual Imparement",
                          "Treatment X Specific Learning Disability",
                          "Treatment X Speech and Language Impairment",
                          "Treatment X Developmental Delay",
                          "Treatment X Autism",
                          "Treatment X Other Health Impairment",
                          "Treatment X Visual Imparement"
          )
          ,
          dv.labels = c(
                        "Model 5",
                        "Model 6"
          ),
          CSS = list( 
            css.table = 'font-family: Open Sans'
          ),
          #p.style = "asterisk",
          wrap.labels = 25,
          digits = 2,
          show.aic = T)

rm(list=setdiff(ls(), "orland"))

#### Subgroups ####
# IDEA Groups

# Specific Learning Disability
table(orland$SLD)
table(orland$Treatment, orland$SLD)
table(orland$Pullout_Dummy, orland$SLD)
MAP_M_SLD1  <-  lmer(SPRING_MAP_RIT ~ 
                   FALL_MAP_RIT_mc +
                   Treatment +
                   (1 | School/class_names),
                 data = orland[orland$SLD == 1, ],
                 na.action = na.omit)
summary(MAP_M_SLD1)
tab_model(MAP_M_SLD1)


MAP_M_SLD2  <-  lmer(SPRING_MAP_RIT ~ 
                       FALL_MAP_RIT_mc +
                       Treatment +
                       Grade +
                       (1 | School/class_names),
                     data = orland[orland$SLD == 1, ],
                     na.action = na.omit)
summary(MAP_M_SLD2)
tab_model(MAP_M_SLD2)


MAP_M_SLD3  <-  lmer(SPRING_MAP_RIT ~ 
                       FALL_MAP_RIT_mc +
                       Treatment*Grade +
                       (1 | School/class_names),
                     data = orland[orland$SLD == 1, ],
                     na.action = na.omit)
summary(MAP_M_SLD3)
tab_model(MAP_M_SLD1, MAP_M_SLD2, MAP_M_SLD3)


tab_model(MAP_M_SLD1, MAP_M_SLD2, MAP_M_SLD3,
          title = c("MLM Predicting Spring MAP Scores - Specific Learning Disability"),
          pred.labels = c("Intercept",
                          "Fall MAP RIT (mean centered)",
                          "Treatment: Core5",
                          "Grade 1",
                          "Grade 2",
                          "Grade 3",
                          "Grade 4",
                          "Grade 5",
                          "Treatment X Grade 1",
                          "Treatment X Grade 2",
                          "Treatment X Grade 3",
                          "Treatment X Grade 4"
          )
          ,
          dv.labels = c(
                        "Model 1",
                        "Model 2",
                        "Model 3"
          ),
          CSS = list( 
            css.table = 'font-family: Open Sans'
          ),
          #p.style = "asterisk",
          wrap.labels = 25,
          digits = 2,
          show.aic = T)

Clone this wiki locally