Multi level with R
邱浩恩
2016年11月24日
隨機截距模型
載入套件，進行不同截距的層級分析
library(lme4)## Loading required package: Matrixlibrary(ggplot2)
library(lattice)
ML<- read.table("D:/104/ML_R/ML.csv",
header=TRUE, sep=",")先將斜率視為相同–>
y=截距j+斜率*xij+誤差ij
j=不同間學校,i=不位學生
y=(固定效果+隨機效果j)+斜率*xij+誤差ij
和lm迴歸寫法相似，(1 | 學校)代表截距在不同學校間的隨機效果
整體模型:
summary(m0 <- lmer(reading_3 ~ reading_1_center + gender + reading_1_center:gender +
( 1 | school_3 ) + ( 1 | senson ), data = ML ) ) ## Linear mixed model fit by REML ['lmerMod']
## Formula:
## reading_3 ~ reading_1_center + gender + reading_1_center:gender +
## (1 | school_3) + (1 | senson)
## Data: ML
##
## REML criterion at convergence: 1415.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5075 -0.5793 -0.0332 0.6254 2.8518
##
## Random effects:
## Groups Name Variance Std.Dev.
## school_3 (Intercept) 0.3448 0.5872
## senson (Intercept) 0.0000 0.0000
## Residual 88.7883 9.4228
## Number of obs: 193, groups: school_3, 33; senson, 2
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 499.87180 1.01387 493.0
## reading_1_center 0.78066 0.04775 16.3
## gendermale -1.08637 1.36946 -0.8
## reading_1_center:gendermale -0.01785 0.06414 -0.3
##
## Correlation of Fixed Effects:
## (Intr) rdn_1_ gndrml
## rdng_1_cntr -0.099
## gendermale -0.726 0.076
## rdng_1_cnt: 0.075 -0.744 -0.016去除senson的隨機效果，並看看去除隨機效果是否顯著
m1 <- update(m0, . ~ . - ( 1 | senson ) )
anova(m0, m1)## refitting model(s) with ML (instead of REML)## Data: ML
## Models:
## m1: reading_3 ~ reading_1_center + gender + (1 | school_3) + reading_1_center:gender
## m0: reading_3 ~ reading_1_center + gender + reading_1_center:gender +
## m0: (1 | school_3) + (1 | senson)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 6 1422.2 1441.8 -705.11 1410.2
## m0 7 1424.2 1447.0 -705.11 1410.2 0 1 1去除交互作用項
drop1(m0, test = 'Chisq')## Single term deletions
##
## Model:
## reading_3 ~ reading_1_center + gender + reading_1_center:gender +
## (1 | school_3) + (1 | senson)
## Df AIC LRT Pr(Chi)
## <none> 1424.2
## reading_1_center:gender 1 1422.3 0.06964 0.7919summary(m0)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## reading_3 ~ reading_1_center + gender + reading_1_center:gender +
## (1 | school_3) + (1 | senson)
## Data: ML
##
## REML criterion at convergence: 1415.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5075 -0.5793 -0.0332 0.6254 2.8518
##
## Random effects:
## Groups Name Variance Std.Dev.
## school_3 (Intercept) 0.3448 0.5872
## senson (Intercept) 0.0000 0.0000
## Residual 88.7883 9.4228
## Number of obs: 193, groups: school_3, 33; senson, 2
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 499.87180 1.01387 493.0
## reading_1_center 0.78066 0.04775 16.3
## gendermale -1.08637 1.36946 -0.8
## reading_1_center:gendermale -0.01785 0.06414 -0.3
##
## Correlation of Fixed Effects:
## (Intr) rdn_1_ gndrml
## rdng_1_cntr -0.099
## gendermale -0.726 0.076
## rdng_1_cnt: 0.075 -0.744 -0.016變異成份
抽取變異成分，計算學校與區域可以解釋部分（ICCs）
print(vc <- VarCorr(m0), comp = 'Variance' )## Groups Name Variance
## school_3 (Intercept) 0.34484
## senson (Intercept) 0.00000
## Residual 88.78830vc <- as.data.frame(vc)
head(vc)## grp var1 var2 vcov sdcor
## 1 school_3 (Intercept) <NA> 0.3448412 0.5872318
## 2 senson (Intercept) <NA> 0.0000000 0.0000000
## 3 Residual <NA> <NA> 88.7883012 9.4227544school_3變異(school_3隨機截距可以佔總變異量的百分比)
cbind(vc$grp,vc$vcov/sum(vc$vcov))## [,1] [,2]
## [1,] "school_3" "0.00386883258608399"
## [2,] "senson" "0"
## [3,] "Residual" "0.996131167413916"vc[vc$grp=='school_3', 'vcov']/ sum(vc$vcov)## [1] 0.003868833senson變異(senson隨機截距可以佔總變異量的百分比)
vc[vc$grp=='senson', 'vcov']/ sum(vc$vcov)## [1] 0總變異
1 - (vc[vc$grp=='Residual', 'vcov']/ sum(vc$vcov))## [1] 0.003868833QQ plots
學校與區域層次的隨機效果QQ圖，檢查隨機效果是否呈常態
qq_rsch <- qqmath(~ ranef(m0, varCond = T)$school_3, type = c('p', 'g', 'r'), pch = '.',
xlab = '常態分位數', ylab = 'school_3隨機效果分位數',scales=list(cex=1.25))
qq_rsch 
qq_rsen <- qqmath(~ ranef(m0, varCond = T)$senson, type = c('p', 'g', 'r'), pch = '.',
xlab = '常態分位數', ylab = 'senson隨機效果分位數',scales=list(cex=1.25))
qq_rsen
殘差的QQ圖，檢驗殘差是否呈常態
qq_rr <- qqmath(~ resid(m0, scale = T), type = c('p', 'g', 'r'), pch = '.',
xlab = '標準化常態分位數', ylab = '標準化殘差',scales=list(cex=1.25))
qq_rr
預測值對殘差圖
m0_f <- fortify(m0)
r_m0 <- ggplot(data = m0_f, aes(x = .fitted, y = .scresid)) +
geom_point(pch = '.') +
stat_smooth(method = 'loess', se = F) +
labs( x = '預測值', y = '標準化殘差')
r_m0