PH6205 R Tutorial 5
Haoran Xue
(Last Updated March 17, 2025)
1 Reading Achievement Example
We analyze the Reading Achievement data. First we load the data into R:
Achieve = read.csv("C:/Users/haoraxue/Dropbox/PH6205/Slides5/ReadAchievement.csv")
head(Achieve)
## school senroll geread gevocab gender age
## 1 767 4.63 3.5 3.1 2 8.666667
## 2 767 4.63 1.2 2.8 2 8.833333
## 3 767 4.63 2.1 1.7 2 9.333333
## 4 767 4.63 1.6 2.1 2 9.083333
## 5 767 4.63 3.7 2.4 1 8.916667
## 6 767 4.63 2.4 2.4 1 9.000000
1.1 Model with Random Intercept, No Covariate
We use the function lmer() from R package lme4 to fit our models. We first fit the model with random
intercept, and no covariate is used.
library(lme4)
## Warning: package 'lme4' was built under R version 4.3.2
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 4.3.2
Model0 = lmer(geread~1 +(1|school), data=Achieve)
summary(Model0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: geread ~ 1 + (1 | school)
## Data: Achieve
##
## REML criterion at convergence: 46268.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3229 -0.6378 -0.2138 0.2850 3.8812
##
## Random effects:
## Groups Name Variance Std.Dev.
## school (Intercept) 0.3915 0.6257
## Residual 5.0450 2.2461
1
## Number of obs: 10320, groups: school, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.30675 0.05498 78.34
1.2 Model with Random Intercept, Fixed Slope
Next we use random intercept, and add the covariate gevocab with fixed slope.
Model1 = lmer(geread~1 + gevocab + (1|school), data=Achieve)
summary(Model1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: geread ~ 1 + gevocab + (1 | school)
## Data: Achieve
##
## REML criterion at convergence: 43137.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0823 -0.5735 -0.2103 0.3207 4.4334
##
## Random effects:
## Groups Name Variance Std.Dev.
## school (Intercept) 0.09978 0.3159
## Residual 3.76647 1.9407
## Number of obs: 10320, groups: school, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.023356 0.049309 41.03
## gevocab 0.512898 0.008373 61.26
##
## Correlation of Fixed Effects:
## (Intr)
## gevocab -0.758
1.3 Model with Random Intercept, Random Slope
Next we use random intercept, and also use random slope for gevocab.
Model2 =
lmer(geread~1+gevocab + (1+gevocab|school),
data=Achieve,
control = lmerControl(optimizer= "nlminbwrap"))
summary(Model2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: geread ~ 1 + gevocab + (1 + gevocab | school)
## Data: Achieve
## Control: lmerControl(optimizer = "nlminbwrap")
##
## REML criterion at convergence: 42992.9
##
## Scaled residuals:
2
## Min 1Q Median 3Q Max
## -3.7102 -0.5674 -0.2074 0.3176 4.6775
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## school (Intercept) 0.2827 0.5317
## gevocab 0.0193 0.1389 -0.86
## Residual 3.6659 1.9147
## Number of obs: 10320, groups: school, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.00570 0.06109 32.83
## gevocab 0.52036 0.01442 36.10
##
## Correlation of Fixed Effects:
## (Intr)
## gevocab -0.867
Note that we add the argument control = lmerControl(optimizer= "nlminbwrap"), this is for the con-
vergence of the underlying optimization procedure.
1.4 Model with Random Intercept, Random Slope, and Level 2 Covariate
Next we use random intercept, fixed slope for age, random slope for gevocab, and add the Level 2 covariate
senroll.
Model3 =
lmer(geread ~ 1 + age + senroll + gevocab + senroll:gevocab + (1+gevocab|school),
data = Achieve,
control = lmerControl(optimizer= "nlminbwrap"))
summary(Model3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: geread ~ 1 + age + senroll + gevocab + senroll:gevocab + (1 +
## gevocab | school)
## Data: Achieve
## Control: lmerControl(optimizer = "nlminbwrap")
##
## REML criterion at convergence: 43004.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6875 -0.5680 -0.2064 0.3176 4.6752
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## school (Intercept) 0.27885 0.5281
## gevocab 0.01924 0.1387 -0.86
## Residual 3.66514 1.9145
## Number of obs: 10320, groups: school, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.881479 0.459986 6.264
3
## age -0.107987 0.045564 -2.370
## senroll 0.019101 0.039252 0.487
## gevocab 0.542388 0.048250 11.241
## senroll:gevocab -0.004612 0.009173 -0.503
##
## Correlation of Fixed Effects:
## (Intr) age senrll gevocb
## age -0.891
## senroll -0.432 -0.002
## gevocab -0.400 0.007 0.834
## senrll:gvcb 0.371 0.003 -0.868 -0.954
