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

library(Hmisc)
suppressPackageStartupMessages(library(psych))
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(lme4))

Load and Explore the Data


projectCHPS <- readr::read_csv("ColonialHeightsMAPReadingScores2020.csv")


-- Column specification --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
cols(
  .default = col_double(),
  Virtual = col_character(),
  ESLStatus = col_character(),
  FedEthnicity = col_character(),
  Gender = col_character()
)
i Use `spec()` for the full column specifications.

glimpse(projectCHPS)

Rows: 1,416
Columns: 20
$ StudentNumber           <dbl> 2000, 1585, 2134, 1905, 2333, 1656, 1750, 2240, 2017, 1424, 2157, 1494, 2084, 2299, 2188, 2010, 2215, 2245, 1883, 2131, 2325, 2222, 1926, 2200, 2252, 2150, 2197, 2393, 2192, 2207, 2...
$ ReadingMAPScore         <dbl> 140, 146, 183, 182, 168, 158, 148, 156, 151, 157, 142, 154, 150, 138, 144, 142, 151, 150, 139, 144, 141, 208, 146, 178, 153, 134, 185, 140, 149, 137, 152, 164, 166, 141, 157, 197, 1...
$ Virtual                 <chr> "Virtual", "Virtual", "Virtual", "Virtual", "Virtual", "Virtual", "Traditional", "Traditional", "Traditional", "Traditional", "Traditional", "Traditional", "Traditional", "Tradition...
$ Grade                   <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
$ Grade0                  <dbl> 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, 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, 0...
$ DisabilityStatus        <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 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, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ ESLStatus               <chr> "Non-ESL", "Non-ESL", "Non-ESL", "Non-ESL", "Non-ESL", "Non-ESL", "Non-ESL", "Non-ESL", "Non-ESL", "Non-ESL", "Non-ESL", "Non-ESL", "Non-ESL", "Non-ESL", "ESL", "ESL", "ESL", "ESL",...
$ SES                     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
$ FedEthnicity            <chr> "Non-Hispanic", "Non-Hispanic", "Non-Hispanic", "Non-Hispanic", "Non-Hispanic", "Non-Hispanic", "Non-Hispanic", "Non-Hispanic", "Non-Hispanic", "Non-Hispanic", "Non-Hispanic", "Non-...
$ Gender                  <chr> "M", "M", "M", "M", "M", "F", "F", "M", "F", "M", "M", "F", "M", "M", "M", "M", "M", "M", "M", "F", "M", "M", "M", "M", "M", "M", "F", "M", "M", "M", "F", "F", "F", "F", "M", "F", "...
$ TeacherID               <dbl> 113, 173, 173, 173, 173, 176, 110, 110, 110, 122, 127, 144, 163, 170, 110, 160, 160, 160, 167, 167, 170, 102, 102, 159, 159, 176, 176, 122, 122, 122, 122, 144, 160, 167, 102, 102, 1...
$ School                  <dbl> 60, 50, 30, 30, 50, 60, 50, 50, 50, 50, 30, 60, 60, 60, 50, 50, 50, 50, 30, 30, 60, 50, 30, 50, 50, 30, 50, 50, 50, 50, 50, 60, 50, 30, 60, 60, 60, 30, 30, 60, 60, 60, 50, 60, 60, 3...
$ SchoolID                <dbl> 70, 70, 70, 70, 70, 70, 50, 50, 50, 50, 30, 60, 60, 60, 50, 50, 50, 50, 30, 30, 60, 70, 70, 70, 70, 70, 70, 50, 50, 50, 50, 60, 50, 30, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 7...
$ AmericanIndian          <dbl> 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, 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, 0...
$ Asian                   <dbl> 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, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ Black                   <dbl> 0, 0, 0, 1, 1, 0, 0, 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, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1...
$ White                   <dbl> 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0...
$ PacificIslander         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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...
$ Multirace               <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0...
$ StandardReadingMAPScore <dbl> -1.21446848, -0.83609485, 1.49720921, 1.43414694, 0.55127514, -0.07934758, -0.70997030, -0.20547213, -0.52078349, -0.14240986, -1.08834393, -0.33159667, -0.58384576, -1.34059302, -0...

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.

Create a Clean Dataset Ready for Modeling

chps.clean <- projectCHPS %>%
  mutate(.,
         
         Virtual.fac = as_factor(Virtual),
         DisabilityStatus.fac = as_factor(DisabilityStatus),
         ESLStatus.fac = as_factor(ESLStatus),
         SES.fac = as.factor(SES),
         FedEthnicity.fac = as.factor(FedEthnicity),
         Gender.fac = as.factor(Gender),
         Asian.fac = as.factor(Asian),
         Black.fac = as.factor(Black),
         White.fac = as.factor (White),
         Multirace.fac = as.factor(Multirace)
         
         ) %>%
  group_by(TeacherID) %>% # Create a new variable, which is the average SWD by classroom:
  mutate(.,
            class_swd = mean(DisabilityStatus, na.rm = TRUE),
            class_ses = mean(SES, na.rm = TRUE)
         ) %>%
  ungroup() %>%
  group_by(SchoolID) %>% # Create a new variable, which is the average SWD by school:
  mutate(.,
            school_swd = mean(DisabilityStatus, na.rm = TRUE),
            school_ses = mean(SES, na.rm = TRUE)
         ) %>%
  ungroup() %>%
  select(.,
         SchoolID,
         TeacherID,
         Grade,
         Grade0,
         ReadingMAPScore,
         StandardReadingMAPScore,
         Virtual.fac,
         DisabilityStatus.fac,
         ESLStatus.fac,
         SES.fac,
         FedEthnicity.fac,
         Gender.fac,
         Asian.fac,
         Black.fac,
         White.fac,
         Multirace.fac, 
         class_swd,
         school_swd, 
         class_ses,
         school_ses)

glimpse(chps.clean)

Rows: 1,416
Columns: 20
$ SchoolID                <dbl> 70, 70, 70, 70, 70, 70, 50, 50, 50, 50, 30, 60, 60, 60, 50, 50, 50, 50, 30, 30, 60, 70, 70, 70, 70, 70, 70, 50, 50, 50, 50, 60, 50, 30, 70, 70, 70, 70, 70, 70, 70, 70,...
$ TeacherID               <dbl> 113, 173, 173, 173, 173, 176, 110, 110, 110, 122, 127, 144, 163, 170, 110, 160, 160, 160, 167, 167, 170, 102, 102, 159, 159, 176, 176, 122, 122, 122, 122, 144, 160, 16...
$ Grade                   <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
$ Grade0                  <dbl> 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, 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,...
$ ReadingMAPScore         <dbl> 140, 146, 183, 182, 168, 158, 148, 156, 151, 157, 142, 154, 150, 138, 144, 142, 151, 150, 139, 144, 141, 208, 146, 178, 153, 134, 185, 140, 149, 137, 152, 164, 166, 14...
$ StandardReadingMAPScore <dbl> -1.21446848, -0.83609485, 1.49720921, 1.43414694, 0.55127514, -0.07934758, -0.70997030, -0.20547213, -0.52078349, -0.14240986, -1.08834393, -0.33159667, -0.58384576, -...
$ Virtual.fac             <fct> Virtual, Virtual, Virtual, Virtual, Virtual, Virtual, Traditional, Traditional, Traditional, Traditional, Traditional, Traditional, Traditional, Traditional, Tradition...
$ DisabilityStatus.fac    <fct> 1, 1, 1, 1, 1, 1, 1, 1, 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, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ ESLStatus.fac           <fct> Non-ESL, Non-ESL, Non-ESL, Non-ESL, Non-ESL, Non-ESL, Non-ESL, Non-ESL, Non-ESL, Non-ESL, Non-ESL, Non-ESL, Non-ESL, Non-ESL, ESL, ESL, ESL, ESL, ESL, ESL, ESL, Non-ES...
$ SES.fac                 <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
$ FedEthnicity.fac        <fct> Non-Hispanic, Non-Hispanic, Non-Hispanic, Non-Hispanic, Non-Hispanic, Non-Hispanic, Non-Hispanic, Non-Hispanic, Non-Hispanic, Non-Hispanic, Non-Hispanic, Non-Hispanic,...
$ Gender.fac              <fct> M, M, M, M, M, F, F, M, F, M, M, F, M, M, M, M, M, M, M, F, M, M, M, M, M, M, F, M, M, M, F, F, F, F, M, F, F, F, F, M, F, M, F, M, M, M, F, M, F, F, F, M, F, F, M, F,...
$ Asian.fac               <fct> 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, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ Black.fac               <fct> 0, 0, 0, 1, 1, 0, 0, 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, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0,...
$ White.fac               <fct> 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1,...
$ Multirace.fac           <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0,...
$ class_swd               <dbl> 0.16666667, 0.26315789, 0.26315789, 0.26315789, 0.26315789, 0.09090909, 0.15000000, 0.15000000, 0.15000000, 0.09090909, 0.16666667, 0.23076923, 0.36363636, 0.08333333,...
$ school_swd              <dbl> 0.1518325, 0.1518325, 0.1518325, 0.1518325, 0.1518325, 0.1518325, 0.1736527, 0.1736527, 0.1736527, 0.1736527, 0.1278195, 0.1639344, 0.1639344, 0.1639344, 0.1736527, 0....
$ class_ses               <dbl> 0.6666667, 0.7894737, 0.7894737, 0.7894737, 0.7894737, 0.9090909, 0.7000000, 0.7000000, 0.7000000, 0.7272727, 0.5000000, 0.4615385, 0.3636364, 0.6666667, 0.7000000, 0....
$ school_ses              <dbl> 0.7382199, 0.7382199, 0.7382199, 0.7382199, 0.7382199, 0.7382199, 0.6107784, 0.6107784, 0.6107784, 0.6107784, 0.7368421, 0.4508197, 0.4508197, 0.4508197, 0.6107784, 0....

describe(chps.clean)

Part Two: The Null 3-Level Model

Run a null 3-level model for reading scores


model.null <- lmer(StandardReadingMAPScore ~ (1|SchoolID) + (1|TeacherID), REML = FALSE, data = chps.clean)
summary(model.null)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: StandardReadingMAPScore ~ (1 | SchoolID) + (1 | TeacherID)
   Data: chps.clean

     AIC      BIC   logLik deviance df.resid 
  3735.8   3756.8  -1863.9   3727.8     1412 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5353 -0.5848  0.0756  0.6573  4.8124 

Random effects:
 Groups    Name        Variance Std.Dev.
 TeacherID (Intercept) 0.310827 0.55752 
 SchoolID  (Intercept) 0.003911 0.06254 
 Residual              0.730001 0.85440 
Number of obs: 1416, groups:  TeacherID, 75; SchoolID, 6

Fixed effects:
            Estimate Std. Error t value
(Intercept) -0.02857    0.07525   -0.38

Calculate L2 ICC (Classroom)


ICC.class <- .310827/(0.310827+0.003911+0.730001)
ICC.class

[1] 0.2975164

Calculate L3 ICC (School)


ICC.school <- 0.003911/(0.310827+0.003911+0.730001)
ICC.school

[1] 0.003743519


model.0 <- lmer(StandardReadingMAPScore ~  (1|TeacherID), REML = FALSE, data = chps.clean)
summary(model.0)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: StandardReadingMAPScore ~ (1 | TeacherID)
   Data: chps.clean

     AIC      BIC   logLik deviance df.resid 
  3734.7   3750.5  -1864.4   3728.7     1413 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5860 -0.5926  0.0901  0.6594  4.8119 

Random effects:
 Groups    Name        Variance Std.Dev.
 TeacherID (Intercept) 0.3135   0.5599  
 Residual              0.7311   0.8550  
Number of obs: 1416, groups:  TeacherID, 75

Fixed effects:
            Estimate Std. Error t value
(Intercept) -0.02764    0.06986  -0.396

Calculate L2 ICC (Classroom)


ICC.class <- .3135/(.3135+.7311)
ICC.class

[1] 0.3001149

lmerTest::rand(model.0)

ANOVA-like table for random-effects: Single term deletions

Model:
StandardReadingMAPScore ~ (1 | TeacherID)
                npar  logLik    AIC    LRT Df Pr(>Chisq)    
<none>             3 -1864.4 3734.7                         
(1 | TeacherID)    2 -2009.2 4022.4 289.72  1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Part Three: Remove School Level, Add a student-level predictor, Disability Status


model.1 <- lmer(StandardReadingMAPScore ~ DisabilityStatus.fac + (1|TeacherID), REML= FALSE, data = chps.clean)
summary(model.1)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: StandardReadingMAPScore ~ DisabilityStatus.fac + (1 | TeacherID)
   Data: chps.clean

     AIC      BIC   logLik deviance df.resid 
  3702.1   3723.1  -1847.0   3694.1     1412 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.7585 -0.5853  0.0671  0.6601  4.7264 

Random effects:
 Groups    Name        Variance Std.Dev.
 TeacherID (Intercept) 0.2442   0.4942  
 Residual              0.7209   0.8491  
Number of obs: 1416, groups:  TeacherID, 75

Fixed effects:
                      Estimate Std. Error t value
(Intercept)            0.04325    0.06380   0.678
DisabilityStatus.fac1 -0.42964    0.07124  -6.031

Correlation of Fixed Effects:
            (Intr)
DsbltyStt.1 -0.177

Add student-level predictor, SES status

model.2 <- lmer(StandardReadingMAPScore ~ DisabilityStatus.fac + SES.fac + (1|TeacherID), REML= FALSE, data = chps.clean)
summary(model.2)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: StandardReadingMAPScore ~ DisabilityStatus.fac + SES.fac + (1 |      TeacherID)
   Data: chps.clean

     AIC      BIC   logLik deviance df.resid 
  3690.4   3716.6  -1840.2   3680.4     1411 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9232 -0.5859  0.0671  0.6446  4.7818 

Random effects:
 Groups    Name        Variance Std.Dev.
 TeacherID (Intercept) 0.2404   0.4903  
 Residual              0.7142   0.8451  
Number of obs: 1416, groups:  TeacherID, 75

Fixed effects:
                      Estimate Std. Error t value
(Intercept)            0.15383    0.06997   2.198
DisabilityStatus.fac1 -0.42488    0.07091  -5.992
SES.fac1              -0.18147    0.04887  -3.714

Correlation of Fixed Effects:
            (Intr) DsbS.1
DsbltyStt.1 -0.152       
SES.fac1    -0.425 -0.020

Add a student-level predictor, black


model.3 <- lmer(StandardReadingMAPScore ~ DisabilityStatus.fac + SES.fac + Black.fac + (1|TeacherID), REML= FALSE, data = chps.clean)
summary(model.3)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: StandardReadingMAPScore ~ DisabilityStatus.fac + SES.fac + Black.fac +      (1 | TeacherID)
   Data: chps.clean

     AIC      BIC   logLik deviance df.resid 
  3677.3   3708.9  -1832.7   3665.3     1410 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0025 -0.5925  0.0752  0.6290  4.7005 

Random effects:
 Groups    Name        Variance Std.Dev.
 TeacherID (Intercept) 0.2462   0.4961  
 Residual              0.7055   0.8400  
Number of obs: 1416, groups:  TeacherID, 75

Fixed effects:
                      Estimate Std. Error t value
(Intercept)            0.18639    0.07091   2.629
DisabilityStatus.fac1 -0.41316    0.07058  -5.854
SES.fac1              -0.16345    0.04880  -3.349
Black.fac1            -0.22344    0.05741  -3.892

Correlation of Fixed Effects:
            (Intr) DsbS.1 SES.f1
DsbltyStt.1 -0.145              
SES.fac1    -0.404 -0.016       
Black.fac1  -0.121 -0.034 -0.094

Add classroom-level predictor, class_swd (calculated cluster mean)


model.4 <- lmer(StandardReadingMAPScore ~ DisabilityStatus.fac + SES.fac + Black.fac + class_swd + (1|TeacherID), REML= FALSE, data = chps.clean)
summary(model.4)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: StandardReadingMAPScore ~ DisabilityStatus.fac + SES.fac + Black.fac +      class_swd + (1 | TeacherID)
   Data: chps.clean

     AIC      BIC   logLik deviance df.resid 
  3649.6   3686.3  -1817.8   3635.6     1409 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9527 -0.6091  0.0728  0.6340  4.7615 

Random effects:
 Groups    Name        Variance Std.Dev.
 TeacherID (Intercept) 0.1625   0.4031  
 Residual              0.7030   0.8384  
Number of obs: 1416, groups:  TeacherID, 75

Fixed effects:
                      Estimate Std. Error t value
(Intercept)            0.40586    0.07203   5.634
DisabilityStatus.fac1 -0.32044    0.07275  -4.405
SES.fac1              -0.16332    0.04856  -3.363
Black.fac1            -0.21538    0.05707  -3.774
class_swd             -1.46810    0.25289  -5.805

Correlation of Fixed Effects:
            (Intr) DsbS.1 SES.f1 Blck.1
DsbltyStt.1  0.010                     
SES.fac1    -0.392 -0.014              
Black.fac1  -0.111 -0.030 -0.098       
class_swd   -0.501 -0.287 -0.006 -0.013

Tested a random slope for class_swd


model.5 <- lmer(StandardReadingMAPScore ~ DisabilityStatus.fac + SES.fac + Black.fac + class_swd + (class_swd|TeacherID), REML= FALSE, data = chps.clean)
summary(model.5)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: StandardReadingMAPScore ~ DisabilityStatus.fac + SES.fac + Black.fac +      class_swd + (class_swd | TeacherID)
   Data: chps.clean

     AIC      BIC   logLik deviance df.resid 
  3648.4   3695.7  -1815.2   3630.4     1407 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9636 -0.6019  0.0652  0.6340  4.7664 

Random effects:
 Groups    Name        Variance Std.Dev. Corr 
 TeacherID (Intercept) 0.1552   0.3939        
           class_swd   0.9896   0.9948   -0.49
 Residual              0.7014   0.8375        
Number of obs: 1416, groups:  TeacherID, 75

Fixed effects:
                      Estimate Std. Error t value
(Intercept)            0.37723    0.07144   5.280
DisabilityStatus.fac1 -0.32034    0.07266  -4.409
SES.fac1              -0.16392    0.04847  -3.382
Black.fac1            -0.21742    0.05697  -3.816
class_swd             -1.20224    0.31049  -3.872

Correlation of Fixed Effects:
            (Intr) DsbS.1 SES.f1 Blck.1
DsbltyStt.1  0.010                     
SES.fac1    -0.393 -0.014              
Black.fac1  -0.107 -0.030 -0.099       
class_swd   -0.530 -0.233 -0.010 -0.022

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

library(modelsummary)
library(broom.mixed)

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

	Model 1	Model 2	Model 3	Model 4	Model 5	Model 6	Model 7
(Intercept)	-0.029	-0.028	0.043	0.154	0.186	0.406	0.377
	(0.075)	(0.070)	(0.064)	(0.070)	(0.071)	(0.072)	(0.071)
TeacherID sd__(Intercept)	0.558
SchoolID sd__(Intercept)	0.063
Residual sd__Observation	0.854
sd__(Intercept)		0.560	0.494	0.490	0.496	0.403	0.394
sd__Observation		0.855	0.849	0.845	0.840	0.838	0.838
DisabilityStatus.fac1			-0.430	-0.425	-0.413	-0.320	-0.320
			(0.071)	(0.071)	(0.071)	(0.073)	(0.073)
SES.fac1				-0.181	-0.163	-0.163	-0.164
				(0.049)	(0.049)	(0.049)	(0.048)
Black.fac1					-0.223	-0.215	-0.217
					(0.057)	(0.057)	(0.057)
class_swd						-1.468	-1.202
						(0.253)	(0.310)
cor__(Intercept).class_swd							-0.489
sd__class_swd							0.995
AIC	3735.8	3734.7	3702.1	3690.4	3677.3	3649.6	3648.4
BIC	3756.8	3750.5	3723.1	3716.6	3708.9	3686.3	3695.7
Log.Lik.	-1863.895	-1864.358	-1847.047	-1840.186	-1832.665	-1817.779	-1815.195

HTML Version That You Can Open in Word:

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

LS0tDQp0aXRsZTogJ011bHRpbGV2ZWwgTW9kZWxpbmcsIEZpbmFsIFByZXNlbnRhdGlvbicNCmF1dGhvcjogJ0pha2UgUmV5bm9sZHMgLSBOb3ZlbWJlciAyNSwgMjAyMCcNCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQoNCiMgTG9hZCBTb21lIFBhY2thZ2VzIHRvIEhlbHAgd2l0aCB0aGUgQW5hbHlzaXMgYW5kIERhdGEgTWFuYWdlbWVudDoNCmBgYHtyfQ0KbGlicmFyeShIbWlzYykNCnN1cHByZXNzUGFja2FnZVN0YXJ0dXBNZXNzYWdlcyhsaWJyYXJ5KHBzeWNoKSkNCnN1cHByZXNzUGFja2FnZVN0YXJ0dXBNZXNzYWdlcyhsaWJyYXJ5KHRpZHl2ZXJzZSkpDQpzdXBwcmVzc1BhY2thZ2VTdGFydHVwTWVzc2FnZXMobGlicmFyeShsbWU0KSkNCmBgYA0KDQojIExvYWQgYW5kIEV4cGxvcmUgdGhlIERhdGENCmBgYHtyfQ0KDQpwcm9qZWN0Q0hQUyA8LSByZWFkcjo6cmVhZF9jc3YoIkNvbG9uaWFsSGVpZ2h0c01BUFJlYWRpbmdTY29yZXMyMDIwLmNzdiIpDQpnbGltcHNlKHByb2plY3RDSFBTKQ0KDQpgYGANCg0KIyBQYXJ0IE9uZTogRGF0YSBFeHBsb3JhdGlvbiBhbmQgTWFuYWdlbWVudA0KDQojIyBVc2UgYHN0cmAgdG8gRmlndXJlIG91dCB0aGUgMi8xIENvZGluZyBTY2hlbWUgYW5kIFZhcmlhYmxlIExhYmVscw0KVGhlIGBzdHJgIGZ1bmN0aW9uIChzaG9ydCBmb3Igc3RydWN0dXJlKSB3aWxsIHNob3cgdXMgYW55IHZhcmlhYmxlcyBsYWJlbHMgb3IgdmFsdWVzIHRoYXQgaGF2ZSBiZWVuIGF0dGFjaGVkIHRvIGEgdmFyaWFibGUuIFlvdSBjYW4gdXNlIHRoaXMgdG8gZmluZCBvdXIgaG93IHZhcmlhYmxlcyBoYXZlIGJlZW4gY29kZWQuDQoNCiMjIENyZWF0ZSBhIENsZWFuIERhdGFzZXQgUmVhZHkgZm9yIE1vZGVsaW5nDQpgYGB7cn0NCmNocHMuY2xlYW4gPC0gcHJvamVjdENIUFMgJT4lDQogIG11dGF0ZSguLA0KICAgICAgICAgDQogICAgICAgICBWaXJ0dWFsLmZhYyA9IGFzX2ZhY3RvcihWaXJ0dWFsKSwNCiAgICAgICAgIERpc2FiaWxpdHlTdGF0dXMuZmFjID0gYXNfZmFjdG9yKERpc2FiaWxpdHlTdGF0dXMpLA0KICAgICAgICAgRVNMU3RhdHVzLmZhYyA9IGFzX2ZhY3RvcihFU0xTdGF0dXMpLA0KICAgICAgICAgU0VTLmZhYyA9IGFzLmZhY3RvcihTRVMpLA0KICAgICAgICAgRmVkRXRobmljaXR5LmZhYyA9IGFzLmZhY3RvcihGZWRFdGhuaWNpdHkpLA0KICAgICAgICAgR2VuZGVyLmZhYyA9IGFzLmZhY3RvcihHZW5kZXIpLA0KICAgICAgICAgQXNpYW4uZmFjID0gYXMuZmFjdG9yKEFzaWFuKSwNCiAgICAgICAgIEJsYWNrLmZhYyA9IGFzLmZhY3RvcihCbGFjayksDQogICAgICAgICBXaGl0ZS5mYWMgPSBhcy5mYWN0b3IgKFdoaXRlKSwNCiAgICAgICAgIE11bHRpcmFjZS5mYWMgPSBhcy5mYWN0b3IoTXVsdGlyYWNlKQ0KICAgICAgICAgDQogICAgICAgICApICU+JQ0KICBncm91cF9ieShUZWFjaGVySUQpICU+JSAjIENyZWF0ZSBhIG5ldyB2YXJpYWJsZSwgd2hpY2ggaXMgdGhlIGF2ZXJhZ2UgU1dEIGJ5IGNsYXNzcm9vbToNCiAgbXV0YXRlKC4sDQogICAgICAgICAgICBjbGFzc19zd2QgPSBtZWFuKERpc2FiaWxpdHlTdGF0dXMsIG5hLnJtID0gVFJVRSksDQogICAgICAgICAgICBjbGFzc19zZXMgPSBtZWFuKFNFUywgbmEucm0gPSBUUlVFKQ0KICAgICAgICAgKSAlPiUNCiAgdW5ncm91cCgpICU+JQ0KICBncm91cF9ieShTY2hvb2xJRCkgJT4lICMgQ3JlYXRlIGEgbmV3IHZhcmlhYmxlLCB3aGljaCBpcyB0aGUgYXZlcmFnZSBTV0QgYnkgc2Nob29sOg0KICBtdXRhdGUoLiwNCiAgICAgICAgICAgIHNjaG9vbF9zd2QgPSBtZWFuKERpc2FiaWxpdHlTdGF0dXMsIG5hLnJtID0gVFJVRSksDQogICAgICAgICAgICBzY2hvb2xfc2VzID0gbWVhbihTRVMsIG5hLnJtID0gVFJVRSkNCiAgICAgICAgICkgJT4lDQogIHVuZ3JvdXAoKSAlPiUNCiAgc2VsZWN0KC4sDQogICAgICAgICBTY2hvb2xJRCwNCiAgICAgICAgIFRlYWNoZXJJRCwNCiAgICAgICAgIEdyYWRlLA0KICAgICAgICAgR3JhZGUwLA0KICAgICAgICAgUmVhZGluZ01BUFNjb3JlLA0KICAgICAgICAgU3RhbmRhcmRSZWFkaW5nTUFQU2NvcmUsDQogICAgICAgICBWaXJ0dWFsLmZhYywNCiAgICAgICAgIERpc2FiaWxpdHlTdGF0dXMuZmFjLA0KICAgICAgICAgRVNMU3RhdHVzLmZhYywNCiAgICAgICAgIFNFUy5mYWMsDQogICAgICAgICBGZWRFdGhuaWNpdHkuZmFjLA0KICAgICAgICAgR2VuZGVyLmZhYywNCiAgICAgICAgIEFzaWFuLmZhYywNCiAgICAgICAgIEJsYWNrLmZhYywNCiAgICAgICAgIFdoaXRlLmZhYywNCiAgICAgICAgIE11bHRpcmFjZS5mYWMsIA0KICAgICAgICAgY2xhc3Nfc3dkLA0KICAgICAgICAgc2Nob29sX3N3ZCwgDQogICAgICAgICBjbGFzc19zZXMsDQogICAgICAgICBzY2hvb2xfc2VzKQ0KDQpnbGltcHNlKGNocHMuY2xlYW4pICANCmBgYA0KYGBge3J9DQpkZXNjcmliZShjaHBzLmNsZWFuKQ0KYGBgDQoNCiMgUGFydCBUd286IFRoZSBOdWxsIDMtTGV2ZWwgTW9kZWwNCg0KIyMgUnVuIGEgbnVsbCAzLWxldmVsIG1vZGVsIGZvciByZWFkaW5nIHNjb3Jlcw0KDQpgYGB7cn0NCg0KbW9kZWwubnVsbCA8LSBsbWVyKFN0YW5kYXJkUmVhZGluZ01BUFNjb3JlIH4gKDF8U2Nob29sSUQpICsgKDF8VGVhY2hlcklEKSwgUkVNTCA9IEZBTFNFLCBkYXRhID0gY2hwcy5jbGVhbikNCnN1bW1hcnkobW9kZWwubnVsbCkNCg0KYGBgDQoNCiMjIENhbGN1bGF0ZSBMMiBJQ0MgKENsYXNzcm9vbSkNCmBgYHtyfQ0KDQpJQ0MuY2xhc3MgPC0gLjMxMDgyNy8oMC4zMTA4MjcrMC4wMDM5MTErMC43MzAwMDEpDQpJQ0MuY2xhc3MNCg0KYGBgDQoNCiMjIENhbGN1bGF0ZSBMMyBJQ0MgKFNjaG9vbCkNCmBgYHtyfQ0KDQpJQ0Muc2Nob29sIDwtIDAuMDAzOTExLygwLjMxMDgyNyswLjAwMzkxMSswLjczMDAwMSkNCklDQy5zY2hvb2wNCmBgYA0KDQpgYGB7cn0NCg0KbW9kZWwuMCA8LSBsbWVyKFN0YW5kYXJkUmVhZGluZ01BUFNjb3JlIH4gICgxfFRlYWNoZXJJRCksIFJFTUwgPSBGQUxTRSwgZGF0YSA9IGNocHMuY2xlYW4pDQpzdW1tYXJ5KG1vZGVsLjApDQoNCmBgYA0KDQojIyBDYWxjdWxhdGUgTDIgSUNDIChDbGFzc3Jvb20pDQpgYGB7cn0NCg0KSUNDLmNsYXNzIDwtIC4zMTM1LyguMzEzNSsuNzMxMSkNCklDQy5jbGFzcw0KDQpgYGANCmBgYHtyfQ0KbG1lclRlc3Q6OnJhbmQobW9kZWwuMCApDQpgYGANCg0KDQoNCiMgUGFydCBUaHJlZTogUmVtb3ZlIFNjaG9vbCBMZXZlbCwgQWRkIGEgc3R1ZGVudC1sZXZlbCBwcmVkaWN0b3IsIERpc2FiaWxpdHkgU3RhdHVzIA0KYGBge3J9DQoNCm1vZGVsLjEgPC0gbG1lcihTdGFuZGFyZFJlYWRpbmdNQVBTY29yZSB+IERpc2FiaWxpdHlTdGF0dXMuZmFjICsgKDF8VGVhY2hlcklEKSwgUkVNTD0gRkFMU0UsIGRhdGEgPSBjaHBzLmNsZWFuKQ0Kc3VtbWFyeShtb2RlbC4xKQ0KDQpgYGANCg0KIyMgQWRkIHN0dWRlbnQtbGV2ZWwgcHJlZGljdG9yLCBTRVMgc3RhdHVzIA0KDQpgYGB7cn0NCm1vZGVsLjIgPC0gbG1lcihTdGFuZGFyZFJlYWRpbmdNQVBTY29yZSB+IERpc2FiaWxpdHlTdGF0dXMuZmFjICsgU0VTLmZhYyArICgxfFRlYWNoZXJJRCksIFJFTUw9IEZBTFNFLCBkYXRhID0gY2hwcy5jbGVhbikNCnN1bW1hcnkobW9kZWwuMikNCg0KYGBgDQoNCiMjIEFkZCBhIHN0dWRlbnQtbGV2ZWwgcHJlZGljdG9yLCBibGFjayANCg0KYGBge3J9DQoNCm1vZGVsLjMgPC0gbG1lcihTdGFuZGFyZFJlYWRpbmdNQVBTY29yZSB+IERpc2FiaWxpdHlTdGF0dXMuZmFjICsgU0VTLmZhYyArIEJsYWNrLmZhYyArICgxfFRlYWNoZXJJRCksIFJFTUw9IEZBTFNFLCBkYXRhID0gY2hwcy5jbGVhbikNCnN1bW1hcnkobW9kZWwuMykNCg0KYGBgDQoNCiMjIEFkZCBjbGFzc3Jvb20tbGV2ZWwgcHJlZGljdG9yLCBjbGFzc19zd2QgKGNhbGN1bGF0ZWQgY2x1c3RlciBtZWFuKQ0KDQpgYGB7cn0NCg0KbW9kZWwuNCA8LSBsbWVyKFN0YW5kYXJkUmVhZGluZ01BUFNjb3JlIH4gRGlzYWJpbGl0eVN0YXR1cy5mYWMgKyBTRVMuZmFjICsgQmxhY2suZmFjICsgY2xhc3Nfc3dkICsgKDF8VGVhY2hlcklEKSwgUkVNTD0gRkFMU0UsIGRhdGEgPSBjaHBzLmNsZWFuKQ0Kc3VtbWFyeShtb2RlbC40KQ0KDQpgYGANCg0KDQojIyBUZXN0ZWQgYSByYW5kb20gc2xvcGUgZm9yIGNsYXNzX3N3ZA0KYGBge3J9DQoNCm1vZGVsLjUgPC0gbG1lcihTdGFuZGFyZFJlYWRpbmdNQVBTY29yZSB+IERpc2FiaWxpdHlTdGF0dXMuZmFjICsgU0VTLmZhYyArIEJsYWNrLmZhYyArIGNsYXNzX3N3ZCArIChjbGFzc19zd2R8VGVhY2hlcklEKSwgUkVNTD0gRkFMU0UsIGRhdGEgPSBjaHBzLmNsZWFuKQ0Kc3VtbWFyeShtb2RlbC41KQ0KDQpgYGANCg0KIyBVc2luZyB0aGUgYG1vZGVsc3VtbWFyeWAgYW5kIGBicm9vbS5taXhlZGAgUGFja2FnZXMgdG8gT3JnYW5pemUgWW91ciBSZXN1bHRzOg0KYGBge3J9DQpsaWJyYXJ5KG1vZGVsc3VtbWFyeSkNCmxpYnJhcnkoYnJvb20ubWl4ZWQpDQoNCm1vZGVscyA8LSBsaXN0KG1vZGVsLm51bGwsIG1vZGVsLjAsIG1vZGVsLjEsIG1vZGVsLjIsIG1vZGVsLjMsIG1vZGVsLjQsIG1vZGVsLjUpDQptb2RlbHN1bW1hcnkobW9kZWxzLCBvdXRwdXQgPSAiZGVmYXVsdCIpDQpgYGANCg0KIyBIVE1MIFZlcnNpb24gVGhhdCBZb3UgQ2FuIE9wZW4gaW4gV29yZDoNCmBgYHtyfQ0KbW9kZWxzdW1tYXJ5KG1vZGVscywgb3V0cHV0ID0gJ21zdW0uaHRtbCcsIHRpdGxlID0gJ01MTSBFc3RpbWF0ZXMnKQ0KDQpgYGANCg0K

Multilevel Modeling, Final Presentation

Jake Reynolds - November 25, 2020

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)

Calculate L2 ICC (Classroom)

Part Three: Remove School Level, Add a student-level predictor, Disability Status

Add student-level predictor, SES status

Add a student-level predictor, black

Add classroom-level predictor, class_swd (calculated cluster mean)

Tested a random slope for class_swd

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

HTML Version That You Can Open in Word:

Multilevel Modeling, Final Presentation

Jake Reynolds - November 25, 2020

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)

Calculate L2 ICC (Classroom)

Part Three: Remove School Level, Add a student-level predictor, Disability Status

Add student-level predictor, SES status

Add a student-level predictor, black

Add classroom-level predictor, class_swd (calculated cluster mean)

Tested a random slope for class_swd

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

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