Read the dataset.

dat <- read_csv("Assignment_data.csv")
Rows: 2310 Columns: 9
── Column specification ───────────────────────────────
Delimiter: ","
dbl (9): stdid, schid, p7vrq, p7read, dadunemp, dad...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(dat)
summary(dat)
     stdid            schid           p7vrq         
 Min.   :   1.0   Min.   : 0.00   Min.   :-27.0280  
 1st Qu.: 578.2   1st Qu.: 5.00   1st Qu.: -7.0280  
 Median :1155.5   Median : 9.00   Median : -0.0280  
 Mean   :1155.5   Mean   :10.01   Mean   :  0.5058  
 3rd Qu.:1732.8   3rd Qu.:16.00   3rd Qu.:  7.9720  
 Max.   :2310.0   Max.   :20.00   Max.   : 42.9720  
     p7read             dadunemp          daded       
 Min.   :-31.86600   Min.   :0.0000   Min.   :0.0000  
 1st Qu.: -9.86600   1st Qu.:0.0000   1st Qu.:0.0000  
 Median : -0.86600   Median :0.0000   Median :0.0000  
 Mean   : -0.04435   Mean   :0.1091   Mean   :0.2152  
 3rd Qu.:  9.13400   3rd Qu.:0.0000   3rd Qu.:0.0000  
 Max.   : 28.13400   Max.   :1.0000   Max.   :1.0000  
     momed             male          highattain    
 Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
 1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
 Median :0.0000   Median :0.0000   Median :0.0000  
 Mean   :0.2485   Mean   :0.4801   Mean   :0.1649  
 3rd Qu.:0.0000   3rd Qu.:1.0000   3rd Qu.:0.0000  
 Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  

Q1 Random Intercept Model and Calculate ICC

mod0 <- lmer(p7read ~ (1|schid), data = dat)
summary(mod0)
Linear mixed model fit by REML ['lmerMod']
Formula: p7read ~ (1 | schid)
   Data: dat

REML criterion at convergence: 18550.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.76533 -0.68850 -0.00339  0.68667  2.56831 

Random effects:
 Groups   Name        Variance Std.Dev.
 schid    (Intercept)  20.08    4.481  
 Residual             176.64   13.290  
Number of obs: 2310, groups:  schid, 17

Fixed effects:
            Estimate Std. Error t value
(Intercept)   -0.296      1.130  -0.262
tab_model(mod0)

var_sch <- as.numeric(VarCorr(mod0)$schid)
var_resid <- attr(VarCorr(mod0), "sc")^2
icc <- var_sch / (var_sch + var_resid)
icc
[1] 0.1020834

The fixed intercept shows the average reading score across schools.

The random intercept variance shows variability between schools.

The ICC measures how much of the total variance is due to school differences.

Q2 Compare ICC and R-squared

mod1 <- lmer(p7read ~ highattain + (1|schid), data = dat)
summary(mod1)
Linear mixed model fit by REML ['lmerMod']
Formula: p7read ~ highattain + (1 | schid)
   Data: dat

REML criterion at convergence: 18047

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.79101 -0.67564  0.00624  0.69202  2.97652 

Random effects:
 Groups   Name        Variance Std.Dev.
 schid    (Intercept)  10.93    3.307  
 Residual             142.49   11.937  
Number of obs: 2310, groups:  schid, 17

Fixed effects:
            Estimate Std. Error t value
(Intercept)  -2.9227     0.8557  -3.416
highattain   16.1563     0.6810  23.723

Correlation of Fixed Effects:
           (Intr)
highattain -0.131
tab_model(mod0, mod1)

var_sch_1 <- as.numeric(VarCorr(mod1)$schid)
var_resid_1 <- attr(VarCorr(mod1), "sc")^2
icc_1 <- var_sch_1 / (var_sch_1 + var_resid_1)
icc_1
[1] 0.0712681

The slope coefficient for high attainment shows its impact on reading scores.

The ICC and R-squared shows the explained variance increase.

Q3 Random Slope Model and ICC

mod2 <- lmer(p7read ~ highattain + (highattain|schid), data = dat)
boundary (singular) fit: see help('isSingular')
summary(mod2)
Linear mixed model fit by REML ['lmerMod']
Formula: p7read ~ highattain + (highattain | schid)
   Data: dat

REML criterion at convergence: 18046.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.80824 -0.66539  0.01059  0.68236  2.98821 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 schid    (Intercept)  11.7750  3.4315       
          highattain    0.5102  0.7143  -1.00
 Residual             142.4261 11.9342       
Number of obs: 2310, groups:  schid, 17

Fixed effects:
            Estimate Std. Error t value
(Intercept)  -2.9131     0.8846  -3.293
highattain   16.3554     0.7002  23.359

Correlation of Fixed Effects:
           (Intr)
highattain -0.361
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
tab_model(mod1, mod2)

var_sch_2 <- as.numeric(VarCorr(mod2)$schid)[1]
var_slope <- as.numeric(VarCorr(mod2)$schid)[2]
var_resid_2 <- attr(VarCorr(mod2), "sc")^2
icc_2 <- var_sch_2 / (var_sch_2 + var_resid_2)
icc_2
[1] 0.07636134

The random slope shows the variability in high attainment effects across schools.

The correlation between random effects shows further insights.

Random Intercept Model (Q1)

ICC = 0.1021 (10.2%)

Shows a moderate school-level effect on reading scores.

High Attainment (Q2)

#ICC = 0.0713 (7.1%) # The ICC decreased, shows that high attainment explains part of the school-level variance in reading scores. # Some variance attributed to schools is explained by high attainment differences.

Random Slope Model (Q3)

ICC = 0.0753 (7.5%)

The variance shows at the school level slightly increases compared to ICC, but remains lower than ICC.

Shows that allowing high attainment to vary across schools allows more variance.

LS0tCnRpdGxlOiAiRURVIDcwOSBDWVJPIEhXIDEiCm91dHB1dDogaHRtbF9ub3RlYm9vawpkYXRlOiAyLzExLzIwMjUKLS0tCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0KbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkobG1lNCkKbGlicmFyeShzalBsb3QpCmtuaXRyOjpvcHRzX2NodW5rJHNldChlY2hvID0gVFJVRSwgbWVzc2FnZSA9IEZBTFNFLCB3YXJuaW5nID0gRkFMU0UsIGVycm9yID0gRkFMU0UpCmBgYAoKCgoKUmVhZCB0aGUgZGF0YXNldC4gCmBgYHtyfQpkYXQgPC0gcmVhZF9jc3YoIkFzc2lnbm1lbnRfZGF0YS5jc3YiKQpoZWFkKGRhdCkKc3VtbWFyeShkYXQpCmBgYAoKIyBRMSBSYW5kb20gSW50ZXJjZXB0IE1vZGVsIGFuZCBDYWxjdWxhdGUgSUNDCmBgYHtyfQptb2QwIDwtIGxtZXIocDdyZWFkIH4gKDF8c2NoaWQpLCBkYXRhID0gZGF0KQpzdW1tYXJ5KG1vZDApCnRhYl9tb2RlbChtb2QwKQoKdmFyX3NjaCA8LSBhcy5udW1lcmljKFZhckNvcnIobW9kMCkkc2NoaWQpCnZhcl9yZXNpZCA8LSBhdHRyKFZhckNvcnIobW9kMCksICJzYyIpXjIKaWNjIDwtIHZhcl9zY2ggLyAodmFyX3NjaCArIHZhcl9yZXNpZCkKaWNjCmBgYAoKIyBUaGUgZml4ZWQgaW50ZXJjZXB0IHNob3dzIHRoZSBhdmVyYWdlIHJlYWRpbmcgc2NvcmUgYWNyb3NzIHNjaG9vbHMuCiMgVGhlIHJhbmRvbSBpbnRlcmNlcHQgdmFyaWFuY2Ugc2hvd3MgdmFyaWFiaWxpdHkgYmV0d2VlbiBzY2hvb2xzLgojIFRoZSBJQ0MgbWVhc3VyZXMgaG93IG11Y2ggb2YgdGhlIHRvdGFsIHZhcmlhbmNlIGlzIGR1ZSB0byBzY2hvb2wgZGlmZmVyZW5jZXMuCgojIFEyIENvbXBhcmUgSUNDIGFuZCBSLXNxdWFyZWQKYGBge3J9Cm1vZDEgPC0gbG1lcihwN3JlYWQgfiBoaWdoYXR0YWluICsgKDF8c2NoaWQpLCBkYXRhID0gZGF0KQpzdW1tYXJ5KG1vZDEpCnRhYl9tb2RlbChtb2QwLCBtb2QxKQoKdmFyX3NjaF8xIDwtIGFzLm51bWVyaWMoVmFyQ29ycihtb2QxKSRzY2hpZCkKdmFyX3Jlc2lkXzEgPC0gYXR0cihWYXJDb3JyKG1vZDEpLCAic2MiKV4yCmljY18xIDwtIHZhcl9zY2hfMSAvICh2YXJfc2NoXzEgKyB2YXJfcmVzaWRfMSkKaWNjXzEKYGBgCgojIFRoZSBzbG9wZSBjb2VmZmljaWVudCBmb3IgaGlnaCBhdHRhaW5tZW50IHNob3dzIGl0cyBpbXBhY3Qgb24gcmVhZGluZyBzY29yZXMuCiMgVGhlIElDQyBhbmQgUi1zcXVhcmVkIHNob3dzIHRoZSBleHBsYWluZWQgdmFyaWFuY2UgaW5jcmVhc2UuCgoKIyBRMyBSYW5kb20gU2xvcGUgTW9kZWwgYW5kIElDQwpgYGB7cn0KbW9kMiA8LSBsbWVyKHA3cmVhZCB+IGhpZ2hhdHRhaW4gKyAoaGlnaGF0dGFpbnxzY2hpZCksIGRhdGEgPSBkYXQpCnN1bW1hcnkobW9kMikKdGFiX21vZGVsKG1vZDEsIG1vZDIpCgp2YXJfc2NoXzIgPC0gYXMubnVtZXJpYyhWYXJDb3JyKG1vZDIpJHNjaGlkKVsxXQp2YXJfc2xvcGUgPC0gYXMubnVtZXJpYyhWYXJDb3JyKG1vZDIpJHNjaGlkKVsyXQp2YXJfcmVzaWRfMiA8LSBhdHRyKFZhckNvcnIobW9kMiksICJzYyIpXjIKaWNjXzIgPC0gdmFyX3NjaF8yIC8gKHZhcl9zY2hfMiArIHZhcl9yZXNpZF8yKQppY2NfMgpgYGAKCiMgVGhlIHJhbmRvbSBzbG9wZSBzaG93cyB0aGUgdmFyaWFiaWxpdHkgaW4gaGlnaCBhdHRhaW5tZW50IGVmZmVjdHMgYWNyb3NzIHNjaG9vbHMuCiMgVGhlIGNvcnJlbGF0aW9uIGJldHdlZW4gcmFuZG9tIGVmZmVjdHMgc2hvd3MgZnVydGhlciBpbnNpZ2h0cy4KCgojIFJhbmRvbSBJbnRlcmNlcHQgTW9kZWwgKFExKQoKIyBJQ0MgPSAwLjEwMjEgKDEwLjIlKQojIFNob3dzIGEgbW9kZXJhdGUgc2Nob29sLWxldmVsIGVmZmVjdCBvbiByZWFkaW5nIHNjb3Jlcy4KCiMgSGlnaCBBdHRhaW5tZW50IChRMikKCiNJQ0MgPSAwLjA3MTMgKDcuMSUpCiMgVGhlIElDQyBkZWNyZWFzZWQsIHNob3dzIHRoYXQgaGlnaCBhdHRhaW5tZW50IGV4cGxhaW5zIHBhcnQgb2YgdGhlIHNjaG9vbC1sZXZlbCB2YXJpYW5jZSBpbiByZWFkaW5nIHNjb3Jlcy4KIyBTb21lIHZhcmlhbmNlIGF0dHJpYnV0ZWQgdG8gc2Nob29scyBpcyBleHBsYWluZWQgYnkgaGlnaCBhdHRhaW5tZW50IGRpZmZlcmVuY2VzLgoKIyBSYW5kb20gU2xvcGUgTW9kZWwgKFEzKQoKIyBJQ0MgPSAwLjA3NTMgKDcuNSUpCiMgVGhlIHZhcmlhbmNlIHNob3dzIGF0IHRoZSBzY2hvb2wgbGV2ZWwgc2xpZ2h0bHkgaW5jcmVhc2VzIGNvbXBhcmVkIHRvIElDQywgYnV0IHJlbWFpbnMgbG93ZXIgdGhhbiBJQ0MuCiMgU2hvd3MgdGhhdCBhbGxvd2luZyBoaWdoIGF0dGFpbm1lbnQgdG8gdmFyeSBhY3Jvc3Mgc2Nob29scyBhbGxvd3MgbW9yZSB2YXJpYW5jZS4K