This week, we are going to use the Kindergarten data from Project STAR, the Tennessee Class Size Experiment conducted from 1985-1990.

Load Some Packages to Help with the Analysis and Data Management:

suppressPackageStartupMessages(library(tidyverse))

package ‘tibble’ was built under R version 3.6.2package ‘dplyr’ was built under R version 3.6.2

suppressPackageStartupMessages(library(lme4))

Load and Explore the Data


projectSTAR <- haven::read_dta("projectSTAR.dta")
glimpse(projectSTAR)

Rows: 6,325
Columns: 27
$ stdntid       [3m[38;5;246m<dbl>[39m[23m 10001, 10133, 10246, 10263, 10266, 10275, 10281, 10282, 10285, 10286, 10287…
$ race          [3m[38;5;246m<dbl+lbl>[39m[23m 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1,…
$ gender        [3m[38;5;246m<dbl+lbl>[39m[23m 1, 1, 2, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 2, 2, 2, 2, 1, 1, 1, 1,…
$ FLAGSGK       [3m[38;5;246m<dbl+lbl>[39m[23m 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ flaggk        [3m[38;5;246m<dbl+lbl>[39m[23m 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1,…
$ gkclasstype   [3m[38;5;246m<dbl+lbl>[39m[23m 3, 3, 3, 1, 2, 3, 1, 3, 1, 2, 3, 3, 2, 3, 3, 3, 1, 2, 3, 3, 2, 1, 1, 3,…
$ gkschid       [3m[38;5;246m<dbl>[39m[23m 169229, 169280, 218562, 205492, 257899, 161176, 189382, 189382, 201449, 230…
$ gksurban      [3m[38;5;246m<dbl+lbl>[39m[23m 2, 2, 4, 2, 3, 3, 2, 2, 3, 3, 3, 2, 2, 2, 2, 3, 3, 3, 3, 3, 2, 1, 2, 3,…
$ gktchid       [3m[38;5;246m<dbl>[39m[23m 16922904, 16928003, 21856202, 20549204, 25789904, 16117602, 18938204, 18938…
$ gktgen        [3m[38;5;246m<dbl+lbl>[39m[23m 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,…
$ gktrace       [3m[38;5;246m<dbl+lbl>[39m[23m 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1,…
$ gkthighdegree [3m[38;5;246m<dbl+lbl>[39m[23m 2, 2, 3, 2, 3, 2, 3, 3, 3, 3, 2, 2, 3, 2, 2, 2, 2, 2, 2, 3, 2, 2, 3, 2,…
$ gktcareer     [3m[38;5;246m<dbl+lbl>[39m[23m 4, 4, 4, 4, 4, NA, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, NA, 2, 4, …
$ gktyears      [3m[38;5;246m<dbl>[39m[23m 5, 7, 8, 3, 12, 2, 7, 14, 4, 6, 11, 16, 12, 5, 17, 10, 6, 10, 13, 9, 18, 1,…
$ gkclasssize   [3m[38;5;246m<dbl>[39m[23m 24, 22, 17, 17, 24, 24, 13, 24, 14, 23, 23, 22, 20, 24, 23, 27, 17, 24, 22,…
$ gkfreelunch   [3m[38;5;246m<dbl+lbl>[39m[23m 2, 2, 2, 1, 2, 2, 2, 1, 1, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 1, 2, 2, 2, 2,…
$ gkrepeat      [3m[38;5;246m<dbl+lbl>[39m[23m 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ gkspeced      [3m[38;5;246m<dbl+lbl>[39m[23m 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,…
$ gkspecin      [3m[38;5;246m<dbl+lbl>[39m[23m 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1,…
$ gkpresent     [3m[38;5;246m<dbl>[39m[23m 161, 175, 132, 178, 170, 94, 160, 154, 172, 95, 163, 172, 180, 149, 173, NA…
$ gkabsent      [3m[38;5;246m<dbl>[39m[23m 19, 5, 28, 2, 10, 3, 2, 7, 8, 2, 17, 8, 0, 31, 7, NA, 20, 7, 2, 19, 24, 0, …
$ gktreadss     [3m[38;5;246m<dbl>[39m[23m NA, 427, 450, 483, 456, 411, 443, 448, 463, 472, 428, 545, 408, 422, 472, N…
$ gktmathss     [3m[38;5;246m<dbl>[39m[23m NA, 478, 494, 513, 513, 468, 473, 449, 520, 536, 484, 626, 454, 439, 528, N…
$ gktlistss     [3m[38;5;246m<dbl>[39m[23m NA, 509, 549, 554, 520, 571, 595, 540, 565, 595, NA, 622, 474, 536, 578, NA…
$ gkwordskillss [3m[38;5;246m<dbl>[39m[23m NA, 418, 444, 431, 468, 396, 444, 444, 480, 486, 423, 524, 410, 423, 458, N…
$ gkmotivraw    [3m[38;5;246m<dbl>[39m[23m 23, 24, 28, 27, 25, 24, NA, NA, 26, 27, 24, 24, 23, 28, 24, NA, 26, 25, NA,…
$ gkselfconcraw [3m[38;5;246m<dbl>[39m[23m 52, 53, 56, 61, 54, 55, NA, NA, 52, 61, 55, 49, 49, 59, 50, NA, 58, 45, NA,…

Part One: Data Exploration and Management

Use `str` to Figure out the 2/1 Coding Scheme and Variable Labels

The str function (short for structure) will show us any variables labels or values that have been attached to a variable. You can use this to find our how variables have been coded.

str(projectSTAR)

tibble [6,325 × 27] (S3: tbl_df/tbl/data.frame)
 $ stdntid      : num [1:6325] 10001 10133 10246 10263 10266 ...
  ..- attr(*, "label")= chr "STUDENT ID"
  ..- attr(*, "format.stata")= chr "%10.0g"
 $ race         : 'haven_labelled' num [1:6325] 1 2 1 2 1 1 1 1 1 1 ...
  ..- attr(*, "label")= chr "STUDENT RACE/ETHNICITY"
  ..- attr(*, "format.stata")= chr "%10.0g"
  ..- attr(*, "labels")= Named num [1:6] 1 2 3 4 5 6
  .. ..- attr(*, "names")= chr [1:6] "WHITE" "BLACK" "ASIAN" "HISPANIC" ...
 $ gender       : 'haven_labelled' num [1:6325] 1 1 2 2 1 1 1 1 2 1 ...
  ..- attr(*, "label")= chr "STUDENT GENDER"
  ..- attr(*, "format.stata")= chr "%10.0g"
  ..- attr(*, "labels")= Named num [1:2] 1 2
  .. ..- attr(*, "names")= chr [1:2] "MALE" "FEMALE"
 $ FLAGSGK      : 'haven_labelled' num [1:6325] 1 1 1 1 1 1 1 1 1 1 ...
  ..- attr(*, "label")= chr "IN STAR IN KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
  ..- attr(*, "labels")= Named num [1:2] 0 1
  .. ..- attr(*, "names")= chr [1:2] "NO" "YES"
 $ flaggk       : 'haven_labelled' num [1:6325] 0 1 1 1 1 1 1 1 1 1 ...
  ..- attr(*, "label")= chr "ACHIEVEMENT DATA AVAILABLE GRADE K"
  ..- attr(*, "format.stata")= chr "%10.0g"
  ..- attr(*, "labels")= Named num [1:2] 0 1
  .. ..- attr(*, "names")= chr [1:2] "NO" "YES"
 $ gkclasstype  : 'haven_labelled' num [1:6325] 3 3 3 1 2 3 1 3 1 2 ...
  ..- attr(*, "label")= chr "CLASSROOM TYPE KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
  ..- attr(*, "labels")= Named num [1:3] 1 2 3
  .. ..- attr(*, "names")= chr [1:3] "SMALL CLASS" "REGULAR CLASS" "REGULAR + AIDE CLASS"
 $ gkschid      : num [1:6325] 169229 169280 218562 205492 257899 ...
  ..- attr(*, "label")= chr "KINDERGARTEN SCHOOL ID"
  ..- attr(*, "format.stata")= chr "%10.0g"
 $ gksurban     : 'haven_labelled' num [1:6325] 2 2 4 2 3 3 2 2 3 3 ...
  ..- attr(*, "label")= chr "SCHOOL URBANICITY KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
  ..- attr(*, "labels")= Named num [1:4] 1 2 3 4
  .. ..- attr(*, "names")= chr [1:4] "INNER CITY" "SUBURBAN" "RURAL" "URBAN"
 $ gktchid      : num [1:6325] 16922904 16928003 21856202 20549204 25789904 ...
  ..- attr(*, "label")= chr "KINDERGARTEN TEACHER ID"
  ..- attr(*, "format.stata")= chr "%10.0g"
 $ gktgen       : 'haven_labelled' num [1:6325] 2 2 2 2 2 2 2 2 2 2 ...
  ..- attr(*, "label")= chr "TEACHER GENDER KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
  ..- attr(*, "labels")= Named num [1:2] 1 2
  .. ..- attr(*, "names")= chr [1:2] "MALE" "FEMALE"
 $ gktrace      : 'haven_labelled' num [1:6325] 2 1 1 1 1 1 1 1 1 1 ...
  ..- attr(*, "label")= chr "TEACHER RACE/ETHNICITY KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
  ..- attr(*, "labels")= Named num [1:6] 1 2 3 4 5 6
  .. ..- attr(*, "names")= chr [1:6] "WHITE" "BLACK" "ASIAN" "HISPANIC" ...
 $ gkthighdegree: 'haven_labelled' num [1:6325] 2 2 3 2 3 2 3 3 3 3 ...
  ..- attr(*, "label")= chr "TEACHER HIGHEST DEGREE KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
  ..- attr(*, "labels")= Named num [1:6] 1 2 3 4 5 6
  .. ..- attr(*, "names")= chr [1:6] "ASSOCIATES" "BACHELORS" "MASTERS" "MASTERS +" ...
 $ gktcareer    : 'haven_labelled' num [1:6325] 4 4 4 4 4 NA 4 4 4 4 ...
  ..- attr(*, "label")= chr "TEACHER CAREER LADDER LEVEL KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
  ..- attr(*, "labels")= Named num [1:7] 1 2 3 4 5 6 7
  .. ..- attr(*, "names")= chr [1:7] "CHOSE NO TO BE ON CAREER LADDER" "APPRENTICE" "PROBATION" "LADDER LEVEL 1" ...
 $ gktyears     : num [1:6325] 5 7 8 3 12 2 7 14 4 6 ...
  ..- attr(*, "label")= chr "YEARS OF TOTAL TEACHING EXPERIENCE KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
 $ gkclasssize  : num [1:6325] 24 22 17 17 24 24 13 24 14 23 ...
  ..- attr(*, "label")= chr "CLASS SIZE KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
 $ gkfreelunch  : 'haven_labelled' num [1:6325] 2 2 2 1 2 2 2 1 1 2 ...
  ..- attr(*, "label")= chr "FREE/REDUCED LUNCH STATUS KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
  ..- attr(*, "labels")= Named num [1:2] 1 2
  .. ..- attr(*, "names")= chr [1:2] "FREE LUNCH" "NON-FREE LUNCH"
 $ gkrepeat     : 'haven_labelled' num [1:6325] 1 1 1 1 1 1 1 1 1 1 ...
  ..- attr(*, "label")= chr "REPEATING KINDERGARTEN IN 1985-1986 SCHOOL YEAR"
  ..- attr(*, "format.stata")= chr "%10.0g"
  ..- attr(*, "labels")= Named num [1:2] 1 2
  .. ..- attr(*, "names")= chr [1:2] "YES, PROMOTION RECOMMENDED" "NO, PROMOTION NOT RECOMMENDED"
 $ gkspeced     : 'haven_labelled' num [1:6325] 2 2 2 2 2 2 2 2 2 2 ...
  ..- attr(*, "label")= chr "SPECIAL EDUCATION STATUS KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
  ..- attr(*, "labels")= Named num [1:2] 1 2
  .. ..- attr(*, "names")= chr [1:2] "YES" "NO"
 $ gkspecin     : 'haven_labelled' num [1:6325] 2 2 2 2 2 2 2 2 2 2 ...
  ..- attr(*, "label")= chr "PULLED OUT FOR SPECIAL INSTRUCTION KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
  ..- attr(*, "labels")= Named num [1:2] 1 2
  .. ..- attr(*, "names")= chr [1:2] "YES" "NO"
 $ gkpresent    : num [1:6325] 161 175 132 178 170 94 160 154 172 95 ...
  ..- attr(*, "label")= chr "DAYS PRESENT AT SCHOOL KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
 $ gkabsent     : num [1:6325] 19 5 28 2 10 3 2 7 8 2 ...
  ..- attr(*, "label")= chr "DAYS ABSENT FROM SCHOOL KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
 $ gktreadss    : num [1:6325] NA 427 450 483 456 411 443 448 463 472 ...
  ..- attr(*, "label")= chr "TOTAL READING SCALE SCORE SAT KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
 $ gktmathss    : num [1:6325] NA 478 494 513 513 468 473 449 520 536 ...
  ..- attr(*, "label")= chr "TOTAL MATH SCALE SCORE SAT KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
 $ gktlistss    : num [1:6325] NA 509 549 554 520 571 595 540 565 595 ...
  ..- attr(*, "label")= chr "TOTAL LISTENING SCALE SCORE SAT KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
 $ gkwordskillss: num [1:6325] NA 418 444 431 468 396 444 444 480 486 ...
  ..- attr(*, "label")= chr "WORD STUDY SKILLS SCALE SCORE SAT KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
 $ gkmotivraw   : num [1:6325] 23 24 28 27 25 24 NA NA 26 27 ...
  ..- attr(*, "label")= chr "MOTIVATION RAW SCORE SCAMIN KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"
 $ gkselfconcraw: num [1:6325] 52 53 56 61 54 55 NA NA 52 61 ...
  ..- attr(*, "label")= chr "SELF-CONCEPT RAW SCORE SCAMIN KINDERGARTEN"
  ..- attr(*, "format.stata")= chr "%10.0g"

Create a Clean Dataset Ready for Modeling

star.clean <- projectSTAR %>%
  mutate(.,
         schoolid = gkschid,
         classid = gktchid,
         read = gktreadss,
         classtype = as_factor(gkclasstype),
         years_exp = gktyears,
         urbanicity = as_factor(gksurban),
         stu_frl = as_factor(gkfreelunch)
         ) %>%
  group_by(classid) %>% # Create a new variable, which is the average FRL by classroom:
  mutate(.,
            class_frl = mean(gkfreelunch, na.rm = TRUE) - 1
         ) %>%
  ungroup() %>%
  group_by(schoolid) %>% # Create a new variable, which is the average FRL by school:
  mutate(.,
            school_frl = mean(gkfreelunch, na.rm = TRUE) - 1
         ) %>%
  ungroup() %>%
  select(.,
         schoolid,
         classid,
         read,
         classtype,
         years_exp,
         urbanicity,
         stu_frl,
         class_frl,
         school_frl)

glimpse(star.clean)

Rows: 6,325
Columns: 9
$ schoolid   [3m[38;5;246m<dbl>[39m[23m 169229, 169280, 218562, 205492, 257899, 161176, 189382, 189382, 201449, 230612…
$ classid    [3m[38;5;246m<dbl>[39m[23m 16922904, 16928003, 21856202, 20549204, 25789904, 16117602, 18938204, 18938203…
$ read       [3m[38;5;246m<dbl>[39m[23m NA, 427, 450, 483, 456, 411, 443, 448, 463, 472, 428, 545, 408, 422, 472, NA, …
$ classtype  [3m[38;5;246m<fct>[39m[23m REGULAR + AIDE CLASS, REGULAR + AIDE CLASS, REGULAR + AIDE CLASS, SMALL CLASS,…
$ years_exp  [3m[38;5;246m<dbl>[39m[23m 5, 7, 8, 3, 12, 2, 7, 14, 4, 6, 11, 16, 12, 5, 17, 10, 6, 10, 13, 9, 18, 1, 7,…
$ urbanicity [3m[38;5;246m<fct>[39m[23m SUBURBAN, SUBURBAN, URBAN, SUBURBAN, RURAL, RURAL, SUBURBAN, SUBURBAN, RURAL, …
$ stu_frl    [3m[38;5;246m<fct>[39m[23m NON-FREE LUNCH, NON-FREE LUNCH, NON-FREE LUNCH, FREE LUNCH, NON-FREE LUNCH, NO…
$ class_frl  [3m[38;5;246m<dbl>[39m[23m 0.8750000, 0.9545455, 0.7647059, 0.2352941, 0.4583333, 0.6086957, 0.8461538, 0…
$ school_frl [3m[38;5;246m<dbl>[39m[23m 0.86206897, 0.86206897, 0.71666667, 0.15189873, 0.62857143, 0.63291139, 0.8918…

This is a great example of when you might “aggregate up” an individual measure to see if it is a significant predictor at level 2 or level 3. We did this in the past with the task significance measure (TSIG) and the lq2002 dataset.

The above command has two parts. First, we use the group_by command from our friend, the dplyr package, to group all responses in their groups determined by classid. This is a first step to declare that all the operations that follow are to be performed at the cluster (classroom) level. Second, we use the mutate command, which we have used before, to create new variables that are equal to the mean of the original variables. Since we just grouped by classid, this means we get the mean of the variables by classroom.

We have to get just a little fancier - since the values of gkfreelunch are 1 for free lunch and 2 for non-free lunch, first we take the mean of that variable. Then, we subtract 1 from that value to give us the proportion of students who are non-free lunch. This gives us the new variable class_frl.

Finally, we do the same thing at the school level, grouping by school_id, and then performing the same operation on gkfreelunch to give us school_frl, which is the proportion of students at each school who are non- free lunch.

Other than that, we are just tidying and renaming first, and then at the bottom selecting only the variables we need for the models we will run.

Part Two: The Null 3-Level Model

Run a null 3-level model for reading scores


model.0 <- lmer(read ~ (1|schoolid) + (1|classid), data = star.clean)
summary(model.0)

Linear mixed model fit by REML ['lmerMod']
Formula: read ~ (1 | schoolid) + (1 | classid)
   Data: star.clean

REML criterion at convergence: 55139.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.6315 -0.6013 -0.0966  0.4466  7.3775 

Random effects:
 Groups   Name        Variance Std.Dev.
 classid  (Intercept) 115.0    10.72   
 schoolid (Intercept) 169.8    13.03   
 Residual             728.2    26.99   
Number of obs: 5789, groups:  classid, 325; schoolid, 79

Fixed effects:
            Estimate Std. Error t value
(Intercept)  437.049      1.633   267.6

Calculate L2 ICC (Classroom)


ICC.class <- 114.9/(728.2 + 169.8 + 114.9)
ICC.class

[1] 0.1134367

Calculate L3 ICC (School)


ICC.school <- 169.8/(728.2 + 169.8 + 114.9)
ICC.school

[1] 0.1676375

Part Three: Add Predictors at Each Level

Add a student-level predictor, FRL


model.1 <- lmer(read ~ stu_frl + (1|schoolid) + (1|classid), data = star.clean)
summary(model.1)

Linear mixed model fit by REML ['lmerMod']
Formula: read ~ stu_frl + (1 | schoolid) + (1 | classid)
   Data: star.clean

REML criterion at convergence: 54635.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8204 -0.5915 -0.0926  0.4431  7.3227 

Random effects:
 Groups   Name        Variance Std.Dev.
 classid  (Intercept) 113.4    10.65   
 schoolid (Intercept) 142.4    11.93   
 Residual             687.5    26.22   
Number of obs: 5771, groups:  classid, 325; schoolid, 79

Fixed effects:
                      Estimate Std. Error t value
(Intercept)           428.8394     1.5808  271.28
stu_frlNON-FREE LUNCH  15.6106     0.8338   18.72

Correlation of Fixed Effects:
            (Intr)
s_NON-FREEL -0.278

Add teacher-level predictor, yrs. of experience and class type


model.2 <- lmer(read ~ stu_frl + years_exp + classtype + (1|schoolid) + (1|classid), data = star.clean)
summary(model.2)

Linear mixed model fit by REML ['lmerMod']
Formula: read ~ stu_frl + years_exp + classtype + (1 | schoolid) + (1 |      classid)
   Data: star.clean

REML criterion at convergence: 54432.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8496 -0.5895 -0.0972  0.4409  7.2503 

Random effects:
 Groups   Name        Variance Std.Dev.
 classid  (Intercept) 102.4    10.12   
 schoolid (Intercept) 141.0    11.88   
 Residual             689.1    26.25   
Number of obs: 5751, groups:  classid, 324; schoolid, 79

Fixed effects:
                              Estimate Std. Error t value
(Intercept)                   429.6765     2.1265 202.063
stu_frlNON-FREE LUNCH          15.5597     0.8362  18.607
years_exp                       0.3125     0.1262   2.476
classtypeREGULAR CLASS         -6.2966     1.6306  -3.861
classtypeREGULAR + AIDE CLASS  -5.3194     1.6309  -3.262

Correlation of Fixed Effects:
            (Intr) s_NONL yrs_xp cREGUC
s_NON-FREEL -0.200                     
years_exp   -0.532 -0.017              
cREGULARCLA -0.360  0.000  0.003       
cREGULAR+AC -0.324  0.009 -0.065  0.464

Add a school-level predictor, urbanicity of school


model.3 <- lmer(read ~ stu_frl + years_exp + classtype + urbanicity + (1|schoolid) + (1|classid), data = star.clean)
summary(model.3)

Linear mixed model fit by REML ['lmerMod']
Formula: read ~ stu_frl + years_exp + classtype + urbanicity + (1 | schoolid) +      (1 | classid)
   Data: star.clean

REML criterion at convergence: 54417.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8510 -0.5896 -0.0971  0.4411  7.2542 

Random effects:
 Groups   Name        Variance Std.Dev.
 classid  (Intercept) 102.4    10.12   
 schoolid (Intercept) 146.1    12.09   
 Residual             689.0    26.25   
Number of obs: 5751, groups:  classid, 324; schoolid, 79

Fixed effects:
                              Estimate Std. Error t value
(Intercept)                   427.2909     3.6298 117.716
stu_frlNON-FREE LUNCH          15.4736     0.8422  18.374
years_exp                       0.3052     0.1265   2.412
classtypeREGULAR CLASS         -6.2813     1.6312  -3.851
classtypeREGULAR + AIDE CLASS  -5.2977     1.6314  -3.247
urbanicitySUBURBAN              3.8798     4.6802   0.829
urbanicityRURAL                 2.5357     4.0529   0.626
urbanicityURBAN                 4.3470     6.1330   0.709

Correlation of Fixed Effects:
            (Intr) s_NONL yrs_xp cREGUC cRE+AC uSUBUR uRURAL
s_NON-FREEL -0.026                                          
years_exp   -0.266 -0.011                                   
cREGULARCLA -0.214 -0.001  0.003                            
cREGULAR+AC -0.203  0.008 -0.066  0.464                     
urbSUBURBAN -0.657 -0.109 -0.048  0.010  0.013              
urbnctRURAL -0.758 -0.098 -0.053 -0.003  0.013  0.614       
urbnctURBAN -0.505 -0.063 -0.030  0.012  0.013  0.405  0.466

Add in FRL at all three levels (student, teacher, school)


model.4 <- lmer(read ~ stu_frl +  class_frl + school_frl +  years_exp + classtype + urbanicity + (1|schoolid) + (1|classid), data = star.clean)
summary(model.4)

Linear mixed model fit by REML ['lmerMod']
Formula: read ~ stu_frl + class_frl + school_frl + years_exp + classtype +  
    urbanicity + (1 | schoolid) + (1 | classid)
   Data: star.clean

REML criterion at convergence: 54405.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8519 -0.5895 -0.0974  0.4423  7.2565 

Random effects:
 Groups   Name        Variance Std.Dev.
 classid  (Intercept) 103.0    10.15   
 schoolid (Intercept) 148.5    12.18   
 Residual             689.0    26.25   
Number of obs: 5751, groups:  classid, 324; schoolid, 79

Fixed effects:
                              Estimate Std. Error t value
(Intercept)                   427.2476     3.8148 111.996
stu_frlNON-FREE LUNCH          15.4340     0.8516  18.124
class_frl                       2.5997     7.1910   0.362
school_frl                     -2.1258    11.5341  -0.184
years_exp                       0.3018     0.1270   2.376
classtypeREGULAR CLASS         -6.2803     1.6345  -3.842
classtypeREGULAR + AIDE CLASS  -5.2600     1.6379  -3.211
urbanicitySUBURBAN              3.6239     7.2776   0.498
urbanicityRURAL                 2.3449     5.9095   0.397
urbanicityURBAN                 4.1521     7.4252   0.559

Correlation of Fixed Effects:
            (Intr) s_NONL clss_f schl_f yrs_xp cREGUC cRE+AC uSUBUR uRURAL
s_NON-FREEL  0.001                                                        
class_frl    0.000 -0.118                                                 
school_frl  -0.230  0.001 -0.606                                          
years_exp   -0.252 -0.003 -0.061  0.034                                   
cREGULARCLA -0.200  0.000  0.005 -0.013  0.003                            
cREGULAR+AC -0.193  0.000  0.065 -0.038 -0.070  0.464                     
urbSUBURBAN -0.188 -0.001 -0.006 -0.603 -0.027  0.016  0.007              
urbnctRURAL -0.293 -0.002 -0.003 -0.574 -0.033  0.007  0.008  0.826       
urbnctURBAN -0.242 -0.002 -0.005 -0.439 -0.022  0.016  0.010  0.641  0.669

Part Four: Try Adding a Random Slope

Add a random slope for years experience at level 3 (school)


model.5 <- lmer(read ~ stu_frl + years_exp + classtype + urbanicity + (1|classid) + (years_exp|schoolid), data = star.clean)

Model failed to converge with max|grad| = 0.0371295 (tol = 0.002, component 1)

summary(model.5)

Linear mixed model fit by REML ['lmerMod']
Formula: read ~ stu_frl + years_exp + classtype + urbanicity + (1 | classid) +  
    (years_exp | schoolid)
   Data: star.clean

REML criterion at convergence: 54411.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8570 -0.5877 -0.1013  0.4444  7.3116 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 classid  (Intercept)  95.3880  9.7667      
 schoolid (Intercept)  85.4100  9.2418      
          years_exp     0.1933  0.4396  0.56
 Residual             689.2198 26.2530      
Number of obs: 5751, groups:  classid, 324; schoolid, 79

Fixed effects:
                              Estimate Std. Error t value
(Intercept)                   426.7155     3.2799 130.102
stu_frlNON-FREE LUNCH          15.4281     0.8416  18.331
years_exp                       0.3063     0.1372   2.233
classtypeREGULAR CLASS         -6.5457     1.6088  -4.069
classtypeREGULAR + AIDE CLASS  -5.7739     1.6060  -3.595
urbanicitySUBURBAN              5.2602     4.3975   1.196
urbanicityRURAL                 3.0632     3.8078   0.804
urbanicityURBAN                 5.2028     5.8291   0.893

Correlation of Fixed Effects:
            (Intr) s_NONL yrs_xp cREGUC cRE+AC uSUBUR uRURAL
s_NON-FREEL -0.028                                          
years_exp   -0.184 -0.016                                   
cREGULARCLA -0.224 -0.002  0.000                            
cREGULAR+AC -0.203  0.008 -0.059  0.464                     
urbSUBURBAN -0.660 -0.111 -0.042  0.009 -0.009              
urbnctRURAL -0.753 -0.103 -0.090 -0.008 -0.001  0.593       
urbnctURBAN -0.491 -0.069 -0.046 -0.021 -0.006  0.387  0.448
convergence code: 0
Model failed to converge with max|grad| = 0.0371295 (tol = 0.002, component 1)

Should We Include That Random Slope? Comparing `model.4` with `model.5`:

anova(model.4, model.5)

refitting model(s) with ML (instead of REML)

Data: star.clean
Models:
model.4: read ~ stu_frl + class_frl + school_frl + years_exp + classtype + 
model.4:     urbanicity + (1 | schoolid) + (1 | classid)
model.5: read ~ stu_frl + years_exp + classtype + urbanicity + (1 | classid) + 
model.5:     (years_exp | schoolid)
        Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
model.4 13 54465 54551 -27219    54439                             
model.5 13 54458 54545 -27216    54432 6.4427      0  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Using the `modelsummary` and `broom.mixed` Packages to Organize Your Results:

library(modelsummary)

package ‘modelsummary’ was built under R version 3.6.2

library(broom.mixed)

Registered S3 methods overwritten by 'broom.mixed':
  method         from 
  augment.lme    broom
  augment.merMod broom
  glance.lme     broom
  glance.merMod  broom
  glance.stanreg broom
  tidy.brmsfit   broom
  tidy.gamlss    broom
  tidy.lme       broom
  tidy.merMod    broom
  tidy.rjags     broom
  tidy.stanfit   broom
  tidy.stanreg   broom

models <- list(model.0, model.1, model.2, model.3, model.4, model.5)
modelsummary(models, output = "default")

`data_frame()` is deprecated as of tibble 1.1.0.
Please use `tibble()` instead.
[90mThis warning is displayed once every 8 hours.[39m
[90mCall `lifecycle::last_warnings()` to see where this warning was generated.[39m

	Model 1	Model 2	Model 3	Model 4	Model 5	Model 6
(Intercept)	437.049	428.839	429.677	427.291	427.248	426.715
	(1.633)	(1.581)	(2.126)	(3.630)	(3.815)	(3.280)
sd__(Intercept)	10.722	10.651	10.119	10.120	10.147	9.767
sd__(Intercept)	10.722	10.651	10.119	10.120	10.147	9.242
sd__(Intercept)	10.722	10.651	10.119	10.120	12.184	9.767
sd__(Intercept)	10.722	10.651	10.119	10.120	12.184	9.242
sd__(Intercept)	10.722	10.651	10.119	12.088	10.147	9.767
sd__(Intercept)	10.722	10.651	10.119	12.088	10.147	9.242
sd__(Intercept)	10.722	10.651	10.119	12.088	12.184	9.767
sd__(Intercept)	10.722	10.651	10.119	12.088	12.184	9.242
sd__(Intercept)	10.722	10.651	11.876	10.120	10.147	9.767
sd__(Intercept)	10.722	10.651	11.876	10.120	10.147	9.242
sd__(Intercept)	10.722	10.651	11.876	10.120	12.184	9.767
sd__(Intercept)	10.722	10.651	11.876	10.120	12.184	9.242
sd__(Intercept)	10.722	10.651	11.876	12.088	10.147	9.767
sd__(Intercept)	10.722	10.651	11.876	12.088	10.147	9.242
sd__(Intercept)	10.722	10.651	11.876	12.088	12.184	9.767
sd__(Intercept)	10.722	10.651	11.876	12.088	12.184	9.242
sd__(Intercept)	10.722	11.934	10.119	10.120	10.147	9.767
sd__(Intercept)	10.722	11.934	10.119	10.120	10.147	9.242
sd__(Intercept)	10.722	11.934	10.119	10.120	12.184	9.767
sd__(Intercept)	10.722	11.934	10.119	10.120	12.184	9.242
sd__(Intercept)	10.722	11.934	10.119	12.088	10.147	9.767
sd__(Intercept)	10.722	11.934	10.119	12.088	10.147	9.242
sd__(Intercept)	10.722	11.934	10.119	12.088	12.184	9.767
sd__(Intercept)	10.722	11.934	10.119	12.088	12.184	9.242
sd__(Intercept)	10.722	11.934	11.876	10.120	10.147	9.767
sd__(Intercept)	10.722	11.934	11.876	10.120	10.147	9.242
sd__(Intercept)	10.722	11.934	11.876	10.120	12.184	9.767
sd__(Intercept)	10.722	11.934	11.876	10.120	12.184	9.242
sd__(Intercept)	10.722	11.934	11.876	12.088	10.147	9.767
sd__(Intercept)	10.722	11.934	11.876	12.088	10.147	9.242
sd__(Intercept)	10.722	11.934	11.876	12.088	12.184	9.767
sd__(Intercept)	10.722	11.934	11.876	12.088	12.184	9.242
sd__(Intercept)	13.032	10.651	10.119	10.120	10.147	9.767
sd__(Intercept)	13.032	10.651	10.119	10.120	10.147	9.242
sd__(Intercept)	13.032	10.651	10.119	10.120	12.184	9.767
sd__(Intercept)	13.032	10.651	10.119	10.120	12.184	9.242
sd__(Intercept)	13.032	10.651	10.119	12.088	10.147	9.767
sd__(Intercept)	13.032	10.651	10.119	12.088	10.147	9.242
sd__(Intercept)	13.032	10.651	10.119	12.088	12.184	9.767
sd__(Intercept)	13.032	10.651	10.119	12.088	12.184	9.242
sd__(Intercept)	13.032	10.651	11.876	10.120	10.147	9.767
sd__(Intercept)	13.032	10.651	11.876	10.120	10.147	9.242
sd__(Intercept)	13.032	10.651	11.876	10.120	12.184	9.767
sd__(Intercept)	13.032	10.651	11.876	10.120	12.184	9.242
sd__(Intercept)	13.032	10.651	11.876	12.088	10.147	9.767
sd__(Intercept)	13.032	10.651	11.876	12.088	10.147	9.242
sd__(Intercept)	13.032	10.651	11.876	12.088	12.184	9.767
sd__(Intercept)	13.032	10.651	11.876	12.088	12.184	9.242
sd__(Intercept)	13.032	11.934	10.119	10.120	10.147	9.767
sd__(Intercept)	13.032	11.934	10.119	10.120	10.147	9.242
sd__(Intercept)	13.032	11.934	10.119	10.120	12.184	9.767
sd__(Intercept)	13.032	11.934	10.119	10.120	12.184	9.242
sd__(Intercept)	13.032	11.934	10.119	12.088	10.147	9.767
sd__(Intercept)	13.032	11.934	10.119	12.088	10.147	9.242
sd__(Intercept)	13.032	11.934	10.119	12.088	12.184	9.767
sd__(Intercept)	13.032	11.934	10.119	12.088	12.184	9.242
sd__(Intercept)	13.032	11.934	11.876	10.120	10.147	9.767
sd__(Intercept)	13.032	11.934	11.876	10.120	10.147	9.242
sd__(Intercept)	13.032	11.934	11.876	10.120	12.184	9.767
sd__(Intercept)	13.032	11.934	11.876	10.120	12.184	9.242
sd__(Intercept)	13.032	11.934	11.876	12.088	10.147	9.767
sd__(Intercept)	13.032	11.934	11.876	12.088	10.147	9.242
sd__(Intercept)	13.032	11.934	11.876	12.088	12.184	9.767
sd__(Intercept)	13.032	11.934	11.876	12.088	12.184	9.242
































































sd__Observation	26.986	26.221	26.250	26.250	26.249	26.253

stu_frlNON-FREE LUNCH		15.611	15.560	15.474	15.434	15.428
		(0.834)	(0.836)	(0.842)	(0.852)	(0.842)
classtypeREGULAR + AIDE CLASS			-5.319	-5.298	-5.260	-5.774
			(1.631)	(1.631)	(1.638)	(1.606)
classtypeREGULAR CLASS			-6.297	-6.281	-6.280	-6.546
			(1.631)	(1.631)	(1.634)	(1.609)
years_exp			0.313	0.305	0.302	0.306
			(0.126)	(0.127)	(0.127)	(0.137)
urbanicityRURAL				2.536	2.345	3.063
				(4.053)	(5.910)	(3.808)
urbanicitySUBURBAN				3.880	3.624	5.260
				(4.680)	(7.278)	(4.398)
urbanicityURBAN				4.347	4.152	5.203
				(6.133)	(7.425)	(5.829)
class_frl					2.600
					(7.191)
school_frl					-2.126
					(11.534)
cor__(Intercept).years_exp						0.565

sd__years_exp						0.440

AIC	55147.6	54645.9	54448.7	54439.6	54431.5	54437.3
BIC	55174.3	54679.2	54502.0	54512.9	54518.0	54523.9
Log.Lik.	-27569.798	-27317.936	-27216.369	-27208.818	-27202.728	-27205.669
REMLcrit	55139.596	54635.872	54432.739	54417.635	54405.456	54411.337

HTML Version That You Can Open in Word:

modelsummary(models, output = 'msum.html', title = 'MLM Estimates')

[WARNING] This document format requires a nonempty <title> element.
  Please specify either 'title' or 'pagetitle' in the metadata.
  Falling back to 'msum'

LS0tCnRpdGxlOiAnTXVsdGlsZXZlbCBNb2RlbGluZywgTW9kdWxlIDc6IFRocmVlLUxldmVsIE1vZGVscycKYXV0aG9yOiAnRHIuIEJyb2RhJwpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpUaGlzIHdlZWssIHdlIGFyZSBnb2luZyB0byB1c2UgdGhlIEtpbmRlcmdhcnRlbiBkYXRhIGZyb20gUHJvamVjdCBTVEFSLCB0aGUgVGVubmVzc2VlIENsYXNzIFNpemUgRXhwZXJpbWVudCBjb25kdWN0ZWQgZnJvbSAxOTg1LTE5OTAuCgojIExvYWQgU29tZSBQYWNrYWdlcyB0byBIZWxwIHdpdGggdGhlIEFuYWx5c2lzIGFuZCBEYXRhIE1hbmFnZW1lbnQ6CmBgYHtyfQpzdXBwcmVzc1BhY2thZ2VTdGFydHVwTWVzc2FnZXMobGlicmFyeSh0aWR5dmVyc2UpKQpzdXBwcmVzc1BhY2thZ2VTdGFydHVwTWVzc2FnZXMobGlicmFyeShsbWU0KSkKYGBgCgojIExvYWQgYW5kIEV4cGxvcmUgdGhlIERhdGEKYGBge3J9Cgpwcm9qZWN0U1RBUiA8LSBoYXZlbjo6cmVhZF9kdGEoInByb2plY3RTVEFSLmR0YSIpCmdsaW1wc2UocHJvamVjdFNUQVIpCgpgYGAKCiMgUGFydCBPbmU6IERhdGEgRXhwbG9yYXRpb24gYW5kIE1hbmFnZW1lbnQKCiMjIFVzZSBgc3RyYCB0byBGaWd1cmUgb3V0IHRoZSAyLzEgQ29kaW5nIFNjaGVtZSBhbmQgVmFyaWFibGUgTGFiZWxzClRoZSBgc3RyYCBmdW5jdGlvbiAoc2hvcnQgZm9yIHN0cnVjdHVyZSkgd2lsbCBzaG93IHVzIGFueSB2YXJpYWJsZXMgbGFiZWxzIG9yIHZhbHVlcyB0aGF0IGhhdmUgYmVlbiBhdHRhY2hlZCB0byBhIHZhcmlhYmxlLiBZb3UgY2FuIHVzZSB0aGlzIHRvIGZpbmQgb3VyIGhvdyB2YXJpYWJsZXMgaGF2ZSBiZWVuIGNvZGVkLgpgYGB7cn0Kc3RyKHByb2plY3RTVEFSKQoKYGBgCgojIyBDcmVhdGUgYSBDbGVhbiBEYXRhc2V0IFJlYWR5IGZvciBNb2RlbGluZwpgYGB7cn0Kc3Rhci5jbGVhbiA8LSBwcm9qZWN0U1RBUiAlPiUKICBtdXRhdGUoLiwKICAgICAgICAgc2Nob29saWQgPSBna3NjaGlkLAogICAgICAgICBjbGFzc2lkID0gZ2t0Y2hpZCwKICAgICAgICAgcmVhZCA9IGdrdHJlYWRzcywKICAgICAgICAgY2xhc3N0eXBlID0gYXNfZmFjdG9yKGdrY2xhc3N0eXBlKSwKICAgICAgICAgeWVhcnNfZXhwID0gZ2t0eWVhcnMsCiAgICAgICAgIHVyYmFuaWNpdHkgPSBhc19mYWN0b3IoZ2tzdXJiYW4pLAogICAgICAgICBzdHVfZnJsID0gYXNfZmFjdG9yKGdrZnJlZWx1bmNoKQogICAgICAgICApICU+JQogIGdyb3VwX2J5KGNsYXNzaWQpICU+JSAjIENyZWF0ZSBhIG5ldyB2YXJpYWJsZSwgd2hpY2ggaXMgdGhlIGF2ZXJhZ2UgRlJMIGJ5IGNsYXNzcm9vbToKICBtdXRhdGUoLiwKICAgICAgICAgICAgY2xhc3NfZnJsID0gbWVhbihna2ZyZWVsdW5jaCwgbmEucm0gPSBUUlVFKSAtIDEKICAgICAgICAgKSAlPiUKICB1bmdyb3VwKCkgJT4lCiAgZ3JvdXBfYnkoc2Nob29saWQpICU+JSAjIENyZWF0ZSBhIG5ldyB2YXJpYWJsZSwgd2hpY2ggaXMgdGhlIGF2ZXJhZ2UgRlJMIGJ5IHNjaG9vbDoKICBtdXRhdGUoLiwKICAgICAgICAgICAgc2Nob29sX2ZybCA9IG1lYW4oZ2tmcmVlbHVuY2gsIG5hLnJtID0gVFJVRSkgLSAxCiAgICAgICAgICkgJT4lCiAgdW5ncm91cCgpICU+JQogIHNlbGVjdCguLAogICAgICAgICBzY2hvb2xpZCwKICAgICAgICAgY2xhc3NpZCwKICAgICAgICAgcmVhZCwKICAgICAgICAgY2xhc3N0eXBlLAogICAgICAgICB5ZWFyc19leHAsCiAgICAgICAgIHVyYmFuaWNpdHksCiAgICAgICAgIHN0dV9mcmwsCiAgICAgICAgIGNsYXNzX2ZybCwKICAgICAgICAgc2Nob29sX2ZybCkKCmdsaW1wc2Uoc3Rhci5jbGVhbikgIApgYGAKClRoaXMgaXMgYSBncmVhdCBleGFtcGxlIG9mIHdoZW4geW91IG1pZ2h0ICJhZ2dyZWdhdGUgdXAiIGFuIGluZGl2aWR1YWwgbWVhc3VyZSB0byBzZWUgaWYgaXQgaXMgYSBzaWduaWZpY2FudCBwcmVkaWN0b3IgYXQgbGV2ZWwgMiBvciBsZXZlbCAzLiBXZSBkaWQgdGhpcyBpbiB0aGUgcGFzdCB3aXRoIHRoZSB0YXNrIHNpZ25pZmljYW5jZSBtZWFzdXJlIChgVFNJR2ApIGFuZCB0aGUgYGxxMjAwMmAgZGF0YXNldC4KClRoZSBhYm92ZSBjb21tYW5kIGhhcyB0d28gcGFydHMuIEZpcnN0LCB3ZSB1c2UgdGhlIGBncm91cF9ieWAgY29tbWFuZCBmcm9tIG91ciBmcmllbmQsIHRoZSBgZHBseXJgIHBhY2thZ2UsIHRvIGdyb3VwIGFsbCByZXNwb25zZXMgaW4gdGhlaXIgZ3JvdXBzIGRldGVybWluZWQgYnkgIGBjbGFzc2lkYC4gVGhpcyBpcyBhIGZpcnN0IHN0ZXAgdG8gZGVjbGFyZSB0aGF0IGFsbCB0aGUgb3BlcmF0aW9ucyB0aGF0IGZvbGxvdyBhcmUgdG8gYmUgcGVyZm9ybWVkIGF0IHRoZSBjbHVzdGVyIChjbGFzc3Jvb20pIGxldmVsLiBTZWNvbmQsIHdlIHVzZSB0aGUgYG11dGF0ZWAgY29tbWFuZCwgd2hpY2ggd2UgaGF2ZSB1c2VkIGJlZm9yZSwgdG8gY3JlYXRlIG5ldyB2YXJpYWJsZXMgdGhhdCBhcmUgZXF1YWwgdG8gdGhlIG1lYW4gb2YgdGhlIG9yaWdpbmFsIHZhcmlhYmxlcy4gU2luY2Ugd2UganVzdCBncm91cGVkIGJ5IGBjbGFzc2lkYCwgdGhpcyBtZWFucyB3ZSBnZXQgdGhlIG1lYW4gb2YgdGhlIHZhcmlhYmxlcyBieSBjbGFzc3Jvb20uCgpXZSBoYXZlIHRvIGdldCBqdXN0IGEgbGl0dGxlIGZhbmNpZXIgLSBzaW5jZSB0aGUgdmFsdWVzIG9mIGBna2ZyZWVsdW5jaGAgYXJlIDEgZm9yIGZyZWUgbHVuY2ggYW5kIDIgZm9yIG5vbi1mcmVlIGx1bmNoLCBmaXJzdCB3ZSB0YWtlIHRoZSBtZWFuIG9mIHRoYXQgdmFyaWFibGUuIFRoZW4sIHdlIHN1YnRyYWN0IDEgZnJvbSB0aGF0IHZhbHVlIHRvIGdpdmUgdXMgdGhlIHByb3BvcnRpb24gb2Ygc3R1ZGVudHMgd2hvIGFyZSBub24tZnJlZSBsdW5jaC4gVGhpcyBnaXZlcyB1cyB0aGUgbmV3IHZhcmlhYmxlIGBjbGFzc19mcmxgLgoKRmluYWxseSwgd2UgZG8gdGhlIHNhbWUgdGhpbmcgYXQgdGhlIHNjaG9vbCBsZXZlbCwgZ3JvdXBpbmcgYnkgYHNjaG9vbF9pZGAsIGFuZCB0aGVuIHBlcmZvcm1pbmcgdGhlIHNhbWUgb3BlcmF0aW9uIG9uIGBna2ZyZWVsdW5jaGAgdG8gZ2l2ZSB1cyBgc2Nob29sX2ZybGAsIHdoaWNoIGlzIHRoZSBwcm9wb3J0aW9uIG9mIHN0dWRlbnRzIGF0IGVhY2ggc2Nob29sIHdobyBhcmUgbm9uLSBmcmVlIGx1bmNoLgoKT3RoZXIgdGhhbiB0aGF0LCB3ZSBhcmUganVzdCB0aWR5aW5nIGFuZCByZW5hbWluZyBmaXJzdCwgYW5kIHRoZW4gYXQgdGhlIGJvdHRvbSBzZWxlY3Rpbmcgb25seSB0aGUgdmFyaWFibGVzIHdlIG5lZWQgZm9yIHRoZSBtb2RlbHMgd2Ugd2lsbCBydW4uCgojIFBhcnQgVHdvOiBUaGUgTnVsbCAzLUxldmVsIE1vZGVsCgojIyBSdW4gYSBudWxsIDMtbGV2ZWwgbW9kZWwgZm9yIHJlYWRpbmcgc2NvcmVzCgpgYGB7cn0KCm1vZGVsLjAgPC0gbG1lcihyZWFkIH4gKDF8c2Nob29saWQpICsgKDF8Y2xhc3NpZCksIGRhdGEgPSBzdGFyLmNsZWFuKQpzdW1tYXJ5KG1vZGVsLjApCgpgYGAKCiMjIENhbGN1bGF0ZSBMMiBJQ0MgKENsYXNzcm9vbSkKYGBge3J9CgpJQ0MuY2xhc3MgPC0gMTE0LjkvKDcyOC4yICsgMTY5LjggKyAxMTQuOSkKSUNDLmNsYXNzCgpgYGAKCiMjIENhbGN1bGF0ZSBMMyBJQ0MgKFNjaG9vbCkKYGBge3J9CgpJQ0Muc2Nob29sIDwtIDE2OS44Lyg3MjguMiArIDE2OS44ICsgMTE0LjkpCklDQy5zY2hvb2wKCmBgYAoKIyBQYXJ0IFRocmVlOiBBZGQgUHJlZGljdG9ycyBhdCBFYWNoIExldmVsCgojIyBBZGQgYSBzdHVkZW50LWxldmVsIHByZWRpY3RvciwgRlJMCgpgYGB7cn0KCm1vZGVsLjEgPC0gbG1lcihyZWFkIH4gc3R1X2ZybCArICgxfHNjaG9vbGlkKSArICgxfGNsYXNzaWQpLCBkYXRhID0gc3Rhci5jbGVhbikKc3VtbWFyeShtb2RlbC4xKQoKYGBgCgojIyBBZGQgdGVhY2hlci1sZXZlbCBwcmVkaWN0b3IsIHlycy4gb2YgZXhwZXJpZW5jZSBhbmQgY2xhc3MgdHlwZQoKYGBge3J9Cgptb2RlbC4yIDwtIGxtZXIocmVhZCB+IHN0dV9mcmwgKyB5ZWFyc19leHAgKyBjbGFzc3R5cGUgKyAoMXxzY2hvb2xpZCkgKyAoMXxjbGFzc2lkKSwgZGF0YSA9IHN0YXIuY2xlYW4pCnN1bW1hcnkobW9kZWwuMikKCmBgYAoKIyMgQWRkIGEgc2Nob29sLWxldmVsIHByZWRpY3RvciwgdXJiYW5pY2l0eSBvZiBzY2hvb2wKCmBgYHtyfQoKbW9kZWwuMyA8LSBsbWVyKHJlYWQgfiBzdHVfZnJsICsgeWVhcnNfZXhwICsgY2xhc3N0eXBlICsgdXJiYW5pY2l0eSArICgxfHNjaG9vbGlkKSArICgxfGNsYXNzaWQpLCBkYXRhID0gc3Rhci5jbGVhbikKc3VtbWFyeShtb2RlbC4zKQoKYGBgCgojIyBBZGQgaW4gRlJMIGF0IGFsbCB0aHJlZSBsZXZlbHMgKHN0dWRlbnQsIHRlYWNoZXIsIHNjaG9vbCkKCmBgYHtyfQoKbW9kZWwuNCA8LSBsbWVyKHJlYWQgfiBzdHVfZnJsICsgIGNsYXNzX2ZybCArIHNjaG9vbF9mcmwgKyAgeWVhcnNfZXhwICsgY2xhc3N0eXBlICsgdXJiYW5pY2l0eSArICgxfHNjaG9vbGlkKSArICgxfGNsYXNzaWQpLCBkYXRhID0gc3Rhci5jbGVhbikKc3VtbWFyeShtb2RlbC40KQoKYGBgCgojIFBhcnQgRm91cjogVHJ5IEFkZGluZyBhIFJhbmRvbSBTbG9wZQoKIyMgQWRkIGEgcmFuZG9tIHNsb3BlIGZvciB5ZWFycyBleHBlcmllbmNlIGF0IGxldmVsIDMgKHNjaG9vbCkKYGBge3J9Cgptb2RlbC41IDwtIGxtZXIocmVhZCB+IHN0dV9mcmwgKyB5ZWFyc19leHAgKyBjbGFzc3R5cGUgKyB1cmJhbmljaXR5ICsgKDF8Y2xhc3NpZCkgKyAoeWVhcnNfZXhwfHNjaG9vbGlkKSwgZGF0YSA9IHN0YXIuY2xlYW4pCnN1bW1hcnkobW9kZWwuNSkKCmBgYAoKIyBTaG91bGQgV2UgSW5jbHVkZSBUaGF0IFJhbmRvbSBTbG9wZT8gQ29tcGFyaW5nIGBtb2RlbC40YCB3aXRoIGBtb2RlbC41YDoKYGBge3J9CmFub3ZhKG1vZGVsLjQsIG1vZGVsLjUpCmBgYAoKIyBVc2luZyB0aGUgYG1vZGVsc3VtbWFyeWAgYW5kIGBicm9vbS5taXhlZGAgUGFja2FnZXMgdG8gT3JnYW5pemUgWW91ciBSZXN1bHRzOgpgYGB7cn0KbGlicmFyeShtb2RlbHN1bW1hcnkpCmxpYnJhcnkoYnJvb20ubWl4ZWQpCgptb2RlbHMgPC0gbGlzdChtb2RlbC4wLCBtb2RlbC4xLCBtb2RlbC4yLCBtb2RlbC4zLCBtb2RlbC40LCBtb2RlbC41KQptb2RlbHN1bW1hcnkobW9kZWxzLCBvdXRwdXQgPSAiZGVmYXVsdCIpCmBgYAoKIyBIVE1MIFZlcnNpb24gVGhhdCBZb3UgQ2FuIE9wZW4gaW4gV29yZDoKYGBge3J9Cm1vZGVsc3VtbWFyeShtb2RlbHMsIG91dHB1dCA9ICdtc3VtLmh0bWwnLCB0aXRsZSA9ICdNTE0gRXN0aW1hdGVzJykKCmBgYAoK

Multilevel Modeling, Module 7: Three-Level Models

Dr. Broda

Load Some Packages to Help with the Analysis and Data Management:

Load and Explore the Data

Part One: Data Exploration and Management

Use `str` to Figure out the 2/1 Coding Scheme and Variable Labels

Create a Clean Dataset Ready for Modeling

Part Two: The Null 3-Level Model

Run a null 3-level model for reading scores

Calculate L2 ICC (Classroom)

Calculate L3 ICC (School)

Part Three: Add Predictors at Each Level

Add a student-level predictor, FRL

Add teacher-level predictor, yrs. of experience and class type

Add a school-level predictor, urbanicity of school

Add in FRL at all three levels (student, teacher, school)

Part Four: Try Adding a Random Slope

Add a random slope for years experience at level 3 (school)

Should We Include That Random Slope? Comparing `model.4` with `model.5`:

Using the `modelsummary` and `broom.mixed` Packages to Organize Your Results:

HTML Version That You Can Open in Word:

Multilevel Modeling, Module 7: Three-Level Models

Dr. Broda

Load Some Packages to Help with the Analysis and Data Management:

Load and Explore the Data

Part One: Data Exploration and Management

Use str to Figure out the 2/1 Coding Scheme and Variable Labels

Create a Clean Dataset Ready for Modeling

Part Two: The Null 3-Level Model

Run a null 3-level model for reading scores

Calculate L2 ICC (Classroom)

Calculate L3 ICC (School)

Part Three: Add Predictors at Each Level

Add a student-level predictor, FRL

Add teacher-level predictor, yrs. of experience and class type

Add a school-level predictor, urbanicity of school

Add in FRL at all three levels (student, teacher, school)

Part Four: Try Adding a Random Slope

Add a random slope for years experience at level 3 (school)

Should We Include That Random Slope? Comparing model.4 with model.5:

Using the modelsummary and broom.mixed Packages to Organize Your Results:

HTML Version That You Can Open in Word:

Use `str` to Figure out the 2/1 Coding Scheme and Variable Labels

Should We Include That Random Slope? Comparing `model.4` with `model.5`:

Using the `modelsummary` and `broom.mixed` Packages to Organize Your Results: