R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

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 dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(lme4)

Loading required package: Matrix

Attaching package: 㤼㸱Matrix㤼㸲

The following objects are masked from 㤼㸱package:tidyr㤼㸲:

    expand, pack, unpack

library(psych)

Attaching package: 㤼㸱psych㤼㸲

The following objects are masked from 㤼㸱package:ggplot2㤼㸲:

    %+%, alpha

hsbmerged <- read_dta("EDUS 651/Module 3/hsbmerged.dta")
glimpse(hsbmerged)

Rows: 7,185
Columns: 12
$ schoolid <chr> "1224", "1224", "1224", "1224", "1224", "1224", "1224"...
$ minority <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, ...
$ female   <dbl> 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, ...
$ ses      <dbl> -1.528, -0.588, -0.528, -0.668, -0.158, 0.022, -0.618,...
$ mathach  <dbl> 5.876, 19.708, 20.349, 8.781, 17.898, 4.583, -2.832, 0...
$ size     <dbl> 842, 842, 842, 842, 842, 842, 842, 842, 842, 842, 842,...
$ sector   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
$ pracad   <dbl> 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, ...
$ disclim  <dbl> 1.597, 1.597, 1.597, 1.597, 1.597, 1.597, 1.597, 1.597...
$ himinty  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
$ meanses  <dbl> -0.428, -0.428, -0.428, -0.428, -0.428, -0.428, -0.428...
$ readach  <dbl> 51.96463, 39.28524, 43.81590, 52.01301, 38.99202, 39.0...

hsbmerged.1 <- hsbmerged %>%
  mutate(., minority.factor = as_factor(minority),
         female.factor = as_factor(female),
         sector.factor = as_factor(sector))
glimpse(hsbmerged.1)

Rows: 7,185
Columns: 15
$ schoolid        <chr> "1224", "1224", "1224", "1224", "1224", "1224",...
$ minority        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ female          <dbl> 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1,...
$ ses             <dbl> -1.528, -0.588, -0.528, -0.668, -0.158, 0.022, ...
$ mathach         <dbl> 5.876, 19.708, 20.349, 8.781, 17.898, 4.583, -2...
$ size            <dbl> 842, 842, 842, 842, 842, 842, 842, 842, 842, 84...
$ sector          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ pracad          <dbl> 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35,...
$ disclim         <dbl> 1.597, 1.597, 1.597, 1.597, 1.597, 1.597, 1.597...
$ himinty         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ meanses         <dbl> -0.428, -0.428, -0.428, -0.428, -0.428, -0.428,...
$ readach         <dbl> 51.96463, 39.28524, 43.81590, 52.01301, 38.9920...
$ minority.factor <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ female.factor   <fct> 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1,...
$ sector.factor   <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...

1. Run a null model with mathach as the DV, and then run a conditional random intercept model with mathach as the DV and female, ses, size, and sector as IVs.

model.1 <- lmer(mathach ~ (1|schoolid), REML = FALSE, data = hsbmerged.1)
summary(model.1)

Linear mixed model fit by maximum
  likelihood [lmerMod]
Formula: mathach ~ (1 | schoolid)
   Data: hsbmerged.1

     AIC      BIC   logLik deviance 
 47121.8  47142.4 -23557.9  47115.8 
df.resid 
    7182 

Scaled residuals: 
     Min       1Q   Median       3Q 
-3.06262 -0.75365  0.02676  0.76070 
     Max 
 2.74184 

Random effects:
 Groups   Name        Variance Std.Dev.
 schoolid (Intercept)  8.553   2.925   
 Residual             39.148   6.257   
Number of obs: 7185, groups:  schoolid, 160

Fixed effects:
            Estimate Std. Error t value
(Intercept)  12.6371     0.2436   51.87

null.ICC <- 8.553/(8.553 + 39.148)

The ICC is .18 which means 18% of the variance in math achievement occurs between schools. We should run MLM.

model.2 <- lmer(mathach ~ female.factor +  ses + size + sector.factor + (1|schoolid), REML = FALSE, data = hsbmerged.1)
summary(model.2)

Linear mixed model fit by maximum
  likelihood [lmerMod]
Formula: 
mathach ~ female.factor + ses + size + sector.factor + (1 | schoolid)
   Data: hsbmerged.1

     AIC      BIC   logLik deviance 
 46564.2  46612.3 -23275.1  46550.2 
df.resid 
    7178 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2032 -0.7404  0.0323  0.7584  2.8130 

Random effects:
 Groups   Name        Variance Std.Dev.
 schoolid (Intercept)  3.27    1.808   
 Residual             36.81    6.067   
Number of obs: 7185, groups:  schoolid, 160

Fixed effects:
                 Estimate Std. Error
(Intercept)    11.6791431  0.4543239
female.factor1 -1.2010906  0.1641199
ses             2.3451770  0.1050492
size            0.0004972  0.0002894
sector.factor1  2.3777505  0.3633391
               t value
(Intercept)     25.707
female.factor1  -7.318
ses             22.325
size             1.718
sector.factor1   6.544

Correlation of Fixed Effects:
            (Intr) fml.f1 ses    size  
femal.fctr1 -0.202                     
ses          0.028  0.057              
size        -0.858  0.016 -0.008       
sectr.fctr1 -0.671  0.006 -0.089  0.447

2. Find the appropriate variance estimates and calculate Multilevel R Squared, BG-PRV, and WG-PRV for this comparison of two models. Interpret the effect sizes.

Using the speadsheet: The multilevel R Squared is .16, which is a small effect. The target model reduced the within-group variance by 5.97% and the between-group variance by 61.77%. This model is much better at reducing between group variance than the null model.

3. Now, using the target model from Part 1 (mathach as the DV and female, ses, size, and sector as IVs) as your starting model, test an additional model where you allow the slope for ses to vary randomly by school. Conduct a likelihood ratio test comparing these two models.

model.3 <- lmer(mathach ~ female.factor +  ses + size + sector.factor + (ses|schoolid), REML = FALSE, data = hsbmerged.1)
summary(model.3)

Linear mixed model fit by maximum
  likelihood [lmerMod]
Formula: 
mathach ~ female.factor + ses + size + sector.factor + (ses |  
    schoolid)
   Data: hsbmerged.1

     AIC      BIC   logLik deviance 
 46560.0  46621.9 -23271.0  46542.0 
df.resid 
    7176 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2053 -0.7384  0.0281  0.7570  2.8322 

Random effects:
 Groups   Name        Variance Std.Dev.
 schoolid (Intercept)  3.5435  1.8824  
          ses          0.3548  0.5956  
 Residual             36.6025  6.0500  
 Corr
     
 0.60
     
Number of obs: 7185, groups:  schoolid, 160

Fixed effects:
                 Estimate Std. Error
(Intercept)    11.4758217  0.4505476
female.factor1 -1.1940658  0.1641275
ses             2.3527077  0.1151975
size            0.0004656  0.0002850
sector.factor1  2.7891409  0.3651071
               t value
(Intercept)     25.471
female.factor1  -7.275
ses             20.423
size             1.634
sector.factor1   7.639

Correlation of Fixed Effects:
            (Intr) fml.f1 ses    size  
femal.fctr1 -0.207                     
ses          0.109  0.055              
size        -0.850  0.019 -0.006       
sectr.fctr1 -0.660  0.009 -0.078  0.436

lmerTest::rand(model.3)

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

Model:
mathach ~ female.factor + ses + size + sector.factor + (ses | 
    schoolid)
                        npar logLik   AIC    LRT Df Pr(>Chisq)  
<none>                     9 -23271 46560                       
ses in (ses | schoolid)    7 -23275 46564 8.1673  2    0.01685 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

When we add the random slope for ses the AIC goes down and this is significant (p<.05). This means that the model with random slope for ses is a better fit than.

4. Compare the results of the LR test with the AIC and BIC values for the models you ran in #3. Is the evidence consistent as to which model to use? Why or why not? The AIC values for the models is consistent with the LR test. However, the BIC values are not consistent. When we add the random slope for ses the BIC actually gets bigger meaning the model is not a better fit than without the random slope for ses. BIC adjusts for the number of observations while AIC does not. This could be because BIC favors more simple models