Multilevel Models: In-class exercises
Hao-Lun Fu
2020-10-12
1 In-class exercises 1:
Replicate the results of analysis reported the high school and beyond - null 2 3 4 examples.
1.1 Load data file
pacman::p_load(mlmRev, tidyverse, lme4)
# data from mlmRev package
data(Hsb82, "mlmRev")
#
str(Hsb82)'data.frame': 7185 obs. of 8 variables:
$ school : Ord.factor w/ 160 levels "8367"<"8854"<..: 59 59 59 59 59 59 59 59 59 59 ...
$ minrty : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
$ sx : Factor w/ 2 levels "Male","Female": 2 2 1 1 1 1 2 1 2 1 ...
$ ses : num -1.528 -0.588 -0.528 -0.668 -0.158 ...
$ mAch : num 5.88 19.71 20.35 8.78 17.9 ...
$ meanses: num -0.434 -0.434 -0.434 -0.434 -0.434 ...
$ sector : Factor w/ 2 levels "Public","Catholic": 1 1 1 1 1 1 1 1 1 1 ...
$ cses : num -1.0936 -0.1536 -0.0936 -0.2336 0.2764 ...# make a copy of data to the dta object
dta <- Hsb82 # use sensible variable names
names(dta) <- c("School", "Minority", "Gender", "SES", "Math",
"Ave_SES", "Sector", "C_SES")
#
head(dta) School Minority Gender SES Math Ave_SES Sector C_SES
1 1224 No Female -1.528 5.876 -0.43438 Public -1.093617
2 1224 No Female -0.588 19.708 -0.43438 Public -0.153617
3 1224 No Male -0.528 20.349 -0.43438 Public -0.093617
4 1224 No Male -0.668 8.781 -0.43438 Public -0.233617
5 1224 No Male -0.158 17.898 -0.43438 Public 0.276383
6 1224 No Male 0.022 4.583 -0.43438 Public 0.4563831.2 Histograms
ot <- theme_set(theme_bw())
# freedman-diaconis rule for binwidth
bw <- 2*IQR(dta$Math)/(dim(dta)[1]^(1/3))
# histograms
ggplot(dta, aes(Math)) +
stat_bin(aes(fill = ..count..), binwidth=bw) +
labs(x="Math achievement scores")
1.3 Box plots for a sample of 60 schools
## seed the random number to reproduce the sample
set.seed(20160923)
# pick 60 out of 160
n <- sample(160, 60)
# box plots for a sample of 60 schools
ggplot(filter(dta, School %in% levels(School)[n]),
aes(reorder(School, Math, mean), Math, color=Sector)) +
geom_boxplot() +
geom_point(size= rel(0.6)) +
labs(x="School ID", y="Math achievement score")+
theme(axis.text.x=element_text(size=6))+
ggtitle("A sample of 60 schools")
1.4
# math mean by school
dta_mmath <- dta %>%
group_by(School) %>%
mutate( ave_math = mean(Math)) %>%
select( School, ave_math ) %>%
unique %>%
as.data.frame
summary(dta_mmath$ave_math) Min. 1st Qu. Median Mean 3rd Qu. Max.
4.24 10.47 12.90 12.62 14.65 19.72 sd(dta_mmath$ave_math)[1] 3.1177# grand mean
summary(dta$Math) Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.83 7.28 13.13 12.75 18.32 24.99 # sd for overall Math scores
sd(dta$Math)[1] 6.8782# freedman-diaconis rule for binwidth
bw <- 2*IQR(dta_mmath$ave_math)/(160^(1/3))
#
ggplot(dta_mmath, aes(ave_math)) +
geom_histogram(aes(y= ..density.., alpha= .2), binwidth=bw) +
geom_density(fill="NA") +
geom_vline(xintercept=mean(dta_mmath$ave_math), linetype="longdash") +
geom_text(aes(7, 0.17, label= "ave = 12.62, sd = 3.12")) +
labs(x="Math achievement school means", y="Density") +
theme(legend.position="NONE")
1.5 baseline model
summary(m0 <- lmer(Math ~ (1 | School), data = dta))Linear mixed model fit by REML ['lmerMod']
Formula: Math ~ (1 | School)
Data: dta
REML criterion at convergence: 47117
Scaled residuals:
Min 1Q Median 3Q Max
-3.0631 -0.7539 0.0267 0.7606 2.7426
Random effects:
Groups Name Variance Std.Dev.
School (Intercept) 8.61 2.93
Residual 39.15 6.26
Number of obs: 7185, groups: School, 160
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.637 0.244 51.7The variance of school is 8.61 with 2.93 SD.
1.6 CIs for model parameters
confint(m0) 2.5 % 97.5 %
.sig01 2.5947 3.3159
.sigma 6.1548 6.3618
(Intercept) 12.1563 13.11711.7 default residual plot
plot(m0, xlab="Fitted values", ylab="Pearson residuals",
pch=20, cex=.5, type = c("p","g"))
pacman::p_load(mlmRev, tidyverse, lme4)
# data from mlmRev package
data(Hsb82)
# first 6 lines
head(Hsb82) school minrty sx ses mAch meanses sector cses
1 1224 No Female -1.528 5.876 -0.43438 Public -1.093617
2 1224 No Female -0.588 19.708 -0.43438 Public -0.153617
3 1224 No Male -0.528 20.349 -0.43438 Public -0.093617
4 1224 No Male -0.668 8.781 -0.43438 Public -0.233617
5 1224 No Male -0.158 17.898 -0.43438 Public 0.276383
6 1224 No Male 0.022 4.583 -0.43438 Public 0.456383# make a copy of data to the dta object
dta <- Hsb82
# names for human
names(dta) <- c("School", "Minority", "Gender", "SES", "Math",
"Ave_SES", "Sector", "C_SES")
#
str(dta)'data.frame': 7185 obs. of 8 variables:
$ School : Ord.factor w/ 160 levels "8367"<"8854"<..: 59 59 59 59 59 59 59 59 59 59 ...
$ Minority: Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
$ Gender : Factor w/ 2 levels "Male","Female": 2 2 1 1 1 1 2 1 2 1 ...
$ SES : num -1.528 -0.588 -0.528 -0.668 -0.158 ...
$ Math : num 5.88 19.71 20.35 8.78 17.9 ...
$ Ave_SES : num -0.434 -0.434 -0.434 -0.434 -0.434 ...
$ Sector : Factor w/ 2 levels "Public","Catholic": 1 1 1 1 1 1 1 1 1 1 ...
$ C_SES : num -1.0936 -0.1536 -0.0936 -0.2336 0.2764 ...1.8 random effects anova - intercept only model
summary(m1 <- lmer(Math ~ (1 | School), data = dta))Linear mixed model fit by REML ['lmerMod']
Formula: Math ~ (1 | School)
Data: dta
REML criterion at convergence: 47117
Scaled residuals:
Min 1Q Median 3Q Max
-3.0631 -0.7539 0.0267 0.7606 2.7426
Random effects:
Groups Name Variance Std.Dev.
School (Intercept) 8.61 2.93
Residual 39.15 6.26
Number of obs: 7185, groups: School, 160
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.637 0.244 51.7The variance of school is 8.61 with 2.93 SD.
1.9 random intercept; fixed-effect is MeanSES
summary(m2 <- lmer(Math ~ Ave_SES + (1 | School), data = dta))Linear mixed model fit by REML ['lmerMod']
Formula: Math ~ Ave_SES + (1 | School)
Data: dta
REML criterion at convergence: 46961
Scaled residuals:
Min 1Q Median 3Q Max
-3.1349 -0.7525 0.0241 0.7677 2.7852
Random effects:
Groups Name Variance Std.Dev.
School (Intercept) 2.64 1.62
Residual 39.16 6.26
Number of obs: 7185, groups: School, 160
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.685 0.149 85.0
Ave_SES 5.864 0.361 16.2
Correlation of Fixed Effects:
(Intr)
Ave_SES 0.010 The variance of school is 2.64 with 1.62 SD and the estimate of AVE_SES is 12.685 and SE is 0.149.
1.10 Residual plot
plot(m2, xlab = "Fitted values", ylab = "Pearson residuals",
main = "Model 2", pch = 20, cex = .5, type = c("p", "g"))
1.11 one regression for all
summary(m0a <- lm(Math ~ Ave_SES, data = dta))
Call:
lm(formula = Math ~ Ave_SES, data = dta)
Residuals:
Min 1Q Median 3Q Max
-19.452 -4.833 0.174 5.112 16.279
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.7470 0.0762 167 <2e-16
Ave_SES 5.7169 0.1843 31 <2e-16
Residual standard error: 6.46 on 7183 degrees of freedom
Multiple R-squared: 0.118, Adjusted R-squared: 0.118
F-statistic: 962 on 1 and 7183 DF, p-value: <2e-16The estimate of AVE_SES is 5.7169 and SE is 0.1843.
1.12 individual math scores vs mean school ses
ggplot(dta, aes(Ave_SES, Math)) +
geom_smooth(method = "lm") +
geom_point(alpha = .2) +
labs(x = "School mean SES", y = "Math achievement score")
1.13 math mean by school
dta_mmath <- dta %>%
group_by(School) %>%
mutate( ave_math = mean(Math)) %>%
select( School, Ave_SES, ave_math ) %>%
unique %>%
as.data.frame
# regression of mean math over mean ses by school
ggplot(dta_mmath, aes(Ave_SES, ave_math)) +
geom_smooth(method = "lm") +
geom_point(alpha = .2)+
labs(x = "Mean School SES", y = "Mean School Math Scores")
# see the random number generator
set.seed(20160923)
# take a sample of size 31 from total number of schools
n <- sample(length(levels(dta$School)), 31)
#
ggplot(data = filter(dta, School %in% levels(School)[n]),
aes(C_SES, Math, color = Sector)) +
geom_smooth(method = "lm", se = F, lwd = .6) +
geom_point(size = rel(0.5)) +
facet_wrap(~ School, nrow = 2) +
labs(x = "Centered SES", y="Math achievement score")+
theme(axis.text.x = element_text(size = 6))
1.14 one regression line per school
summary(m1b <- lm(Math ~ School / C_SES - 1, data = dta))
Call:
lm(formula = Math ~ School/C_SES - 1, data = dta)
Residuals:
Min 1Q Median 3Q Max
-19.266 -4.324 0.116 4.443 18.362
Coefficients:
Estimate Std. Error t value Pr(>|t|)
School8367 4.5528 1.6195 2.81 0.00495
School8854 4.2398 1.0712 3.96 7.6e-05
School4458 5.8114 0.8746 6.64 3.3e-11
School5762 4.3249 0.9962 4.34 1.4e-05
School6990 5.9768 0.8324 7.18 7.7e-13
School5815 7.2714 1.2119 6.00 2.1e-09
School7172 8.0668 0.9135 8.83 < 2e-16
School4868 12.3102 1.0392 11.85 < 2e-16
School7341 9.7942 0.8485 11.54 < 2e-16
School1358 11.2062 1.1063 10.13 < 2e-16
School4383 11.4657 1.2119 9.46 < 2e-16
School2305 11.1378 0.7403 15.04 < 2e-16
School8800 7.3359 1.0712 6.85 8.1e-12
School3088 9.1458 0.9703 9.43 < 2e-16
School8775 9.4670 0.8746 10.82 < 2e-16
School7890 8.3411 0.8485 9.83 < 2e-16
School6144 8.5451 0.9241 9.25 < 2e-16
School6443 9.4755 1.1063 8.56 < 2e-16
School5192 10.4095 1.1452 9.09 < 2e-16
School6808 9.2861 0.9135 10.17 < 2e-16
School2818 13.8728 0.9350 14.84 < 2e-16
School9340 11.1786 1.1253 9.93 < 2e-16
School4523 8.3517 0.8839 9.45 < 2e-16
School6816 14.5382 0.8171 17.79 < 2e-16
School2277 9.2976 0.7759 11.98 < 2e-16
School8009 14.0847 0.8839 15.93 < 2e-16
School5783 13.1503 1.1253 11.69 < 2e-16
School3013 12.6108 0.8324 15.15 < 2e-16
School7101 11.8499 1.1452 10.35 < 2e-16
School4530 9.0557 0.7635 11.86 < 2e-16
School9021 14.6967 0.8098 18.15 < 2e-16
School4511 13.4090 0.7957 16.85 < 2e-16
School2639 6.6155 0.9350 7.08 1.6e-12
School3377 9.1867 0.9033 10.17 < 2e-16
School6578 11.9940 0.8098 14.81 < 2e-16
School9347 13.5388 0.8026 16.87 < 2e-16
School3705 10.3317 0.9033 11.44 < 2e-16
School3533 10.4090 0.8746 11.90 < 2e-16
School1296 7.6360 0.8746 8.73 < 2e-16
School4350 11.8554 1.0549 11.24 < 2e-16
School9397 10.3555 0.8839 11.72 < 2e-16
School4253 9.4129 0.7957 11.83 < 2e-16
School2655 12.3452 0.8403 14.69 < 2e-16
School7342 11.1664 0.7957 14.03 < 2e-16
School9292 10.2796 1.3902 7.39 1.6e-13
School3499 13.2765 0.9830 13.51 < 2e-16
School7364 14.1721 0.9135 15.51 < 2e-16
School8983 10.9920 0.8485 12.95 < 2e-16
School5650 14.2735 0.9033 15.80 < 2e-16
School2658 13.3962 0.9033 14.83 < 2e-16
School8188 12.7410 1.1063 11.52 < 2e-16
School4410 13.4730 0.9464 14.24 < 2e-16
School9508 13.5747 1.0243 13.25 < 2e-16
School8707 12.8839 0.8746 14.73 < 2e-16
School1499 7.6604 0.8324 9.20 < 2e-16
School8477 12.5222 0.9962 12.57 < 2e-16
School1288 13.5108 1.2119 11.15 < 2e-16
School6291 10.1073 1.0243 9.87 < 2e-16
School1224 9.7154 0.8839 10.99 < 2e-16
School4292 12.8644 0.7516 17.12 < 2e-16
School8857 15.2969 0.7575 20.19 < 2e-16
School3967 12.0351 0.8403 14.32 < 2e-16
School6415 11.8602 0.8246 14.38 < 2e-16
School1317 13.1777 0.8746 15.07 < 2e-16
School2629 14.9078 0.8026 18.57 < 2e-16
School4223 14.6226 0.9033 16.19 < 2e-16
School1462 10.4956 0.8026 13.08 < 2e-16
School9550 11.0891 1.1253 9.85 < 2e-16
School6464 7.0916 1.1253 6.30 3.1e-10
School4931 13.7908 0.7957 17.33 < 2e-16
School5937 16.7760 1.1253 14.91 < 2e-16
School7919 14.8500 0.9962 14.91 < 2e-16
School3716 10.3687 0.9464 10.96 < 2e-16
School1909 14.4233 1.1452 12.59 < 2e-16
School2651 11.0843 0.9830 11.28 < 2e-16
School2467 10.1475 0.8403 12.08 < 2e-16
School1374 9.7285 1.1452 8.50 < 2e-16
School6600 11.7039 0.8098 14.45 < 2e-16
School5667 13.7782 0.7759 17.76 < 2e-16
School5720 14.2823 0.8324 17.16 < 2e-16
School3498 16.3905 0.8324 19.69 < 2e-16
School3881 11.9492 0.9464 12.63 < 2e-16
School2995 9.5461 0.8935 10.68 < 2e-16
School5838 13.6896 1.0884 12.58 < 2e-16
School3688 14.6563 0.9241 15.86 < 2e-16
School9158 8.5452 0.8324 10.27 < 2e-16
School8946 10.3751 0.7957 13.04 < 2e-16
School7232 12.5426 0.8403 14.93 < 2e-16
School2917 7.9790 0.9241 8.63 < 2e-16
School6170 14.1810 1.3223 10.72 < 2e-16
School8165 16.4512 0.8657 19.00 < 2e-16
School9104 16.8321 0.8171 20.60 < 2e-16
School2030 12.0782 0.8839 13.66 < 2e-16
School8150 14.8524 0.9135 16.26 < 2e-16
School4042 14.3154 0.7575 18.90 < 2e-16
School8357 14.3819 1.1662 12.33 < 2e-16
School8531 13.5287 0.9464 14.30 < 2e-16
School6074 13.7791 0.8098 17.02 < 2e-16
School4420 13.8742 1.0712 12.95 < 2e-16
School1906 15.9832 0.8324 19.20 < 2e-16
School3992 14.6452 0.8324 17.59 < 2e-16
School3999 10.9440 0.8935 12.25 < 2e-16
School4173 12.7247 0.9135 13.93 < 2e-16
School4325 13.2400 0.8324 15.91 < 2e-16
School5761 11.1381 0.8403 13.25 < 2e-16
School6484 12.9124 1.0243 12.61 < 2e-16
School6897 15.0976 0.8657 17.44 < 2e-16
School7635 15.0655 0.8485 17.75 < 2e-16
School7734 10.5596 1.2919 8.17 3.5e-16
School8175 11.6981 1.0549 11.09 < 2e-16
School8874 12.0550 1.0100 11.94 < 2e-16
School9225 14.6673 1.0100 14.52 < 2e-16
School2458 13.9857 0.8026 17.42 < 2e-16
School3610 15.3550 0.7575 20.27 < 2e-16
School5640 13.1601 0.8026 16.40 < 2e-16
School3838 16.0628 0.8246 19.48 < 2e-16
School9359 15.2706 0.8324 18.35 < 2e-16
School2208 15.4047 0.7823 19.69 < 2e-16
School6089 15.5696 1.0549 14.76 < 2e-16
School1477 14.2285 0.7696 18.49 < 2e-16
School2768 10.8869 1.2119 8.98 < 2e-16
School3039 16.9639 1.3223 12.83 < 2e-16
School5819 12.1389 0.8570 14.16 < 2e-16
School6397 12.7961 0.7823 16.36 < 2e-16
School1308 16.2555 1.3550 12.00 < 2e-16
School1433 19.7191 1.0243 19.25 < 2e-16
School1436 18.1116 0.9135 19.83 < 2e-16
School1461 16.8426 1.0549 15.97 < 2e-16
School1637 7.0241 1.1662 6.02 1.8e-09
School1942 18.1109 1.1253 16.09 < 2e-16
School1946 12.9084 0.9703 13.30 < 2e-16
School2336 16.5177 0.8839 18.69 < 2e-16
School2526 17.0530 0.8026 21.25 < 2e-16
School2626 13.3966 0.9830 13.63 < 2e-16
School2755 16.4765 0.8839 18.64 < 2e-16
School2771 11.8441 0.8171 14.50 < 2e-16
School2990 18.4479 0.8746 21.09 < 2e-16
School3020 14.3953 0.7889 18.25 < 2e-16
School3152 13.2090 0.8403 15.72 < 2e-16
School3332 14.2782 0.9830 14.52 < 2e-16
School3351 11.4652 0.9703 11.82 < 2e-16
School3427 19.7156 0.8657 22.77 < 2e-16
School3657 9.5212 0.8485 11.22 < 2e-16
School4642 14.5990 0.7759 18.82 < 2e-16
School5404 15.4150 0.8026 19.21 < 2e-16
School5619 15.4162 0.7459 20.67 < 2e-16
School6366 15.6564 0.7957 19.68 < 2e-16
School6469 18.4557 0.8026 22.99 < 2e-16
School7011 13.8136 1.0549 13.10 < 2e-16
School7276 12.6794 0.8324 15.23 < 2e-16
School7332 14.6361 0.8746 16.73 < 2e-16
School7345 11.3386 0.8098 14.00 < 2e-16
School7688 18.4223 0.8246 22.34 < 2e-16
School7697 15.7218 1.0712 14.68 < 2e-16
School8193 16.2323 0.9241 17.57 < 2e-16
School8202 11.7124 1.0243 11.43 < 2e-16
School8627 10.8837 0.8324 13.08 < 2e-16
School8628 16.5284 0.7759 21.30 < 2e-16
School9198 19.0923 1.0884 17.54 < 2e-16
School9586 14.8637 0.7889 18.84 < 2e-16
School8367:C_SES 0.2504 2.2249 0.11 0.91040
School8854:C_SES 1.9388 1.3543 1.43 0.15229
School4458:C_SES 1.1318 1.3679 0.83 0.40803
School5762:C_SES -1.0141 1.9595 -0.52 0.60480
School6990:C_SES 0.9477 0.9817 0.97 0.33439
School5815:C_SES 3.0180 2.1648 1.39 0.16333
School7172:C_SES 0.9945 1.3661 0.73 0.46666
School4868:C_SES 1.2865 1.4857 0.87 0.38658
School7341:C_SES 1.7037 1.0053 1.69 0.09017
School1358:C_SES 5.0680 1.7117 2.96 0.00308
School4383:C_SES 6.1802 2.1160 2.92 0.00350
School2305:C_SES -0.7821 1.1376 -0.69 0.49180
School8800:C_SES 2.5681 1.4267 1.80 0.07189
School3088:C_SES 1.7913 1.2314 1.45 0.14581
School8775:C_SES 1.0015 1.3201 0.76 0.44811
School7890:C_SES -0.6560 1.4446 -0.45 0.64978
School6144:C_SES 2.7703 1.3016 2.13 0.03334
School6443:C_SES -0.7433 1.4603 -0.51 0.61073
School5192:C_SES 1.6035 1.4610 1.10 0.27245
School6808:C_SES 2.2761 1.1528 1.97 0.04838
School2818:C_SES 2.8024 1.5235 1.84 0.06589
School9340:C_SES 3.3095 2.0019 1.65 0.09833
School4523:C_SES 2.3808 1.2455 1.91 0.05597
School6816:C_SES 1.3527 1.1678 1.16 0.24677
School2277:C_SES -2.0150 1.1775 -1.71 0.08706
School8009:C_SES 1.5569 1.5867 0.98 0.32654
School5783:C_SES 3.1141 1.4940 2.08 0.03717
School3013:C_SES 3.8398 1.7509 2.19 0.02834
School7101:C_SES 1.2955 1.4284 0.91 0.36448
School4530:C_SES 1.6474 1.2393 1.33 0.18377
School9021:C_SES 2.5242 1.3948 1.81 0.07038
School4511:C_SES 0.0425 1.3807 0.03 0.97544
School2639:C_SES -0.6301 1.5297 -0.41 0.68042
School3377:C_SES -0.7469 1.3196 -0.57 0.57142
School6578:C_SES 2.3905 1.3073 1.83 0.06749
School9347:C_SES 2.6860 1.1805 2.28 0.02292
School3705:C_SES 1.1585 1.3305 0.87 0.38393
School3533:C_SES -0.3118 1.6447 -0.19 0.84966
School1296:C_SES 1.0760 1.3661 0.79 0.43095
School4350:C_SES 4.3719 1.5946 2.74 0.00613
School9397:C_SES 2.4464 1.4758 1.66 0.09743
School4253:C_SES -0.3995 1.1533 -0.35 0.72902
School2655:C_SES 5.2270 1.3536 3.86 0.00011
School7342:C_SES 1.0125 1.4210 0.71 0.47617
School9292:C_SES 0.7587 2.7481 0.28 0.78250
School3499:C_SES 0.9924 1.5818 0.63 0.53045
School7364:C_SES 0.2595 1.8283 0.14 0.88714
School8983:C_SES 1.3859 1.3430 1.03 0.30212
School5650:C_SES 0.6806 1.1746 0.58 0.56231
School2658:C_SES 2.6299 1.4268 1.84 0.06533
School8188:C_SES 4.3976 1.6556 2.66 0.00792
School4410:C_SES 2.7602 1.5570 1.77 0.07631
School9508:C_SES 3.9538 1.8218 2.17 0.03002
School8707:C_SES 3.3915 1.0990 3.09 0.00204
School1499:C_SES 3.6347 1.1643 3.12 0.00180
School8477:C_SES 3.8122 1.2555 3.04 0.00240
School1288:C_SES 3.2554 1.8482 1.76 0.07821
School6291:C_SES 3.9809 1.6935 2.35 0.01877
School1224:C_SES 2.5086 1.4243 1.76 0.07824
School4292:C_SES -0.1606 1.1633 -0.14 0.89019
School8857:C_SES 0.8022 1.4345 0.56 0.57602
School3967:C_SES 3.3111 1.3507 2.45 0.01426
School6415:C_SES 3.5301 1.2152 2.90 0.00369
School1317:C_SES 1.2739 1.5893 0.80 0.42284
School2629:C_SES 0.2223 1.1465 0.19 0.84622
School4223:C_SES 2.4866 1.5077 1.65 0.09914
School1462:C_SES -0.8288 1.3984 -0.59 0.55341
School9550:C_SES 3.8919 1.4594 2.67 0.00767
School6464:C_SES 1.0035 1.8246 0.55 0.58235
School4931:C_SES 0.9118 1.1760 0.78 0.43813
School5937:C_SES 1.0396 2.0053 0.52 0.60417
School7919:C_SES 3.9894 1.8818 2.12 0.03404
School3716:C_SES 5.8638 1.1809 4.97 7.0e-07
School1909:C_SES 2.8548 1.6372 1.74 0.08126
School2651:C_SES 4.8991 1.5563 3.15 0.00165
School2467:C_SES 3.1371 1.6583 1.89 0.05856
School1374:C_SES 3.8543 1.6684 2.31 0.02090
School6600:C_SES 4.7043 1.1649 4.04 5.4e-05
School5667:C_SES 3.5230 1.2407 2.84 0.00453
School5720:C_SES 2.4663 1.2652 1.95 0.05130
School3498:C_SES -0.1311 1.4571 -0.09 0.92832
School3881:C_SES 2.3907 1.7927 1.33 0.18239
School2995:C_SES 1.4323 1.0180 1.41 0.15949
School5838:C_SES 1.8531 1.7799 1.04 0.29787
School3688:C_SES 1.5367 1.7204 0.89 0.37175
School9158:C_SES 3.8612 1.1253 3.43 0.00060
School8946:C_SES 1.6905 1.0675 1.58 0.11333
School7232:C_SES 5.0016 1.4774 3.39 0.00071
School2917:C_SES 1.1359 1.0298 1.10 0.27008
School6170:C_SES 4.8118 1.6496 2.92 0.00355
School8165:C_SES 1.8022 1.3848 1.30 0.19315
School9104:C_SES 1.4940 1.8344 0.81 0.41543
School2030:C_SES 1.4120 1.3142 1.07 0.28268
School8150:C_SES -0.1857 1.5106 -0.12 0.90216
School4042:C_SES 1.6936 1.2050 1.41 0.15992
School8357:C_SES 2.6758 1.1572 2.31 0.02079
School8531:C_SES 3.3182 1.4029 2.37 0.01804
School6074:C_SES 1.5291 1.3029 1.17 0.24060
School4420:C_SES 2.9587 1.5888 1.86 0.06261
School1906:C_SES 2.1455 1.3695 1.57 0.11726
School3992:C_SES 0.5379 1.3880 0.39 0.69839
School3999:C_SES 3.5670 0.9884 3.61 0.00031
School4173:C_SES 3.3657 1.4596 2.31 0.02114
School4325:C_SES 2.7560 1.0327 2.67 0.00763
School5761:C_SES 3.1080 1.1914 2.61 0.00911
School6484:C_SES 0.6057 1.4935 0.41 0.68509
School6897:C_SES 3.5805 1.1748 3.05 0.00231
School7635:C_SES 2.4485 1.2693 1.93 0.05377
School7734:C_SES 6.0352 1.4847 4.06 4.9e-05
School8175:C_SES 1.6124 1.5433 1.04 0.29619
School8874:C_SES 4.0963 1.4351 2.85 0.00433
School9225:C_SES 2.8859 1.3964 2.07 0.03881
School2458:C_SES 2.9567 1.2299 2.40 0.01624
School3610:C_SES 2.9558 1.2087 2.45 0.01449
School5640:C_SES 3.8277 1.3889 2.76 0.00587
School3838:C_SES 0.5979 1.2339 0.48 0.62801
School9359:C_SES -0.8335 1.4668 -0.57 0.56990
School2208:C_SES 2.6366 1.3190 2.00 0.04565
School6089:C_SES 1.6925 1.5193 1.11 0.26532
School1477:C_SES 1.2306 1.1487 1.07 0.28408
School2768:C_SES 3.5123 1.3752 2.55 0.01067
School3039:C_SES 2.9557 1.9205 1.54 0.12384
School5819:C_SES 1.9725 1.4395 1.37 0.17066
School6397:C_SES 2.7590 1.0385 2.66 0.00791
School1308:C_SES 0.1260 2.8974 0.04 0.96531
School1433:C_SES 1.8543 1.7480 1.06 0.28881
School1436:C_SES 1.6006 1.6165 0.99 0.32215
School1461:C_SES 6.2665 1.5468 4.05 5.2e-05
School1637:C_SES 3.1168 1.5623 2.00 0.04608
School1942:C_SES 0.0894 2.1776 0.04 0.96726
School1946:C_SES 3.5858 1.4078 2.55 0.01088
School2336:C_SES 1.9050 1.3961 1.36 0.17246
School2526:C_SES 0.1595 1.3261 0.12 0.90426
School2626:C_SES 4.0997 1.7786 2.30 0.02120
School2755:C_SES 0.5605 1.4194 0.39 0.69295
School2771:C_SES 4.2682 1.6053 2.66 0.00786
School2990:C_SES 1.3245 1.3233 1.00 0.31688
School3020:C_SES 1.6537 1.2180 1.36 0.17459
School3152:C_SES 2.7682 1.0016 2.76 0.00573
School3332:C_SES 2.0310 1.7177 1.18 0.23711
School3351:C_SES 2.4550 1.1697 2.10 0.03587
School3427:C_SES -0.4882 1.3179 -0.37 0.71109
School3657:C_SES 3.7359 1.1005 3.39 0.00069
School4642:C_SES 3.2724 1.1734 2.79 0.00531
School5404:C_SES 1.2142 1.6071 0.76 0.44996
School5619:C_SES 5.2575 1.2584 4.18 3.0e-05
School6366:C_SES 1.5175 1.2153 1.25 0.21181
School6469:C_SES 1.7553 1.2927 1.36 0.17456
School7011:C_SES 5.0747 1.7259 2.94 0.00329
School7276:C_SES 3.7734 1.0636 3.55 0.00039
School7332:C_SES 2.4632 1.3228 1.86 0.06263
School7345:C_SES 4.2119 0.9895 4.26 2.1e-05
School7688:C_SES 0.1163 1.4747 0.08 0.93712
School7697:C_SES 3.1362 1.7744 1.77 0.07719
School8193:C_SES 2.3352 1.6248 1.44 0.15070
School8202:C_SES 3.7059 1.4761 2.51 0.01207
School8627:C_SES 1.8696 1.1874 1.57 0.11541
School8628:C_SES 1.2314 1.3489 0.91 0.36134
School9198:C_SES 2.6105 1.6990 1.54 0.12446
School9586:C_SES 1.6721 1.3373 1.25 0.21122
Residual standard error: 6.06 on 6865 degrees of freedom
Multiple R-squared: 0.833, Adjusted R-squared: 0.825
F-statistic: 107 on 320 and 6865 DF, p-value: <2e-16Almost of schools showed the effect on math, but the effect did not interact with C_SES.
1.15 collect regression estimates from individual regression by school
m1b_betas <- m1b %>%
coef %>%
matrix(ncol = 2) %>%
as.data.frame %>%
rename( b0 = V1, b1 = V2)
#
colMeans(m1b_betas) b0 b1
12.6208 2.2016 1.16 SD
apply(m1b_betas, 2, sd) b0 b1
3.1177 1.6310 1.17
cor(m1b_betas) b0 b1
b0 1.000000 0.013467
b1 0.013467 1.0000001.18 scatter plot of regression coefficient estimates
ggplot(m1b_betas, aes(b0, b1)) +
geom_point(pch = 20) +
stat_smooth(method = "lm", se = F) +
stat_ellipse() +
annotate("text", x = 5, y = 5, label = "r = 0.013") +
labs(x = "Intercept estimates", y = "Slope estimates")
1.19 random intercept only
summary(m3a <- lmer(Math ~ C_SES + (1 | School), data = dta))Linear mixed model fit by REML ['lmerMod']
Formula: Math ~ C_SES + (1 | School)
Data: dta
REML criterion at convergence: 46724
Scaled residuals:
Min 1Q Median 3Q Max
-3.0969 -0.7322 0.0194 0.7572 2.9147
Random effects:
Groups Name Variance Std.Dev.
School (Intercept) 8.67 2.94
Residual 37.01 6.08
Number of obs: 7185, groups: School, 160
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.636 0.244 51.7
C_SES 2.191 0.109 20.2
Correlation of Fixed Effects:
(Intr)
C_SES 0.000 The variance of school is 8.67 with 2.94 SD and the estimate of C_SES is 2.191 and SE is 0.109
1.22 default residual plot
plot(m3c, xlab = "Fitted values", ylab = "Pearson residuals",
pch = 20, cex = .5, type = c("p", "g"))
1.23 model comparisons
anova(m3a, m3b, m3c)Data: dta
Models:
m3a: Math ~ C_SES + (1 | School)
m3c: Math ~ C_SES + (1 | School) + (C_SES - 1 | School)
m3b: Math ~ C_SES + (C_SES | School)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
m3a 4 46728 46756 -23360 46720
m3c 5 46721 46755 -23355 46711 9.42 1 0.0021
m3b 6 46723 46764 -23355 46711 0.01 1 0.9080Model 3c is significantly different from m3a but no significant difference with m3b.
1.24 CIs for model parameters
confint(m3c, method="boot") 2.5 % 97.5 %
.sig01 2.6100 3.2699
.sig02 0.3848 1.1419
.sigma 5.9425 6.1504
(Intercept) 12.1305 13.0889
C_SES 1.9525 2.43952 In-class exercises 2:
The following R script examines the effect of including predictors at different levels in the school-teacher-pupil example. Compare and explain the variance accounted for in each one of the fitted models.
2.1 Load data file
library(pacman)
# load them
pacman::p_load(mlmRev, tidyverse, lme4, merTools)
# input data
dta <- read.csv("data/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.00002.2 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))2.3 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 )2.4 Model - Null (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.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.7The variance of teacher is 0.137 with 0.370 SD and the variance of school is 0.128 with 0.357 SD
2.5 Model - add a school-level predictor
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
convergence code: 0
boundary (singular) fit: see ?isSingularThe variance of teacher is 0.123 with 0.351 SD and the variance of school is 0 with 0 SD. The effect of mean read scores by school is significant on Read with 1.00 beta and 0.18 SE.
2.6 Model - add a teacher-level predictor
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
convergence code: 0
boundary (singular) fit: see ?isSingularThe variance of teacher is 0.000 with 0.000 SD and the variance of school is 0.04 with 0.2 SD. The effect of mean read scores by teacher is significant on Read with 1.07 beta and 0.11 SE.
2.7 Model - add a teacher-level predictor away from respective school means
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
convergence code: 0
boundary (singular) fit: see ?isSingularThe variance of teacher is 0.000 with 0.000 SD and the variance of school is 0.176 with 0.42 SD. The effect of centered mean read scores by teacher is significant on Read with 1.17 beta and 0.158 SE.
3 In-class exercises 3:
Replicate the results of analysis reported the classroom data example.
3.1 Load data file
library(WWGbook)
data(classroom, "WWGbook")
dta <- classroom
head(dta) sex minority mathkind mathgain ses yearstea mathknow housepov mathprep
1 1 1 448 32 0.46 1 NA 0.082 2.00
2 0 1 460 109 -0.27 1 NA 0.082 2.00
3 1 1 511 56 -0.03 1 NA 0.082 2.00
4 0 1 449 83 -0.38 2 -0.11 0.082 3.25
5 0 1 425 53 -0.03 2 -0.11 0.082 3.25
6 1 1 450 65 0.76 2 -0.11 0.082 3.25
classid schoolid childid
1 160 1 1
2 160 1 2
3 160 1 3
4 217 1 4
5 217 1 5
6 217 1 63.2 Model 1
summary(m1 <- lmer(mathgain ~ mathkind + sex + minority + ses +
(1 | schoolid) + (1 | classid:schoolid), data = dta))Linear mixed model fit by REML ['lmerMod']
Formula: mathgain ~ mathkind + sex + minority + ses + (1 | schoolid) +
(1 | classid:schoolid)
Data: dta
REML criterion at convergence: 11386
Scaled residuals:
Min 1Q Median 3Q Max
-5.826 -0.611 -0.034 0.554 4.268
Random effects:
Groups Name Variance Std.Dev.
classid:schoolid (Intercept) 83.3 9.13
schoolid (Intercept) 75.2 8.67
Residual 734.6 27.10
Number of obs: 1190, groups: classid:schoolid, 312; schoolid, 107
Fixed effects:
Estimate Std. Error t value
(Intercept) 282.7903 10.8532 26.06
mathkind -0.4698 0.0223 -21.10
sex -1.2512 1.6577 -0.75
minority -8.2621 2.3401 -3.53
ses 5.3464 1.2411 4.31
Correlation of Fixed Effects:
(Intr) mthknd sex minrty
mathkind -0.978
sex -0.044 -0.031
minority -0.307 0.163 -0.018
ses 0.140 -0.168 0.019 0.163The variance of school is 75.2 with 8.67 SD and the variance of classid and schoolid is 83.3 with 9.13 SD. The effects of math from kindergarten (with 0.47 beta and 0.02 SE), minority (with 8.26 beta and 2.34 SE), and ses (with 5.34 beta and 1.24 SE) are significant on math.
3.3 Model 2
summary(m2 <- lmer(mathgain ~ (1 | schoolid) + (1 | classid:schoolid), data = dta))Linear mixed model fit by REML ['lmerMod']
Formula: mathgain ~ (1 | schoolid) + (1 | classid:schoolid)
Data: dta
REML criterion at convergence: 11769
Scaled residuals:
Min 1Q Median 3Q Max
-4.644 -0.598 -0.034 0.533 5.634
Random effects:
Groups Name Variance Std.Dev.
classid:schoolid (Intercept) 99.2 9.96
schoolid (Intercept) 77.5 8.80
Residual 1028.2 32.07
Number of obs: 1190, groups: classid:schoolid, 312; schoolid, 107
Fixed effects:
Estimate Std. Error t value
(Intercept) 57.43 1.44 39.8The variance of school is 77.5 with 8.80 SD and the variance of classid and schoolid is 99.2 with 9.96 SD.
3.4 Model 3
summary(m3 <- lmer(mathgain ~ mathkind + sex + minority + ses +
(1 | schoolid) + (1 | classid:schoolid), data = dta))Linear mixed model fit by REML ['lmerMod']
Formula: mathgain ~ mathkind + sex + minority + ses + (1 | schoolid) +
(1 | classid:schoolid)
Data: dta
REML criterion at convergence: 11386
Scaled residuals:
Min 1Q Median 3Q Max
-5.826 -0.611 -0.034 0.554 4.268
Random effects:
Groups Name Variance Std.Dev.
classid:schoolid (Intercept) 83.3 9.13
schoolid (Intercept) 75.2 8.67
Residual 734.6 27.10
Number of obs: 1190, groups: classid:schoolid, 312; schoolid, 107
Fixed effects:
Estimate Std. Error t value
(Intercept) 282.7903 10.8532 26.06
mathkind -0.4698 0.0223 -21.10
sex -1.2512 1.6577 -0.75
minority -8.2621 2.3401 -3.53
ses 5.3464 1.2411 4.31
Correlation of Fixed Effects:
(Intr) mthknd sex minrty
mathkind -0.978
sex -0.044 -0.031
minority -0.307 0.163 -0.018
ses 0.140 -0.168 0.019 0.163The variance of school is 75.2 with 8.67 SD and the variance of classid and schoolid is 83.3 with 9.13 SD. The effects of math from kindergarten (with 0.47 beta and 0.02 SE), minority (with 8.26 beta and 2.34 SE), and ses (with 5.34 beta and 1.24 SE) are significant on math.