# path to working directory
setwd("C:/Users/a0905/Documents/multi level")
# package management
library(pacman)
# load them
pacman::p_load(mlmRev, tidyverse, lme4, merTools, dplyr)
# input data
dta <- read.csv("C:/Users/a0905/Documents/multi level/pts.csv")
#
head(dta)
## School Teacher Pupil Read_0 Read_1
## 1 1 11 18 -0.16407000 4.2426
## 2 1 11 19 -0.98126000 1.0000
## 3 1 11 20 -1.25370000 2.2361
## 4 1 11 15 -0.87230000 3.1623
## 5 1 11 17 -0.00063104 3.4641
## 6 1 11 16 -0.92678000 4.0000
# coerce variables to factor type and
# compute mean Read_0 by school
dta <- dta %>%
mutate(School = factor(School), Teacher = factor(Teacher),
Pupil = factor(Pupil)) %>%
group_by( School ) %>%
mutate(msRead_0 = mean(Read_0))
#
# compute mean Read_0 by teacher and
# centered teacher mean of Read_0 (from respective school means)
#
dta <- dta %>%
group_by( Teacher ) %>%
mutate(mtRead_0 = mean(Read_0), ctRead_0 = mean(Read_0) - msRead_0 )
##
# no predictor
summary(m0 <- lmer(Read_1 ~ (1 | School) + (1 | Teacher), data = dta))
## Linear mixed model fit by REML ['lmerMod']
## Formula: Read_1 ~ (1 | School) + (1 | Teacher)
## Data: dta
##
## REML criterion at convergence: 2486.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.05438 -0.65686 0.06497 0.62523 2.47233
##
## Random effects:
## Groups Name Variance Std.Dev.
## Teacher (Intercept) 0.1370 0.3701
## School (Intercept) 0.1277 0.3574
## Residual 1.3250 1.1511
## Number of obs: 777, groups: Teacher, 46; School, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.3755 0.1065 31.69
# add a school-level predictor
summary(m0_s <- update(m0, . ~ . + msRead_0))
## boundary (singular) fit: see help('isSingular')
## Linear mixed model fit by REML ['lmerMod']
## Formula: Read_1 ~ (1 | School) + (1 | Teacher) + msRead_0
## Data: dta
##
## REML criterion at convergence: 2467.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.06110 -0.65597 0.06822 0.62772 2.45391
##
## Random effects:
## Groups Name Variance Std.Dev.
## Teacher (Intercept) 0.1234 0.3513
## School (Intercept) 0.0000 0.0000
## Residual 1.3224 1.1500
## Number of obs: 777, groups: Teacher, 46; School, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.39822 0.06808 49.918
## msRead_0 1.00159 0.17807 5.625
##
## Correlation of Fixed Effects:
## (Intr)
## msRead_0 0.086
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# add a teacher-level predictor
summary(m0_t <- update(m0, . ~ . + mtRead_0))
## boundary (singular) fit: see help('isSingular')
## Linear mixed model fit by REML ['lmerMod']
## Formula: Read_1 ~ (1 | School) + (1 | Teacher) + mtRead_0
## Data: dta
##
## REML criterion at convergence: 2434.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.07539 -0.65979 0.06707 0.63259 2.49973
##
## Random effects:
## Groups Name Variance Std.Dev.
## Teacher (Intercept) 3.296e-15 5.741e-08
## School (Intercept) 4.471e-02 2.114e-01
## Residual 1.309e+00 1.144e+00
## Number of obs: 777, groups: Teacher, 46; School, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.40649 0.06307 54.010
## mtRead_0 1.07319 0.11495 9.336
##
## Correlation of Fixed Effects:
## (Intr)
## mtRead_0 0.024
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# add a teacher-level predictor away from respective school means
summary(m0_ct <- update(m0, . ~ . + ctRead_0))
## boundary (singular) fit: see help('isSingular')
## Linear mixed model fit by REML ['lmerMod']
## Formula: Read_1 ~ (1 | School) + (1 | Teacher) + ctRead_0
## Data: dta
##
## REML criterion at convergence: 2454.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1804 -0.6461 0.0791 0.6420 2.5272
##
## Random effects:
## Groups Name Variance Std.Dev.
## Teacher (Intercept) 2.462e-10 1.569e-05
## School (Intercept) 1.764e-01 4.200e-01
## Residual 1.311e+00 1.145e+00
## Number of obs: 777, groups: Teacher, 46; School, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.3896 0.1029 32.932
## ctRead_0 1.1742 0.1581 7.427
##
## Correlation of Fixed Effects:
## (Intr)
## ctRead_0 0.000
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
anova(m0, m0_s, m0_t, m0_ct)
## refitting model(s) with ML (instead of REML)
## Data: dta
## Models:
## m0: Read_1 ~ (1 | School) + (1 | Teacher)
## m0_s: Read_1 ~ (1 | School) + (1 | Teacher) + msRead_0
## m0_t: Read_1 ~ (1 | School) + (1 | Teacher) + mtRead_0
## m0_ct: Read_1 ~ (1 | School) + (1 | Teacher) + ctRead_0
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 4 2492.1 2510.7 -1242.0 2484.1
## m0_s 5 2472.2 2495.4 -1231.1 2462.2 21.927 1 2.833e-06 ***
## m0_t 5 2438.7 2461.9 -1214.3 2428.7 33.495 0
## m0_ct 5 2459.5 2482.8 -1224.7 2449.5 0.000 0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# result: m0_s fits best
###