This week, we continue with the data from the US Sustaining Effects Study (Carter 1984), which attempted to assess the sustaining impacts of compensatory schooling for students over time. We will dig into nonlinear growth patterns and heterogeneous variance models, two advanced topics in growth modeling.

Load in Our MVP Packages

suppressPackageStartupMessages(library(tidyverse))
library(lme4)

Loading required package: Matrix

Attaching package: ‘Matrix’

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

    expand, pack, unpack

Part 1: Loading in the Data, Reshaping Data from Long to Wide, Create New Variables

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, 361…
$ childid  <dbl> 289376791, 288278731, 292789461, 275898121, 300246721, 24…
$ year     <dbl> 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, …
$ math     <dbl> -3.887, -4.055, -3.525, -3.592, -3.525, -3.068, -2.596, -…
$ retained <dbl> 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, …
$ black    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, …
$ hispanic <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, …
$ size     <dbl> 641, 1062, 777, 724, 185, 526, 710, 532, 1062, 1062, 1062…
$ lowinc   <dbl> 84.8, 97.8, 96.6, 45.3, 94.6, 100.0, 66.2, 66.4, 97.8, 97…
$ mobility <dbl> 38.2, 26.8, 44.4, 18.7, 47.4, 55.4, 25.4, 20.9, 26.8, 26.…

Super Handy Trick: Reshaping Data

Using `pivot_wider` to reshape from long to wide….

egmerged.wide <- egmerged %>%
  pivot_wider(.,
              id_cols = c("childid", "schoolid", "female", "black", "hispanic", "size", "lowinc", "mobility"),
              names_from = year,
              values_from = math,
              names_prefix = "math_year",
              )

glimpse(egmerged.wide)

Rows: 1,721
Columns: 14
$ childid    <dbl> 289376791, 288278731, 292789461, 275898121, 300246721, …
$ schoolid   <dbl> 3430, 3610, 2180, 3040, 3221, 3360, 4450, 4390, 3610, 3…
$ female     <dbl> 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1…
$ black      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1…
$ hispanic   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0…
$ size       <dbl> 641, 1062, 777, 724, 185, 526, 710, 532, 1062, 1062, 10…
$ lowinc     <dbl> 84.8, 97.8, 96.6, 45.3, 94.6, 100.0, 66.2, 66.4, 97.8, …
$ mobility   <dbl> 38.2, 26.8, 44.4, 18.7, 47.4, 55.4, 25.4, 20.9, 26.8, 2…
$ math_year1 <dbl> -3.887, -4.055, -3.525, -3.592, -3.525, -3.068, -2.596,…
$ math_year2 <dbl> -2.839, -3.395, -0.365, -1.130, -2.097, -1.538, -1.218,…
$ math_year3 <dbl> -2.596, -3.395, -0.574, -2.009, -1.314, -1.757, -1.120,…
$ math_year4 <dbl> -1.236, -1.641, NA, NA, NA, -2.270, -0.419, NA, -0.522,…
$ math_year5 <dbl> -0.456, -2.263, NA, NA, NA, -1.484, -0.403, NA, -0.945,…
$ math_year6 <dbl> 0.601, -1.552, NA, NA, NA, NA, NA, NA, -0.170, -1.763, …

Using `pivot_longer` to reshape from wide to long….

egmerged.long <- egmerged.wide %>%
  pivot_longer(.,
              cols = starts_with("math"),
              names_to = "year",
              values_to = "math",
              names_prefix = "math_year",
              )

glimpse(egmerged.long)

Rows: 10,326
Columns: 10
$ childid  <dbl> 289376791, 289376791, 289376791, 289376791, 289376791, 28…
$ schoolid <dbl> 3430, 3430, 3430, 3430, 3430, 3430, 3610, 3610, 3610, 361…
$ female   <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, …
$ black    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ hispanic <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ size     <dbl> 641, 641, 641, 641, 641, 641, 1062, 1062, 1062, 1062, 106…
$ lowinc   <dbl> 84.8, 84.8, 84.8, 84.8, 84.8, 84.8, 97.8, 97.8, 97.8, 97.…
$ mobility <dbl> 38.2, 38.2, 38.2, 38.2, 38.2, 38.2, 26.8, 26.8, 26.8, 26.…
$ year     <chr> "1", "2", "3", "4", "5", "6", "1", "2", "3", "4", "5", "6…
$ math     <dbl> -3.887, -2.839, -2.596, -1.236, -0.456, 0.601, -4.055, -3…

A little more cleanup needed - change `year` to numeric, get ridding of observations with missing math scores…

egmerged.long.clean <- egmerged.long %>%
  mutate(.,
         year  = as.numeric(year)) %>%
  na.omit()

glimpse(egmerged.long.clean)

Rows: 7,230
Columns: 10
$ childid  <dbl> 289376791, 289376791, 289376791, 289376791, 289376791, 28…
$ schoolid <dbl> 3430, 3430, 3430, 3430, 3430, 3430, 3610, 3610, 3610, 361…
$ female   <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, …
$ black    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ hispanic <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ size     <dbl> 641, 641, 641, 641, 641, 641, 1062, 1062, 1062, 1062, 106…
$ lowinc   <dbl> 84.8, 84.8, 84.8, 84.8, 84.8, 84.8, 97.8, 97.8, 97.8, 97.…
$ mobility <dbl> 38.2, 38.2, 38.2, 38.2, 38.2, 38.2, 26.8, 26.8, 26.8, 26.…
$ year     <dbl> 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 1, 2, 3, 1, …
$ math     <dbl> -3.887, -2.839, -2.596, -1.236, -0.456, 0.601, -4.055, -3…

A little housekeeping - creating squared and cubic terms for time…

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,
         year_sq = year0^2,
         year_cube = year0^3)

glimpse(egmerged.clean)

Rows: 7,230
Columns: 19
$ schoolid     <dbl> 3430, 3610, 2180, 3040, 3221, 3360, 4450, 4390, 3610,…
$ childid      <dbl> 289376791, 288278731, 292789461, 275898121, 300246721…
$ year         <dbl> 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,…
$ math         <dbl> -3.887, -4.055, -3.525, -3.592, -3.525, -3.068, -2.59…
$ retained     <dbl> 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,…
$ black        <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,…
$ hispanic     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,…
$ size         <dbl> 641, 1062, 777, 724, 185, 526, 710, 532, 1062, 1062, …
$ lowinc       <dbl> 84.8, 97.8, 96.6, 45.3, 94.6, 100.0, 66.2, 66.4, 97.8…
$ mobility     <dbl> 38.2, 26.8, 44.4, 18.7, 47.4, 55.4, 25.4, 20.9, 26.8,…
$ retained.fac <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ female.fac   <fct> 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0,…
$ black.fac    <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,…
$ hispanic.fac <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,…
$ year0        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ year_sq      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ year_cube    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…

Part 2: Nonlinear Growth Models

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 nonlinear growth terms (square)

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

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

REML criterion at convergence: 16951.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.6833 -0.5749 -0.0259  0.5607  6.2588 

Random effects:
 Groups   Name        Variance Std.Dev.
 childid  (Intercept) 0.8669   0.9311  
 Residual             0.3409   0.5839  
Number of obs: 7230, groups:  childid, 1721

Fixed effects:
            Estimate Std. Error t value
(Intercept) -2.96054    0.03746  -79.03
year0        0.97320    0.02276   42.77
year_sq     -0.03890    0.00381  -10.21

Correlation of Fixed Effects:
        (Intr) year0 
year0   -0.740       
year_sq  0.662 -0.972

Add nonlinear growth terms (cubic)

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

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

REML criterion at convergence: 16980.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.6878 -0.5733 -0.0233  0.5630  6.2327 

Random effects:
 Groups   Name        Variance Std.Dev.
 childid  (Intercept) 0.8677   0.9315  
 Residual             0.3423   0.5851  
Number of obs: 7230, groups:  childid, 1721

Fixed effects:
              Estimate Std. Error t value
(Intercept) -2.8770188  0.0339696  -84.69
year0        0.8617959  0.0139383   61.83
year_cube   -0.0039513  0.0004445   -8.89

Correlation of Fixed Effects:
          (Intr) year0 
year0     -0.693       
year_cube  0.562 -0.923

Add nonlinear growth terms (square AND cubic)

model.both <- lmer(math ~ year0 + year_sq + year_cube + (1|childid), REML = FALSE, data = egmerged.clean)
summary(model.both)

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

     AIC      BIC   logLik deviance df.resid 
 16882.5  16923.8  -8435.3  16870.5     7224 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.7115 -0.5724 -0.0298  0.5588  6.2789 

Random effects:
 Groups   Name        Variance Std.Dev.
 childid  (Intercept) 0.8630   0.9290  
 Residual             0.3378   0.5812  
Number of obs: 7230, groups:  childid, 1721

Fixed effects:
             Estimate Std. Error t value
(Intercept) -3.223851   0.051032 -63.173
year0        1.408171   0.061794  22.788
year_sq     -0.220628   0.024317  -9.073
year_cube    0.021417   0.002831   7.566

Correlation of Fixed Effects:
          (Intr) year0  yer_sq
year0     -0.833              
year_sq    0.749 -0.975       
year_cube -0.682  0.930 -0.988

Should we keep that cubic term time (year)? Use `anova` to test…

anova(model.square, model.both)

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

Data: egmerged.clean
Models:
model.square: math ~ year0 + year_sq + (1 | childid)
model.both: math ~ year0 + year_sq + year_cube + (1 | childid)
             npar   AIC   BIC  logLik deviance  Chisq Df Pr(>Chisq)    
model.square    5 16938 16972 -8463.8    16928                         
model.both      6 16882 16924 -8435.3    16870 57.008  1   4.34e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Same thing, PLUS adding a random slope for year0

model.both.rand <- lmer(math ~ year0 + year_sq + year_cube + (year0|childid), REML = FALSE, data = egmerged.clean)
summary(model.both.rand)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: math ~ year0 + year_sq + year_cube + (year0 | childid)
   Data: egmerged.clean

     AIC      BIC   logLik deviance df.resid 
 16529.3  16584.4  -8256.6  16513.3     7222 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0747 -0.5525 -0.0342  0.5293  5.6148 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 childid  (Intercept) 0.61268  0.7827       
          year0       0.02458  0.1568   0.09
 Residual             0.28497  0.5338       
Number of obs: 7230, groups:  childid, 1721

Fixed effects:
             Estimate Std. Error t value
(Intercept) -3.322138   0.048116 -69.044
year0        1.541861   0.059322  25.991
year_sq     -0.269158   0.023175 -11.614
year_cube    0.026471   0.002685   9.857

Correlation of Fixed Effects:
          (Intr) year0  yer_sq
year0     -0.852              
year_sq    0.773 -0.974       
year_cube -0.706  0.930 -0.988

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

lmerTest::rand(model.both.rand)

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

Model:
math ~ year0 + year_sq + year_cube + (year0 | childid)
                           npar  logLik   AIC    LRT Df Pr(>Chisq)    
<none>                        8 -8256.6 16529                         
year0 in (year0 | childid)    6 -8435.3 16882 357.22  2  < 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:

	Model 1	Model 2	Model 3	Model 4	Model 5
(Intercept)	-2.707	-2.961	-2.877	-3.224	-3.322
	(0.028)	(0.037)	(0.034)	(0.051)	(0.048)
year0	0.747	0.973	0.862	1.408	1.542
	(0.005)	(0.023)	(0.014)	(0.062)	(0.059)
sd__(Intercept)	0.932	0.931	0.931	0.929	0.783
sd__Observation	0.589	0.584	0.585	0.581	0.534
year_sq		-0.039		-0.221	-0.269
		(0.004)		(0.024)	(0.023)
year_cube			-0.004	0.021	0.026
			(0.000)	(0.003)	(0.003)
cor__(Intercept).year0					0.090
sd__year0					0.157
AIC	17053.1	16961.1	16990.2	16882.5	16529.3
BIC	17080.7	16995.5	17024.6	16923.8	16584.4
Log.Lik.	-8522.569	-8475.554	-8490.104	-8435.256	-8256.647
REMLcrit	17045.139	16951.108	16980.208

Visualize Our Model by Predicting Scores…

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

glimpse(predicted.vals)

Rows: 7,230
Columns: 16
$ math      <dbl> -3.887, -4.055, -3.525, -3.592, -3.525, -3.068, -2.596, …
$ year0     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ year_sq   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ year_cube <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ childid   <dbl> 289376791, 288278731, 292789461, 275898121, 300246721, 2…
$ .fitted   <dbl> -4.058049, -4.418839, -2.812483, -3.365132, -3.456950, -…
$ .resid    <dbl> 0.17104854, 0.36383866, -0.71251725, -0.22686787, -0.068…
$ .hat      <dbl> 0.3167563, 0.3167563, 0.3360663, 0.3360663, 0.3360663, 0…
$ .cooksd   <dbl> 0.017416159, 0.078800926, 0.339551761, 0.034423991, 0.00…
$ .fixed    <dbl> -3.322138, -3.322138, -3.322138, -3.322138, -3.322138, -…
$ .mu       <dbl> -4.058049, -4.418839, -2.812483, -3.365132, -3.456950, -…
$ .offset   <dbl> 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,…
$ .sqrtrwt  <dbl> 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,…
$ .wtres    <dbl> 0.17104854, 0.36383866, -0.71251725, -0.22686787, -0.068…

… and Then Plotting Them…

ggplot(data = predicted.vals, mapping = aes(x = year0, y = .fitted)) + 
  geom_jitter(alpha = .40, color = "seagreen1") + 
  geom_smooth(method = "loess", formula = 'y ~ x', color = "seagreen") +
  labs(title = "Predicted Change in Math Scores Over Time",
      subtitle = "With Square and Cubic Effects for Nonlinear Growth",
      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_point(alpha = .50) + 
  geom_line() +
  labs(title = "Predicted Change in Math Scores Over Time",
      subtitle = "With Square and Cubic Effects for Nonlinear Growth",
      caption = "Notes. Subset of N = 50 students in the Sustaining Effects Study (Carter, 1984). Year 0 = Kindergarten.") +
  ylab("Predicted Math Score") +
  xlab("Year") +
  theme(legend.position = "none")

Part 2: Heterogeneous Variance Repeated Measures Model

library(nlme)


Attaching package: ‘nlme’

The following object is masked from ‘package:lme4’:

    lmList

The following object is masked from ‘package:dplyr’:

    collapse

model.null.growth.homog <- lme(math ~ year0, 
                         random = ~ year0|childid, 
                         data = egmerged.clean,
                         method = "ML")
summary(model.null.growth.homog)

Linear mixed-effects model fit by maximum likelihood
 Data: egmerged.clean


Random effects:
 Formula: ~year0 | childid
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 0.7922145 (Intr)
year0       0.1461291 0.113 
Residual    0.5488263       

Fixed effects: math ~ year0

 Correlation: 
      (Intr)
year0 -0.442

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.07338195 -0.56120525 -0.03801274  0.52482971  5.46800933 

Number of Observations: 7230
Number of Groups: 1721

model.null.growth.het <- lme(math ~ year0, 
                             random =  ~ year0|childid, 
                             data = egmerged.clean, 
                             method = "ML", 
                             weights = varIdent(form = ~ 1 | year0)
)
summary(model.null.growth.het)

Linear mixed-effects model fit by maximum likelihood
 Data: egmerged.clean


Random effects:
 Formula: ~year0 | childid
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 0.7909930 (Intr)
year0       0.1526271 0.069 
Residual    0.5643307       

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | year0 
 Parameter estimates:
        2         3         4         5         1         0 
1.0000000 1.0900648 0.8026221 0.6839841 1.1024842 1.0356196 
Fixed effects: math ~ year0

 Correlation: 
      (Intr)
year0 -0.455

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.62560022 -0.55340209 -0.03040342  0.53323108  5.05558475 

Number of Observations: 7230
Number of Groups: 1721

LS0tCnRpdGxlOiAnTXVsdGlsZXZlbCBNb2RlbGluZywgTW9kdWxlIDgsIFBhcnQgMjogQWR2YW5jZWQgR3Jvd3RoIE1vZGVscycKYXV0aG9yOiAnRHIuIEJyb2RhJwpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpUaGlzIHdlZWssIHdlIGNvbnRpbnVlIHdpdGggdGhlIGRhdGEgZnJvbSB0aGUgVVMgU3VzdGFpbmluZyBFZmZlY3RzIFN0dWR5IChDYXJ0ZXIgMTk4NCksIHdoaWNoIGF0dGVtcHRlZCB0byBhc3Nlc3MgdGhlIHN1c3RhaW5pbmcgaW1wYWN0cyBvZiBjb21wZW5zYXRvcnkgc2Nob29saW5nIGZvciBzdHVkZW50cyBvdmVyIHRpbWUuIFdlIHdpbGwgZGlnIGludG8gbm9ubGluZWFyIGdyb3d0aCBwYXR0ZXJucyBhbmQgaGV0ZXJvZ2VuZW91cyB2YXJpYW5jZSBtb2RlbHMsIHR3byBhZHZhbmNlZCB0b3BpY3MgaW4gZ3Jvd3RoIG1vZGVsaW5nLgoKIyBMb2FkIGluIE91ciBNVlAgUGFja2FnZXMKYGBge3J9CnN1cHByZXNzUGFja2FnZVN0YXJ0dXBNZXNzYWdlcyhsaWJyYXJ5KHRpZHl2ZXJzZSkpCmxpYnJhcnkobG1lNCkKYGBgCgojIFBhcnQgMTogTG9hZGluZyBpbiB0aGUgRGF0YSwgUmVzaGFwaW5nIERhdGEgZnJvbSBMb25nIHRvIFdpZGUsIENyZWF0ZSBOZXcgVmFyaWFibGVzCiMjIExvYWQgaW4gdGhlIERhdGEKYGBge3J9CgplZ21lcmdlZCA8LSBoYXZlbjo6cmVhZF9kdGEoImVnbWVyZ2VkLmR0YSIpCgpnbGltcHNlKGVnbWVyZ2VkKQpgYGAKCiMjIFN1cGVyIEhhbmR5IFRyaWNrOiBSZXNoYXBpbmcgRGF0YQojIyBVc2luZyBgcGl2b3Rfd2lkZXJgIHRvIHJlc2hhcGUgZnJvbSBsb25nIHRvIHdpZGUuLi4uCmBgYHtyfQplZ21lcmdlZC53aWRlIDwtIGVnbWVyZ2VkICU+JQogIHBpdm90X3dpZGVyKC4sCiAgICAgICAgICAgICAgaWRfY29scyA9IGMoImNoaWxkaWQiLCAic2Nob29saWQiLCAiZmVtYWxlIiwgImJsYWNrIiwgImhpc3BhbmljIiwgInNpemUiLCAibG93aW5jIiwgIm1vYmlsaXR5IiksCiAgICAgICAgICAgICAgbmFtZXNfZnJvbSA9IHllYXIsCiAgICAgICAgICAgICAgdmFsdWVzX2Zyb20gPSBtYXRoLAogICAgICAgICAgICAgIG5hbWVzX3ByZWZpeCA9ICJtYXRoX3llYXIiLAogICAgICAgICAgICAgICkKCmdsaW1wc2UoZWdtZXJnZWQud2lkZSkKYGBgCgojIyBVc2luZyBgcGl2b3RfbG9uZ2VyYCB0byByZXNoYXBlIGZyb20gd2lkZSB0byBsb25nLi4uLgpgYGB7cn0KZWdtZXJnZWQubG9uZyA8LSBlZ21lcmdlZC53aWRlICU+JQogIHBpdm90X2xvbmdlciguLAogICAgICAgICAgICAgIGNvbHMgPSBzdGFydHNfd2l0aCgibWF0aCIpLAogICAgICAgICAgICAgIG5hbWVzX3RvID0gInllYXIiLAogICAgICAgICAgICAgIHZhbHVlc190byA9ICJtYXRoIiwKICAgICAgICAgICAgICBuYW1lc19wcmVmaXggPSAibWF0aF95ZWFyIiwKICAgICAgICAgICAgICApCgpnbGltcHNlKGVnbWVyZ2VkLmxvbmcpCmBgYAoKIyMgQSBsaXR0bGUgbW9yZSBjbGVhbnVwIG5lZWRlZCAtIGNoYW5nZSBgeWVhcmAgdG8gbnVtZXJpYywgZ2V0IHJpZGRpbmcgb2Ygb2JzZXJ2YXRpb25zIHdpdGggbWlzc2luZyBtYXRoIHNjb3Jlcy4uLgpgYGB7cn0KZWdtZXJnZWQubG9uZy5jbGVhbiA8LSBlZ21lcmdlZC5sb25nICU+JQogIG11dGF0ZSguLAogICAgICAgICB5ZWFyICA9IGFzLm51bWVyaWMoeWVhcikpICU+JQogIG5hLm9taXQoKQoKZ2xpbXBzZShlZ21lcmdlZC5sb25nLmNsZWFuKQpgYGAKCiMjIEEgbGl0dGxlIGhvdXNla2VlcGluZyAtIGNyZWF0aW5nIHNxdWFyZWQgYW5kIGN1YmljIHRlcm1zIGZvciB0aW1lLi4uCmBgYHtyfQplZ21lcmdlZC5jbGVhbiA8LSBlZ21lcmdlZCAlPiUKICBtdXRhdGUoLiwKICAgICAgICAgcmV0YWluZWQuZmFjID0gYXNfZmFjdG9yKHJldGFpbmVkKSwKICAgICAgICAgZmVtYWxlLmZhYyA9IGFzX2ZhY3RvcihmZW1hbGUpLAogICAgICAgICBibGFjay5mYWMgPSBhc19mYWN0b3IoYmxhY2spLAogICAgICAgICBoaXNwYW5pYy5mYWMgPSBhc19mYWN0b3IoaGlzcGFuaWMpLAogICAgICAgICB5ZWFyMCA9IHllYXIgLSAxLAogICAgICAgICB5ZWFyX3NxID0geWVhcjBeMiwKICAgICAgICAgeWVhcl9jdWJlID0geWVhcjBeMykKCmdsaW1wc2UoZWdtZXJnZWQuY2xlYW4pCmBgYAoKIyBQYXJ0IDI6IE5vbmxpbmVhciBHcm93dGggTW9kZWxzCiMjIE5vdywgcnVuIGEgR1JPV1RIIG51bGwgbW9kZWwgZm9yIG1hdGggc2NvcmVzICh3aXRoIHllYXIgaW5jbHVkZWQpCmBgYHtyfQptb2RlbC5udWxsLmdyb3d0aCA8LSBsbWVyKG1hdGggfiB5ZWFyMCArICgxfGNoaWxkaWQpLCBkYXRhID0gZWdtZXJnZWQuY2xlYW4pCnN1bW1hcnkobW9kZWwubnVsbC5ncm93dGgpCmBgYAoKIyMjIENhbGN1bGF0ZSBHUk9XVEggSUNDCmBgYHtyfQpudWxsLmdyb3d0aC5pY2MgPC0gMC44NjgzLygwLjg2ODMgKyAwLjM0NzApCm51bGwuZ3Jvd3RoLmljYwpgYGAKCiMjIEFkZCBub25saW5lYXIgZ3Jvd3RoIHRlcm1zIChzcXVhcmUpCmBgYHtyfQptb2RlbC5zcXVhcmUgPC0gbG1lcihtYXRoIH4geWVhcjAgKyB5ZWFyX3NxICsgKDF8Y2hpbGRpZCksIGRhdGEgPSBlZ21lcmdlZC5jbGVhbikKc3VtbWFyeShtb2RlbC5zcXVhcmUpCmBgYAoKIyMgQWRkIG5vbmxpbmVhciBncm93dGggdGVybXMgKGN1YmljKQpgYGB7cn0KbW9kZWwuY3ViZSA8LSBsbWVyKG1hdGggfiB5ZWFyMCArIHllYXJfY3ViZSArICgxfGNoaWxkaWQpLCBkYXRhID0gZWdtZXJnZWQuY2xlYW4pCnN1bW1hcnkobW9kZWwuY3ViZSkKYGBgCgojIyBBZGQgbm9ubGluZWFyIGdyb3d0aCB0ZXJtcyAoc3F1YXJlIEFORCBjdWJpYykKYGBge3J9Cm1vZGVsLmJvdGggPC0gbG1lcihtYXRoIH4geWVhcjAgKyB5ZWFyX3NxICsgeWVhcl9jdWJlICsgKDF8Y2hpbGRpZCksIFJFTUwgPSBGQUxTRSwgZGF0YSA9IGVnbWVyZ2VkLmNsZWFuKQpzdW1tYXJ5KG1vZGVsLmJvdGgpCmBgYAoKIyMgU2hvdWxkIHdlIGtlZXAgdGhhdCBjdWJpYyB0ZXJtIHRpbWUgKHllYXIpPyBVc2UgYGFub3ZhYCB0byB0ZXN0Li4uCmBgYHtyfQphbm92YShtb2RlbC5zcXVhcmUsIG1vZGVsLmJvdGgpCmBgYAoKIyMgU2FtZSB0aGluZywgUExVUyBhZGRpbmcgYSByYW5kb20gc2xvcGUgZm9yIHllYXIwCmBgYHtyfQptb2RlbC5ib3RoLnJhbmQgPC0gbG1lcihtYXRoIH4geWVhcjAgKyB5ZWFyX3NxICsgeWVhcl9jdWJlICsgKHllYXIwfGNoaWxkaWQpLCBSRU1MID0gRkFMU0UsIGRhdGEgPSBlZ21lcmdlZC5jbGVhbikKc3VtbWFyeShtb2RlbC5ib3RoLnJhbmQpCmBgYAoKIyMgU2hvdWxkIHdlIGtlZXAgdGhhdCByYW5kb20gc2xvcGUgZm9yIHRpbWUgKHllYXIpPyBVc2UgYHJhbmRgIHRvIHRlc3QuLi4KYGBge3J9CmxtZXJUZXN0OjpyYW5kKG1vZGVsLmJvdGgucmFuZCkKYGBgCgojIFVzaW5nIHRoZSBgbW9kZWxzdW1tYXJ5YCBhbmQgYGJyb29tLm1peGVkYCBQYWNrYWdlcyB0byBPcmdhbml6ZSBZb3VyIFJlc3VsdHM6CmBgYHtyfQpsaWJyYXJ5KG1vZGVsc3VtbWFyeSkKbGlicmFyeShicm9vbS5taXhlZCkKCm1vZGVsczEgPC0gbGlzdChtb2RlbC5udWxsLmdyb3d0aCwgbW9kZWwuc3F1YXJlLCBtb2RlbC5jdWJlLCBtb2RlbC5ib3RoLCBtb2RlbC5ib3RoLnJhbmQpCm1vZGVsc3VtbWFyeShtb2RlbHMxLCBvdXRwdXQgPSAiZGVmYXVsdCIpCmBgYAojIFZpc3VhbGl6ZSBPdXIgTW9kZWwgYnkgUHJlZGljdGluZyBTY29yZXMuLi4KYGBge3J9CnByZWRpY3RlZC52YWxzIDwtIGJyb29tLm1peGVkOjphdWdtZW50KG1vZGVsLmJvdGgucmFuZCwgZWZmZWN0cyA9IGMoInJhbl9wYXJzIiwgImZpeGVkIiksIGNvbmYuaW50ID0gVFJVRSkKCmdsaW1wc2UocHJlZGljdGVkLnZhbHMpCmBgYAoKIyAuLi4gYW5kIFRoZW4gUGxvdHRpbmcgVGhlbS4uLgpgYGB7cn0KZ2dwbG90KGRhdGEgPSBwcmVkaWN0ZWQudmFscywgbWFwcGluZyA9IGFlcyh4ID0geWVhcjAsIHkgPSAuZml0dGVkKSkgKyAKICBnZW9tX2ppdHRlcihhbHBoYSA9IC40MCwgY29sb3IgPSAic2VhZ3JlZW4xIikgKyAKICBnZW9tX3Ntb290aChtZXRob2QgPSAibG9lc3MiLCBmb3JtdWxhID0gJ3kgfiB4JywgY29sb3IgPSAic2VhZ3JlZW4iKSArCiAgbGFicyh0aXRsZSA9ICJQcmVkaWN0ZWQgQ2hhbmdlIGluIE1hdGggU2NvcmVzIE92ZXIgVGltZSIsCiAgICAgIHN1YnRpdGxlID0gIldpdGggU3F1YXJlIGFuZCBDdWJpYyBFZmZlY3RzIGZvciBOb25saW5lYXIgR3Jvd3RoIiwKICAgICAgY2FwdGlvbiA9ICJOb3Rlcy4gTiA9IDEsNzMxIHN0dWRlbnRzIGluIHRoZSBTdXN0YWluaW5nIEVmZmVjdHMgU3R1ZHkgKENhcnRlciwgMTk4NCkuIFllYXIgMCA9IEtpbmRlcmdhcnRlbi4iKSArCiAgeWxhYigiUHJlZGljdGVkIE1hdGggU2NvcmUiKSArCiAgeGxhYigiWWVhciIpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAibm9uZSIpCmBgYAoKIyBSYW5kb20gU2xvcGVzIGZvciBhIFNhbXBsZSBvZiA1MCBTdHVkZW50cy4uLlRvIEdldCBhIFNlbnNlIG9mIHRoZSBWYXJpYWJpbGl0eSBpbiBTbG9wZXMuLi4KYGBge3J9CgpwcmVkaWN0ZWQudmFsczUwIDwtIHByZWRpY3RlZC52YWxzICU+JQogIGFycmFuZ2UoY2hpbGRpZCwgeWVhcjApICU+JQogIHNsaWNlKC4sIDE6MjUyKQoKZ2dwbG90KGRhdGEgPSBwcmVkaWN0ZWQudmFsczUwLCBtYXBwaW5nID0gYWVzKHggPSB5ZWFyMCwgeSA9IC5maXR0ZWQsIGNvbG9yID0gYXNfZmFjdG9yKGNoaWxkaWQpKSkgKyAKICBnZW9tX3BvaW50KGFscGhhID0gLjUwKSArIAogIGdlb21fbGluZSgpICsKICBsYWJzKHRpdGxlID0gIlByZWRpY3RlZCBDaGFuZ2UgaW4gTWF0aCBTY29yZXMgT3ZlciBUaW1lIiwKICAgICAgc3VidGl0bGUgPSAiV2l0aCBTcXVhcmUgYW5kIEN1YmljIEVmZmVjdHMgZm9yIE5vbmxpbmVhciBHcm93dGgiLAogICAgICBjYXB0aW9uID0gIk5vdGVzLiBTdWJzZXQgb2YgTiA9IDUwIHN0dWRlbnRzIGluIHRoZSBTdXN0YWluaW5nIEVmZmVjdHMgU3R1ZHkgKENhcnRlciwgMTk4NCkuIFllYXIgMCA9IEtpbmRlcmdhcnRlbi4iKSArCiAgeWxhYigiUHJlZGljdGVkIE1hdGggU2NvcmUiKSArCiAgeGxhYigiWWVhciIpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAibm9uZSIpCgpgYGAKCiMgUGFydCAyOiBIZXRlcm9nZW5lb3VzIFZhcmlhbmNlIFJlcGVhdGVkIE1lYXN1cmVzIE1vZGVsCmBgYHtyfQpsaWJyYXJ5KG5sbWUpCgptb2RlbC5udWxsLmdyb3d0aC5ob21vZyA8LSBsbWUobWF0aCB+IHllYXIwLCAKICAgICAgICAgICAgICAgICAgICAgICAgIHJhbmRvbSA9IH4geWVhcjB8Y2hpbGRpZCwgCiAgICAgICAgICAgICAgICAgICAgICAgICBkYXRhID0gZWdtZXJnZWQuY2xlYW4sCiAgICAgICAgICAgICAgICAgICAgICAgICBtZXRob2QgPSAiTUwiKQpzdW1tYXJ5KG1vZGVsLm51bGwuZ3Jvd3RoLmhvbW9nKQpgYGAKCmBgYHtyfQptb2RlbC5udWxsLmdyb3d0aC5oZXQgPC0gbG1lKG1hdGggfiB5ZWFyMCwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcmFuZG9tID0gIH4geWVhcjB8Y2hpbGRpZCwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZGF0YSA9IGVnbWVyZ2VkLmNsZWFuLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtZXRob2QgPSAiTUwiLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICB3ZWlnaHRzID0gdmFySWRlbnQoZm9ybSA9IH4gMSB8IHllYXIwKQopCnN1bW1hcnkobW9kZWwubnVsbC5ncm93dGguaGV0KQpgYGA=

Multilevel Modeling, Module 8, Part 2: Advanced Growth Models

Dr. Broda

Load in Our MVP Packages

Part 1: Loading in the Data, Reshaping Data from Long to Wide, Create New Variables

Load in the Data

Super Handy Trick: Reshaping Data

Using `pivot_wider` to reshape from long to wide….

Using `pivot_longer` to reshape from wide to long….

A little more cleanup needed - change `year` to numeric, get ridding of observations with missing math scores…

A little housekeeping - creating squared and cubic terms for time…

Part 2: Nonlinear Growth Models

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

Calculate GROWTH ICC

Add nonlinear growth terms (square)

Add nonlinear growth terms (cubic)

Add nonlinear growth terms (square AND cubic)

Should we keep that cubic term time (year)? Use `anova` to test…

Same thing, PLUS adding a random slope for year0

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

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

Visualize Our Model by Predicting Scores…

… and Then Plotting Them…

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

Part 2: Heterogeneous Variance Repeated Measures Model

Multilevel Modeling, Module 8, Part 2: Advanced Growth Models

Dr. Broda

Load in Our MVP Packages

Part 1: Loading in the Data, Reshaping Data from Long to Wide, Create New Variables

Load in the Data

Super Handy Trick: Reshaping Data

Using pivot_wider to reshape from long to wide….

Using pivot_longer to reshape from wide to long….

A little more cleanup needed - change year to numeric, get ridding of observations with missing math scores…

A little housekeeping - creating squared and cubic terms for time…

Part 2: Nonlinear Growth Models

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

Calculate GROWTH ICC

Add nonlinear growth terms (square)

Add nonlinear growth terms (cubic)

Add nonlinear growth terms (square AND cubic)

Should we keep that cubic term time (year)? Use anova to test…

Same thing, PLUS adding a random slope for year0

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

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

Visualize Our Model by Predicting Scores…

… and Then Plotting Them…

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

Part 2: Heterogeneous Variance Repeated Measures Model

Using `pivot_wider` to reshape from long to wide….

Using `pivot_longer` to reshape from wide to long….

A little more cleanup needed - change `year` to numeric, get ridding of observations with missing math scores…

Should we keep that cubic term time (year)? Use `anova` to test…

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

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