分析流程
- 建議分析流程三步驟：探索資料、正式分析、診斷
- 大致上順序是這樣，但有時發現問題，程序會來來回回
- 不過我猜大家急著看正式分析，沒耐心看探索資料，所以先讓大家看正式分析
- 還是建議大家要記得探索與診斷，可以跟正式分析互補
讀進資料
並將 ID 視為類別變項（factor）
dta <- read.csv('Chemistry.csv',header=T)
head(dta)
## area schID stuID chem gender sciscore
## 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
dta$schID <- as.factor(dta$schID)
dta$area <- as.factor(dta$area)
載入 lme4 套件，用來分析多層次資料
library(pacman)
## Warning: 套件 'pacman' 是用 R 版本 4.2.3 來建造的
p_load(lme4)
形成學校層次的資料
dta_m <- aggregate(gender=="F" ~ schID,
data = dta, mean)
names(dta_m) <- c('schID', 'girlprop')
head(dta_m)
## schID girlprop
## 1 1 1.0000000
## 2 2 0.0000000
## 3 3 0.6666667
## 4 4 0.1428571
## 5 5 0.4375000
## 6 6 0.7857143
包含不同層次獨變項的多層次模型
- 先要將不同層次資料併在一起
- 可以看看，隨機斜率效果、學校層次的預測變項以及跨層次交互作用如何設定
dta <- merge(dta,dta_m, by="schID")
rst <- lmer(sciscore ~ gender + girlprop + girlprop:gender + (1|schID) + (1|area), data=dta,REML = FALSE)
summary(rst)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sciscore ~ gender + girlprop + girlprop:gender + (1 | schID) +
## (1 | area)
## Data: dta
##
## AIC BIC logLik deviance df.resid
## 72293.3 72351.7 -36139.6 72279.3 31015
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.2937 -0.6150 0.0564 0.6796 3.4018
##
## Random effects:
## Groups Name Variance Std.Dev.
## schID (Intercept) 0.18588 0.4311
## area (Intercept) 0.02498 0.1580
## Residual 0.53656 0.7325
## Number of obs: 31022, groups: schID, 2410; area, 131
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.15730 0.03177 193.790
## genderM -0.10375 0.02931 -3.540
## girlprop 0.29848 0.04414 6.762
## genderM:girlprop -0.59695 0.06253 -9.547
##
## Correlation of Fixed Effects:
## (Intr) gendrM grlprp
## genderM -0.616
## girlprop -0.779 0.631
## gndrM:grlpr 0.537 -0.937 -0.587
話說從頭，我們先從學校層級性別隨機斜率與隨機截距開始
rst1 <- lmer(sciscore ~ gender + (1+gender|schID) , data=dta, REML = FALSE)
summary(rst1)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sciscore ~ gender + (1 + gender | schID)
## Data: dta
##
## AIC BIC logLik deviance df.resid
## 72385.6 72435.6 -36186.8 72373.6 31016
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.2752 -0.6142 0.0566 0.6764 3.4032
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## schID (Intercept) 0.19036 0.4363
## genderM 0.03632 0.1906 0.19
## Residual 0.53026 0.7282
## Number of obs: 31022, groups: schID, 2410
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.36859 0.01195 532.84
## genderM -0.37955 0.01136 -33.43
##
## Correlation of Fixed Effects:
## (Intr)
## genderM -0.374
- 為何性別的隨機效果（看隨機斜率變異數）很大？
- 這顯示同樣是性別，但在不同學校有不同效果。為什麼？
- 跟男校、女校、混合學校有關嗎？把性別比率加進去看看。
我們在想，性別在不同學校有效果，可能是因為各校校風不同？
rst2 <- lmer(sciscore ~ gender + girlprop + gender:girlprop + (1+gender|schID) , data=dta, REML = FALSE)
summary(rst2)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sciscore ~ gender + girlprop + gender:girlprop + (1 + gender |
## schID)
## Data: dta
##
## AIC BIC logLik deviance df.resid
## 72313.4 72380.1 -36148.7 72297.4 31014
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.2477 -0.6129 0.0544 0.6785 3.4092
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## schID (Intercept) 0.18237 0.4271
## genderM 0.03555 0.1885 0.13
## Residual 0.53085 0.7286
## Number of obs: 31022, groups: schID, 2410
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.18595 0.02748 225.105
## genderM -0.12673 0.03087 -4.105
## girlprop 0.29520 0.04367 6.760
## genderM:girlprop -0.56993 0.06674 -8.539
##
## Correlation of Fixed Effects:
## (Intr) gendrM grlprp
## genderM -0.679
## girlprop -0.896 0.600
## gndrM:grlpr 0.581 -0.923 -0.545
底下是第三部份，診斷（檢查假設）
- 診斷就是要看，資料是否與分析模型的假設一致
- 如果一致，我們對分析結果會比較有信心
- 如果不一致，需要修正模型、或是說明分析結果時要小心
學校與區域層次的隨機效果QQ圖，檢查隨機效果是否呈常態
p_load(ggplot2,lattice,gridExtra)
qq_r1 <- qqmath(~ ranef(rst, condVar = T)$schID, type = c('p', 'g', 'r'), pch = '.',xlab = '常態分位數', ylab = '學校隨機截距分位數')
qq_r2 <- qqmath(~ ranef(rst, condVar = T)$area, type = c('p', 'g', 'r'), pch = '.',xlab = '常態分位數', ylab = '區域隨機截距分位數')
殘差的QQ圖，檢驗殘差是否呈常態
qq_r0 <- qqmath(~ resid(rst, scale = T), type = c('p', 'g', 'r'), pch = '.',xlab = '標準化常態分位數', ylab = '標準化殘差')
預測值對殘差圖，檢查殘差與獨變項無關的假設
rst_f <- fortify.merMod(rst)
r_m0 <- ggplot(data = rst_f, aes(x = .fitted, y = .scresid)) +
geom_point(pch = '.') +
stat_smooth(method = 'loess', se = F) +
labs( x = '預測值', y = '標準化殘差')
把圖放在一起呈現
grid.arrange(qq_r1, qq_r2,qq_r0, r_m0, nrow = 2, ncol = 2)
## `geom_smooth()` using formula = 'y ~ x'
補充部分：交互作用顯著時，以圖顯示
p_load(effects)
m <- mean(dta$girlprop)
s <- sd(dta$girlprop)
ef <- effect(term="gender*girlprop", xlevels= list(girlprop=c(m-s,m,m+s)), mod=rst)
efdta<-as.data.frame(ef)
efdta$girlpropn <- round((efdta$girlprop-m)/s,0)
efdta
## gender girlprop fit se lower upper girlpropn
## 1 F 0.1576881 6.204371 0.02670554 6.152027 6.256715 -1
## 2 M 0.1576881 6.006489 0.02224123 5.962896 6.050083 -1
## 3 F 0.4435562 6.289698 0.02056768 6.249384 6.330011 0
## 4 M 0.4435562 5.921168 0.02035336 5.881274 5.961061 0
## 5 F 0.7294243 6.375024 0.02124523 6.333382 6.416665 1
## 6 M 0.7294243 5.835846 0.02759087 5.781767 5.889925 1
efdta$girlpropn<-as.factor(efdta$girlpropn)
ggplot(efdta, aes(x=gender, y=fit, color=girlpropn,group=girlpropn)) +
geom_point() +
geom_line(linewidth=1.2) +
geom_ribbon(aes(ymin=fit-se, ymax=fit+se, fill=girlpropn),alpha=0.3) +
labs(title = "Interaction Plot", x= "Gender", y="Science score", color="girlprop", fill="girlprop") +
theme_classic() + theme(text=element_text(size=20))
