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, ...
$ grade      <dbl> 3, 4, 5, 6, 7, 8, 3, 4, 5, 6...
$ race       <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, ...
$ gender     <dbl+lbl> 1, 1, 1, 1, 1, 1, 2, 2, ...
$ ac_mot     <dbl> 50, 50, 50, 50, 50, 50, 50, ...
$ ever_lunch <dbl+lbl> 2, 2, 2, 2, 2, 2, 1, 1, ...
$ science    <dbl> 632, 756, 729, 762, 778, 783...
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,...
$ grade      <dbl> 3, 4, 5, 6, 7, 8, 3, 4, 5, ...
$ race       <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1,...
$ gender     <dbl+lbl> 1, 1, 1, 1, 1, 1, 2, 2,...
$ ac_mot     <dbl> 50, 50, 50, 50, 50, 50, 50,...
$ ever_lunch <dbl+lbl> 2, 2, 2, 2, 2, 2, 1, 1,...
$ science    <dbl> 632, 756, 729, 762, 778, 78...
$ race.fac   <fct> WHITE, WHITE, WHITE, WHITE,...
$ gender.fac <fct> MALE, MALE, MALE, MALE, MAL...
$ lunch.fac  <fct> NON-FREE LUNCH, NON-FREE LU...
$ grade_0    <dbl> 0, 1, 2, 3, 4, 5, 0, 1, 2, ...

#null model

model.null <- lmer(science ~ (1|stdntid), data = starlong.clean)
summary(model.null)
Linear mixed model fit by REML. t-tests use
  Satterthwaite's method [lmerModLmerTest]
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        df t value
(Intercept)  730.0741     0.9363 1469.9997   779.7
            Pr(>|t|)    
(Intercept)   <2e-16 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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), data = starlong.clean)
summary(model.null.growth)
Linear mixed model fit by REML. t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: science ~ grade_0 + (1 | stdntid)
   Data: starlong.clean

REML criterion at convergence: 90950.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0250 -0.6361  0.0105  0.6366  4.0222 

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

Fixed effects:
             Estimate Std. Error        df t value
(Intercept)  669.3920     1.0919 2650.0494     613
grade_0       24.2729     0.2247 7354.0000     108
            Pr(>|t|)    
(Intercept)   <2e-16 ***
grade_0       <2e-16 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
        (Intr)
grade_0 -0.514

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), data = starlong.clean)
summary(model.1)
Linear mixed model fit by REML. t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: 
science ~ grade_0 + gender.fac + race.fac + ac_mot + (1 | stdntid)
   Data: starlong.clean

REML criterion at convergence: 90762.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0869 -0.6308  0.0111  0.6372  3.9901 

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

Fixed effects:
                  Estimate Std. Error        df
(Intercept)       691.3021    10.9427 1471.7467
grade_0            24.2729     0.2247 7354.0000
gender.facFEMALE   -6.4660     1.8286 1464.0000
race.facBLACK     -28.7016     2.3604 1463.9999
race.facASIAN      22.4758    15.3035 1463.9999
race.facHISPANIC   -6.8712    19.7052 1463.9999
race.facOTHER      27.0196    24.1222 1463.9999
ac_mot             -0.2766     0.2250 1464.0015
                 t value Pr(>|t|)    
(Intercept)       63.175  < 2e-16 ***
grade_0          108.021  < 2e-16 ***
gender.facFEMALE  -3.536 0.000419 ***
race.facBLACK    -12.160  < 2e-16 ***
race.facASIAN      1.469 0.142136    
race.facHISPANIC  -0.349 0.727364    
race.facOTHER      1.120 0.262849    
ac_mot            -1.229 0.219235    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) grad_0 g.FEMA r.BLAC r.ASIA
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
ac_mot      -0.991  0.000 -0.217 -0.043  0.049
            r.HISP r.OTHE
grade_0                  
gndr.FEMALE              
rac.fcBLACK              
rac.fcASIAN              
rc.HISPANIC              
rac.fcOTHER  0.002       
ac_mot      -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. t-tests use
  Satterthwaite's method [lmerModLmerTest]
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         df
(Intercept)       6.921e+02  1.289e+01  2.751e+03
grade_0           2.397e+01  2.735e+00  7.353e+03
gender.facFEMALE -6.466e+00  1.829e+00  1.464e+03
race.facBLACK    -2.870e+01  2.360e+00  1.464e+03
race.facASIAN     2.248e+01  1.530e+01  1.464e+03
race.facHISPANIC -6.871e+00  1.971e+01  1.464e+03
race.facOTHER     2.702e+01  2.412e+01  1.464e+03
ac_mot           -2.919e-01  2.642e-01  2.707e+03
grade_0:ac_mot    6.121e-03  5.542e-02  7.353e+03
                 t value Pr(>|t|)    
(Intercept)       53.683  < 2e-16 ***
grade_0            8.763  < 2e-16 ***
gender.facFEMALE  -3.536 0.000419 ***
race.facBLACK    -12.160  < 2e-16 ***
race.facASIAN      1.469 0.142136    
race.facHISPANIC  -0.349 0.727364    
race.facOTHER      1.120 0.262849    
ac_mot            -1.105 0.269466    
grade_0:ac_mot     0.110 0.912047    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) grad_0 g.FEMA r.BLAC r.ASIA
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
ac_mot      -0.993  0.523 -0.185 -0.036  0.042
grad_0:c_mt  0.529 -0.997  0.000  0.000  0.000
            r.HISP r.OTHE ac_mot
grade_0                         
gndr.FEMALE                     
rac.fcBLACK                     
rac.fcASIAN                     
rc.HISPANIC                     
rac.fcOTHER  0.002              
ac_mot      -0.009 -0.015       
grad_0:c_mt  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. t-tests use
  Satterthwaite's method [lmerModLmerTest]
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        df
(Intercept)       691.3784    10.9397 1471.4844
grade_0            24.2729     0.2247 7324.5700
gender.facFEMALE   -6.4023     1.8282 1464.1854
race.facBLACK     -28.6254     2.3599 1464.1854
race.facASIAN      22.2725    15.3006 1464.1854
race.facHISPANIC   -6.5058    19.7015 1464.1854
race.facOTHER      26.3064    24.1177 1464.1854
ac_mot             -0.2790     0.2250 1464.1858
                 t value Pr(>|t|)    
(Intercept)       63.199  < 2e-16 ***
grade_0          108.002  < 2e-16 ***
gender.facFEMALE  -3.502 0.000476 ***
race.facBLACK    -12.130  < 2e-16 ***
race.facASIAN      1.456 0.145699    
race.facHISPANIC  -0.330 0.741280    
race.facOTHER      1.091 0.275561    
ac_mot            -1.240 0.215024    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) grad_0 g.FEMA r.BLAC r.ASIA
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
ac_mot      -0.991  0.000 -0.217 -0.043  0.049
            r.HISP r.OTHE
grade_0                  
gndr.FEMALE              
rac.fcBLACK              
rac.fcASIAN              
rc.HISPANIC              
rac.fcOTHER  0.002       
ac_mot      -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
<none>                           12 -45381 90786
grade_0 in (grade_0 | stdntid)   10 -45381 90782
                                 LRT Df Pr(>Chisq)
<none>                                            
grade_0 in (grade_0 | stdntid) 0.607  2     0.7382
LS0tDQp0aXRsZTogJ01vZHVsZSA4LCBQYXJ0IDE6IEludHJvIHRvIEdyb3d0aCBNb2RlbHMnDQphdXRob3I6ICdKYWtlIFJleW5vbGRzIC0gT2N0b2JlciAxMiwgMjAyMCcNCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQoNCiMgTG9hZCBpbiBPdXIgTVZQIFBhY2thZ2VzDQpgYGB7cn0NCnN1cHByZXNzUGFja2FnZVN0YXJ0dXBNZXNzYWdlcyhsaWJyYXJ5KHRpZHl2ZXJzZSkpDQpzdXBwcmVzc1BhY2thZ2VTdGFydHVwTWVzc2FnZXMobGlicmFyeShIbWlzYykpDQpzdXBwcmVzc1BhY2thZ2VTdGFydHVwTWVzc2FnZXMobGlicmFyeShsbWU0KSkNCmBgYA0KDQojIExvYWQgaW4gdGhlIERhdGENCmBgYHtyfQ0Kc3RhcmxvbmcgPC0gaGF2ZW46OnJlYWRfZHRhKCJTVEFSX2xvbmcuZHRhIikNCg0KZ2xpbXBzZShzdGFybG9uZykNCmBgYA0KDQpgYGB7cn0NCnN0YXJsb25nLmNsZWFuIDwtIHN0YXJsb25nICU+JQ0KICBtdXRhdGUoLiwNCiAgICAgICAgIHJhY2UuZmFjID0gYXNfZmFjdG9yKHJhY2UpLA0KICAgICAgICAgZ2VuZGVyLmZhYyA9IGFzX2ZhY3RvcihnZW5kZXIpLA0KICAgICAgICAgbHVuY2guZmFjID0gYXNfZmFjdG9yKGV2ZXJfbHVuY2gpLA0KICAgICAgICAgZ3JhZGVfMCA9IGdyYWRlIC0gMykNCg0KZ2xpbXBzZShzdGFybG9uZy5jbGVhbikNCmBgYA0KDQojbnVsbCBtb2RlbCANCmBgYHtyfQ0KbW9kZWwubnVsbCA8LSBsbWVyKHNjaWVuY2UgfiAoMXxzdGRudGlkKSwgZGF0YSA9IHN0YXJsb25nLmNsZWFuKQ0Kc3VtbWFyeShtb2RlbC5udWxsKQ0KYGBgDQoNCiMjIENhbGN1bGF0ZSBJQ0MNCmBgYHtyfQ0KbnVsbC5pY2MgPC0gNzI5LjQvKDcyOS40ICsgMzM2MS43KQ0KbnVsbC5pY2MNCmBgYA0KDQojIEdST1dUSCBudWxsIG1vZGVsIA0KYGBge3J9DQptb2RlbC5udWxsLmdyb3d0aCA8LSBsbWVyKHNjaWVuY2UgfiBncmFkZV8wICsgKDF8c3RkbnRpZCksIGRhdGEgPSBzdGFybG9uZy5jbGVhbikNCnN1bW1hcnkobW9kZWwubnVsbC5ncm93dGgpDQpgYGANCg0KIyMgQ2FsY3VsYXRlIEdST1dUSCBJQ0MNCmBgYHtyfQ0KbnVsbC5ncm93dGguaWNjIDwtIDY2OS4zOTIwLyg2NjkuMzkyMCArIDI0LjI3MjkpDQpudWxsLmdyb3d0aC5pY2MNCmBgYA0KDQojIEFkZCBhIHN0dWRlbnQtbGV2ZWwgcHJlZGljdG9ycw0KYGBge3J9DQptb2RlbC4xIDwtIGxtZXIoc2NpZW5jZSB+IGdyYWRlXzAgKyBnZW5kZXIuZmFjICsgcmFjZS5mYWMgKyBhY19tb3QgKyAoMXxzdGRudGlkKSwgZGF0YSA9IHN0YXJsb25nLmNsZWFuKQ0Kc3VtbWFyeShtb2RlbC4xKQ0KYGBgDQojIFRoaXMgaXMgYW4gaW50ZXJhY3Rpb24gZWZmZWN0Lg0KYGBge3J9DQptb2RlbC4yIDwtIGxtZXIoc2NpZW5jZSB+IGdyYWRlXzAgKyBnZW5kZXIuZmFjICsgcmFjZS5mYWMgKyBhY19tb3QgKyArYWNfbW90OmdyYWRlXzAgKyAoMXxzdGRudGlkKSwgZGF0YSA9IHN0YXJsb25nLmNsZWFuKQ0Kc3VtbWFyeShtb2RlbC4yKQ0KYGBgDQoNCiNyYW5kb20gc2xvcGUNCmBgYHtyfQ0KbW9kZWwuMyA8LSBsbWVyKHNjaWVuY2UgfiBncmFkZV8wICsgZ2VuZGVyLmZhYyArIHJhY2UuZmFjICsgYWNfbW90ICsgKGdyYWRlXzB8c3RkbnRpZCksIGRhdGEgPSBzdGFybG9uZy5jbGVhbikNCnN1bW1hcnkobW9kZWwuMykNCmBgYA0KDQojIFNob3VsZCB3ZSBrZWVwIHRoYXQgcmFuZG9tIHNsb3BlIGZvciB0aW1lICh5ZWFyKT8gDQpgYGB7cn0NCmxtZXJUZXN0OjpyYW5kKG1vZGVsLjMpDQpgYGANCg==