# 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.0000Column 1: School ID (1-35)
Column 2: Teacher ID (11-352)
Column 3: Pupil ID (1-1948)
Column 4: Reading attainment at the end
of reception (mean = 0, SD = 1)
Column 5: Reading attainment at the
end of year one
# 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 )# DV: read_1; random intercept: School, Teacher
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.054 -0.657 0.065 0.625 2.472
Random effects:
Groups Name Variance Std.Dev.
Teacher (Intercept) 0.137 0.370
School (Intercept) 0.128 0.357
Residual 1.325 1.151
Number of obs: 777, groups: Teacher, 46; School, 20
Fixed effects:
Estimate Std. Error t value
(Intercept) 3.375 0.107 31.7結果顯示：
# add IV: msRead_0
summary(m0_s <- update(m0, . ~ . + msRead_0))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.0611 -0.6560 0.0682 0.6277 2.4539
Random effects:
Groups Name Variance Std.Dev.
Teacher (Intercept) 0.123 0.351
School (Intercept) 0.000 0.000
Residual 1.322 1.150
Number of obs: 777, groups: Teacher, 46; School, 20
Fixed effects:
Estimate Std. Error t value
(Intercept) 3.3982 0.0681 49.92
msRead_0 1.0016 0.1781 5.62
Correlation of Fixed Effects:
(Intr)
msRead_0 0.086
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')結果顯示：
# add IV: mtRead_0
summary(m0_t <- update(m0, . ~ . + mtRead_0))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.0754 -0.6598 0.0671 0.6326 2.4997
Random effects:
Groups Name Variance Std.Dev.
Teacher (Intercept) 3.30e-15 5.74e-08
School (Intercept) 4.47e-02 2.11e-01
Residual 1.31e+00 1.14e+00
Number of obs: 777, groups: Teacher, 46; School, 20
Fixed effects:
Estimate Std. Error t value
(Intercept) 3.4065 0.0631 54.01
mtRead_0 1.0732 0.1149 9.34
Correlation of Fixed Effects:
(Intr)
mtRead_0 0.024
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')sjPlot::tab_model(m0_t, show.p=FALSE, show.r2=FALSE)| Read_1 | ||
|---|---|---|
| Predictors | Estimates | CI |
| (Intercept) | 3.41 | 3.28 – 3.53 |
| mtRead 0 | 1.07 | 0.85 – 1.30 |
| Random Effects | ||
| σ2 | 1.31 | |
| τ00 Teacher | 0.00 | |
| τ00 School | 0.04 | |
| N School | 20 | |
| N Teacher | 46 | |
| Observations | 777 | |
結果顯示：
# add IV: ctRead_0
summary(m0_ct <- update(m0, . ~ . + ctRead_0))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.180 -0.646 0.079 0.642 2.527
Random effects:
Groups Name Variance Std.Dev.
Teacher (Intercept) 2.46e-10 1.57e-05
School (Intercept) 1.76e-01 4.20e-01
Residual 1.31e+00 1.14e+00
Number of obs: 777, groups: Teacher, 46; School, 20
Fixed effects:
Estimate Std. Error t value
(Intercept) 3.390 0.103 32.93
ctRead_0 1.174 0.158 7.43
Correlation of Fixed Effects:
(Intr)
ctRead_0 0.000
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')結果顯示：