R Notebook

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages --------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.3.2     v purrr   0.3.4
v tibble  3.0.3     v dplyr   1.0.2
v tidyr   1.1.1     v stringr 1.4.0
v readr   1.3.1     v forcats 0.5.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x tidyr::expand() masks Matrix::expand()
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
x tidyr::pack()   masks Matrix::pack()
x tidyr::unpack() masks Matrix::unpack()

library(lme4)
library(dplyr)

glimpse(STAR_wide)

Rows: 1,471
Columns: 11
$ stdntid    <dbl> 10023, 10042, 10056, 10082, 10097, 10103, 10111, ...
$ race       <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, ...
$ gender     <dbl+lbl> 1, 2, 2, 1, 1, 1, 2, 1, 1, 2, 2, 1, 2, 1, 2, ...
$ ac_mot     <dbl> 50, 50, 52, 49, 51, 50, 51, 57, 40, 48, 53, 28, 5...
$ ever_lunch <dbl+lbl> 2, 1, 2, 2, 1, 2, 1, 2, 2, 2, 1, 1, 2, 1, 1, ...
$ science3   <dbl> 632, 618, 618, 649, 577, 655, 662, 627, 573, 627,...
$ science4   <dbl> 756, 710, 719, 722, 688, 719, 767, 709, 668, 694,...
$ science5   <dbl> 729, 738, 760, 729, 681, 738, 749, 754, 586, 731,...
$ science6   <dbl> 762, 740, 776, 753, 720, 710, 782, 755, 687, 750,...
$ science7   <dbl> 778, 734, 782, 809, 698, 766, 763, 729, 634, 744,...
$ science8   <dbl> 783, 771, 772, 806, 745, 774, 804, 793, 765, 769,...

STAR.long <- STAR_wide %>%
  pivot_longer(.,
              cols = starts_with("science"),
              names_to = "year",
              values_to = "science",
              names_prefix = "science"
              )

glimpse(STAR.long)

Rows: 8,826
Columns: 7
$ stdntid    <dbl> 10023, 10023, 10023, 10023, 10023, 10023, 10042, ...
$ race       <dbl+lbl> 1, 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, 2, ...
$ ac_mot     <dbl> 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 5...
$ ever_lunch <dbl+lbl> 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, ...
$ year       <chr> "3", "4", "5", "6", "7", "8", "3", "4", "5", "6",...
$ science    <dbl> 632, 756, 729, 762, 778, 783, 618, 710, 738, 740,...

STAR.long.clean <- STAR.long %>%
  mutate(.,
         year  = as.numeric(year),
         grade_0 = year - 3,
         grade_sq = grade_0^2,
         grade_cube = grade_0^3)
glimpse(STAR.long.clean)

Rows: 8,826
Columns: 10
$ stdntid    <dbl> 10023, 10023, 10023, 10023, 10023, 10023, 10042, ...
$ race       <dbl+lbl> 1, 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, 2, ...
$ ac_mot     <dbl> 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 5...
$ ever_lunch <dbl+lbl> 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, ...
$ year       <dbl> 3, 4, 5, 6, 7, 8, 3, 4, 5, 6, 7, 8, 3, 4, 5, 6, 7...
$ science    <dbl> 632, 756, 729, 762, 778, 783, 618, 710, 738, 740,...
$ grade_0    <dbl> 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4...
$ grade_sq   <dbl> 0, 1, 4, 9, 16, 25, 0, 1, 4, 9, 16, 25, 0, 1, 4, ...
$ grade_cube <dbl> 0, 1, 8, 27, 64, 125, 0, 1, 8, 27, 64, 125, 0, 1,...

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

Linear mixed model fit by maximum likelihood . t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: science ~ grade_0 + (1 | stdntid)
   Data: STAR.long.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        df t value Pr(>|t|)    
(Intercept)  669.3920     1.0916 2652.5068   613.2   <2e-16 ***
grade_0       24.2729     0.2247 7355.0000   108.0   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
        (Intr)
grade_0 -0.515

model.square <- lmer(science ~ grade_0 + grade_sq + (1|stdntid), REML=FALSE, data = STAR.long.clean)
summary(model.square)

Linear mixed model fit by maximum likelihood . t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: science ~ grade_0 + grade_sq + (1 | stdntid)
   Data: STAR.long.clean

     AIC      BIC   logLik deviance df.resid 
 88570.1  88605.5 -44280.0  88560.1     8821 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.6568 -0.5429 -0.0099  0.5539  4.5442 

Random effects:
 Groups   Name        Variance Std.Dev.
 stdntid  (Intercept) 1132.3   33.65   
 Residual              938.9   30.64   
Number of obs: 8826, groups:  stdntid, 1471

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)  646.2231     1.1376 3069.3181  568.08   <2e-16 ***
grade_0       59.0262     0.6811 7354.9999   86.66   <2e-16 ***
grade_sq      -6.9507     0.1308 7354.9999  -53.16   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
         (Intr) grad_0
grade_0  -0.485       
grade_sq  0.383 -0.960

model.both <- lmer(science ~ grade_0 + grade_sq + grade_cube + (1|stdntid), REML=FALSE, data = STAR.long.clean)
summary(model.both)

Linear mixed model fit by maximum likelihood . t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: science ~ grade_0 + grade_sq + grade_cube + (1 | stdntid)
   Data: STAR.long.clean

     AIC      BIC   logLik deviance df.resid 
 87223.4  87265.9 -43605.7  87211.4     8820 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.2663 -0.4985  0.0027  0.5236  4.4394 

Random effects:
 Groups   Name        Variance Std.Dev.
 stdntid  (Intercept) 1158.5   34.04   
 Residual              781.6   27.96   
Number of obs: 8826, groups:  stdntid, 1471

Fixed effects:
              Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)  635.77161    1.13923 3084.84335  558.07   <2e-16 ***
grade_0      106.75456    1.38750 7355.00518   76.94   <2e-16 ***
grade_sq     -33.07934    0.68954 7355.00521  -47.97   <2e-16 ***
grade_cube     3.48382    0.09055 7355.00481   38.47   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
           (Intr) grad_0 grd_sq
grade_0    -0.411              
grade_sq    0.295 -0.955       
grade_cube -0.238  0.894 -0.985

model.both.random <- lmer(science ~ grade_0 + grade_sq + grade_cube + (grade_0|stdntid), REML=FALSE, data = STAR.long.clean)

boundary (singular) fit: see ?isSingular

summary(model.both.random)

Linear mixed model fit by maximum likelihood . t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: science ~ grade_0 + grade_sq + grade_cube + (grade_0 | stdntid)
   Data: STAR.long.clean

     AIC      BIC   logLik deviance df.resid 
 87223.8  87280.5 -43603.9  87207.8     8818 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.2775 -0.4967  0.0019  0.5228  4.4707 

Random effects:
 Groups   Name        Variance  Std.Dev. Corr
 stdntid  (Intercept) 1100.2324 33.1698      
          grade_0        0.1208  0.3475  1.00
 Residual              781.1644 27.9493      
Number of obs: 8826, groups:  stdntid, 1471

Fixed effects:
              Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)  635.77161    1.12157 2221.90567  566.86   <2e-16 ***
grade_0      106.75456    1.38715 7355.64725   76.96   <2e-16 ***
grade_sq     -33.07934    0.68935 7355.00321  -47.99   <2e-16 ***
grade_cube     3.48382    0.09053 7355.00314   38.48   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
           (Intr) grad_0 grd_sq
grade_0    -0.413              
grade_sq    0.300 -0.955       
grade_cube -0.242  0.894 -0.985
convergence code: 0
boundary (singular) fit: see ?isSingular

library(nlme)

Attaching package: 㤼㸱nlme㤼㸲

The following object is masked from 㤼㸱package:dplyr㤼㸲:

    collapse

The following object is masked from 㤼㸱package:lme4㤼㸲:

    lmList

model.null.growth.homog <- lme(science ~ grade_0, 
                         random = ~ 1|stdntid, 
                         data = STAR.long.clean,
                         method = "ML")
summary(model.null.growth.homog)

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

Random effects:
 Formula: ~1 | stdntid
        (Intercept) Residual
StdDev:    32.74456 36.05019

Fixed effects: science ~ grade_0 

 Correlation: 
        (Intr)
grade_0 -0.515

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-4.02535175 -0.63612595  0.01049814  0.63666093  4.02256978 

Number of Observations: 8826
Number of Groups: 1471 

model.null.growth.het <- lme(science ~ grade_0, 
                             random =  ~ 1|stdntid, 
                             data = STAR.long.clean, 
                             method = "ML", 
                             weights = varIdent(form = ~ 1 | grade_0)
)
summary(model.null.growth.het)

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

Random effects:
 Formula: ~1 | stdntid
        (Intercept) Residual
StdDev:    34.32642 67.40924

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | grade_0 
 Parameter estimates:
        0         1         2         3         4         5 
1.0000000 0.4987552 0.4549923 0.4852018 0.3421510 0.3497103 
Fixed effects: science ~ grade_0 

 Correlation: 
        (Intr)
grade_0 -0.593

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-4.91503937 -0.71613675 -0.05974243  0.51662541  4.30180514 

Number of Observations: 8826
Number of Groups: 1471