Module 8, Part 1: Intro to Growth Models

1. Create a new variable, grade_0, which goes from 0-5 instead of 3-8 like our current grade variable does. Why is using a time variable that starts at 0 helpful in interpreting our results?

It is helpful to use a time variable that starts with 0 because it is easier to interpret the intercept at 0, than 3. As stated in our text, “the time variable may be grand-mean centered, which means at 0, so that the intercept is the value of the outcome half-way through the time period under study” (p. 39). If you scale the year to 0, it can be controled at the starting point.

2. Run and interpret a null model, AND a null growth model, with science scores science as the DV, with scores clustered within students (stdntid). How much variation in science scores is at the student level for each null model? Why does this value change when moving from a regular null model to a null growth model?

The average science score in the entire sample is 730.The random variance component at the student-level is 729.4. The random variance component at the observation level is 3361.7. Most of the variability is within student, rather than between student.17.8% of the variability is at the student-level. However, the element of time is not accounted for. When we run a growth model, we add a predictor of grade_0 to the equation. The student-level variability in the growth model is 1073the residual variability went down to 1300. The intercept also changed to 669. The intercept is the predicted science score at grade_0 (grade 3). The growth ICC increased significantly because we introduced a growth parameter that explained a lot of variability within the student.

Part 2: Adding Predictors, Interactions, and Testing Random Slopes

3. Now, run a random intercept growth model with science scores (science) as the DV, and grade (grade0), gender (gender), race/ethnicity (race), and academic motivation (ac_mot) as predictors. Interpret the results and evaluate model fit using AIC/BIC and your choice of effect size.

The AIC of the Growth Null Model is 90959.4, BIC is 90987.7. The AIC and BIC of the new model is 90811.3 and 90882.1 respectively. This decrease in AIC and BIC indicate that the model is improving, compared to the null model. Grade, gender, and black all appear to be significant predictors. This is not talking about slope, but rather starting point.

4. Use an interaction effect between ac_mot and grade0 to test whether students with different levels of academic motivation have different growth in science scores. Interpret the results.

The interaction between grade_0 and ac_mot does not show to be significant. It does not appear to impact the growth measures in change over time.

5. Try adding a random slope for time (grade) at the student level. What do these random slopes “mean”? Use lrtest (Stata) or the rand function (R) to evaluate whether the slope should be treated as random.

Random slopes allow variability for each student. That is, the model accounts for variation of the slope (grade) over time. However, for this model, when you include the random slope, it does not improve the model significantly. The Chi-squared statistic is 0.7382.

Load in Our MVP Packages

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

Load in the Data

starlong <- haven::read_dta("STAR_long.dta")

glimpse(starlong)

Rows: 8,826
Columns: 7
$ stdntid    <dbl> 10023, 10023, 10023, 10023, 10023, 10023, 1004...
$ grade      <dbl> 3, 4, 5, 6, 7, 8, 3, 4, 5, 6, 7, 8, 3, 4, 5, 6...
$ race       <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
$ gender     <dbl+lbl> 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, ...
$ ac_mot     <dbl> 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50...
$ ever_lunch <dbl+lbl> 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, ...
$ science    <dbl> 632, 756, 729, 762, 778, 783, 618, 710, 738, 7...

starlong.clean <- starlong %>%
  mutate(.,
         race.fac = as_factor(race),
         gender.fac = as_factor(gender),
         lunch.fac = as_factor(ever_lunch),
         grade_0 = grade - 3)

glimpse(starlong.clean)

Rows: 8,826
Columns: 11
$ stdntid    <dbl> 10023, 10023, 10023, 10023, 10023, 10023, 1004...
$ grade      <dbl> 3, 4, 5, 6, 7, 8, 3, 4, 5, 6, 7, 8, 3, 4, 5, 6...
$ race       <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
$ gender     <dbl+lbl> 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, ...
$ ac_mot     <dbl> 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50...
$ ever_lunch <dbl+lbl> 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, ...
$ science    <dbl> 632, 756, 729, 762, 778, 783, 618, 710, 738, 7...
$ race.fac   <fct> WHITE, WHITE, WHITE, WHITE, WHITE, WHITE, WHIT...
$ gender.fac <fct> MALE, MALE, MALE, MALE, MALE, MALE, FEMALE, FE...
$ lunch.fac  <fct> NON-FREE LUNCH, NON-FREE LUNCH, NON-FREE LUNCH...
$ grade_0    <dbl> 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3...

#null model

model.null <- lmer(science ~ (1|stdntid), data = starlong.clean)
summary(model.null)

Linear mixed model fit by REML ['lmerMod']
Formula: science ~ (1 | stdntid)
   Data: starlong.clean

REML criterion at convergence: 97939.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1898 -0.3885  0.2381  0.6370  2.7353 

Random effects:
 Groups   Name        Variance Std.Dev.
 stdntid  (Intercept)  729.4   27.01   
 Residual             3361.7   57.98   
Number of obs: 8826, groups:  stdntid, 1471

Fixed effects:
            Estimate Std. Error t value
(Intercept) 730.0741     0.9363   779.7

Calculate ICC

null.icc <- 729.4/(729.4 + 3361.7)
null.icc

[1] 0.1782895

GROWTH null model

model.null.growth <- lmer(science ~ grade_0 + (1|stdntid), REML = FALSE, data = starlong.clean)
summary(model.null.growth)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: science ~ grade_0 + (1 | stdntid)
   Data: starlong.clean

     AIC      BIC   logLik deviance df.resid 
 90959.4  90987.7 -45475.7  90951.4     8822 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0254 -0.6361  0.0105  0.6367  4.0226 

Random effects:
 Groups   Name        Variance Std.Dev.
 stdntid  (Intercept) 1072     32.74   
 Residual             1300     36.05   
Number of obs: 8826, groups:  stdntid, 1471

Fixed effects:
            Estimate Std. Error t value
(Intercept) 669.3920     1.0916   613.2
grade_0      24.2729     0.2247   108.0

Correlation of Fixed Effects:
        (Intr)
grade_0 -0.515

Calculate GROWTH ICC

null.growth.icc <- 669.3920/(669.3920 + 24.2729)
null.growth.icc

[1] 0.9650077

Add a student-level predictors

model.1 <- lmer(science ~ grade_0 + gender.fac + race.fac + ac_mot + (1|stdntid), REML=FALSE, data = starlong.clean)
summary(model.1)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: 
science ~ grade_0 + gender.fac + race.fac + ac_mot + (1 | stdntid)
   Data: starlong.clean

     AIC      BIC   logLik deviance df.resid 
 90811.3  90882.1 -45395.6  90791.3     8816 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0880 -0.6310  0.0112  0.6372  3.9908 

Random effects:
 Groups   Name        Variance Std.Dev.
 stdntid  (Intercept)  939.3   30.65   
 Residual             1299.6   36.05   
Number of obs: 8826, groups:  stdntid, 1471

Fixed effects:
                 Estimate Std. Error t value
(Intercept)      691.3021    10.9167  63.325
grade_0           24.2729     0.2247 108.029
gender.facFEMALE  -6.4660     1.8242  -3.545
race.facBLACK    -28.7016     2.3547 -12.189
race.facASIAN     22.4758    15.2670   1.472
race.facHISPANIC  -6.8712    19.6583  -0.350
race.facOTHER     27.0196    24.0648   1.123
ac_mot            -0.2766     0.2245  -1.232

Correlation of Fixed Effects:
            (Intr) grad_0 g.FEMA r.BLAC r.ASIA r.HISP r.OTHE
grade_0     -0.051                                          
gndr.FEMALE  0.128  0.000                                   
rac.fcBLACK  0.002  0.000  0.049                            
rac.fcASIAN -0.058  0.000  0.029  0.026                     
rc.HISPANIC  0.004  0.000  0.021  0.022  0.003              
rac.fcOTHER  0.014  0.000  0.007  0.018  0.002  0.002       
ac_mot      -0.991  0.000 -0.217 -0.043  0.049 -0.010 -0.018

This is an interaction effect.

model.2 <- lmer(science ~ grade_0 + gender.fac + race.fac + ac_mot + +ac_mot:grade_0 + (1|stdntid), data = starlong.clean)
summary(model.2)

Linear mixed model fit by REML ['lmerMod']
Formula: 
science ~ grade_0 + gender.fac + race.fac + ac_mot + +ac_mot:grade_0 +  
    (1 | stdntid)
   Data: starlong.clean

REML criterion at convergence: 90766.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0845 -0.6315  0.0109  0.6369  3.9896 

Random effects:
 Groups   Name        Variance Std.Dev.
 stdntid  (Intercept)  944.8   30.74   
 Residual             1300.0   36.06   
Number of obs: 8826, groups:  stdntid, 1471

Fixed effects:
                   Estimate Std. Error t value
(Intercept)      692.054985  12.891599  53.683
grade_0           23.971719   2.735455   8.763
gender.facFEMALE  -6.465992   1.828568  -3.536
race.facBLACK    -28.701558   2.360372 -12.160
race.facASIAN     22.475804  15.303458   1.469
race.facHISPANIC  -6.871236  19.705237  -0.349
race.facOTHER     27.019612  24.122216   1.120
ac_mot            -0.291857   0.264239  -1.105
grade_0:ac_mot     0.006121   0.055418   0.110

Correlation of Fixed Effects:
            (Intr) grad_0 g.FEMA r.BLAC r.ASIA r.HISP r.OTHE ac_mot
grade_0     -0.530                                                 
gndr.FEMALE  0.108  0.000                                          
rac.fcBLACK  0.001  0.000  0.049                                   
rac.fcASIAN -0.049  0.000  0.029  0.026                            
rc.HISPANIC  0.003  0.000  0.021  0.022  0.003                     
rac.fcOTHER  0.012  0.000  0.007  0.018  0.002  0.002              
ac_mot      -0.993  0.523 -0.185 -0.036  0.042 -0.009 -0.015       
grad_0:c_mt  0.529 -0.997  0.000  0.000  0.000  0.000  0.000 -0.524

#random slope

model.3 <- lmer(science ~ grade_0 + gender.fac + race.fac + ac_mot + (grade_0|stdntid), data = starlong.clean)

boundary (singular) fit: see ?isSingular

summary(model.3)

Linear mixed model fit by REML ['lmerMod']
Formula: 
science ~ grade_0 + gender.fac + race.fac + ac_mot + (grade_0 |  
    stdntid)
   Data: starlong.clean

REML criterion at convergence: 90761.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0652 -0.6306  0.0104  0.6359  4.0013 

Random effects:
 Groups   Name        Variance  Std.Dev. Corr
 stdntid  (Intercept) 9.166e+02 30.2751      
          grade_0     3.431e-02  0.1852  1.00
 Residual             1.300e+03 36.0508      
Number of obs: 8826, groups:  stdntid, 1471

Fixed effects:
                 Estimate Std. Error t value
(Intercept)      691.3784    10.9397  63.199
grade_0           24.2729     0.2247 108.002
gender.facFEMALE  -6.4023     1.8282  -3.502
race.facBLACK    -28.6254     2.3599 -12.130
race.facASIAN     22.2725    15.3006   1.456
race.facHISPANIC  -6.5058    19.7015  -0.330
race.facOTHER     26.3064    24.1177   1.091
ac_mot            -0.2790     0.2250  -1.240

Correlation of Fixed Effects:
            (Intr) grad_0 g.FEMA r.BLAC r.ASIA r.HISP r.OTHE
grade_0     -0.050                                          
gndr.FEMALE  0.128  0.000                                   
rac.fcBLACK  0.002  0.000  0.049                            
rac.fcASIAN -0.058  0.000  0.029  0.026                     
rc.HISPANIC  0.004  0.000  0.021  0.022  0.003              
rac.fcOTHER  0.014  0.000  0.007  0.018  0.002  0.002       
ac_mot      -0.991  0.000 -0.217 -0.043  0.049 -0.010 -0.018
convergence code: 0
boundary (singular) fit: see ?isSingular

Should we keep that random slope for time (year)?

lmerTest::rand(model.3)

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

Model:
science ~ grade_0 + gender.fac + race.fac + ac_mot + (grade_0 | 
    stdntid)
                               npar logLik   AIC   LRT Df Pr(>Chisq)
<none>                           12 -45381 90786                    
grade_0 in (grade_0 | stdntid)   10 -45381 90782 0.607  2     0.7382