R Notebook

library(tidyverse)
library(Hmisc)

Loading required package: lattice
Loading required package: survival
Loading required package: Formula
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

Attaching package: 㤼㸱Hmisc㤼㸲

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

    src, summarize

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

    format.pval, units

library(lme4)

scot.clean <- scotland %>%
  mutate(., 
         dadunemp.fac = as_factor(dadunemp),
         daded.fac = as_factor(daded),
         momed.fac = as_factor(momed),
         male.fac = as_factor(male))

levels(scot.clean$dadunemp.fac) = c("No","Yes")
levels(scot.clean$daded.fac) = c("No","Yes")
levels(scot.clean$momed.fac) = c("No","Yes")
levels(scot.clean$male.fac) = c("No","Yes")

glimpse(scot.clean)

Rows: 2,310
Columns: 15
$ neighid      <dbl> 26, 26, 26, 26, 26, 27, 29, 29, 29, 29, 29, 29,...
$ schid        <dbl> 20, 20, 20, 20, 20, 20, 18, 20, 20, 20, 20, 20,...
$ attain       <dbl> 1.5177, -1.3276, 0.5610, 1.5177, -1.3276, -0.13...
$ p7vrq        <dbl> 17.972, -10.028, 2.972, 1.972, -1.028, 3.972, 8...
$ p7read       <dbl> 17.134, -27.866, 6.134, 11.134, -0.866, -0.866,...
$ dadocc       <dbl> 16.196, -3.454, 2.316, -9.094, -3.454, -3.454, ...
$ dadunemp     <dbl> 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,...
$ momed        <dbl> 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
$ male         <dbl> 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1,...
$ deprive      <dbl> -0.551, -0.551, -0.551, -0.551, -0.551, 0.147, ...
$ dadunemp.fac <fct> No, No, No, No, No, No, No, No, No, No, No, No,...
$ 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, N...
$ male.fac     <fct> Yes, No, No, No, Yes, No, Yes, No, Yes, No, No,...

model.null.school <- lmer(p7read ~ (1|schid), REML = FALSE, data = scot.clean)
summary(model.null.school)

Linear mixed model fit by maximum likelihood . t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: p7read ~ (1 | schid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
 18558.8  18576.0  -9276.4  18552.8     2307 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7611 -0.6895 -0.0027  0.6854  2.5664 

Random effects:
 Groups   Name        Variance Std.Dev.
 schid    (Intercept)  18.75    4.33   
 Residual             176.64   13.29   
Number of obs: 2310, groups:  schid, 17

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)
(Intercept)  -0.2912     1.0951 16.1869  -0.266    0.794

icc.school <- 18.75/(18.75+176.64)
icc.school

[1] 0.09596192

model.null.neigh <- lmer(p7read ~ (1|neighid), REML = FALSE, data = scot.clean)
summary(model.null.neigh)

Linear mixed model fit by maximum likelihood . t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: p7read ~ (1 | neighid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
 18606.3  18623.5  -9300.1  18600.3     2307 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.65885 -0.67581 -0.01035  0.67042  2.51256 

Random effects:
 Groups   Name        Variance Std.Dev.
 neighid  (Intercept)  32.89    5.735  
 Residual             160.37   12.664  
Number of obs: 2310, groups:  neighid, 524

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)
(Intercept)  -0.2407     0.3816 454.9417  -0.631    0.529

icc.neighboorhood <- 32.89 / (32.89 +  160.37)
icc.neighboorhood

[1] 0.1701852

model.null.crossed <- lmer(p7read ~ (1|schid) + (1|neighid), REML = FALSE, data = scot.clean)
summary(model.null.crossed)

Linear mixed model fit by maximum likelihood . t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: p7read ~ (1 | schid) + (1 | neighid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
 18524.9  18547.8  -9258.4  18516.9     2306 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.61239 -0.66099  0.01158  0.66532  2.55698 

Random effects:
 Groups   Name        Variance Std.Dev.
 neighid  (Intercept)  17.21    4.148  
 schid    (Intercept)  16.59    4.074  
 Residual             160.36   12.663  
Number of obs: 2310, groups:  neighid, 524; schid, 17

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)
(Intercept)  -0.3701     1.0522 16.2037  -0.352     0.73

icc.school.crossed <- 16.59 / (16.59 + 17.21 + 160.36)
icc.school.crossed

[1] 0.08544499

icc.neighborhood.crossed <- 17.21 / (17.21 + 16.59 + 160.36)
icc.neighborhood.crossed

[1] 0.08863824

model.1 <- lmer(p7read ~ daded.fac + momed.fac + (1|schid) + (1|neighid), REML = FALSE, data = scot.clean)
summary(model.1)

Linear mixed model fit by maximum likelihood . t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: p7read ~ daded.fac + momed.fac + (1 | schid) + (1 | neighid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
 18424.4  18458.9  -9206.2  18412.4     2304 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.99961 -0.66409  0.01792  0.66336  2.63998 

Random effects:
 Groups   Name        Variance Std.Dev.
 neighid  (Intercept)  13.24    3.639  
 schid    (Intercept)  13.06    3.614  
 Residual             155.66   12.476  
Number of obs: 2310, groups:  neighid, 524; schid, 17

Fixed effects:
              Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)    -2.1132     0.9580   17.2968  -2.206 0.041199 *  
daded.facYes    5.3704     0.7504 2286.9432   7.157 1.11e-12 ***
momed.facYes    2.4043     0.7069 2272.7606   3.401 0.000683 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) ddd.fY
daded.facYs -0.083       
momed.facYs -0.110 -0.460

model.2 <- lmer(p7read ~ daded.fac + momed.fac + dadocc + (1|schid) + (1|neighid), REML = FALSE, data = scot.clean)
summary(model.2)

Linear mixed model fit by maximum likelihood . t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: p7read ~ daded.fac + momed.fac + dadocc + (1 | schid) + (1 |  
    neighid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
 18310.2  18350.5  -9148.1  18296.2     2303 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2651 -0.6765  0.0008  0.6506  2.6851 

Random effects:
 Groups   Name        Variance Std.Dev.
 neighid  (Intercept)   7.921   2.814  
 schid    (Intercept)   9.900   3.146  
 Residual             151.718  12.317  
Number of obs: 2310, groups:  neighid, 524; schid, 17

Fixed effects:
               Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)    -1.38112    0.84752   17.86409  -1.630  0.12069    
daded.facYes    3.47317    0.75559 2291.80157   4.597 4.53e-06 ***
momed.facYes    1.83651    0.69334 2283.63704   2.649  0.00813 ** 
dadocc          0.26808    0.02425 2279.12078  11.057  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) ddd.fY mmd.fY
daded.facYs -0.107              
momed.facYs -0.127 -0.428       
dadocc       0.078 -0.242 -0.079

model.3 <- lmer(p7read ~ daded.fac + momed.fac + dadocc + deprive + (1|schid) + (1|neighid), REML = FALSE, data = scot.clean)
summary(model.3)

Linear mixed model fit by maximum likelihood . t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: 
p7read ~ daded.fac + momed.fac + dadocc + deprive + (1 | schid) +  
    (1 | neighid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
 18253.9  18299.9  -9119.0  18237.9     2302 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.3799 -0.6834  0.0091  0.6557  3.1614 

Random effects:
 Groups   Name        Variance Std.Dev.
 neighid  (Intercept)   3.716   1.928  
 schid    (Intercept)   6.381   2.526  
 Residual             151.646  12.314  
Number of obs: 2310, groups:  neighid, 524; schid, 17

Fixed effects:
               Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)    -1.15892    0.70581   17.75224  -1.642   0.1182    
daded.facYes    3.31275    0.74800 2300.75393   4.429 9.92e-06 ***
momed.facYes    1.66801    0.68663 2294.65484   2.429   0.0152 *  
dadocc          0.23752    0.02435 2284.76701   9.752  < 2e-16 ***
deprive        -3.84492    0.48529  408.72238  -7.923 2.22e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) ddd.fY mmd.fY dadocc
daded.facYs -0.128                     
momed.facYs -0.151 -0.427              
dadocc       0.085 -0.232 -0.075       
deprive     -0.033  0.039  0.032  0.197

model.4 <- lmer(p7read ~ daded.fac + momed.fac + dadocc + deprive + momed.fac:daded.fac + (1|schid) + (1|neighid), REML = FALSE, data = scot.clean)
summary(model.4)

Linear mixed model fit by maximum likelihood . t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: 
p7read ~ daded.fac + momed.fac + dadocc + deprive + momed.fac:daded.fac +  
    (1 | schid) + (1 | neighid)
   Data: scot.clean

     AIC      BIC   logLik deviance df.resid 
 18247.3  18299.0  -9114.6  18229.3     2301 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.3143 -0.6715  0.0026  0.6539  3.1707 

Random effects:
 Groups   Name        Variance Std.Dev.
 neighid  (Intercept)   3.513   1.874  
 schid    (Intercept)   6.506   2.551  
 Residual             151.222  12.297  
Number of obs: 2310, groups:  neighid, 524; schid, 17

Fixed effects:
                            Estimate Std. Error         df t value
(Intercept)                 -1.37684    0.71429   18.13613  -1.928
daded.facYes                 5.34180    1.01514 2296.98834   5.262
momed.facYes                 3.12965    0.84626 2292.70713   3.698
dadocc                       0.23530    0.02432 2284.28187   9.676
deprive                     -3.75830    0.48407  410.53776  -7.764
daded.facYes:momed.facYes   -4.24203    1.43861 2272.83105  -2.949
                          Pr(>|t|)    
(Intercept)               0.069719 .  
daded.facYes              1.56e-07 ***
momed.facYes              0.000222 ***
dadocc                     < 2e-16 ***
deprive                   6.62e-14 ***
daded.facYes:momed.facYes 0.003224 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) ddd.fY mmd.fY dadocc depriv
daded.facYs -0.163                            
momed.facYs -0.181  0.143                     
dadocc       0.087 -0.192 -0.079              
deprive     -0.039  0.067  0.059  0.195       
ddd.fcYs:.Y  0.103 -0.678 -0.587  0.032 -0.057

interplot::interplot(model.4, var1 = "momed", var2 = "daded")

Error in interplot.lmerMod(model.4, var1 = "momed", var2 = "daded") : 
  The function does not support interactions between two factors.