if (!require("pacman")) install.packages("pacman", repos="cloud-r.project.org")## 載入需要的套件：pacmanpacman::p_load(mlmRev, gridExtra, tidyverse, coda, reshape,lme4)## 將程式套件安載入 'C:/Users/a0905/AppData/Local/R/win-library/4.2'
## (因為 'lib' 沒有被指定)## Warning: 無法在貯藏處 http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.2 中讀寫索引:
## 無法開啟網址 'http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.2/PACKAGES'## 程式套件 'coda' 開啟成功，MD5 和檢查也透過
##
## 下載的二進位程式套件在
## C:\Users\a0905\AppData\Local\Temp\Rtmp8qvcqD\downloaded_packages 裡##
## coda installed## 將程式套件安載入 'C:/Users/a0905/AppData/Local/R/win-library/4.2'
## (因為 'lib' 沒有被指定)## Warning: 無法在貯藏處 http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.2 中讀寫索引:
## 無法開啟網址 'http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.2/PACKAGES'## 程式套件 'reshape' 開啟成功，MD5 和檢查也透過
##
## 下載的二進位程式套件在
## C:\Users\a0905\AppData\Local\Temp\Rtmp8qvcqD\downloaded_packages 裡##
## reshape installed資料管理
#載入mlmRev套件
library(mlmRev)#載進資料
data(Chem97)#檢視套件提供的資訊
?Chem97## 開啟 httpd 求助伺服器… 好了#看資料結構
str(Chem97)## 'data.frame': 31022 obs. of 8 variables:
## $ lea : Factor w/ 131 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ school : Factor w/ 2410 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ student : Factor w/ 31022 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ score : num 4 10 10 10 8 10 6 8 4 10 ...
## $ gender : Factor w/ 2 levels "M","F": 2 2 2 2 2 2 2 2 2 2 ...
## $ age : num 3 -3 -4 -2 -1 4 1 4 3 0 ...
## $ gcsescore: num 6.62 7.62 7.25 7.5 6.44 ...
## $ gcsecnt : num 0.339 1.339 0.964 1.214 0.158 ...#程式報表5.1
head(Chem97)## lea school student score gender age gcsescore gcsecnt
## 1 1 1 1 4 F 3 6.625 0.3393157
## 2 1 1 2 10 F -3 7.625 1.3393157
## 3 1 1 3 10 F -4 7.250 0.9643157
## 4 1 1 4 10 F -2 7.500 1.2143157
## 5 1 1 5 8 F -1 6.444 0.1583157
## 6 1 1 6 10 F 4 7.750 1.4643157#去掉第六與第八個變項，只將需要的資料留下，並重新命名
#原始文章分析是以 GCSE 預測 Chem，我們預計倒過來做
dta <- Chem97[, -c(6, 8)]
names(dta) <- c("區域", "學校", "學生", "化學", "性別", "總成績")#程式報表5.1
head(dta)## 區域 學校 學生 化學 性別 總成績
## 1 1 1 1 4 F 6.625
## 2 1 1 2 10 F 7.625
## 3 1 1 3 10 F 7.250
## 4 1 1 4 10 F 7.500
## 5 1 1 5 8 F 6.444
## 6 1 1 6 10 F 7.750描述統計
#程式報表5.2
sapply(dta[, c('化學', '總成績')], summary)## 化學 總成績
## Min. 0.000000 0.000000
## 1st Qu. 4.000000 5.750000
## Median 6.000000 6.375000
## Mean 5.812585 6.285684
## 3rd Qu. 8.000000 6.900000
## Max. 10.000000 8.000000sapply(dta[, c('化學', '總成績')], sd)## 化學 總成績
## 3.3191673 0.8734512#計算各校化學與總成績平均分數，記錄成資料檔，留待後用
dta_m <- aggregate(cbind(化學, 總成績, as.numeric(性別) - 1) ~ 學校,
data = dta, mean)
names(dta_m) <- c('學校', '化學平均', '總成績平均', '女性比率')#程式報表5.3
head(dta_m)## 學校 化學平均 總成績平均 女性比率
## 1 1 8.307692 7.224462 1.0000000
## 2 2 8.714286 7.023893 0.0000000
## 3 3 4.666667 5.310000 0.6666667
## 4 4 8.285714 6.232143 0.1428571
## 5 5 5.500000 6.174313 0.4375000
## 6 6 4.714286 6.705000 0.7857143繪圖
#載進 lattice，準備畫圖
library(lattice)#以總成績與化學原始分數，繪製直方圖
p1 <- histogram(~ 總成績, data = dta, type = 'density',
xlab = '總成績', ylab = '機率密度')
p2 <- histogram(~ 化學, data = dta, type = 'density',
xlab = '化學分數', ylab = '機率密度')#以學校平均分數，繪製直方圖
p3 <- histogram(~ 總成績平均, data = dta_m, type = 'density',
xlab = '各校總成績平均', ylab = '機率密度')
p4 <- histogram(~ 化學平均, data = dta_m, type = 'density',
xlab = '各校化學平均分數', ylab = '機率密度')#計算學校標準差，存檔留後用
dta_sd <- aggregate(cbind(化學, 總成績) ~ 學校, data = dta, sd)
names(dta_sd)[2:3] <- c('化學標準差', '總成績標準差')#以學校標準差，繪製直方圖。
p5 <- histogram(~ 總成績標準差, data = dta_sd, type = 'density',
xlab = '各校總成績標準差', ylab = '機率密度')
p6 <- histogram(~ 化學標準差, data = dta_sd, type = 'density',
xlab = '各校化學分數標準差', ylab = '機率密度')#載入gridExtra，把六張放一起
#圖 5.1
library(gridExtra)
grid.arrange(p6, p5, p4, p3, p2, p1, as.table = T)
#看看化學與總成績相關，以學生層次分數計算是 .662
cor(dta[, c('化學', '總成績')])## 化學 總成績
## 化學 1.0000000 0.6622476
## 總成績 0.6622476 1.0000000#以學校層次分數計算生態相關（ecological correlation），為 .698
##程式報表5.4
cor(dta_m[, -1])## 化學平均 總成績平均 女性比率
## 化學平均 1.00000000 0.6980651 0.06625497
## 總成績平均 0.69806507 1.0000000 0.24520916
## 女性比率 0.06625497 0.2452092 1.00000000#畫看兩者間關聯
#圖5.2
plot(dta[, 4], dta[, 6], type = 'n', xlab ='化學分數', ylab = '總成績',
asp = 1)
grid()
#學生
points(dta[, 4], dta[, 6], pch = '.', cex = 2)
abline(lm(dta[,6] ~ dta[, 4]))
#學校
points(dta_m[, 2], dta_m[, 3], cex = 0.5, col = 'grey')
abline(lm(dta_m[,3] ~ dta_m[,2]), col = 'grey')
#對角線
abline(0, 1, lty = 3, col = 'blue')
#看看各校與各區域以化學預測總分時的截距與斜率
#請注意原始文章的預測變項與解釋變項跟此所作是相反的。
#載入 ggplot2 套件，準備畫圖。
library(ggplot2)# 記錄下原始配色
old <- theme_set(theme_bw())#各校以化學預測總分時的截距與斜率
#圖5.3
ggplot(data = dta, aes(x = 化學, y = 總成績, group = 學校))+
stat_smooth(method = 'lm', formula = y ~ x, se = F, color = 'lightgray') +
geom_point(size = 1) +
stat_smooth(aes(group = 1), method= 'lm', se = F, color = 'black') +
labs(x = '化學分數', y = '總成績', title = '學校')## `geom_smooth()` using formula 'y ~ x'
#各區域以化學預測總分時的截距與斜率
#圖5.4
ggplot(data = dta, aes(x = 化學, y = 總成績, group = 區域))+
stat_smooth(method = 'lm', formula = y ~ x, se = F, color = 'lightgray') +
geom_point(size = 1) +
stat_smooth(aes(group = 1), method = 'lm', se = F, color = 'black') +
labs(x = '化學分數', y = '總成績' , title = '區域') ## `geom_smooth()` using formula 'y ~ x'
# 將變項以總分均置中，亦即，減去總平均。
dta$化學置中 <- scale(dta$化學, scale = F)#選取25個學校與區域
set.seed(1225)
ns25 <- sample(levels(dta$學校), 25)
set.seed(1225)
nr25 <- sample(levels(dta$區域), 25)#重看一次各校以化學預測總分時的截距與斜率
#圖5.5
ggplot(data = dta[dta$學校 %in% ns25, ], aes(x = 化學置中, y = 總成績, color = 性別))+
geom_point(size = 1) +
stat_smooth(method = 'lm', formula = y ~ x, se = F) +
facet_wrap( ~ 學校 )
#重看一次各區域以化學預測總分時的截距與斜率
#圖5.6
ggplot(data = dta[dta$區域 %in% nr25, ], aes(x = 化學置中, y = 總成績, color = 性別))+
geom_point(size = 1) +
stat_smooth(method = 'lm', formula = y ~ x, se = F) +
facet_wrap( ~ 區域 )
多層次分析
#載入 lme4 套件，用來分析多層次資料
library(lme4)# 將變項以總分均置中，亦即，減去總平均。
dta$化學置中 <- scale(dta$化學, scale = F)#先以完整模型嘗試
#程式報表5.5
summary(m0 <- lmer(總成績 ~ 化學置中 + 性別 + 化學置中:性別 +
( 1 | 學校 ) + ( 1 | 區域 ), data = dta ) ) ## Linear mixed model fit by REML ['lmerMod']
## Formula: 總成績 ~ 化學置中 + 性別 + 化學置中:性別 + (1 |
## 學校) + (1 | 區域)
## Data: dta
##
## REML criterion at convergence: 55708.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.2306 -0.5671 0.0638 0.6527 4.0788
##
## Random effects:
## Groups Name Variance Std.Dev.
## 學校 (Intercept) 0.08763 0.2960
## 區域 (Intercept) 0.01064 0.1031
## Residual 0.31811 0.5640
## Number of obs: 31022, groups: 學校, 2410; 區域, 131
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.080254 0.013263 458.434
## 化學置中 0.168078 0.001408 119.333
## 性別F 0.334091 0.007518 44.440
## 化學置中:性別F -0.012266 0.002091 -5.866
##
## Correlation of Fixed Effects:## Warning in abbreviate(rn, minlength = 6): abbreviate 不適用於非 ASCII 字元## (Intr) 化學置中 性別F
## 化學置中 0.051
## 性別F -0.265 -0.065
## 化學置中:性別F -0.023 -0.630 0.067# 試著去除區域的隨機效果，並看看去除隨機效果是否顯著。
m1 <- update(m0, . ~ . - ( 1 | 區域 ) )#程式報表5.6
anova(m0, m1)## refitting model(s) with ML (instead of REML)## Data: dta
## Models:
## m1: 總成績 ~ 化學置中 + 性別 + (1 | 學校) + 化學置中:性別
## m0: 總成績 ~ 化學置中 + 性別 + 化學置中:性別 + (1 | 學校) + (1 | 區域)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 6 55743 55793 -27866 55731
## m0 7 55685 55744 -27836 55671 59.744 1 1.08e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#去除交互作用項
#程式報表5.7
drop1(m0, test = 'Chisq')## Single term deletions
##
## Model:
## 總成績 ~ 化學置中 + 性別 + 化學置中:性別 + (1 |
## 學校) + (1 | 區域)
## npar AIC LRT Pr(Chi)
## <none> 55685
## 化學置中:性別 1 55718 34.388 4.515e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m0)## Linear mixed model fit by REML ['lmerMod']
## Formula: 總成績 ~ 化學置中 + 性別 + 化學置中:性別 + (1 |
## 學校) + (1 | 區域)
## Data: dta
##
## REML criterion at convergence: 55708.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.2306 -0.5671 0.0638 0.6527 4.0788
##
## Random effects:
## Groups Name Variance Std.Dev.
## 學校 (Intercept) 0.08763 0.2960
## 區域 (Intercept) 0.01064 0.1031
## Residual 0.31811 0.5640
## Number of obs: 31022, groups: 學校, 2410; 區域, 131
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.080254 0.013263 458.434
## 化學置中 0.168078 0.001408 119.333
## 性別F 0.334091 0.007518 44.440
## 化學置中:性別F -0.012266 0.002091 -5.866
##
## Correlation of Fixed Effects:## Warning in abbreviate(rn, minlength = 6): abbreviate 不適用於非 ASCII 字元## (Intr) 化學置中 性別F
## 化學置中 0.051
## 性別F -0.265 -0.065
## 化學置中:性別F -0.023 -0.630 0.067固定效果
#載入 coefplot2 套件，繪製固定效果
if (!require("coefplot2")) install.packages("coefplot2",
repos="http://www.math.mcmaster.ca/bolker/R")## 載入需要的套件：coefplot2## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : 不存在叫 'coefplot2' 這個名稱的套件## 將程式套件安載入 'C:/Users/a0905/AppData/Local/R/win-library/4.2'
## (因為 'lib' 沒有被指定)## Warning: 無法在貯藏處 http://www.math.mcmaster.ca/bolker/R/bin/windows/contrib/4.2 中讀寫索引:
## 無法開啟網址 'http://www.math.mcmaster.ca/bolker/R/bin/windows/contrib/4.2/PACKAGES'## 安裝原始碼套件 'coefplot2'library(coefplot2)
coefplot2(m0)
#抽取變異成分，計算學校與區域可以解釋部分（ICCs）
print(vc <- lme4::VarCorr(m0), comp = 'Variance' )## Groups Name Variance
## 學校 (Intercept) 0.087635
## 區域 (Intercept) 0.010637
## Residual 0.318108vc <- as.data.frame(vc)
vc[vc$grp=='學校', 'vcov']/ sum(vc$vcov)## [1] 0.2104676vc[vc$grp=='區域', 'vcov']/ sum(vc$vcov)## [1] 0.025547331 - (vc[vc$grp=='Residual', 'vcov']/ sum(vc$vcov))## [1] 0.2360149檢驗模型
#學校與區域層次的隨機效果QQ圖，檢查隨機效果是否呈常態
qq_r21 <- qqmath(~ ranef(m0)$學校, type = c('p', 'g', 'r'), pch = '.',
xlab = '常態分位數', ylab = '學校隨機效果分位數')qq_r22 <- qqmath(~ ranef(m0)$區域, type = c('p', 'g', 'r'), pch = '.',
xlab = '常態分位數', ylab = '區域隨機效果分位數')#殘差的QQ圖，檢驗殘差是否呈常態
qq_r0 <- qqmath(~ resid(m0, scale = T), type = c('p', 'g', 'r'), pch = '.',
xlab = '標準化常態分位數', ylab = '標準化殘差')#預測值對殘差圖
m0_f <- data.frame(yhat=fitted(m0), zres=resid(m0, scaled=TRUE))
r_m0 <- ggplot(data = m0_f, aes(x = yhat, y = zres)) +
geom_point(pch = '.') +
stat_smooth(method = 'gam', formula = y ~ s(x, bs='cs'),
se = FALSE) +
labs( x = '預測值', y = '標準化殘差')#把圖放在一起呈現
#圖5.7
library(gridExtra)
grid.arrange(qq_r21, qq_r22, qq_r0, r_m0, nrow = 2, ncol = 2)
複雜的多層次模型
#顯示複雜的多層次模型，包括隨機斜率效果、學校層次的預測變項以及跨層次交互作用
#先要將不同層次資料併在一起
dta <- merge(dta, dta_m, by="學校") rslt <- lmer(總成績 ~ 化學置中 + 性別 + 化學置中:性別 +
女性比率 + 女性比率:性別+
(1+化學置中+性別|學校) + (1|區域), data=dta)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00683933 (tol = 0.002, component 1)summary(rslt)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## 總成績 ~ 化學置中 + 性別 + 化學置中:性別 + 女性比率 +
## 女性比率:性別 + (1 + 化學置中 + 性別 | 學校) +
## (1 | 區域)
## Data: dta
##
## REML criterion at convergence: 55393
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.3496 -0.5612 0.0682 0.6500 4.1231
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## 學校 (Intercept) 0.0981171 0.31324
## 化學置中 0.0005485 0.02342 -0.56
## 性別F 0.0180773 0.13445 -0.63 0.07
## 區域 (Intercept) 0.0106268 0.10309
## Residual 0.3103917 0.55713
## Number of obs: 31022, groups: 學校, 2410; 區域, 131
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.140239 0.019252 318.944
## 化學置中 0.169861 0.001557 109.110
## 性別F 0.146319 0.022258 6.574
## 女性比率 -0.222862 0.040019 -5.569
## 化學置中:性別F -0.011096 0.002156 -5.147
## 性別F:女性比率 0.430826 0.048234 8.932
##
## Correlation of Fixed Effects:## Warning in abbreviate(rn, minlength = 6): abbreviate 不適用於非 ASCII 字元## (Intr) 化學置中 性別F 女性比率 化中:性別F
## 化學置中 -0.095
## 性別F -0.482 -0.004
## 女性比率 -0.705 0.069 0.645
## 化學置中:性別F 0.066 -0.608 0.004 -0.101
## 性別F:女性比率 0.494 -0.026 -0.915 -0.797 0.030
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00683933 (tol = 0.002, component 1)結束
顯示演練單元信息
sessionInfo()## R version 4.2.1 (2022-06-23 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22000)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Chinese (Traditional)_Taiwan.utf8
## [2] LC_CTYPE=Chinese (Traditional)_Taiwan.utf8
## [3] LC_MONETARY=Chinese (Traditional)_Taiwan.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=Chinese (Traditional)_Taiwan.utf8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] coefplot2_0.1.3.2 lattice_0.20-45 reshape_0.8.9 coda_0.19-4
## [5] forcats_0.5.2 stringr_1.4.1 dplyr_1.0.10 purrr_0.3.4
## [9] readr_2.1.2 tidyr_1.2.1 tibble_3.1.8 ggplot2_3.3.6
## [13] tidyverse_1.3.2 gridExtra_2.3 mlmRev_1.0-8 lme4_1.1-30
## [17] Matrix_1.5-1 pacman_0.5.1
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.4 sass_0.4.2 jsonlite_1.8.0
## [4] splines_4.2.1 modelr_0.1.9 bslib_0.4.0
## [7] assertthat_0.2.1 highr_0.9 googlesheets4_1.0.1
## [10] cellranger_1.1.0 pillar_1.8.1 backports_1.4.1
## [13] glue_1.6.2 digest_0.6.29 rvest_1.0.3
## [16] minqa_1.2.4 colorspace_2.0-3 plyr_1.8.7
## [19] htmltools_0.5.3 pkgconfig_2.0.3 broom_1.0.1
## [22] haven_2.5.1 scales_1.2.1 tzdb_0.3.0
## [25] googledrive_2.0.0 mgcv_1.8-40 farver_2.1.1
## [28] generics_0.1.3 ellipsis_0.3.2 cachem_1.0.6
## [31] withr_2.5.0 cli_3.4.0 crayon_1.5.1
## [34] magrittr_2.0.3 readxl_1.4.1 evaluate_0.16
## [37] fs_1.5.2 fansi_1.0.3 nlme_3.1-157
## [40] MASS_7.3-57 xml2_1.3.3 tools_4.2.1
## [43] hms_1.1.2 gargle_1.2.1 lifecycle_1.0.2
## [46] munsell_0.5.0 reprex_2.0.2 compiler_4.2.1
## [49] jquerylib_0.1.4 rlang_1.0.5 grid_4.2.1
## [52] nloptr_2.0.3 rstudioapi_0.14 labeling_0.4.2
## [55] rmarkdown_2.16 boot_1.3-28 gtable_0.3.1
## [58] DBI_1.1.3 R6_2.5.1 lubridate_1.8.0
## [61] knitr_1.40 fastmap_1.1.0 utf8_1.2.2
## [64] stringi_1.7.8 Rcpp_1.0.9 vctrs_0.4.1
## [67] dbplyr_2.2.1 tidyselect_1.1.2 xfun_0.33