Module 10

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

glimpse(projectSTAR_1_)

Rows: 6,325
Columns: 27
$ stdntid       <dbl> 10001, 10133, 10246, 10263, 10266, 10275, 1028...
$ race          <dbl+lbl> 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, ...
$ gender        <dbl+lbl> 1, 1, 2, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, ...
$ FLAGSGK       <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
$ flaggk        <dbl+lbl> 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
$ gkclasstype   <dbl+lbl> 3, 3, 3, 1, 2, 3, 1, 3, 1, 2, 3, 3, 2, 3, ...
$ gkschid       <dbl> 169229, 169280, 218562, 205492, 257899, 161176...
$ gksurban      <dbl+lbl> 2, 2, 4, 2, 3, 3, 2, 2, 3, 3, 3, 2, 2, 2, ...
$ gktchid       <dbl> 16922904, 16928003, 21856202, 20549204, 257899...
$ gktgen        <dbl+lbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, ...
$ gktrace       <dbl+lbl> 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, ...
$ gkthighdegree <dbl+lbl> 2, 2, 3, 2, 3, 2, 3, 3, 3, 3, 2, 2, 3, 2, ...
$ gktcareer     <dbl+lbl>  4,  4,  4,  4,  4, NA,  4,  4,  4,  4,  4...
$ gktyears      <dbl> 5, 7, 8, 3, 12, 2, 7, 14, 4, 6, 11, 16, 12, 5,...
$ gkclasssize   <dbl> 24, 22, 17, 17, 24, 24, 13, 24, 14, 23, 23, 22...
$ gkfreelunch   <dbl+lbl> 2, 2, 2, 1, 2, 2, 2, 1, 1, 2, 2, 2, 2, 2, ...
$ gkrepeat      <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
$ gkspeced      <dbl+lbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, ...
$ gkspecin      <dbl+lbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, ...
$ gkpresent     <dbl> 161, 175, 132, 178, 170, 94, 160, 154, 172, 95...
$ gkabsent      <dbl> 19, 5, 28, 2, 10, 3, 2, 7, 8, 2, 17, 8, 0, 31,...
$ gktreadss     <dbl> NA, 427, 450, 483, 456, 411, 443, 448, 463, 47...
$ gktmathss     <dbl> NA, 478, 494, 513, 513, 468, 473, 449, 520, 53...
$ gktlistss     <dbl> NA, 509, 549, 554, 520, 571, 595, 540, 565, 59...
$ gkwordskillss <dbl> NA, 418, 444, 431, 468, 396, 444, 444, 480, 48...
$ gkmotivraw    <dbl> 23, 24, 28, 27, 25, 24, NA, NA, 26, 27, 24, 24...
$ gkselfconcraw <dbl> 52, 53, 56, 61, 54, 55, NA, NA, 52, 61, 55, 49...

model.final <- lmer(gktmathss ~ gkselfconcraw + gkclasstype + gkthighdegree + (1|gktchid), REML=FALSE, data = projectSTAR_1_)
summary(model.final)

Linear mixed model fit by maximum likelihood . t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: gktmathss ~ gkselfconcraw + gkclasstype + gkthighdegree + (1 |  
    gktchid)
   Data: projectSTAR_1_

     AIC      BIC   logLik deviance df.resid 
 49055.9  49094.8 -24522.0  49043.9     4750 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8391 -0.6386 -0.0618  0.5814  3.8970 

Random effects:
 Groups   Name        Variance Std.Dev.
 gktchid  (Intercept)  663.4   25.76   
 Residual             1550.2   39.37   
Number of obs: 4756, groups:  gktchid, 301

Fixed effects:
               Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)    464.8157    10.1601  888.2763  45.749  < 2e-16 ***
gkselfconcraw    0.4458     0.1196 4639.9291   3.726 0.000197 ***
gkclasstype     -4.6260     1.9336  300.9683  -2.392 0.017349 *  
gkthighdegree    2.7561     2.7143  304.9732   1.015 0.310722    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) gkslfc gkclss
gkselfcncrw -0.660              
gkclasstype -0.377  0.017       
gkthighdegr -0.630 -0.010 -0.002

library(broom.mixed)

diagnostics <- augment(model.final)

ggplot(data = diagnostics, mapping = aes(x = .resid)) +
  geom_histogram(binwidth = .25) + theme_classic() + 
  labs(title = "Histogram of Residuals for Math Scores Model",
                      x = "Residual Value") +
  geom_vline(xintercept = c(-2.5, 2.5), linetype = "dotted")

shapiro.test(diagnostics$.resid)

    Shapiro-Wilk normality test

data:  diagnostics$.resid
W = 0.99119, p-value < 2.2e-16

ggplot(data = diagnostics, mapping = aes(x = .fitted, y = .resid)) +
  geom_point() + labs(title = "RVF Plot for Math Scores Model",
                      x = "Predicted Value, Math Scores",
                      y = "Residual Value") + theme_classic()

ggplot(data = diagnostics, mapping = aes(x = .fitted, y = .cooksd, label = gktchid)) +
  geom_point() + geom_text(nudge_x = .25) + theme_classic() + 
  labs(title = "Cook's Distance Plot for Math Scores Model",
                      x = "Predicted Value, Math Scores",
                      y = "Cook's Distance") + 
  geom_hline(yintercept = 4/6325, linetype = "dotted")

library(robustlmm)
model.robust <- lmer(gktmathss ~ gkselfconcraw + gkclasstype + gkthighdegree + (1|gktchid), REML=FALSE, data = projectSTAR_1_)
summary(model.robust)

Linear mixed model fit by maximum likelihood . t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: gktmathss ~ gkselfconcraw + gkclasstype + gkthighdegree + (1 |  
    gktchid)
   Data: projectSTAR_1_

     AIC      BIC   logLik deviance df.resid 
 49055.9  49094.8 -24522.0  49043.9     4750 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8391 -0.6386 -0.0618  0.5814  3.8970 

Random effects:
 Groups   Name        Variance Std.Dev.
 gktchid  (Intercept)  663.4   25.76   
 Residual             1550.2   39.37   
Number of obs: 4756, groups:  gktchid, 301

Fixed effects:
               Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)    464.8157    10.1601  888.2763  45.749  < 2e-16 ***
gkselfconcraw    0.4458     0.1196 4639.9291   3.726 0.000197 ***
gkclasstype     -4.6260     1.9336  300.9683  -2.392 0.017349 *  
gkthighdegree    2.7561     2.7143  304.9732   1.015 0.310722    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) gkslfc gkclss
gkselfcncrw -0.660              
gkclasstype -0.377  0.017       
gkthighdegr -0.630 -0.010 -0.002

STAR.trimmed <- diagnostics %>%
  filter(., .cooksd < .10)
model.trimmed <- lmer(gktmathss ~ gkselfconcraw + gkclasstype + gkthighdegree + (1|gktchid), REML=FALSE, data = STAR.trimmed)
summary(model.trimmed)

Linear mixed model fit by maximum likelihood . t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: gktmathss ~ gkselfconcraw + gkclasstype + gkthighdegree + (1 |  
    gktchid)
   Data: STAR.trimmed

     AIC      BIC   logLik deviance df.resid 
 47346.0  47384.7 -23667.0  47334.0     4657 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.05922 -0.65534 -0.04224  0.61962  3.08934 

Random effects:
 Groups   Name        Variance Std.Dev.
 gktchid  (Intercept)  704.3   26.54   
 Residual             1300.7   36.06   
Number of obs: 4663, groups:  gktchid, 301

Fixed effects:
               Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)    461.5853     9.9617  782.0993  46.336  < 2e-16 ***
gkselfconcraw    0.4567     0.1113 4516.5362   4.102 4.16e-05 ***
gkclasstype     -4.0598     1.9670  300.5117  -2.064   0.0399 *  
gkthighdegree    2.8253     2.7594  303.7239   1.024   0.3067    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) gkslfc gkclss
gkselfcncrw -0.627              
gkclasstype -0.390  0.017       
gkthighdegr -0.655 -0.009 -0.002