Exercise 2: Read and replicate the results in Chapter 5 of Cheng, C.-P., & Sheu, C.-F. (2015). Behavioral data analysis with R (Chinese). Yeh Yeh Book Gallery.

1 資料管理

# 載入 mlmRev 套件
library(mlmRev)

# 載進資料 
data(Chem97)

# 檢視套件提供的資訊
?Chem97

# 看資料結構
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 ...

1.1 程式報表5.1

# 顯示前六筆資料
head(Chem97)

  lea school student score gender age gcsescore  gcsecnt
1   1      1       1     4      F   3     6.625 0.339316
2   1      1       2    10      F  -3     7.625 1.339316
3   1      1       3    10      F  -4     7.250 0.964316
4   1      1       4    10      F  -2     7.500 1.214316
5   1      1       5     8      F  -1     6.444 0.158316
6   1      1       6    10      F   4     7.750 1.464316

# 去掉第六與第八個變項，並將變數重新命名
dta <- Chem97[, -c(6, 8)]
names(dta) <- c("區域", "學校", "學生", "化學", "性別", "總成績")
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

2 描述性統計

2.1 程式報表5.2

# 顯示化學與總成績的描述性統計
sapply(dta[, c('化學', '總成績')], summary)

            化學  總成績
Min.     0.00000 0.00000
1st Qu.  4.00000 5.75000
Median   6.00000 6.37500
Mean     5.81258 6.28568
3rd Qu.  8.00000 6.90000
Max.    10.00000 8.00000

sapply(dta[, c('化學', '總成績')], sd)

    化學   總成績 
3.319167 0.873451

2.2 程式報表5.3

# 計算各校化學和總成績平均分數及女性比率
dta_m <- aggregate(cbind(化學, 總成績, as.numeric(性別) - 1) ~ 學校, 
                   data = dta, mean)
names(dta_m) <- c('學校', '化學平均', '總成績平均', '女性比率')
head(dta_m)

  學校 化學平均 總成績平均 女性比率
1    1  8.30769    7.22446 1.000000
2    2  8.71429    7.02389 0.000000
3    3  4.66667    5.31000 0.666667
4    4  8.28571    6.23214 0.142857
5    5  5.50000    6.17431 0.437500
6    6  4.71429    6.70500 0.785714

3 繪製直方圖

# 載入 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 = '機率密度')

3.1 圖 5.1

# 載入 gridExtra ，把六張直方圖放一起
library(gridExtra)
grid.arrange(p6, p5, p4, p3, p2, p1, as.table = T)

4 相關分析

4.1 學生層次

# 以學生層次計算化學與總成績相關係數為 .662
cor(dta[, c('化學', '總成績')])

           化學   總成績
化學   1.000000 0.662248
總成績 0.662248 1.000000

4.2 學校層次，程式報表5.4

# 以學校層次計算化學與總成績生態相關（ecological correlation）係數為 .698
cor(dta_m[, -1])

           化學平均 總成績平均 女性比率
化學平均   1.000000   0.698065 0.066255
總成績平均 0.698065   1.000000 0.245209
女性比率   0.066255   0.245209 1.000000

5 繪製散佈圖，圖 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')

6 繪製迴歸線(所有迴歸線呈現在同一張圖上)

6.1 以學校為單位，圖 5.3

# 載入 ggplot2 套件，準備畫圖。
library(ggplot2)

# 記錄下原始配色
old <- theme_set(theme_bw())

# 各校以化學預測總分時的截距與斜率
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 = '學校')

6.2 以區域為單位，圖 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 = '區域')

7 繪製迴歸線(依照學校或區域分別繪製)

7.1 以學校為單位，圖 5.5

# 將變項以總分均置中，亦即，減去總平均。
dta$化學置中 <- scale(dta$化學, scale = F)

# 隨機選取25個學校與區域 
set.seed(1225)
ns25 <- sample(levels(dta$學校), 25)
set.seed(1225)
nr25 <- sample(levels(dta$區域), 25)

# 重看一次各校以化學預測總分時的截距與斜率
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( ~ 學校 )

7.2 以區域為單位，圖 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( ~ 區域 )

8 多層次分析

# 載入 lme4 套件，用來分析多層次資料
library(lme4)

# 將變項以總分均置中，亦即，減去總平均。
dta$化學置中 <- scale(dta$化學, scale = F)

8.1 完整模型，程式報表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.231 -0.567  0.064  0.653  4.079 

Random effects:
 Groups   Name        Variance Std.Dev.
 學校     (Intercept) 0.0876   0.296   
 區域     (Intercept) 0.0106   0.103   
 Residual             0.3181   0.564   
Number of obs: 31022, groups:  學校, 2410; 區域, 131

Fixed effects:
               Estimate Std. Error t value
(Intercept)     6.08025    0.01326  458.43
化學置中        0.16808    0.00141  119.33
性別F           0.33409    0.00752   44.44
化學置中:性別F -0.01227    0.00209   -5.87

Correlation of Fixed Effects:
               (Intr) 化學置中 性別F 
化學置中        0.051                
性別F          -0.265 -0.065         
化學置中:性別F -0.023 -0.630    0.067

8.2 去除區域的隨機效果，程式報表5.6

# 試著去除區域的隨機效果，並看看去除隨機效果是否顯著。
m1 <- update(m0, . ~ . - ( 1 | 區域 ) )
anova(m0, m1)

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.74  1   1.08e-14

8.3 去除交互作用項，程式報表5.7

drop1(m0, test = 'Chisq')

Single term deletions

Model:
總成績 ~ 化學置中 + 性別 + 化學置中:性別 + (1 | 
    學校) + (1 | 區域)
              npar   AIC   LRT  Pr(Chi)
<none>             55685               
化學置中:性別    1 55718 34.39 4.52e-09

8.4 固定效果

# 載入 coefplot2 套件，繪製固定效果
library(coefplot2)
coefplot2(m0)

# 抽取變異成分，計算學校與區域可以解釋部分（ICCs）
print(vc <- lme4::VarCorr(m0), comp = 'Variance' )

 Groups   Name        Variance
 學校     (Intercept) 0.08763 
 區域     (Intercept) 0.01064 
 Residual             0.31811

vc <- as.data.frame(vc)
vc[vc$grp=='學校', 'vcov']/ sum(vc$vcov)

[1] 0.210468

vc[vc$grp=='區域', 'vcov']/ sum(vc$vcov)

[1] 0.0255473

1 - (vc[vc$grp=='Residual', 'vcov']/ sum(vc$vcov))

[1] 0.236015

8.5 常態假設檢驗與繪圖，圖5.7

# 學校與區域層次的隨機效果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 = '標準化殘差')

# 把圖放在一起呈現
library(gridExtra)
grid.arrange(qq_r21, qq_r22, qq_r0, r_m0, nrow = 2, ncol = 2)

8.6 複雜的多層次模型

# 顯示複雜的多層次模型，包括隨機斜率效果、學校層次的預測變項以及跨層次交互作用
# 先要將不同層次資料併在一起
dta <- merge(dta, dta_m, by="學校") 
rslt <- lmer(總成績 ~ 化學置中 + 性別 + 化學置中:性別 +
             女性比率 + 女性比率:性別+
             (1+化學置中+性別|學校) + (1|區域), data=dta)
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.350 -0.561  0.068  0.650  4.123 

Random effects:
 Groups   Name        Variance Std.Dev. Corr       
 學校     (Intercept) 0.098117 0.3132              
          化學置中    0.000548 0.0234   -0.56      
          性別F       0.018077 0.1345   -0.63  0.07
 區域     (Intercept) 0.010627 0.1031              
 Residual             0.310392 0.5571              
Number of obs: 31022, groups:  學校, 2410; 區域, 131

Fixed effects:
               Estimate Std. Error t value
(Intercept)     6.14024    0.01925  318.94
化學置中        0.16986    0.00156  109.11
性別F           0.14632    0.02226    6.57
女性比率       -0.22286    0.04002   -5.57
化學置中:性別F -0.01110    0.00216   -5.15
性別F:女性比率  0.43083    0.04823    8.93

Correlation of Fixed Effects:
               (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)

Multilevel Analysis (9/17) Exercise 2

Chin Yuan Chang

2022-09-22

1 資料管理

1.1 程式報表5.1

2 描述性統計

2.1 程式報表5.2

2.2 程式報表5.3

3 繪製直方圖

3.1 圖 5.1

4 相關分析

4.1 學生層次

4.2 學校層次，程式報表5.4

5 繪製散佈圖，圖 5.2

6 繪製迴歸線(所有迴歸線呈現在同一張圖上)

6.1 以學校為單位，圖 5.3

6.2 以區域為單位，圖 5.4

7 繪製迴歸線(依照學校或區域分別繪製)

7.1 以學校為單位，圖 5.5

7.2 以區域為單位，圖 5.6

8 多層次分析

8.1 完整模型，程式報表5.5

8.2 去除區域的隨機效果，程式報表5.6

8.3 去除交互作用項，程式報表5.7

8.4 固定效果

8.5 常態假設檢驗與繪圖，圖5.7

8.6 複雜的多層次模型