# package management
library(pacman)
# load them
pacman::p_load(mlmRev, tidyverse, lme4, merTools)
# input data
dta <- read.csv("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.05439 -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| Read_1 | ||
|---|---|---|
| Predictors | Estimates | CI |
| (Intercept) | 3.38 | 3.17 – 3.58 |
| Random Effects | ||
| σ2 | 1.32 | |
| τ00 Teacher | 0.14 | |
| τ00 School | 0.13 | |
| ICC | 0.17 | |
| N School | 20 | |
| N Teacher | 46 | |
| Observations | 777 | |
The estimate reading attainment score at the end of year is 3.38
Intraclass correlation coefficient is 0.17
## boundary (singular) fit: see ?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
## convergence code: 0
## boundary (singular) fit: see ?isSingular| Read_1 | Read_1 | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | Estimates | std. Error |
| (Intercept) | 3.38 | 0.11 | 3.40 | 0.07 |
| msRead_0 | 1.00 | 0.18 | ||
| Random Effects | ||||
| σ2 | 1.32 | 1.32 | ||
| τ00 | 0.14 Teacher | 0.12 Teacher | ||
| 0.13 School | 0.00 School | |||
| ICC | 0.17 | |||
Add a school-level predictor, we can find τ00 of school become 0.
## boundary (singular) fit: see ?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.65980 0.06707 0.63259 2.49973
##
## Random effects:
## Groups Name Variance Std.Dev.
## Teacher (Intercept) 0.00000 0.0000
## School (Intercept) 0.04471 0.2114
## Residual 1.30860 1.1439
## 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
## convergence code: 0
## boundary (singular) fit: see ?isSingular| Read_1 | Read_1 | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | Estimates | std. Error |
| (Intercept) | 3.38 | 0.11 | 3.41 | 0.06 |
| mtRead_0 | 1.07 | 0.11 | ||
| Random Effects | ||||
| σ2 | 1.32 | 1.31 | ||
| τ00 | 0.14 Teacher | 0.00 Teacher | ||
| 0.13 School | 0.04 School | |||
| ICC | 0.17 | |||
Add a teacher-level predictor, we can find τ00 of teacher become 0.
# add a teacher-level predictor away from respective school means
summary(m0_ct <- update(m0, . ~ . + ctRead_0))## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00376183 (tol = 0.002, component 1)## 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) 7.227e-08 0.0002688
## School (Intercept) 1.763e-01 0.4199281
## Residual 1.311e+00 1.1449203
## Number of obs: 777, groups: Teacher, 46; School, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.3896 0.1029 32.935
## ctRead_0 1.1742 0.1581 7.427
##
## Correlation of Fixed Effects:
## (Intr)
## ctRead_0 0.000
## convergence code: 0
## Model failed to converge with max|grad| = 0.00376183 (tol = 0.002, component 1)| Read_1 | Read_1 | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | Estimates | std. Error |
| (Intercept) | 3.38 | 0.11 | 3.39 | 0.10 |
| ctRead_0 | 1.17 | 0.16 | ||
| Random Effects | ||||
| σ2 | 1.32 | 1.31 | ||
| τ00 | 0.14 Teacher | 0.00 Teacher | ||
| 0.13 School | 0.18 School | |||
| ICC | 0.17 | 0.12 | ||
add a teacher-level predictor away from respective school means,the τ00 of teacher also become 0.