Multilevel models: Exercise & Homework
CHIU, Ming-Tzu
Sun Oct 11 23:29:35 2020
1 Exercise 1
Replicate the results of analysis reported the high school and beyond - null 2 3 4 examples.
1.1 high school and beyond - null
pacman::p_load(mlmRev, tidyverse, lme4)讀取資料
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 ...dta1 <- Hsb82 資料管理好習慣: 重新命名變量
names(dta1) <- c("School", "Minority", "Gender", "SES", "Math",
"Ave_SES", "Sector", "C_SES")head(dta1) 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.456383資料視覺化前置作業
ot <- theme_set(theme_bw())freedman-diaconis rule for binwidth
bw <- 2*IQR(dta1$Math)/(dim(dta1)[1]^(1/3))長條圖
ggplot(dta1, aes(Math)) +
stat_bin(aes(fill = ..count..), binwidth=bw) +
labs(x="Math achievement scores")
隨機取樣設定
set.seed(20160923)n <- sample(160, 60)隨機取樣60間學校的散布圖
ggplot(filter(dta1, School %in% levels(School)[n]), # 挑選符合隨機結果的60間學校(學校為order factor，可用level選取)
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")
以學校分組的數學成績平均
dta1_mmath <- dta1 %>%
group_by(School) %>%
mutate( ave_math = mean(Math)) %>%
select( School, ave_math ) %>%
unique %>%
as.data.frame各組學校數學成績平均的描述統計
summary(dta1_mmath$ave_math) Min. 1st Qu. Median Mean 3rd Qu. Max.
4.24 10.47 12.90 12.62 14.65 19.72 各組學校數學成績平均的標準差
sd(dta1_mmath$ave_math)[1] 3.1177數學成績的描述統計，包含總平均
summary(dta1$Math) Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.83 7.28 13.13 12.75 18.32 24.99 整體數學成績的標準差
sd(dta1$Math)[1] 6.8782以學校分組的資料視覺化
bw <- 2*IQR(dta1_mmath$ave_math)/(160^(1/3))ggplot(dta1_mmath, aes(ave_math)) +
geom_histogram(aes(y= ..density.., alpha= .2), binwidth=bw) + #各學校數學成績平均的密度圖
geom_density(fill="NA") +
geom_vline(xintercept=mean(dta1_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")
baseline model：以學校隨機分組的迴歸模型(含截距)
summary(m0 <- lmer(Math ~ (1 | School), data = dta1))Linear mixed model fit by REML ['lmerMod']
Formula: Math ~ (1 | School)
Data: dta1
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.7baseline model中參數(指示學校的截距)的信賴區間
confint(m0, method="boot") 2.5 % 97.5 %
.sig01 2.5963 3.3335
.sigma 6.1515 6.3605
(Intercept) 12.1289 13.1146殘差圖
plot(m0, xlab="Fitted values", ylab="Pearson residuals",
pch=20, cex=.5, type = c("p","g"))
1.2 high school and beyond - 2
資料副本建立與管理
dta12 <- Hsb82names(dta12) <- c("School", "Minority", "Gender", "SES", "Math",
"Ave_SES", "Sector", "C_SES")str(dta12)'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 ...random effects ANOVA - intercept only model
summary(m1 <- lmer(Math ~ (1 | School), data = dta12))Linear mixed model fit by REML ['lmerMod']
Formula: Math ~ (1 | School)
Data: dta12
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.7# random effect: schoolrandom intercept; fixed-effect is Mean SES (平均社經地位)
summary(m2 <- lmer(Math ~ Ave_SES + (1 | School), data = dta12))Linear mixed model fit by REML ['lmerMod']
Formula: Math ~ Ave_SES + (1 | School)
Data: dta12
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 # random effect: school
# fixed effect: Ave_SES殘差圖
plot(m2, xlab = "Fitted values", ylab = "Pearson residuals",
main = "Model 2", pch = 20, cex = .5, type = c("p", "g"))
整體的迴歸模型
summary(m0a <- lm(Math ~ Ave_SES, data = dta12))
Call:
lm(formula = Math ~ Ave_SES, data = dta12)
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-16資料視覺化
theme_set(theme_bw())學生個體數學成績與各校學生家庭平均社經地位的平均
ggplot(dta12, aes(Ave_SES, Math)) +
geom_smooth(method = "lm") +
geom_point(alpha = .2) +
labs(x = "School mean SES", y = "Math achievement score")
各校數學成績平均
dta12_mmath <- dta12 %>%
group_by(School) %>%
mutate( ave_math = mean(Math)) %>%
select( School, Ave_SES, ave_math ) %>%
unique %>%
as.data.frame各校平均數學成績與平均社經地位散布圖與迴歸線
ggplot(dta12_mmath, aes(Ave_SES, ave_math)) +
geom_smooth(method = "lm") +
geom_point(alpha = .2)+
labs(x = "Mean School SES", y = "Mean School Math Scores")
1.3 high school and beyond - 3
資料副本建立與管理
dta13 <- Hsb82names(dta13) <- c("School", "Minority", "Gender", "SES", "Math",
"Ave_SES", "Sector", "C_SES")隨機分配設置
set.seed(20160923)從所有學校裡隨機選出31所
n <- sample(length(levels(dta13$School)), 31)資料視覺化
theme_set(theme_bw())各校數學成績與各校標準化後的社經地位
ggplot(data = filter(dta13, School %in% levels(School)[n]), # 選取隨機選取的31所學校
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))
根據每一學校建立一不含截距之迴歸模型(math ~ c_ses)
summary(m1b <- lm(Math ~ School / C_SES - 1, data = dta13))
Call:
lm(formula = Math ~ School/C_SES - 1, data = dta13)
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-16收集上述模型裡各個迴歸檢定值(coefficient)
m1b_betas <- m1b %>%
coef %>%
matrix(ncol = 2) %>%
as.data.frame %>%
rename( b0 = V1, b1 = V2)
# b0 : 各校 > 斜率
# b1 : school:C_SES > 截距head(m1b_betas) b0 b1
1 4.5528 0.25037
2 4.2398 1.93884
3 5.8114 1.13184
4 4.3249 -1.01410
5 5.9768 0.94769
6 7.2714 3.01800各欄取平均
colMeans(m1b_betas) b0 b1
12.6208 2.2016 標準差
apply(m1b_betas, 2, sd) b0 b1
3.1177 1.6310 相關係數
cor(m1b_betas) b0 b1
b0 1.000000 0.013467
b1 0.013467 1.000000迴歸共變數檢定值散布圖(斜率與截距)
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")
random effect 為學校
summary(m3a <- lmer(Math ~ C_SES + (1 | School), data = dta13))Linear mixed model fit by REML ['lmerMod']
Formula: Math ~ C_SES + (1 | School)
Data: dta13
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 random effect 為學校，模型截距與各學校的斜率交互作用
summary(m3b <- lmer(Math ~ C_SES + (C_SES | School), data = dta13))Linear mixed model fit by REML ['lmerMod']
Formula: Math ~ C_SES + (C_SES | School)
Data: dta13
REML criterion at convergence: 46714
Scaled residuals:
Min 1Q Median 3Q Max
-3.0968 -0.7319 0.0186 0.7539 2.8992
Random effects:
Groups Name Variance Std.Dev. Corr
School (Intercept) 8.681 2.946
C_SES 0.694 0.833 0.02
Residual 36.700 6.058
Number of obs: 7185, groups: School, 160
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.636 0.245 51.7
C_SES 2.193 0.128 17.1
Correlation of Fixed Effects:
(Intr)
C_SES 0.009 random effect 為學校，模型截距與各學校斜率無相關
summary(m3c <- lmer(Math ~ C_SES + (1 | School) +
(C_SES - 1 | School), data = dta13))Linear mixed model fit by REML ['lmerMod']
Formula: Math ~ C_SES + (1 | School) + (C_SES - 1 | School)
Data: dta13
REML criterion at convergence: 46714
Scaled residuals:
Min 1Q Median 3Q Max
-3.0961 -0.7310 0.0169 0.7539 2.9017
Random effects:
Groups Name Variance Std.Dev.
School (Intercept) 8.681 2.946
School.1 C_SES 0.694 0.833
Residual 36.700 6.058
Number of obs: 7185, groups: School, 160
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.636 0.245 51.7
C_SES 2.193 0.128 17.1
Correlation of Fixed Effects:
(Intr)
C_SES 0.000 殘差圖
plot(m3c, xlab = "Fitted values", ylab = "Pearson residuals",
pch = 20, cex = .5, type = c("p", "g"))
三模型比較
anova(m3a, m3b, m3c)Data: dta13
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.9080信賴區間
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.43951.4 high school and beyond - 4
資料複本與整理
dta14 <- Hsb82names(dta14) <- c("School", "Minority", "Gender", "SES", "Math",
"Ave_SES", "Sector", "C_SES")資料視覺化
theme_set(theme_bw())以 sector 看數學與社經地位(centered at grand mean)
ggplot(dta14, aes(C_SES, Math, color = Sector)) +
geom_smooth(method = "lm") +
geom_jitter(alpha = 0.2) +
labs(x = "Centered SES", y = "Math achievement score") +
theme(legend.position = c(.1, .8))
考慮所有參數的 mixed model
summary(m4a <- lmer(Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES +
Sector:C_SES + (C_SES | School), data = dta14))Linear mixed model fit by REML ['lmerMod']
Formula: Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES + Sector:C_SES +
(C_SES | School)
Data: dta14
REML criterion at convergence: 46504
Scaled residuals:
Min 1Q Median 3Q Max
-3.159 -0.723 0.017 0.754 2.958
Random effects:
Groups Name Variance Std.Dev. Corr
School (Intercept) 2.380 1.543
C_SES 0.101 0.318 0.39
Residual 36.721 6.060
Number of obs: 7185, groups: School, 160
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.128 0.199 60.86
Ave_SES 5.333 0.369 14.45
SectorCatholic 1.227 0.306 4.00
C_SES 2.945 0.156 18.93
Ave_SES:C_SES 1.039 0.299 3.48
SectorCatholic:C_SES -1.643 0.240 -6.85
Correlation of Fixed Effects:
(Intr) Av_SES SctrCt C_SES A_SES:
Ave_SES 0.256
SectorCthlc -0.699 -0.356
C_SES 0.075 0.019 -0.053
A_SES:C_SES 0.019 0.074 -0.026 0.293
SctrC:C_SES -0.052 -0.027 0.077 -0.696 -0.351random effect 只考慮學校
summary(m4b <- lmer(Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES +
Sector:C_SES + (1 | School), data = dta14))Linear mixed model fit by REML ['lmerMod']
Formula: Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES + Sector:C_SES +
(1 | School)
Data: dta14
REML criterion at convergence: 46505
Scaled residuals:
Min 1Q Median 3Q Max
-3.170 -0.725 0.015 0.754 2.966
Random effects:
Groups Name Variance Std.Dev.
School (Intercept) 2.38 1.54
Residual 36.77 6.06
Number of obs: 7185, groups: School, 160
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.128 0.199 60.88
Ave_SES 5.337 0.369 14.46
SectorCatholic 1.225 0.306 4.00
C_SES 2.942 0.151 19.46
Ave_SES:C_SES 1.044 0.291 3.59
SectorCatholic:C_SES -1.642 0.233 -7.05
Correlation of Fixed Effects:
(Intr) Av_SES SctrCt C_SES A_SES:
Ave_SES 0.256
SectorCthlc -0.699 -0.356
C_SES 0.000 0.000 0.000
A_SES:C_SES 0.000 0.000 0.000 0.295
SctrC:C_SES 0.000 0.000 0.000 -0.696 -0.351二個模型的比較
anova(m4a, m4b)Data: dta14
Models:
m4b: Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES + Sector:C_SES +
m4b: (1 | School)
m4a: Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES + Sector:C_SES +
m4a: (C_SES | School)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
m4b 8 46513 46568 -23249 46497
m4a 10 46516 46585 -23248 46496 1 2 0.61h殘差圖
plot(m4b, xlab = "Fitted values", ylab = "Pearson residuals",
pch = 20, cex = .5, type = c("p", "g"))
模型參數共變數與標準差
library(merTools)
fastdisp(m4b)lmer(formula = Math ~ Ave_SES + Sector + C_SES + Ave_SES:C_SES +
Sector:C_SES + (1 | School), data = dta14)
coef.est coef.se
(Intercept) 12.13 0.20
Ave_SES 5.34 0.37
SectorCatholic 1.22 0.31
C_SES 2.94 0.15
Ave_SES:C_SES 1.04 0.29
SectorCatholic:C_SES -1.64 0.23
Error terms:
Groups Name Std.Dev.
School (Intercept) 1.54
Residual 6.06
---
number of obs: 7185, groups: School, 160
AIC = 46520.8利用模擬估計 fixed-effect 的參數
m4b_fe <- FEsim(m4b, 1000)視覺化 fixed-effect
plotFEsim(m4b_fe) +
labs(title = "Coefficient Plot", x = "Median Effect Estimate")
利用模擬估計 randeom-effect 參數
m4b_re <- REsim(m4b)視覺化 randeom effect 標準化後的結果
plotREsim(m4b_re)
The end
2 Exercise 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.
library(pacman)
pacman::p_load(mlmRev, tidyverse, lme4, merTools)資料輸入
dta2 <- read.csv("~/data/pts.csv")head(dta2) 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資料屬性改變，增加以 school 分組計算結束時閱讀能力平均值的平均
dta2 <- dta2 %>%
mutate(School = factor(School), Teacher = factor(Teacher),
Pupil = factor(Pupil)) %>%
group_by( School ) %>%
mutate(msRead_0 = mean(Read_0))以 teacher 為分組增加結束時閱讀能力平均值的平均，並計算與上述以學校分組之平均
dta2 <- dta2 %>%
group_by( Teacher ) %>%
mutate(mtRead_0 = mean(Read_0), ctRead_0 = mean(Read_0) - msRead_0 )無預測子的模型
summary(m0 <- lmer(Read_1 ~ (1 | School) + (1 | Teacher), data = dta2))Linear mixed model fit by REML ['lmerMod']
Formula: Read_1 ~ (1 | School) + (1 | Teacher)
Data: dta2
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加入學校層級的預測子 (加上學校分組曲的平均)
summary(m0_s <- update(m0, . ~ . + msRead_0))Linear mixed model fit by REML ['lmerMod']
Formula: Read_1 ~ (1 | School) + (1 | Teacher) + msRead_0
Data: dta2
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 ?isSingularRandom effect 裡學校因子 variance 變成 0 了，teacher 的 variance 和 sd 都變小。
加上老師分組曲的平均作為預測子
summary(m0_t <- update(m0, . ~ . + mtRead_0))Linear mixed model fit by REML ['lmerMod']
Formula: Read_1 ~ (1 | School) + (1 | Teacher) + mtRead_0
Data: dta2
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 ?isSingularRandom effect 裡teacher 的 variance 和 sd 都變非常小，school 的 variance 和 sd 都變小。
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: dta2
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 end
3 Exercise 3
Replicate the results of analysis reported the classroom data example.
4 Homework 1
Use the models specified in the IQ and language example to replicate the results of analysis.
5 Homework 2
Replicate the results of analysis reported the math scores over years data example.
pacman::p_load(foreign, tidyverse, lme4, merTools)資料讀取
fLoc <- "~/data/eg_hlm.sav"dtah2 <- read.spss(fLoc) %>% as.data.framestr(dtah2)'data.frame': 7230 obs. of 12 variables:
$ size : num 380 380 380 380 380 380 380 380 380 380 ...
$ lowinc : num 40.3 40.3 40.3 40.3 40.3 40.3 40.3 40.3 40.3 40.3 ...
$ mobility: num 12.5 12.5 12.5 12.5 12.5 12.5 12.5 12.5 12.5 12.5 ...
$ female : num 0 0 0 0 0 0 0 0 0 0 ...
$ black : num 0 0 0 0 0 0 0 0 0 0 ...
$ hispanic: num 1 1 1 0 0 0 0 0 1 1 ...
$ year : num 0.5 1.5 2.5 -1.5 -0.5 0.5 1.5 2.5 -1.5 -0.5 ...
$ grade : num 2 3 4 0 1 2 3 4 0 1 ...
$ math : num 1.146 1.134 2.3 -1.303 0.439 ...
$ retained: num 0 0 0 0 0 0 0 0 0 0 ...
$ school : Factor w/ 60 levels " 2020",..: 1 1 1 1 1 1 1 1 1 1 ...
$ cid : Factor w/ 1721 levels " 101480302",..: 244 244 244 248 248 248 248 248 253 253 ...資料視覺化
theme_set(theme_bw())ggplot(dtah2, aes(year, math, group = cid)) +
geom_point(alpha = .1) +
geom_line(alpha = .1) +
stat_smooth(aes(group = school), method = "lm", se = F, lwd = .6) +
stat_smooth(aes(group = 1), method = "lm", se = F, col = "magenta") +
labs(x = "Year of age (centered)", y = "Math score")
藍色線為各學校分組的迴歸線，桃紅色線為整體的迴歸線。
nested
fixed-effect for slope
summary(m1 <- lmer(math ~ year + (1 | school/cid), data = dtah2))Linear mixed model fit by REML ['lmerMod']
Formula: math ~ year + (1 | school/cid)
Data: dtah2
REML criterion at convergence: 16759
Scaled residuals:
Min 1Q Median 3Q Max
-3.809 -0.583 -0.027 0.557 6.087
Random effects:
Groups Name Variance Std.Dev.
cid:school (Intercept) 0.670 0.818
school (Intercept) 0.187 0.432
Residual 0.347 0.589
Number of obs: 7230, groups: cid:school, 1721; school, 60
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.7805 0.0611 -12.8
year 0.7461 0.0054 138.3
Correlation of Fixed Effects:
(Intr)
year -0.031以下相同
summary(lmer(math ~ year + (1 | school) +( 1 | school:cid) , data = dtah2))Linear mixed model fit by REML ['lmerMod']
Formula: math ~ year + (1 | school) + (1 | school:cid)
Data: dtah2
REML criterion at convergence: 16759
Scaled residuals:
Min 1Q Median 3Q Max
-3.809 -0.583 -0.027 0.557 6.087
Random effects:
Groups Name Variance Std.Dev.
school:cid (Intercept) 0.670 0.818
school (Intercept) 0.187 0.432
Residual 0.347 0.589
Number of obs: 7230, groups: school:cid, 1721; school, 60
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.7805 0.0611 -12.8
year 0.7461 0.0054 138.3
Correlation of Fixed Effects:
(Intr)
year -0.031plot(m1, ylim = c(-4, 4),
xlab = "Fitted values", ylab = "Pearson residuals", cex =.5)
random slope
summary(m2 <- lmer(math ~ ( year | school/cid), data = dtah2))Linear mixed model fit by REML ['lmerMod']
Formula: math ~ (year | school/cid)
Data: dtah2
REML criterion at convergence: 16554
Scaled residuals:
Min 1Q Median 3Q Max
-3.223 -0.561 -0.031 0.531 5.158
Random effects:
Groups Name Variance Std.Dev. Corr
cid:school (Intercept) 0.6413 0.801
year 0.0113 0.106 0.55
school (Intercept) 0.9244 0.961
year 0.6079 0.780 0.92
Residual 0.3011 0.549
Number of obs: 7230, groups: cid:school, 1721; school, 60
Fixed effects:
Estimate Std. Error t value
(Intercept) -1.6442 0.0546 -30.1
convergence code: 0
Model failed to converge with max|grad| = 0.00321233 (tol = 0.002, component 1)比較二者
anova(m1, m2)Data: dtah2
Models:
m1: math ~ year + (1 | school/cid)
m2: math ~ (year | school/cid)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
m1 5 16757 16791 -8374 16747
m2 8 16566 16621 -8275 16550 197 3 <2e-16使用 merTools 利用模擬估計 random-effect的參數
m2_re <- REsim(m2)plotREsim(m2_re)
plot(m2, ylim = c(-4, 4),
xlab = "Fitted values", ylab = "Pearson residuals", cex =.5)
The end