This week, we will be using George Leckie’s (2009) scotland dataset, which includes 2,310 students, nested within 17 schools and separately nested within 500+ neighborhoods.

Load in Our MVP Packages

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

Part 1: Loading in the Data, and then Recoding Variables

Load in the Data


scotland <- haven::read_dta("scotland.dta")

glimpse(scotland)
Rows: 2,310
Columns: 11
$ neighid  <dbl> 26, 26, 26, 26, 26, 27, 29, 29, 29, 29, 29, 29, 29, 2…
$ schid    <dbl> 20, 20, 20, 20, 20, 20, 18, 20, 20, 20, 20, 20, 20, 2…
$ attain   <dbl> 1.5177, -1.3276, 0.5610, 1.5177, -1.3276, -0.1325, 0.…
$ p7vrq    <dbl> 17.972, -10.028, 2.972, 1.972, -1.028, 3.972, 8.972, …
$ p7read   <dbl> 17.134, -27.866, 6.134, 11.134, -0.866, -0.866, 6.134…
$ dadocc   <dbl> 16.196, -3.454, 2.316, -9.094, -3.454, -3.454, 16.196…
$ dadunemp <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ daded    <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1,…
$ momed    <dbl> 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,…
$ male     <dbl> 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1,…
$ deprive  <dbl> -0.551, -0.551, -0.551, -0.551, -0.551, 0.147, -0.083…

The describe function from Hmisc is a great codebook…

describe(scotland)
scotland 

 11  Variables      2310  Observations
-------------------------------------------------------------------------
neighid  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    2310        0      524        1    495.3    305.5     61.0    143.9 
     .25      .50      .75      .90      .95 
   240.2    530.0    707.0    808.0    861.0 

lowest :   26   27   29   30   31, highest: 1092 1095 1096 1097 1098
-------------------------------------------------------------------------
schid  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    2310        0       17    0.995    10.01    7.174        0        2 
     .25      .50      .75      .90      .95 
       5        9       16       19       20 

lowest :  0  1  2  3  5, highest: 16 17 18 19 20
                                                                      
Value          0     1     2     3     5     6     7     8     9    10
Frequency    146    22   146   159   155   101   286   112   136   133
Proportion 0.063 0.010 0.063 0.069 0.067 0.044 0.124 0.048 0.059 0.058
                                                    
Value         13    15    16    17    18    19    20
Frequency     92   190   111   154    91   102   174
Proportion 0.040 0.082 0.048 0.067 0.039 0.044 0.075
-------------------------------------------------------------------------
attain  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    2310        0       14    0.987   0.0934    1.124  -1.3276  -1.3276 
     .25      .50      .75      .90      .95 
 -0.5812   0.1582   0.7350   1.5177   1.5177 

lowest : -1.3276 -0.5812 -0.3600 -0.1325  0.0293
highest:  0.7350  0.9127  1.1405  1.5177  2.4151
-------------------------------------------------------------------------
p7vrq  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    2310        0       68    0.999   0.5058    11.95  -17.028  -13.028 
     .25      .50      .75      .90      .95 
  -7.028   -0.028    7.972   13.972   16.972 

lowest : -27.028 -26.028 -25.028 -24.028 -23.028
highest:  35.972  36.972  40.972  41.972  42.972
-------------------------------------------------------------------------
p7read  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    2310        0       61        1 -0.04435     15.8  -23.866  -17.866 
     .25      .50      .75      .90      .95 
  -9.866   -0.866    9.134   19.134   24.134 

lowest : -31.866 -30.866 -29.866 -28.866 -27.866
highest:  24.134  25.134  26.134  27.134  28.134
-------------------------------------------------------------------------
dadocc  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd 
    2310        0        7    0.933  -0.4642    12.33 

lowest : -23.454 -11.494  -9.094  -3.454   2.316
highest:  -9.094  -3.454   2.316  16.196  29.226
                                                                  
Value      -23.454 -11.494  -9.094  -3.454   2.316  16.196  29.226
Frequency       91     285     303     884     242     397     108
Proportion   0.039   0.123   0.131   0.383   0.105   0.172   0.047
-------------------------------------------------------------------------
dadunemp  Format:%12.0g 
       n  missing distinct     Info      Sum     Mean      Gmd 
    2310        0        2    0.292      252   0.1091   0.1945 

-------------------------------------------------------------------------
daded  Format:%12.0g 
       n  missing distinct     Info      Sum     Mean      Gmd 
    2310        0        2    0.507      497   0.2152   0.3379 

-------------------------------------------------------------------------
momed  Format:%12.0g 
       n  missing distinct     Info      Sum     Mean      Gmd 
    2310        0        2     0.56      574   0.2485   0.3736 

-------------------------------------------------------------------------
male  Format:%12.0g 
       n  missing distinct     Info      Sum     Mean      Gmd 
    2310        0        2    0.749     1109   0.4801   0.4994 

-------------------------------------------------------------------------
deprive  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    2310        0      458        1  0.02167   0.6664  -0.8250  -0.6930 
     .25      .50      .75      .90      .95 
 -0.3970  -0.0620   0.2957   0.8410   1.1400 

lowest : -1.082 -1.048 -1.030 -0.983 -0.975, highest:  2.330  2.419  2.438  2.498  2.959
-------------------------------------------------------------------------

Part One: Data Exploration and Coding/ Labelling

Start with variable labelling using the label function

label(scotland$neighid)="Neighborhood ID"
label(scotland$schid)="School ID"
label(scotland$attain)="Student Measure of Educational Attainment"
label(scotland$p7vrq)="Primary 7 Verbal Reasoning Quotient"
label(scotland$p7read)="Primary 7 Reading Test Scores"
label(scotland$dadocc)="School Mean of Dad's Occupation Score on Hope-Goldthorpe Scale"
label(scotland$dadunemp)="Dad Currently Unemployed"
label(scotland$daded)="Dad School After Age 15"
label(scotland$momed)="Mom School After Age 15"
label(scotland$male)="Male"
label(scotland$deprive)="Neighborhood Deprivation Score (Poverty, Health, and Housing)"

glimpse(scot.clean)

Rows: 2,310
Columns: 15
$ neighid      <labelled> 26, 26, 26, 26, 26, 27, 29, 29, 29, 29, 29, …
$ schid        <labelled> 20, 20, 20, 20, 20, 20, 18, 20, 20, 20, 20, …
$ attain       <labelled> 1.5177, -1.3276, 0.5610, 1.5177, -1.3276, -0…
$ p7vrq        <labelled> 17.972, -10.028, 2.972, 1.972, -1.028, 3.972…
$ p7read       <labelled> 17.134, -27.866, 6.134, 11.134, -0.866, -0.8…
$ dadocc       <labelled> 16.196, -3.454, 2.316, -9.094, -3.454, -3.45…
$ dadunemp     <labelled> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ daded        <labelled> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,…
$ momed        <labelled> 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,…
$ male         <labelled> 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1,…
$ deprive      <labelled> -0.551, -0.551, -0.551, -0.551, -0.551, 0.14…
$ dadunemp.fac <fct> No, No, No, No, No, No, No, No, No, No, No, No, N…
$ daded.fac    <fct> No, No, No, No, No, No, Yes, No, No, No, No, No, …
$ momed.fac    <fct> No, No, No, No, No, Yes, Yes, No, No, No, No, No,…
$ male.fac     <fct> Yes, No, No, No, Yes, No, Yes, No, Yes, No, No, N…
describe(scot.clean)
scot.clean 

 15  Variables      2310  Observations
-------------------------------------------------------------------------
neighid : Neighborhood ID  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    2310        0      524        1    495.3    305.5     61.0    143.9 
     .25      .50      .75      .90      .95 
   240.2    530.0    707.0    808.0    861.0 

lowest :   26   27   29   30   31, highest: 1092 1095 1096 1097 1098
-------------------------------------------------------------------------
schid : School ID  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    2310        0       17    0.995    10.01    7.174        0        2 
     .25      .50      .75      .90      .95 
       5        9       16       19       20 

lowest :  0  1  2  3  5, highest: 16 17 18 19 20
                                                                      
Value          0     1     2     3     5     6     7     8     9    10
Frequency    146    22   146   159   155   101   286   112   136   133
Proportion 0.063 0.010 0.063 0.069 0.067 0.044 0.124 0.048 0.059 0.058
                                                    
Value         13    15    16    17    18    19    20
Frequency     92   190   111   154    91   102   174
Proportion 0.040 0.082 0.048 0.067 0.039 0.044 0.075
-------------------------------------------------------------------------
attain : Student Measure of Educational Attainment  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    2310        0       14    0.987   0.0934    1.124  -1.3276  -1.3276 
     .25      .50      .75      .90      .95 
 -0.5812   0.1582   0.7350   1.5177   1.5177 

lowest : -1.3276 -0.5812 -0.3600 -0.1325  0.0293
highest:  0.7350  0.9127  1.1405  1.5177  2.4151
-------------------------------------------------------------------------
p7vrq : Primary 7 Verbal Reasoning Quotient  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    2310        0       68    0.999   0.5058    11.95  -17.028  -13.028 
     .25      .50      .75      .90      .95 
  -7.028   -0.028    7.972   13.972   16.972 

lowest : -27.028 -26.028 -25.028 -24.028 -23.028
highest:  35.972  36.972  40.972  41.972  42.972
-------------------------------------------------------------------------
p7read : Primary 7 Reading Test Scores  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    2310        0       61        1 -0.04435     15.8  -23.866  -17.866 
     .25      .50      .75      .90      .95 
  -9.866   -0.866    9.134   19.134   24.134 

lowest : -31.866 -30.866 -29.866 -28.866 -27.866
highest:  24.134  25.134  26.134  27.134  28.134
-------------------------------------------------------------------------
dadocc : School Mean of Dad's Occupation Score on Hope-Goldthorpe Scale  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd 
    2310        0        7    0.933  -0.4642    12.33 

lowest : -23.454 -11.494  -9.094  -3.454   2.316
highest:  -9.094  -3.454   2.316  16.196  29.226
                                                                  
Value      -23.454 -11.494  -9.094  -3.454   2.316  16.196  29.226
Frequency       91     285     303     884     242     397     108
Proportion   0.039   0.123   0.131   0.383   0.105   0.172   0.047
-------------------------------------------------------------------------
dadunemp : Dad Currently Unemployed  Format:%12.0g 
       n  missing distinct     Info      Sum     Mean      Gmd 
    2310        0        2    0.292      252   0.1091   0.1945 

-------------------------------------------------------------------------
daded : Dad School After Age 15  Format:%12.0g 
       n  missing distinct     Info      Sum     Mean      Gmd 
    2310        0        2    0.507      497   0.2152   0.3379 

-------------------------------------------------------------------------
momed : Mom School After Age 15  Format:%12.0g 
       n  missing distinct     Info      Sum     Mean      Gmd 
    2310        0        2     0.56      574   0.2485   0.3736 

-------------------------------------------------------------------------
male : Male  Format:%12.0g 
       n  missing distinct     Info      Sum     Mean      Gmd 
    2310        0        2    0.749     1109   0.4801   0.4994 

-------------------------------------------------------------------------
deprive : Neighborhood Deprivation Score (Poverty, Health, and Housing)  Format:%12.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    2310        0      458        1  0.02167   0.6664  -0.8250  -0.6930 
     .25      .50      .75      .90      .95 
 -0.3970  -0.0620   0.2957   0.8410   1.1400 

lowest : -1.082 -1.048 -1.030 -0.983 -0.975, highest:  2.330  2.419  2.438  2.498  2.959
-------------------------------------------------------------------------
dadunemp.fac : Dad Currently Unemployed 
       n  missing distinct 
    2310        0        2 
                      
Value         No   Yes
Frequency   2058   252
Proportion 0.891 0.109
-------------------------------------------------------------------------
daded.fac : Dad School After Age 15 
       n  missing distinct 
    2310        0        2 
                      
Value         No   Yes
Frequency   1813   497
Proportion 0.785 0.215
-------------------------------------------------------------------------
momed.fac : Mom School After Age 15 
       n  missing distinct 
    2310        0        2 
                      
Value         No   Yes
Frequency   1736   574
Proportion 0.752 0.248
-------------------------------------------------------------------------
male.fac : Male 
       n  missing distinct 
    2310        0        2 
                    
Value        No  Yes
Frequency  1201 1109
Proportion 0.52 0.48
-------------------------------------------------------------------------

Explore our Outcome, attain

ggplot(data = scot.clean, mapping = aes(x = attain)) +
  geom_histogram(bins = 20) + 
  labs(title = "Distribution of Educational Attainment Scores",
       x = "Standardized Educational Attainment Score") +
  theme_minimal()

Explore a few other predictors

ggplot(data = scot.clean, mapping = aes(x = dadocc)) +
  geom_histogram(bins = 30) + 
  labs(title = "Distribution of (School Mean) of Dad's Occupation Scores",
       x = "Standardized Educational Attainment Score") +
  theme_minimal()

ggplot(data = scot.clean, mapping = aes(x = deprive)) +
  geom_histogram(bins = 50) + 
  labs(title = "Distribution of Neighborhood Context Scores",
       x = "Standardized Neighborhood Deprivation Score") +
  theme_minimal()

Part Two: Modeling

First, run two separate models with school and neighborhood and the L2 clusters…

model.null.school <- lmer(attain ~ (1|schid), REML = FALSE, data = scot.clean)
summary(model.null.school)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: attain ~ (1 | schid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
  6448.2   6465.4  -3221.1   6442.2     2307 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.12306 -0.70227  0.00178  0.60966  2.81884 

Random effects:
 Groups   Name        Variance Std.Dev.
 schid    (Intercept) 0.08874  0.2979  
 Residual             0.93441  0.9666  
Number of obs: 2310, groups:  schid, 17

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.08227    0.07568   1.087

Calculate School ICC

icc.school <- 0.08874/(0.08874 + 0.93441)
icc.school
[1] 0.08673215
model.null.neigh <- lmer(attain ~ (1|neighid), REML = FALSE, data = scot.clean)
summary(model.null.neigh)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: attain ~ (1 | neighid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
  6422.0   6439.2  -3208.0   6416.0     2307 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.33164 -0.65532  0.01513  0.58177  2.96174 

Random effects:
 Groups   Name        Variance Std.Dev.
 neighid  (Intercept) 0.2015   0.4489  
 Residual             0.8044   0.8969  
Number of obs: 2310, groups:  neighid, 524

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.08202    0.02844   2.885

Calculate Neighborhood ICC

icc.neigh <- 0.2015/(0.2015 + 0.8044)
icc.neigh
[1] 0.2003181

Now, a null additive cross-classified model…

model.null.crossed <- lmer(attain ~ (1|schid) + (1|neighid), REML = FALSE, data = scot.clean)
summary(model.null.crossed)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: attain ~ (1 | schid) + (1 | neighid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
  6364.7   6387.7  -3178.4   6356.7     2306 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.36115 -0.69814  0.01345  0.58371  2.92697 

Random effects:
 Groups   Name        Variance Std.Dev.
 neighid  (Intercept) 0.14122  0.3758  
 schid    (Intercept) 0.07545  0.2747  
 Residual             0.79902  0.8939  
Number of obs: 2310, groups:  neighid, 524; schid, 17

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.07535    0.07222   1.043

Calculate Neighborhood ICC in Crossed Model

icc.neigh.crossed <- 0.14122/(0.14122 + 0.07545 + 0.79902)
icc.neigh.crossed
[1] 0.1390385

Calculate School ICC in Crossed Model

icc.school.crossed <- 0.07545/(0.14122 + 0.07545 + 0.79902)
icc.school.crossed
[1] 0.07428448

Add an L1 Predictor, Male…

model.1 <- lmer(attain ~ male.fac + (1|schid) + (1|neighid), REML = FALSE, data = scot.clean)
summary(model.1)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: attain ~ male.fac + (1 | schid) + (1 | neighid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
  6360.6   6389.3  -3175.3   6350.6     2305 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.41396 -0.69992  0.00977  0.58312  2.97855 

Random effects:
 Groups   Name        Variance Std.Dev.
 neighid  (Intercept) 0.14113  0.3757  
 schid    (Intercept) 0.07439  0.2727  
 Residual             0.79681  0.8926  
Number of obs: 2310, groups:  neighid, 524; schid, 17

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.12206    0.07421   1.645
male.facYes -0.09656    0.03898  -2.477

Correlation of Fixed Effects:
            (Intr)
male.facYes -0.254

L1 Plus School-Level Predictor Only (dadocc)…

model.2 <- lmer(attain ~ male.fac + dadocc + (1|schid) + (1|neighid), REML = FALSE, data = scot.clean)
summary(model.2)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: attain ~ male.fac + dadocc + (1 | schid) + (1 | neighid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
  6146.5   6181.0  -3067.2   6134.5     2304 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.48797 -0.69407  0.01097  0.59325  3.10701 

Random effects:
 Groups   Name        Variance Std.Dev.
 neighid  (Intercept) 0.07413  0.2723  
 schid    (Intercept) 0.04277  0.2068  
 Residual             0.76119  0.8725  
Number of obs: 2310, groups:  neighid, 524; schid, 17

Fixed effects:
             Estimate Std. Error t value
(Intercept)  0.139052   0.058591   2.373
male.facYes -0.093390   0.037575  -2.485
dadocc       0.025708   0.001662  15.472

Correlation of Fixed Effects:
            (Intr) ml.fcY
male.facYes -0.309       
dadocc       0.012  0.008

L1 Plus Neighborhood-Level Predictors Only (deprive)…

model.3 <- lmer(attain ~ male.fac + deprive + (1|schid) + (1|neighid), REML = FALSE, data = scot.clean)
summary(model.3)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: attain ~ male.fac + deprive + (1 | schid) + (1 | neighid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
  6238.1   6272.5  -3113.0   6226.1     2304 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.4480 -0.6746  0.0031  0.5900  3.5529 

Random effects:
 Groups   Name        Variance Std.Dev.
 neighid  (Intercept) 0.06053  0.2460  
 schid    (Intercept) 0.03840  0.1960  
 Residual             0.80436  0.8969  
Number of obs: 2310, groups:  neighid, 524; schid, 17

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.14566    0.05638   2.584
male.facYes -0.10269    0.03843  -2.672
deprive     -0.46836    0.03820 -12.260

Correlation of Fixed Effects:
            (Intr) ml.fcY
male.facYes -0.329       
deprive     -0.023  0.009

L1 Plus Both School and Neighborhood Predictors…

model.4 <- lmer(attain ~ male.fac + deprive + dadocc + (1|schid) + (1|neighid), REML = FALSE, data = scot.clean)
summary(model.4)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: 
attain ~ male.fac + deprive + dadocc + (1 | schid) + (1 | neighid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
  6055.3   6095.6  -3020.7   6041.3     2303 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6406 -0.6716  0.0152  0.6010  3.5775 

Random effects:
 Groups   Name        Variance Std.Dev.
 neighid  (Intercept) 0.03005  0.1733  
 schid    (Intercept) 0.02313  0.1521  
 Residual             0.76444  0.8743  
Number of obs: 2310, groups:  neighid, 524; schid, 17

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.15565    0.04635   3.358
male.facYes -0.09836    0.03709  -2.652
deprive     -0.36485    0.03540 -10.306
dadocc       0.02330    0.00166  14.037

Correlation of Fixed Effects:
            (Intr) ml.fcY depriv
male.facYes -0.386              
deprive     -0.022  0.010       
dadocc       0.008  0.011  0.224

Test Random effects for male.fac at neighborhood level and at school level…

model.5 <- lmer(attain ~ male.fac + deprive + dadocc + (male.fac|schid) + (male.fac|neighid), REML = FALSE, data = scot.clean)
Model failed to converge with max|grad| = 0.00226165 (tol = 0.002, component 1)
summary(model.5)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: 
attain ~ male.fac + deprive + dadocc + (male.fac | schid) + (male.fac |  
    neighid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
  6061.9   6125.1  -3020.0   6039.9     2299 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6547 -0.6662  0.0194  0.6131  3.5925 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 neighid  (Intercept) 0.035545 0.18853       
          male.facYes 0.014163 0.11901  -0.37
 schid    (Intercept) 0.030597 0.17492       
          male.facYes 0.008165 0.09036  -0.60
 Residual             0.758356 0.87084       
Number of obs: 2310, groups:  neighid, 524; schid, 17

Fixed effects:
             Estimate Std. Error t value
(Intercept)  0.154424   0.051225   3.015
male.facYes -0.095600   0.043996  -2.173
deprive     -0.364446   0.035435 -10.285
dadocc       0.023301   0.001659  14.042

Correlation of Fixed Effects:
            (Intr) ml.fcY depriv
male.facYes -0.568              
deprive     -0.022  0.011       
dadocc       0.005  0.013  0.225
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.00226165 (tol = 0.002, component 1)
lmerTest::rand(model.5)
ANOVA-like table for random-effects: Single term deletions

Model:
attain ~ male.fac + deprive + dadocc + (male.fac | schid) + (male.fac | neighid)
                                 npar  logLik    AIC     LRT Df
<none>                             11 -3020.0 6061.9           
male.fac in (male.fac | schid)      9 -3020.5 6059.1 1.11528  2
male.fac in (male.fac | neighid)    9 -3020.0 6058.1 0.11791  2
                                 Pr(>Chisq)
<none>                                     
male.fac in (male.fac | schid)       0.5726
male.fac in (male.fac | neighid)     0.9427

Finally, test a model with school-neighborhood interaction…

model.6 <- lmer(attain ~ male.fac + deprive + dadocc + deprive:dadocc + (1|schid) + (1|neighid), REML = FALSE, data = scot.clean)
summary(model.6)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: attain ~ male.fac + deprive + dadocc + deprive:dadocc + (1 |  
    schid) + (1 | neighid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
  6046.5   6092.5  -3015.3   6030.5     2302 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7227 -0.6989  0.0235  0.6097  3.5774 

Random effects:
 Groups   Name        Variance Std.Dev.
 neighid  (Intercept) 0.02719  0.1649  
 schid    (Intercept) 0.02192  0.1481  
 Residual             0.76333  0.8737  
Number of obs: 2310, groups:  neighid, 524; schid, 17

Fixed effects:
                Estimate Std. Error t value
(Intercept)     0.140420   0.045710   3.072
male.facYes    -0.105781   0.037087  -2.852
deprive        -0.374940   0.035199 -10.652
dadocc          0.022321   0.001686  13.243
deprive:dadocc -0.009277   0.002808  -3.303

Correlation of Fixed Effects:
            (Intr) ml.fcY depriv dadocc
male.facYes -0.383                     
deprive     -0.013  0.014              
dadocc       0.027  0.022  0.237       
depriv:ddcc  0.104  0.060  0.082  0.187

Use interplot to visualize the interaction between deprive and dadocc…

interplot::interplot(model.6, var1 = "deprive", var2 = "dadocc")

Step 6: Compare Original “Final” Model, Robust Model, and Trimmed Model for Differences (Sensitivity Analysis)

library(modelsummary)
library(broom.mixed)
models <- list(model.null.neigh, model.null.school, model.null.crossed, model.1, model.2, model.3, model.4, model.5, model.6)
modelsummary(models, output = "markdown")
Model 1 Model 2 Model 3 Model 4 Model 5 Model 6 Model 7 Model 8 Model 9
(Intercept) 0.082 0.082 0.075 0.122 0.139 0.146 0.156 0.154 0.140
(0.028) (0.076) (0.072) (0.074) (0.059) (0.056) (0.046) (0.051) (0.046)
sd__(Intercept) 0.449 0.298
sd__Observation 0.897 0.967
neighid sd__(Intercept) 0.376 0.376 0.272 0.246 0.173 0.189 0.165
schid sd__(Intercept) 0.275 0.273 0.207 0.196 0.152 0.175 0.148
Residual sd__Observation 0.894 0.893 0.872 0.897 0.874 0.871 0.874
male.facYes -0.097 -0.093 -0.103 -0.098 -0.096 -0.106
(0.039) (0.038) (0.038) (0.037) (0.044) (0.037)
dadocc 0.026 0.023 0.023 0.022
(0.002) (0.002) (0.002) (0.002)
deprive -0.468 -0.365 -0.364 -0.375
(0.038) (0.035) (0.035) (0.035)
neighid cor__(Intercept).male.facYes -0.375
neighid sd__male.facYes 0.119
schid cor__(Intercept).male.facYes -0.602
schid sd__male.facYes 0.090
deprive × dadocc -0.009
(0.003)
AIC 6422.0 6448.2 6364.7 6360.6 6146.5 6238.1 6055.3 6061.9 6046.5
BIC 6439.2 6465.4 6387.7 6389.3 6181.0 6272.5 6095.6 6125.1 6092.5
Log.Lik. -3207.985 -3221.082 -3178.356 -3175.291 -3067.244 -3113.035 -3020.674 -3019.968 -3015.257
LS0tCnRpdGxlOiAnTXVsdGlsZXZlbCBNb2RlbGluZywgTW9kdWxlIDExOiBDcm9zcy1DbGFzc2lmaWVkIExpbmVhciBNaXhlZCBNb2RlbHMnCmF1dGhvcjogJ0RyLiBCcm9kYScKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKVGhpcyB3ZWVrLCB3ZSB3aWxsIGJlIHVzaW5nIEdlb3JnZSBMZWNraWUncyAoMjAwOSkgYHNjb3RsYW5kYCBkYXRhc2V0LCB3aGljaCBpbmNsdWRlcyAyLDMxMCBzdHVkZW50cywgbmVzdGVkIHdpdGhpbiAxNyBzY2hvb2xzIGFuZCBzZXBhcmF0ZWx5IG5lc3RlZCB3aXRoaW4gNTAwKyBuZWlnaGJvcmhvb2RzLgoKIyBMb2FkIGluIE91ciBNVlAgUGFja2FnZXMKYGBge3J9CnN1cHByZXNzUGFja2FnZVN0YXJ0dXBNZXNzYWdlcyhsaWJyYXJ5KHRpZHl2ZXJzZSkpCmxpYnJhcnkobG1lNCkKbGlicmFyeShIbWlzYykKYGBgCgojIFBhcnQgMTogTG9hZGluZyBpbiB0aGUgRGF0YSwgYW5kIHRoZW4gUmVjb2RpbmcgVmFyaWFibGVzCiMjIExvYWQgaW4gdGhlIERhdGEKYGBge3J9CgpzY290bGFuZCA8LSBoYXZlbjo6cmVhZF9kdGEoInNjb3RsYW5kLmR0YSIpCgpnbGltcHNlKHNjb3RsYW5kKQpgYGAKIyMgVGhlIGBkZXNjcmliZWAgZnVuY3Rpb24gZnJvbSBgSG1pc2NgIGlzIGEgZ3JlYXQgY29kZWJvb2suLi4KYGBge3J9CmRlc2NyaWJlKHNjb3RsYW5kKQpgYGAKIyBQYXJ0IE9uZTogRGF0YSBFeHBsb3JhdGlvbiBhbmQgQ29kaW5nLyBMYWJlbGxpbmcKIyMgU3RhcnQgd2l0aCB2YXJpYWJsZSBsYWJlbGxpbmcgdXNpbmcgdGhlIGBsYWJlbGAgZnVuY3Rpb24KYGBge3J9CmxhYmVsKHNjb3RsYW5kJG5laWdoaWQpPSJOZWlnaGJvcmhvb2QgSUQiCmxhYmVsKHNjb3RsYW5kJHNjaGlkKT0iU2Nob29sIElEIgpsYWJlbChzY290bGFuZCRhdHRhaW4pPSJTdHVkZW50IE1lYXN1cmUgb2YgRWR1Y2F0aW9uYWwgQXR0YWlubWVudCIKbGFiZWwoc2NvdGxhbmQkcDd2cnEpPSJQcmltYXJ5IDcgVmVyYmFsIFJlYXNvbmluZyBRdW90aWVudCIKbGFiZWwoc2NvdGxhbmQkcDdyZWFkKT0iUHJpbWFyeSA3IFJlYWRpbmcgVGVzdCBTY29yZXMiCmxhYmVsKHNjb3RsYW5kJGRhZG9jYyk9IlNjaG9vbCBNZWFuIG9mIERhZCdzIE9jY3VwYXRpb24gU2NvcmUgb24gSG9wZS1Hb2xkdGhvcnBlIFNjYWxlIgpsYWJlbChzY290bGFuZCRkYWR1bmVtcCk9IkRhZCBDdXJyZW50bHkgVW5lbXBsb3llZCIKbGFiZWwoc2NvdGxhbmQkZGFkZWQpPSJEYWQgU2Nob29sIEFmdGVyIEFnZSAxNSIKbGFiZWwoc2NvdGxhbmQkbW9tZWQpPSJNb20gU2Nob29sIEFmdGVyIEFnZSAxNSIKbGFiZWwoc2NvdGxhbmQkbWFsZSk9Ik1hbGUiCmxhYmVsKHNjb3RsYW5kJGRlcHJpdmUpPSJOZWlnaGJvcmhvb2QgRGVwcml2YXRpb24gU2NvcmUgKFBvdmVydHksIEhlYWx0aCwgYW5kIEhvdXNpbmcpIgoKYGBgCgoKYGBge3J9CnNjb3QuY2xlYW4gPC0gc2NvdGxhbmQgJT4lCiAgbXV0YXRlKC4sIAogICAgICAgICBkYWR1bmVtcC5mYWMgPSBhc19mYWN0b3IoZGFkdW5lbXApLAogICAgICAgICBkYWRlZC5mYWMgPSBhc19mYWN0b3IoZGFkZWQpLAogICAgICAgICBtb21lZC5mYWMgPSBhc19mYWN0b3IobW9tZWQpLAogICAgICAgICBtYWxlLmZhYyA9IGFzX2ZhY3RvcihtYWxlKSkKCmxldmVscyhzY290LmNsZWFuJGRhZHVuZW1wLmZhYykgPSBjKCJObyIsIlllcyIpCmxldmVscyhzY290LmNsZWFuJGRhZGVkLmZhYykgPSBjKCJObyIsIlllcyIpCmxldmVscyhzY290LmNsZWFuJG1vbWVkLmZhYykgPSBjKCJObyIsIlllcyIpCmxldmVscyhzY290LmNsZWFuJG1hbGUuZmFjKSA9IGMoIk5vIiwiWWVzIikKCmdsaW1wc2Uoc2NvdC5jbGVhbikKYGBgCgpgYGB7cn0KZGVzY3JpYmUoc2NvdC5jbGVhbikKYGBgCgojIyBFeHBsb3JlIG91ciBPdXRjb21lLCBgYXR0YWluYApgYGB7cn0KZ2dwbG90KGRhdGEgPSBzY290LmNsZWFuLCBtYXBwaW5nID0gYWVzKHggPSBhdHRhaW4pKSArCiAgZ2VvbV9oaXN0b2dyYW0oYmlucyA9IDIwKSArIAogIGxhYnModGl0bGUgPSAiRGlzdHJpYnV0aW9uIG9mIEVkdWNhdGlvbmFsIEF0dGFpbm1lbnQgU2NvcmVzIiwKICAgICAgIHggPSAiU3RhbmRhcmRpemVkIEVkdWNhdGlvbmFsIEF0dGFpbm1lbnQgU2NvcmUiKSArCiAgdGhlbWVfbWluaW1hbCgpCmBgYAoKIyMgRXhwbG9yZSBhIGZldyBvdGhlciBwcmVkaWN0b3JzCgpgYGB7cn0KZ2dwbG90KGRhdGEgPSBzY290LmNsZWFuLCBtYXBwaW5nID0gYWVzKHggPSBkYWRvY2MpKSArCiAgZ2VvbV9oaXN0b2dyYW0oYmlucyA9IDMwKSArIAogIGxhYnModGl0bGUgPSAiRGlzdHJpYnV0aW9uIG9mIChTY2hvb2wgTWVhbikgb2YgRGFkJ3MgT2NjdXBhdGlvbiBTY29yZXMiLAogICAgICAgeCA9ICJTdGFuZGFyZGl6ZWQgRWR1Y2F0aW9uYWwgQXR0YWlubWVudCBTY29yZSIpICsKICB0aGVtZV9taW5pbWFsKCkKYGBgCgpgYGB7cn0KZ2dwbG90KGRhdGEgPSBzY290LmNsZWFuLCBtYXBwaW5nID0gYWVzKHggPSBkZXByaXZlKSkgKwogIGdlb21faGlzdG9ncmFtKGJpbnMgPSA1MCkgKyAKICBsYWJzKHRpdGxlID0gIkRpc3RyaWJ1dGlvbiBvZiBOZWlnaGJvcmhvb2QgQ29udGV4dCBTY29yZXMiLAogICAgICAgeCA9ICJTdGFuZGFyZGl6ZWQgTmVpZ2hib3Job29kIERlcHJpdmF0aW9uIFNjb3JlIikgKwogIHRoZW1lX21pbmltYWwoKQpgYGAKCiMgUGFydCBUd286IE1vZGVsaW5nCiMjIEZpcnN0LCBydW4gdHdvIHNlcGFyYXRlIG1vZGVscyB3aXRoIHNjaG9vbCBhbmQgbmVpZ2hib3Job29kIGFuZCB0aGUgTDIgY2x1c3RlcnMuLi4KCmBgYHtyfQptb2RlbC5udWxsLnNjaG9vbCA8LSBsbWVyKGF0dGFpbiB+ICgxfHNjaGlkKSwgUkVNTCA9IEZBTFNFLCBkYXRhID0gc2NvdC5jbGVhbikKc3VtbWFyeShtb2RlbC5udWxsLnNjaG9vbCkKYGBgCgojIyBDYWxjdWxhdGUgU2Nob29sIElDQwpgYGB7cn0KaWNjLnNjaG9vbCA8LSAwLjA4ODc0LygwLjA4ODc0ICsgMC45MzQ0MSkKaWNjLnNjaG9vbApgYGAKCmBgYHtyfQptb2RlbC5udWxsLm5laWdoIDwtIGxtZXIoYXR0YWluIH4gKDF8bmVpZ2hpZCksIFJFTUwgPSBGQUxTRSwgZGF0YSA9IHNjb3QuY2xlYW4pCnN1bW1hcnkobW9kZWwubnVsbC5uZWlnaCkKYGBgCgojIyBDYWxjdWxhdGUgTmVpZ2hib3Job29kIElDQwpgYGB7cn0KaWNjLm5laWdoIDwtIDAuMjAxNS8oMC4yMDE1ICsgMC44MDQ0KQppY2MubmVpZ2gKYGBgCgojIyBOb3csIGEgbnVsbCBhZGRpdGl2ZSBjcm9zcy1jbGFzc2lmaWVkIG1vZGVsLi4uCgpgYGB7cn0KbW9kZWwubnVsbC5jcm9zc2VkIDwtIGxtZXIoYXR0YWluIH4gKDF8c2NoaWQpICsgKDF8bmVpZ2hpZCksIFJFTUwgPSBGQUxTRSwgZGF0YSA9IHNjb3QuY2xlYW4pCnN1bW1hcnkobW9kZWwubnVsbC5jcm9zc2VkKQpgYGAKCiMjIENhbGN1bGF0ZSBOZWlnaGJvcmhvb2QgSUNDIGluIENyb3NzZWQgTW9kZWwKYGBge3J9CmljYy5uZWlnaC5jcm9zc2VkIDwtIDAuMTQxMjIvKDAuMTQxMjIgKyAwLjA3NTQ1ICsgMC43OTkwMikKaWNjLm5laWdoLmNyb3NzZWQKYGBgCgojIyBDYWxjdWxhdGUgU2Nob29sIElDQyBpbiBDcm9zc2VkIE1vZGVsCmBgYHtyfQppY2Muc2Nob29sLmNyb3NzZWQgPC0gMC4wNzU0NS8oMC4xNDEyMiArIDAuMDc1NDUgKyAwLjc5OTAyKQppY2Muc2Nob29sLmNyb3NzZWQKYGBgCgojIyBBZGQgYW4gTDEgUHJlZGljdG9yLCBNYWxlLi4uCgpgYGB7cn0KbW9kZWwuMSA8LSBsbWVyKGF0dGFpbiB+IG1hbGUuZmFjICsgKDF8c2NoaWQpICsgKDF8bmVpZ2hpZCksIFJFTUwgPSBGQUxTRSwgZGF0YSA9IHNjb3QuY2xlYW4pCnN1bW1hcnkobW9kZWwuMSkKYGBgCgojIyBMMSBQbHVzIFNjaG9vbC1MZXZlbCBQcmVkaWN0b3IgT25seSAoYGRhZG9jY2ApLi4uCgpgYGB7cn0KbW9kZWwuMiA8LSBsbWVyKGF0dGFpbiB+IG1hbGUuZmFjICsgZGFkb2NjICsgKDF8c2NoaWQpICsgKDF8bmVpZ2hpZCksIFJFTUwgPSBGQUxTRSwgZGF0YSA9IHNjb3QuY2xlYW4pCnN1bW1hcnkobW9kZWwuMikKYGBgCgojIyBMMSBQbHVzIE5laWdoYm9yaG9vZC1MZXZlbCBQcmVkaWN0b3JzIE9ubHkgKGBkZXByaXZlYCkuLi4KCmBgYHtyfQptb2RlbC4zIDwtIGxtZXIoYXR0YWluIH4gbWFsZS5mYWMgKyBkZXByaXZlICsgKDF8c2NoaWQpICsgKDF8bmVpZ2hpZCksIFJFTUwgPSBGQUxTRSwgZGF0YSA9IHNjb3QuY2xlYW4pCnN1bW1hcnkobW9kZWwuMykKYGBgCgojIyBMMSBQbHVzIEJvdGggU2Nob29sIGFuZCBOZWlnaGJvcmhvb2QgUHJlZGljdG9ycy4uLgoKYGBge3J9Cm1vZGVsLjQgPC0gbG1lcihhdHRhaW4gfiBtYWxlLmZhYyArIGRlcHJpdmUgKyBkYWRvY2MgKyAoMXxzY2hpZCkgKyAoMXxuZWlnaGlkKSwgUkVNTCA9IEZBTFNFLCBkYXRhID0gc2NvdC5jbGVhbikKc3VtbWFyeShtb2RlbC40KQpgYGAKCiMjIFRlc3QgUmFuZG9tIGVmZmVjdHMgZm9yIGBtYWxlLmZhY2AgYXQgbmVpZ2hib3Job29kIGxldmVsIGFuZCBhdCBzY2hvb2wgbGV2ZWwuLi4KCmBgYHtyfQptb2RlbC41IDwtIGxtZXIoYXR0YWluIH4gbWFsZS5mYWMgKyBkZXByaXZlICsgZGFkb2NjICsgKG1hbGUuZmFjfHNjaGlkKSArIChtYWxlLmZhY3xuZWlnaGlkKSwgUkVNTCA9IEZBTFNFLCBkYXRhID0gc2NvdC5jbGVhbikKc3VtbWFyeShtb2RlbC41KQpgYGAKCmBgYHtyfQpsbWVyVGVzdDo6cmFuZChtb2RlbC41KQpgYGAKCiMjIEZpbmFsbHksIHRlc3QgYSBtb2RlbCB3aXRoIHNjaG9vbC1uZWlnaGJvcmhvb2QgaW50ZXJhY3Rpb24uLi4KYGBge3J9Cm1vZGVsLjYgPC0gbG1lcihhdHRhaW4gfiBtYWxlLmZhYyArIGRlcHJpdmUgKyBkYWRvY2MgKyBkZXByaXZlOmRhZG9jYyArICgxfHNjaGlkKSArICgxfG5laWdoaWQpLCBSRU1MID0gRkFMU0UsIGRhdGEgPSBzY290LmNsZWFuKQpzdW1tYXJ5KG1vZGVsLjYpCmBgYAojIyBVc2UgYGludGVycGxvdGAgdG8gdmlzdWFsaXplIHRoZSBpbnRlcmFjdGlvbiBiZXR3ZWVuIGBkZXByaXZlYCBhbmQgYGRhZG9jY2AuLi4KYGBge3J9CmludGVycGxvdDo6aW50ZXJwbG90KG1vZGVsLjYsIHZhcjEgPSAiZGVwcml2ZSIsIHZhcjIgPSAiZGFkb2NjIikKYGBgCgojIFN0ZXAgNjogQ29tcGFyZSBPcmlnaW5hbCAiRmluYWwiIE1vZGVsLCBSb2J1c3QgTW9kZWwsIGFuZCBUcmltbWVkIE1vZGVsIGZvciBEaWZmZXJlbmNlcyAoU2Vuc2l0aXZpdHkgQW5hbHlzaXMpCmBgYHtyfQpsaWJyYXJ5KG1vZGVsc3VtbWFyeSkKbGlicmFyeShicm9vbS5taXhlZCkKbW9kZWxzIDwtIGxpc3QobW9kZWwubnVsbC5uZWlnaCwgbW9kZWwubnVsbC5zY2hvb2wsIG1vZGVsLm51bGwuY3Jvc3NlZCwgbW9kZWwuMSwgbW9kZWwuMiwgbW9kZWwuMywgbW9kZWwuNCwgbW9kZWwuNSwgbW9kZWwuNikKbW9kZWxzdW1tYXJ5KG1vZGVscywgb3V0cHV0ID0gIm1hcmtkb3duIikKYGBgCgo=