Multilevel Modeling, Module 8, Part 1: Intro to Growth Models

This week, we go to a dataset collected as part of the US Sustaining Effects Study (Carter 1984), which attemtped to assess the sustaining impacts of compensatory schooling for students over time.

Load in Our MVP Packages

suppressPackageStartupMessages(library(tidyverse))

package ‘ggplot2’ was built under R version 3.6.2package ‘tibble’ was built under R version 3.6.2package ‘tidyr’ was built under R version 3.6.2package ‘purrr’ was built under R version 3.6.2package ‘dplyr’ was built under R version 3.6.2

library(lme4)

package ‘lme4’ was built under R version 3.6.2Loading required package: Matrix

Attaching package: ‘Matrix’

The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack

Load in the Data

egmerged <- haven::read_dta("egmerged.dta")

glimpse(egmerged)

Rows: 7,230
Columns: 12
$ schoolid <dbl> 3430, 3610, 2180, 3040, 3221, 3360, 4450, 4390, 3610, 3610, 3610, 3450,…
$ childid  <dbl> 289376791, 288278731, 292789461, 275898121, 300246721, 244834951, 29685…
$ year     <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,…
$ grade    <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,…
$ math     <dbl> -3.887, -4.055, -3.525, -3.592, -3.525, -3.068, -2.596, -3.721, -3.232,…
$ retained <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,…
$ female   <dbl> 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1,…
$ black    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1,…
$ hispanic <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,…
$ size     <dbl> 641, 1062, 777, 724, 185, 526, 710, 532, 1062, 1062, 1062, 1052, 576, 6…
$ lowinc   <dbl> 84.8, 97.8, 96.6, 45.3, 94.6, 100.0, 66.2, 66.4, 97.8, 97.8, 97.8, 75.7…
$ mobility <dbl> 38.2, 26.8, 44.4, 18.7, 47.4, 55.4, 25.4, 20.9, 26.8, 26.8, 26.8, 62.3,…

Explore the data and ID variables using describe from the Hmisc package…

Hmisc::describe(egmerged)

egmerged 

 12  Variables      7230  Observations
-------------------------------------------------------------------------------------------
schoolid  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10      .25      .50 
    7230        0       60    0.999     3310    778.3     2330     2340     2820     3280 
     .75      .90      .95 
    3850     4280     4390 

lowest : 2020 2040 2180 2330 2340, highest: 4380 4390 4410 4440 4450
-------------------------------------------------------------------------------------------
childid  Format:%12.0g 
        n   missing  distinct      Info      Mean       Gmd       .05       .10       .25 
     7230         0      1721         1 290098439  23285178 227501781 261741391 286008442 
      .50       .75       .90       .95 
297153001 303767301 308611751 310002702 

lowest : 101480302 173559292 174743401 174755092 179986251
highest: 316843741 317345681 317654202 839559420 839578020
                                                                                          
Value      1.0e+08 1.7e+08 1.8e+08 1.9e+08 2.0e+08 2.1e+08 2.2e+08 2.3e+08 2.4e+08 2.5e+08
Frequency        4      11      40      19      38      79      66     162      65     175
Proportion   0.001   0.002   0.006   0.003   0.005   0.011   0.009   0.022   0.009   0.024
                                                                          
Value      2.6e+08 2.7e+08 2.8e+08 2.9e+08 3.0e+08 3.1e+08 3.2e+08 8.4e+08
Frequency      141     310     527    1735    2290    1507      54       7
Proportion   0.020   0.043   0.073   0.240   0.317   0.208   0.007   0.001

For the frequency table, variable is rounded to the nearest 10000000
-------------------------------------------------------------------------------------------
year  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd 
    7230        0        6    0.961     3.88    1.576 

lowest : 1 2 3 4 5, highest: 2 3 4 5 6
                                              
Value          1     2     3     4     5     6
Frequency    131  1346  1520  1672  1387  1174
Proportion 0.018 0.186 0.210 0.231 0.192 0.162
-------------------------------------------------------------------------------------------
grade  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd 
    7230        0        6    0.957    1.812    1.523 

lowest : 0 1 2 3 4, highest: 1 2 3 4 5
                                              
Value          0     1     2     3     4     5
Frequency   1582  1612  1653  1356  1018     9
Proportion 0.219 0.223 0.229 0.188 0.141 0.001
-------------------------------------------------------------------------------------------
math  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10      .25      .50 
    7230        0      523        1  -0.5369    1.743   -3.068   -2.461   -1.631   -0.619 
     .75      .90      .95 
   0.551    1.514    2.034 

lowest : -5.219 -4.591 -4.432 -4.406 -4.258, highest:  4.051  4.319  4.627  4.791  5.766
-------------------------------------------------------------------------------------------
retained  Format:%12.0g 
       n  missing distinct     Info      Sum     Mean      Gmd 
    7230        0        2    0.145      368   0.0509  0.09663 

-------------------------------------------------------------------------------------------
female  Format:%12.0g 
       n  missing distinct     Info      Sum     Mean      Gmd 
    7230        0        2     0.75     3685   0.5097   0.4999 

-------------------------------------------------------------------------------------------
black  Format:%12.0g 
       n  missing distinct     Info      Sum     Mean      Gmd 
    7230        0        2    0.644     4977   0.6884   0.4291 

-------------------------------------------------------------------------------------------
hispanic  Format:%12.0g 
       n  missing distinct     Info      Sum     Mean      Gmd 
    7230        0        2    0.365     1024   0.1416   0.2432 

-------------------------------------------------------------------------------------------
size  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10      .25      .50 
    7230        0       59    0.999    753.4    352.3      339      380      511      724 
     .75      .90      .95 
     944     1203     1424 

lowest :  113  127  185  240  247, highest: 1203 1222 1302 1424 1486
-------------------------------------------------------------------------------------------
lowinc  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10      .25      .50 
    7230        0       50    0.997    77.79     27.6     18.7     36.9     59.3     92.6 
     .75      .90      .95 
    98.1    100.0    100.0 

lowest :   0.0  17.2  18.7  26.7  29.6, highest:  98.7  99.3  99.6  99.7 100.0
-------------------------------------------------------------------------------------------
mobility  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10      .25      .50 
    7230        0       55    0.999    34.17 "print"     10.5     16.4     25.4     31.7 
     .75      .90      .95 
    43.7     52.4     55.4 

lowest :  8.8 10.5 12.5 14.3 16.4, highest: 54.4 55.4 56.4 62.3 67.0
-------------------------------------------------------------------------------------------

A Little Recoding - Create Factor Versions of Categorical Variables

A little housekeeping - our model is easier to interpret if time starts at 0…

egmerged.clean <- egmerged %>%
  mutate(.,
         retained.fac = as_factor(retained),
         female.fac = as_factor(female),
         black.fac = as_factor(black),
         hispanic.fac = as_factor(hispanic),
         year0 = year - 1)

Use arrange function to examine variables that vary by student, year, and school

egmerged %>%
  arrange(., childid, year)

Summarize math score by year, to look for preliminary changes over time

egmerged %>%
  group_by(year) %>%
  summarize(avg_math_score = mean(math, na.rm = TRUE)
            )

`summarise()` ungrouping output (override with `.groups` argument)

Same thing, by grade: only difference should be missing/ retained students

egmerged %>%
  group_by(grade) %>%
  summarize(avg_math_score = mean(math, na.rm = TRUE)
            )

`summarise()` ungrouping output (override with `.groups` argument)

Create Summary Data Frame of Student-Level Variables

student.level <-egmerged.clean %>%
  group_by(childid) %>%
  summarize(.,
            black = mean(black),
            hispanic = mean(hispanic),
            female = mean(female),
            retained = max(retained),
            math = mean(math, na.rm = TRUE))

`summarise()` ungrouping output (override with `.groups` argument)

glimpse(student.level)

Rows: 1,721
Columns: 6
$ childid  <dbl> 101480302, 173559292, 174743401, 174755092, 179986251, 179989511, 17999…
$ black    <dbl> 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0,…
$ hispanic <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,…
$ female   <dbl> 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0,…
$ retained <dbl> 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1,…
$ math     <dbl> -0.4517500, 1.1223333, -0.8530000, 0.5643333, -2.5121667, -2.4475000, -…

Table 1 for Student-Level Variables

library(table1)

Attaching package: ‘table1’

The following objects are masked from ‘package:base’:

    units, units<-

table1(~ math + retained + black + hispanic, data = student.level)

	Overall (N=1721)
math
Mean (SD)	-0.574 (1.12)
Median [Min, Max]	-0.616 [-3.65, 3.14]
retained
Mean (SD)	0.206 (0.405)
Median [Min, Max]	0 [0, 1.00]
black
Mean (SD)	0.694 (0.461)
Median [Min, Max]	1.00 [0, 1.00]
hispanic
Mean (SD)	0.145 (0.352)
Median [Min, Max]	0 [0, 1.00]

NA

Create Summary Data Frame of School-Level Variables

school.level <-egmerged.clean %>%
  group_by(schoolid) %>%
  summarize(.,
            lowinc = mean(lowinc),
            mobility = mean(mobility),
            size = mean(size))

`summarise()` ungrouping output (override with `.groups` argument)

glimpse(school.level)

Rows: 60
Columns: 4
$ schoolid <dbl> 2020, 2040, 2180, 2330, 2340, 2380, 2390, 2440, 2480, 2520, 2540, 2560,…
$ lowinc   <dbl> 40.3, 83.1, 96.6, 78.9, 93.7, 36.9, 100.0, 44.8, 59.3, 76.9, 99.7, 94.3…
$ mobility <dbl> 12.5, 18.6, 44.4, 31.7, 67.0, 39.3, 39.9, 37.1, 25.8, 10.5, 42.7, 28.6,…
$ size     <dbl> 380, 502, 777, 800, 1133, 439, 566, 767, 113, 828, 339, 511, 1222, 364,…

Table 1 for Student-Level Variables

library(table1)
table1(~ lowinc + mobility + size, data = school.level)

	Overall (N=60)
lowinc
Mean (SD)	73.7 (27.3)
Median [Min, Max]	82.8 [0, 100]
mobility
Mean (SD)	34.7 (13.2)
Median [Min, Max]	35.8 [8.80, 67.0]
size
Mean (SD)	643 (317)
Median [Min, Max]	582 [113, 1490]

NA

density(egmerged.clean$math) %>%
  plot(., main = "Distribution of Math Scores")

Getting Fancy: Maybe plot over time?

scatter.smooth is a nice way to get a sense of the trend in the data. But keep in mind, this type of plot does NOT account for observations nested within students, so it should be treated as exploratory.

scatter.smooth(x = egmerged.clean$year, 
               y = egmerged.clean$math,
               main = "Distribution of Math Scores, Over Time")

Run a regular null model for math scores

model.null <- lmer(math ~ (1|childid), data = egmerged.clean)
summary(model.null)

Linear mixed model fit by REML ['lmerMod']
Formula: math ~ (1 | childid)
   Data: egmerged.clean

REML criterion at convergence: 25624.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7703 -0.6468  0.0158  0.6545  3.3731 

Random effects:
 Groups   Name        Variance Std.Dev.
 childid  (Intercept) 0.8483   0.921   
 Residual             1.5247   1.235   
Number of obs: 7230, groups:  childid, 1721

Fixed effects:
            Estimate Std. Error t value
(Intercept) -0.56065    0.02668  -21.02

Calculate ICC

null.icc <- 0.8483/(0.8483 + 1.5247)
null.icc

[1] 0.35748

Now, run a GROWTH null model for math scores (with year included)

model.null.growth <- lmer(math ~ year0 + (1|childid), data = egmerged.clean)
summary(model.null.growth)

Linear mixed model fit by REML ['lmerMod']
Formula: math ~ year0 + (1 | childid)
   Data: egmerged.clean

REML criterion at convergence: 17045.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.7528 -0.5789 -0.0243  0.5598  6.0968 

Random effects:
 Groups   Name        Variance Std.Dev.
 childid  (Intercept) 0.8683   0.9318  
 Residual             0.3470   0.5891  
Number of obs: 7230, groups:  childid, 1721

Fixed effects:
             Estimate Std. Error t value
(Intercept) -2.707304   0.028169  -96.11
year0        0.747452   0.005401  138.40

Correlation of Fixed Effects:
      (Intr)
year0 -0.548

Calculate GROWTH ICC

null.growth.icc <- 0.8683/(0.8683 + 0.3470)
null.growth.icc

[1] 0.7144738

Add a student-level predictor, female

model.1 <- lmer(math ~ year0 + female + (1|childid), data = egmerged.clean)
summary(model.1)

Linear mixed model fit by REML ['lmerMod']
Formula: math ~ year0 + female + (1 | childid)
   Data: egmerged.clean

REML criterion at convergence: 17049.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.7547 -0.5792 -0.0257  0.5598  6.0984 

Random effects:
 Groups   Name        Variance Std.Dev.
 childid  (Intercept) 0.8686   0.9320  
 Residual             0.3470   0.5891  
Number of obs: 7230, groups:  childid, 1721

Fixed effects:
             Estimate Std. Error t value
(Intercept) -2.721185   0.036900 -73.744
year0        0.747427   0.005401 138.387
female       0.027455   0.047131   0.583

Correlation of Fixed Effects:
       (Intr) year0 
year0  -0.414       
female -0.646 -0.007

Add predictors for race/ethnicity, black and hispanic…this is just for the intercept…

model.2 <- lmer(math ~ year0 + female + black + hispanic + (1|childid), data = egmerged.clean)
summary(model.2)

Linear mixed model fit by REML ['lmerMod']
Formula: math ~ year0 + female + black + hispanic + (1 | childid)
   Data: egmerged.clean

REML criterion at convergence: 16894

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8102 -0.5773 -0.0232  0.5564  6.1171 

Random effects:
 Groups   Name        Variance Std.Dev.
 childid  (Intercept) 0.7833   0.8850  
 Residual             0.3470   0.5891  
Number of obs: 7230, groups:  childid, 1721

Fixed effects:
             Estimate Std. Error t value
(Intercept) -2.069134   0.062751 -32.974
year0        0.747614   0.005396 138.561
female       0.017957   0.044996   0.399
black       -0.808716   0.062144 -13.014
hispanic    -0.597271   0.081412  -7.336

Correlation of Fixed Effects:
         (Intr) year0  female black 
year0    -0.242                     
female   -0.380 -0.008              
black    -0.810  0.000  0.018       
hispanic -0.623 -0.006  0.030  0.620

Do students differ in their growth by race/ethnicity? This is an interaction effect.

model.3 <- lmer(math ~ year0 + female + black + hispanic + year0:black + year0:hispanic + (1|childid), data = egmerged.clean)
summary(model.3)

Linear mixed model fit by REML ['lmerMod']
Formula: math ~ year0 + female + black + hispanic + year0:black + year0:hispanic +  
    (1 | childid)
   Data: egmerged.clean

REML criterion at convergence: 16830.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5762 -0.5761 -0.0243  0.5596  6.0517 

Random effects:
 Groups   Name        Variance Std.Dev.
 childid  (Intercept) 0.7854   0.8862  
 Residual             0.3423   0.5850  
Number of obs: 7230, groups:  childid, 1721

Fixed effects:
               Estimate Std. Error t value
(Intercept)    -2.29171    0.07063 -32.449
year0           0.82549    0.01256  65.713
female          0.01953    0.04502   0.434
black          -0.49632    0.07408  -6.700
hispanic       -0.53621    0.09859  -5.439
year0:black    -0.10970    0.01414  -7.758
year0:hispanic -0.02296    0.01920  -1.196

Correlation of Fixed Effects:
            (Intr) year0  female black  hispnc yr0:bl
year0       -0.506                                   
female      -0.339 -0.002                            
black       -0.849  0.483  0.016                     
hispanic    -0.642  0.363  0.024  0.605              
year0:black  0.451 -0.888 -0.003 -0.544 -0.322       
year0:hspnc  0.331 -0.654  0.001 -0.316 -0.563  0.581

#Add a random slope for year at the child level (each child has their own slope)

model.4 <- lmer(math ~ year0 + female + black + hispanic +  (year0|childid), data = egmerged.clean)
summary(model.4)

Linear mixed model fit by REML ['lmerMod']
Formula: math ~ year0 + female + black + hispanic + (year0 | childid)
   Data: egmerged.clean

REML criterion at convergence: 16645.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0852 -0.5547 -0.0399  0.5237  5.4432 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 childid  (Intercept) 0.60360  0.7769       
          year0       0.02137  0.1462   0.02
 Residual             0.30134  0.5489       
Number of obs: 7230, groups:  childid, 1721

Fixed effects:
             Estimate Std. Error t value
(Intercept) -2.144503   0.059617 -35.971
year0        0.748421   0.006352 117.826
female       0.005656   0.043468   0.130
black       -0.697309   0.059775 -11.665
hispanic    -0.566439   0.078695  -7.198

Correlation of Fixed Effects:
         (Intr) year0  female black 
year0    -0.189                     
female   -0.387 -0.008              
black    -0.817 -0.011  0.019       
hispanic -0.625 -0.016  0.032  0.615

Should we keep that random slope for time (year)? Use rand to test…

lmerTest::rand(model.4)

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

Model:
math ~ year0 + female + black + hispanic + (year0 | childid)
                           npar  logLik   AIC   LRT Df Pr(>Chisq)    
<none>                        9 -8322.8 16664                        
year0 in (year0 | childid)    7 -8447.0 16908 248.4  2  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Put it all together- random slope, interactions (not for Latinx students)

model.5 <- lmer(math ~ year0 + female + black + hispanic + year0:black + (year0|childid), data = egmerged.clean)
summary(model.5)

Linear mixed model fit by REML ['lmerMod']
Formula: math ~ year0 + female + black + hispanic + year0:black + (year0 |      childid)
   Data: egmerged.clean

REML criterion at convergence: 16598.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1733 -0.5567 -0.0362  0.5267  5.4746 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 childid  (Intercept) 0.59434  0.7709       
          year0       0.01911  0.1383   0.07
 Residual             0.30155  0.5491       
Number of obs: 7230, groups:  childid, 1721

Fixed effects:
             Estimate Std. Error t value
(Intercept) -2.267271   0.061772 -36.704
year0        0.816435   0.011093  73.601
female       0.005709   0.043453   0.131
black       -0.514994   0.064558  -7.977
hispanic    -0.579556   0.078694  -7.365
year0:black -0.099556   0.013407  -7.426

Correlation of Fixed Effects:
            (Intr) year0  female black  hispnc
year0       -0.319                            
female      -0.374 -0.004                     
black       -0.830  0.307  0.018              
hispanic    -0.597 -0.028  0.031  0.561       
year0:black  0.266 -0.827  0.000 -0.378  0.023

Visualize Our Model by Predicting Scores…

predicted.vals <- broom.mixed::augment(model.5, effects = c("ran_pars", "fixed"), conf.int = TRUE)

Package version inconsistency detected.
TMB was built with Matrix version 1.2.18
Current Matrix version is 1.2.17
Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' packageRegistered S3 method overwritten by 'broom.mixed':
  method      from 
  tidy.gamlss broom

glimpse(predicted.vals)

Rows: 7,230
Columns: 17
$ math     <dbl> -3.887, -4.055, -3.525, -3.592, -3.525, -3.068, -2.596, -3.721, -3.232,…
$ year0    <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,…
$ female   <dbl> 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1,…
$ black    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1,…
$ hispanic <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,…
$ childid  <dbl> 289376791, 288278731, 292789461, 275898121, 300246721, 244834951, 29685…
$ .fitted  <dbl> -3.585709, -4.002690, -2.368814, -2.936436, -3.020023, -3.051652, -2.49…
$ .resid   <dbl> -0.30129073, -0.05231038, -1.15618618, -0.65556371, -0.50497672, -0.016…
$ .hat     <dbl> 0.2930998, 0.2930998, 0.3195197, 0.3195197, 0.3195197, 0.2933761, 0.293…
$ .cooksd  <dbl> 2.942833e-02, 8.870945e-04, 5.098199e-01, 1.639046e-01, 9.725321e-02, 8…
$ .fixed   <dbl> -2.776557, -2.776557, -2.782265, -2.782265, -2.782265, -2.782265, -2.78…
$ .mu      <dbl> -3.585709, -4.002690, -2.368814, -2.936436, -3.020023, -3.051652, -2.49…
$ .offset  <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,…
$ .sqrtXwt <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,…
$ .sqrtrwt <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,…
$ .weights <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,…
$ .wtres   <dbl> -0.30129073, -0.05231038, -1.15618618, -0.65556371, -0.50497672, -0.016…

… and Then Plotting Them…

ggplot(data = predicted.vals, mapping = aes(x = year0, y = .fitted, color = as_factor(black))) + 
  geom_jitter(alpha = .50) + 
  geom_smooth(method = "lm", formula = 'y ~ x') +
  labs(title = "Predicted Change in Math Scores Over Time",
      subtitle = "Black Students (Teal Line) vs. All Others (Orange Line)",
      caption = "Notes. N = 1,731 students in the Sustaining Effects Study (Carter, 1984). Year 0 = Kindergarten.") +
  ylab("Predicted Math Score") +
  xlab("Year") +
  theme(legend.position = "none")

Random Slopes for a Sample of 50 Students…To Get a Sense of the Variability in Slopes…

predicted.vals50 <- predicted.vals %>%
  arrange(childid, year0) %>%
  slice(., 1:252)

ggplot(data = predicted.vals50, mapping = aes(x = year0, y = .fitted, color = as_factor(childid))) + 
  geom_jitter(alpha = .50) + 
  geom_smooth(method = "lm", formula = 'y ~ x') +
  labs(title = "Predicted Change in Math Scores Over Time",
      caption = "Notes. N = 1,731 students in the Sustaining Effects Study (Carter, 1984). Year 0 = Kindergarten.") +
  ylab("Predicted Math Score") +
  xlab("Year") +
  theme(legend.position = "none")

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)

package ‘broom.mixed’ was built under R version 3.6.2

models <- list(model.null, model.null.growth, 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.561	-2.707	-2.721	-2.069	-2.292	-2.145	-2.267
	(0.027)	(0.028)	(0.037)	(0.063)	(0.071)	(0.060)	(0.062)
sd__(Intercept)	0.921	0.932	0.932	0.885	0.886	0.777	0.771
sd__Observation	1.235	0.589	0.589	0.589	0.585	0.549	0.549
year0		0.747	0.747	0.748	0.825	0.748	0.816
		(0.005)	(0.005)	(0.005)	(0.013)	(0.006)	(0.011)
female			0.027	0.018	0.020	0.006	0.006
			(0.047)	(0.045)	(0.045)	(0.043)	(0.043)
black				-0.809	-0.496	-0.697	-0.515
				(0.062)	(0.074)	(0.060)	(0.065)
hispanic				-0.597	-0.536	-0.566	-0.580
				(0.081)	(0.099)	(0.079)	(0.079)
year0 × black					-0.110		-0.100
					(0.014)		(0.013)
year0 × hispanic					-0.023
					(0.019)
cor__(Intercept).year0						0.024	0.069
sd__year0						0.146	0.138
AIC	25630.4	17053.1	17059.1	16908.0	16848.9	16663.6	16618.4
BIC	25651.1	17080.7	17093.5	16956.2	16910.9	16725.6	16687.2
Log.Lik.	-12812.206	-8522.569	-8524.536	-8447.005	-8415.462	-8322.803	-8299.183
REMLcrit	25624.412	17045.139	17049.072	16894.010	16830.924	16645.607	16598.367

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'