Part 1: Calculating Multilevel Effect Sizes

  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.

Null Model AIC: 47121.8 BIC: 47142.4

Model.1 with female, ses, size, and sector AIC: 46564.2 BIC: 46612.3

Difference between the null model and Model.1 AIC and BIC AIC: -557.6 BIC: -530.1

  1. 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.

The R-squared value is 0.16. (Not bad) R-squared explains the predictive power of our MLM. WG-PRV is 5.97% BG-PRV is 61.77%. Target model reduced the within-group variance by 6% and the between-group variance by 62%. The target model is much better at reducing between group variance.

Part 2: Using Goodness-of-Fit Measures

  1. 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 without random slope. AIC: 46564.2 BIC: 46612.3

Model with random slope. AIC: 46560.0
BIC: 46621.9

  1. 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?

I first started with looking at the AIC and BIC. However, they provided different results. AIC went down by 4.2 points. BIC went up by 9.6 points. So, the evidence is not consistent. I then ran the likelihood ratio test for both the model without the random slope and with the random slope. The results indicated that we should not include the random slope, ses, in the model.

#Load Packages & Data

library(tidyverse)
library(lme4)
library(psych)
library(haven)
library(tibble)

hsbme <- haven::read_dta("hsbmergedd.dta")
glimpse(hsbme)
Rows: 7,185
Columns: 12
$ schoolid <chr> "1224", "1224", "1224", "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, 0, 0, 0, 0, 0, 0, 1, 0, ...
$ female   <dbl> 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, ...
$ ses      <dbl> -1.528, -0.588, -0.528, -0.668, -0.158, 0.022, -0.618, -0.998, -0.888, -0.458,...
$ mathach  <dbl> 5.876, 19.708, 20.349, 8.781, 17.898, 4.583, -2.832, 0.523, 1.527, 21.521, 9.4...
$ size     <dbl> 842, 842, 842, 842, 842, 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, 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, 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, 1.597, 1.597, 1.597, 1...
$ himinty  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 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, -0.428, -0.428, -0.428...
$ readach  <dbl> 51.96463, 39.28524, 43.81590, 52.01301, 38.99202, 39.04975, 54.93985, 58.78317...

#Clean Data

hsbme.clean <- hsbme %>%
  mutate(.,
         sector.factor = as.factor(sector), 
         female.factor = as.factor(female),
         minority.factor = as.factor(minority))

levels(hsb.clean$sector.factor)=c("Public", "Parochial")
levels(hsb.clean$female.factor)=c("Male", "Female")
levels(hsb.clean$minority.factor)=c("Non-Minority", "Minority")

str(hsb.clean$sector.factor)
 Factor w/ 2 levels "Public","Parochial": 1 1 1 1 1 1 1 1 1 1 ...
glimpse(hsbme.clean)  
Rows: 7,185
Columns: 15
$ schoolid        <chr> "1224", "1224", "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, 0, 0, 0, 0, 0, 0,...
$ female          <dbl> 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0,...
$ ses             <dbl> -1.528, -0.588, -0.528, -0.668, -0.158, 0.022, -0.618, -0.998, -0.888, ...
$ mathach         <dbl> 5.876, 19.708, 20.349, 8.781, 17.898, 4.583, -2.832, 0.523, 1.527, 21.5...
$ size            <dbl> 842, 842, 842, 842, 842, 842, 842, 842, 842, 842, 842, 842, 842, 842, 8...
$ sector          <dbl> 0, 0, 0, 0, 0, 0, 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, 0.35, 0.35, 0.35,...
$ disclim         <dbl> 1.597, 1.597, 1.597, 1.597, 1.597, 1.597, 1.597, 1.597, 1.597, 1.597, 1...
$ himinty         <dbl> 0, 0, 0, 0, 0, 0, 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, -0.428, -0.428,...
$ readach         <dbl> 51.96463, 39.28524, 43.81590, 52.01301, 38.99202, 39.04975, 54.93985, 5...
$ sector.factor   <fct> 0, 0, 0, 0, 0, 0, 0, 0, 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, 0, 0, 1, 0, 1, 1, 1, 0,...
$ minority.factor <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,...

Null Model

model.null <- lmer(mathach ~  (1|schoolid), REML = FALSE, data = hsbme.clean)
summary(model.null)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: mathach ~ (1 | schoolid)
   Data: hsbme.clean

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

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.06262 -0.75365  0.02676  0.76070  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       df t value Pr(>|t|)    
(Intercept)  12.6371     0.2436 157.6209   51.87   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Question 1-2

lmerTest::rand(model.null)
ANOVA-like table for random-effects: Single term deletions

Model:
mathach ~ (1 | schoolid)
               npar logLik   AIC    LRT Df Pr(>Chisq)    
<none>            3 -23558 47122                         
(1 | schoolid)    2 -24050 48104 983.92  1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
model.1 <- lmer(mathach ~ female.factor + ses + size + sector.factor + (1|schoolid), REML = FALSE, data = hsbme.clean)
summary(model.1)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: mathach ~ female.factor + ses + size + sector.factor + (1 | schoolid)
   Data: hsbme.clean

     AIC      BIC   logLik deviance df.resid 
 46564.2  46612.3 -23275.1  46550.2     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         df t value Pr(>|t|)    
(Intercept)     1.168e+01  4.543e-01  1.613e+02  25.707  < 2e-16 ***
female.factor1 -1.201e+00  1.641e-01  6.464e+03  -7.318 2.81e-13 ***
ses             2.345e+00  1.050e-01  6.643e+03  22.325  < 2e-16 ***
size            4.972e-04  2.894e-04  1.506e+02   1.718   0.0878 .  
sector.factor1  2.378e+00  3.633e-01  1.465e+02   6.544 9.40e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

Question 3

model.2 <- lmer(mathach ~ female.factor + ses + size + sector.factor + (ses|schoolid), REML = FALSE, data = hsbme.clean)
summary(model.2)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method [lmerModLmerTest]
Formula: mathach ~ female.factor + ses + size + sector.factor + (ses |      schoolid)
   Data: hsbme.clean

     AIC      BIC   logLik deviance df.resid 
 46560.0  46621.9 -23271.0  46542.0     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. Corr
 schoolid (Intercept)  3.5435  1.8824       
          ses          0.3548  0.5956   0.60
 Residual             36.6025  6.0500       
Number of obs: 7185, groups:  schoolid, 160

Fixed effects:
                 Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)     1.148e+01  4.505e-01  1.580e+02  25.471  < 2e-16 ***
female.factor1 -1.194e+00  1.641e-01  6.561e+03  -7.275 3.86e-13 ***
ses             2.353e+00  1.152e-01  1.609e+02  20.423  < 2e-16 ***
size            4.656e-04  2.850e-04  1.416e+02   1.634    0.104    
sector.factor1  2.789e+00  3.651e-01  1.496e+02   7.639 2.42e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

Likelihood Ratio Test

lmerTest::rand(model.1)
ANOVA-like table for random-effects: Single term deletions

Model:
mathach ~ female.factor + ses + size + sector.factor + (1 | schoolid)
               npar logLik   AIC    LRT Df Pr(>Chisq)    
<none>            7 -23275 46564                         
(1 | schoolid)    6 -23421 46855 292.62  1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
lmerTest::rand(model.2)
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

```

LS0tDQp0aXRsZTogIk1vZHVsZSA2IC0gTW9kZWwgRml0IGFuZCBFZmZlY3QgU2l6ZSBmb3IgTXVsdGlsZXZlbCBNb2RlbHMiDQphdXRob3I6ICJKYWtlIFJleW5vbGRzIC0gU2VwdGVtYmVyIDMwLCAyMDIwIg0Kb3V0cHV0OiBodG1sX25vdGVib29rDQotLS0NCg0KDQpQYXJ0IDE6IENhbGN1bGF0aW5nIE11bHRpbGV2ZWwgRWZmZWN0IFNpemVzDQoNCjEuIFJ1biBhIG51bGwgbW9kZWwgd2l0aCBtYXRoYWNoIGFzIHRoZSBEViwgYW5kIHRoZW4gcnVuIGEgY29uZGl0aW9uYWwgcmFuZG9tIGludGVyY2VwdCBtb2RlbCB3aXRoIG1hdGhhY2ggYXMgdGhlIERWIGFuZCBmZW1hbGUsIHNlcywgc2l6ZSwgYW5kIHNlY3RvciBhcyBJVnMuDQoNCk51bGwgTW9kZWwgDQpBSUM6IDQ3MTIxLjgNCkJJQzogNDcxNDIuNA0KDQpNb2RlbC4xIHdpdGggZmVtYWxlLCBzZXMsIHNpemUsIGFuZCBzZWN0b3INCkFJQzogNDY1NjQuMg0KQklDOiA0NjYxMi4zDQoNCkRpZmZlcmVuY2UgYmV0d2VlbiB0aGUgbnVsbCBtb2RlbCBhbmQgTW9kZWwuMSBBSUMgYW5kIEJJQw0KQUlDOiAtNTU3LjYNCkJJQzogLTUzMC4xDQoNCg0KMi4gRmluZCB0aGUgYXBwcm9wcmlhdGUgdmFyaWFuY2UgZXN0aW1hdGVzIGFuZCBjYWxjdWxhdGUgTXVsdGlsZXZlbCBSIFNxdWFyZWQsIEJHLVBSViwgYW5kIFdHLVBSViBmb3IgdGhpcyBjb21wYXJpc29uIG9mIHR3byBtb2RlbHMuIEludGVycHJldCB0aGUgZWZmZWN0IHNpemVzLg0KDQpUaGUgUi1zcXVhcmVkIHZhbHVlIGlzIDAuMTYuIChOb3QgYmFkKQ0KUi1zcXVhcmVkIGV4cGxhaW5zIHRoZSBwcmVkaWN0aXZlIHBvd2VyIG9mIG91ciBNTE0uDQpXRy1QUlYgaXMgNS45NyUNCkJHLVBSViBpcyA2MS43NyUuDQpUYXJnZXQgbW9kZWwgcmVkdWNlZCB0aGUgd2l0aGluLWdyb3VwIHZhcmlhbmNlIGJ5IDYlIGFuZCB0aGUgYmV0d2Vlbi1ncm91cCB2YXJpYW5jZSBieSA2MiUuIFRoZSB0YXJnZXQgbW9kZWwgaXMgbXVjaCBiZXR0ZXIgYXQgcmVkdWNpbmcgYmV0d2VlbiBncm91cCB2YXJpYW5jZS4gDQoNCg0KDQpQYXJ0IDI6IFVzaW5nIEdvb2RuZXNzLW9mLUZpdCBNZWFzdXJlcw0KDQozLiAgTm93LCB1c2luZyB0aGUgdGFyZ2V0IG1vZGVsIGZyb20gUGFydCAxIChtYXRoYWNoIGFzIHRoZSBEViBhbmQgZmVtYWxlLCBzZXMsIHNpemUsIGFuZCBzZWN0b3IgYXMgSVZzKSBhcyB5b3VyIHN0YXJ0aW5nIG1vZGVsLCB0ZXN0IGFuIGFkZGl0aW9uYWwgbW9kZWwgd2hlcmUgeW91IGFsbG93IHRoZSBzbG9wZSBmb3Igc2VzIHRvIHZhcnkgcmFuZG9tbHkgYnkgc2Nob29sLiBDb25kdWN0IGEgbGlrZWxpaG9vZCByYXRpbyB0ZXN0IGNvbXBhcmluZyB0aGVzZSB0d28gbW9kZWxzLg0KDQoNCg0KTW9kZWwgd2l0aG91dCByYW5kb20gc2xvcGUuIA0KQUlDOiA0NjU2NC4yDQpCSUM6IDQ2NjEyLjMNCg0KTW9kZWwgd2l0aCByYW5kb20gc2xvcGUuIA0KQUlDOiA0NjU2MC4wICANCkJJQzogNDY2MjEuOQ0KDQo0LiAgQ29tcGFyZSB0aGUgcmVzdWx0cyBvZiB0aGUgTFIgdGVzdCB3aXRoIHRoZSBBSUMgYW5kIEJJQyB2YWx1ZXMgZm9yIHRoZSBtb2RlbHMgeW91IHJhbiBpbiAjMy4gSXMgdGhlIGV2aWRlbmNlIGNvbnNpc3RlbnQgYXMgdG8gd2hpY2ggbW9kZWwgdG8gdXNlPyBXaHkgb3Igd2h5IG5vdD8NCg0KSSBmaXJzdCBzdGFydGVkIHdpdGggbG9va2luZyBhdCB0aGUgQUlDIGFuZCBCSUMuIEhvd2V2ZXIsIHRoZXkgcHJvdmlkZWQgZGlmZmVyZW50IHJlc3VsdHMuIA0KQUlDIHdlbnQgZG93biBieSA0LjIgcG9pbnRzLiBCSUMgd2VudCB1cCBieSA5LjYgcG9pbnRzLiBTbywgdGhlIGV2aWRlbmNlIGlzIG5vdCBjb25zaXN0ZW50LiBJIHRoZW4gcmFuIHRoZSBsaWtlbGlob29kIHJhdGlvIHRlc3QgZm9yIGJvdGggdGhlIG1vZGVsIHdpdGhvdXQgdGhlIHJhbmRvbSBzbG9wZSBhbmQgd2l0aCB0aGUgcmFuZG9tIHNsb3BlLiBUaGUgcmVzdWx0cyBpbmRpY2F0ZWQgdGhhdCB3ZSBzaG91bGQgbm90IGluY2x1ZGUgdGhlIHJhbmRvbSBzbG9wZSwgc2VzLCBpbiB0aGUgbW9kZWwuIA0KDQojTG9hZCBQYWNrYWdlcyAmIERhdGENCmBgYHtyfQ0KbGlicmFyeSh0aWR5dmVyc2UpDQpsaWJyYXJ5KGxtZTQpDQpsaWJyYXJ5KHBzeWNoKQ0KbGlicmFyeShoYXZlbikNCmxpYnJhcnkodGliYmxlKQ0KDQpoc2JtZSA8LSBoYXZlbjo6cmVhZF9kdGEoImhzYm1lcmdlZGQuZHRhIikNCmdsaW1wc2UoaHNibWUpDQoNCmBgYA0KI0NsZWFuIERhdGENCg0KYGBge3J9DQpoc2JtZS5jbGVhbiA8LSBoc2JtZSAlPiUNCiAgbXV0YXRlKC4sDQogICAgICAgICBzZWN0b3IuZmFjdG9yID0gYXMuZmFjdG9yKHNlY3RvciksIA0KICAgICAgICAgZmVtYWxlLmZhY3RvciA9IGFzLmZhY3RvcihmZW1hbGUpLA0KICAgICAgICAgbWlub3JpdHkuZmFjdG9yID0gYXMuZmFjdG9yKG1pbm9yaXR5KSkNCg0KbGV2ZWxzKGhzYi5jbGVhbiRzZWN0b3IuZmFjdG9yKT1jKCJQdWJsaWMiLCAiUGFyb2NoaWFsIikNCmxldmVscyhoc2IuY2xlYW4kZmVtYWxlLmZhY3Rvcik9YygiTWFsZSIsICJGZW1hbGUiKQ0KbGV2ZWxzKGhzYi5jbGVhbiRtaW5vcml0eS5mYWN0b3IpPWMoIk5vbi1NaW5vcml0eSIsICJNaW5vcml0eSIpDQoNCnN0cihoc2IuY2xlYW4kc2VjdG9yLmZhY3RvcikNCg0KZ2xpbXBzZShoc2JtZS5jbGVhbikgIA0KYGBgDQoNCg0KIyBOdWxsIE1vZGVsDQpgYGB7cn0NCm1vZGVsLm51bGwgPC0gbG1lcihtYXRoYWNoIH4gICgxfHNjaG9vbGlkKSwgUkVNTCA9IEZBTFNFLCBkYXRhID0gaHNibWUuY2xlYW4pDQpzdW1tYXJ5KG1vZGVsLm51bGwpDQpgYGAgIA0KI1F1ZXN0aW9uIDEtMg0KDQpgYGB7cn0NCmxtZXJUZXN0OjpyYW5kKG1vZGVsLm51bGwpDQpgYGANCg0KDQpgYGB7cn0NCm1vZGVsLjEgPC0gbG1lcihtYXRoYWNoIH4gZmVtYWxlLmZhY3RvciArIHNlcyArIHNpemUgKyBzZWN0b3IuZmFjdG9yICsgKDF8c2Nob29saWQpLCBSRU1MID0gRkFMU0UsIGRhdGEgPSBoc2JtZS5jbGVhbikNCnN1bW1hcnkobW9kZWwuMSkNCmBgYA0KIA0KIyBRdWVzdGlvbiAzDQoNCmBgYHtyfQ0KbW9kZWwuMiA8LSBsbWVyKG1hdGhhY2ggfiBmZW1hbGUuZmFjdG9yICsgc2VzICsgc2l6ZSArIHNlY3Rvci5mYWN0b3IgKyAoc2VzfHNjaG9vbGlkKSwgUkVNTCA9IEZBTFNFLCBkYXRhID0gaHNibWUuY2xlYW4pDQpzdW1tYXJ5KG1vZGVsLjIpDQpgYGANCiMgTGlrZWxpaG9vZCBSYXRpbyBUZXN0DQoNCmBgYHtyfQ0KbG1lclRlc3Q6OnJhbmQobW9kZWwuMSkNCmBgYA0KDQpgYGB7cn0NCmxtZXJUZXN0OjpyYW5kKG1vZGVsLjIpDQpgYGANCg0KDQoNCmBgYA==