原文鏈接：http://tecdat.cn/?p=3059
介紹
處理分組數據和復雜層次結構的分析師，從嵌入在參與者中的測量，嵌套在州內的縣或嵌套在教室內的學生，經常發現他們需要建模工具來反映他們數據的這種結構。在R中，有兩種主要的方法來擬合多級模型，這些模型考慮了數據中的這種結構。這些教程將向用戶展示如何使用lme4R中的包來擬合線性和非線性混合效果模型，以及如何使用rstan以完全適合貝葉斯多級模型。這里的重點是如何使模型適合R而不是模型背后的理論。有關多級建模的背景知識，請參閱參考資料。
本教程將介紹如何lme4 設置和運行一些基本模型，其中包括：
- 在R中構造變化的截距，變化的斜率以及變化的斜率和截距模型
- 從混合效應模型中生成預測和解釋參數
- 廣義和非線性多層次模型
- 完全貝葉斯多級模型適合
rstan或其他MCMC方法
設置 環境
在R中開始多級建模很簡單。lme4是在R中實現多級模型的規范包，盡管有許多包依賴並增強其功能集，包括貝葉斯擴展。lme4 最近已被重寫以提高速度並整合C ++代碼庫，因此封裝的功能有些不斷變化。
要安裝lme4，我們只需運行：
# Main version
install.packages("lme4")
# Or to install the dev version
library(devtools)
install_github("lme4", user = "lme4")
讀入數據
多級模型適用於特定類型的數據結構，其中單元嵌套在組內（通常為5個以上組），並且我們希望對數據的組結構進行建模。對於我們的介紹性示例，我們將從lme4文檔中的一個簡單示例開始，並解釋模型正在執行的操作。
library(lme4) # load library
library(arm) # convenience functions for regression in R
# summary(lmm.data)
head(lmm.data)
## id extro open agree social class school
## 1 1 63.69 43.43 38.03 75.06 d IV
## 2 2 69.48 46.87 31.49 98.13 a VI
## 3 3 79.74 32.27 40.21 116.34 d VI
## 4 4 62.97 44.41 30.51 90.47 c IV
## 5 5 64.25 36.86 37.44 98.52 d IV
## 6 6 50.97 46.26 38.83 75.22 d I
模型
讓我們首先擬合一個簡單的OLS回歸 。
OLSexamp <- lm(extro ~ open + agree + social, data = lmm.data)
display(OLSexamp)
## lm(formula = extro ~ open + agree + social, data = lmm.data)
## coef.est coef.se
## (Intercept) 57.84 3.15
## open 0.02 0.05
## agree 0.03 0.05
## social 0.01 0.02
## ---
## n = 1200, k = 4
## residual sd = 9.34, R-Squared = 0.00
R模型接口非常簡單，首先指定因變量，然后是 ~符號。右手邊，預測變量，每個都被命名。加法符號表明這些被建模為加性效應。最后，我們指定要計算模型的數據幀。這里我們使用該lm函數執行OLS回歸，但R中還有許多其他選項。
如果我們想要提取諸如AIC之類的度量 。
MLexamp <- glm(extro ~ open + agree + social, data = lmm.data)
display(MLexamp)
## glm(formula = extro ~ open + agree + social, data = lmm.data)
## coef.est coef.se
## (Intercept) 57.84 3.15
## open 0.02 0.05
## agree 0.03 0.05
## social 0.01 0.02
## ---
## n = 1200, k = 4
## residual deviance = 104378.2, null deviance = 104432.7 (difference = 54.5)
## overdispersion parameter = 87.3
## residual sd is sqrt(overdispersion) = 9.34
AIC(MLexamp)
## [1] 8774
這導致模型擬合較差。現在讓我們看一個簡單的變化 模型。
適合不同的 模型
根據學科規范，我們的下一步可能是使用分組變量（如學校或班級）來擬合不同的 模型。
MLexamp.2 <- glm(extro ~ open + agree + social + class, data = lmm.data)
display(MLexamp.2)
## glm(formula = extro ~ open + agree + social + class, data = lmm.data)
## coef.est coef.se
## (Intercept) 56.05 3.09
## open 0.03 0.05
## agree -0.01 0.05
## social 0.01 0.02
## classb 2.06 0.75
## classc 3.70 0.75
## classd 5.67 0.75
## ---
## n = 1200, k = 7
## residual deviance = 99187.7, null deviance = 104432.7 (difference = 5245.0)
## overdispersion parameter = 83.1
## residual sd is sqrt(overdispersion) = 9.12
AIC(MLexamp.2)
## [1] 8719
anova(MLexamp, MLexamp.2, test = "F")
## Analysis of Deviance Table
##
## Model 1: extro ~ open + agree + social
## Model 2: extro ~ open + agree + social + class
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 1196 104378
## 2 1193 99188 3 5190 20.8 3.8e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
這通常被稱為固定效應規范。
MLexamp.3 <- glm(extro ~ open + agree + social + school, data = lmm.data)
display(MLexamp.3)
## glm(formula = extro ~ open + agree + social + school, data = lmm.data)
## coef.est coef.se
## (Intercept) 45.02 0.92
## open 0.01 0.01
## agree 0.03 0.02
## social 0.00 0.00
## schoolII 7.91 0.27
## schoolIII 12.12 0.27
## schoolIV 16.06 0.27
## schoolV 20.43 0.27
## schoolVI 28.05 0.27
## ---
## n = 1200, k = 9
## residual deviance = 8496.2, null deviance = 104432.7 (difference = 95936.5)
## overdispersion parameter = 7.1
## residual sd is sqrt(overdispersion) = 2.67
AIC(MLexamp.3)
## [1] 5774
anova(MLexamp, MLexamp.3, test = "F")
## Analysis of Deviance Table
##
## Model 1: extro ~ open + agree + social
## Model 2: extro ~ open + agree + social + school
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 1196 104378
## 2 1191 8496 5 95882 2688 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
學校效果大大提高了我們的模型擬合度。但是，我們如何解釋這些影響呢？
table(lmm.data$school, lmm.data$class)
##
## a b c d
## I 50 50 50 50
## II 50 50 50 50
## III 50 50 50 50
## IV 50 50 50 50
## V 50 50 50 50
## VI 50 50 50 50
在這里，我們可以看到我們有一個完美平衡的設計，在課堂和學校的每個組合中有50個觀察結果（ 。
讓我們嘗試對這些獨特的細胞進行建模。
MLexamp.4 <- glm(extro ~ open + agree + social + school:class, data = lmm.data)
display(MLexamp.4)
## glm(formula = extro ~ open + agree + social + school:class, data = lmm.data)
## coef.est coef.se
## (Intercept) 80.36 0.37
## open 0.01 0.00
## agree -0.01 0.01
## social 0.00 0.00
## schoolI:classa -40.39 0.20
## schoolII:classa -28.15 0.20
## schoolIII:classa -23.58 0.20
## schoolIV:classa -19.76 0.20
## schoolV:classa -15.50 0.20
## schoolVI:classa -10.46 0.20
## schoolI:classb -34.60 0.20
## schoolII:classb -26.76 0.20
## schoolIII:classb -22.59 0.20
## schoolIV:classb -18.71 0.20
## schoolV:classb -14.31 0.20
## schoolVI:classb -8.54 0.20
## schoolI:classc -31.86 0.20
## schoolII:classc -25.64 0.20
## schoolIII:classc -21.58 0.20
## schoolIV:classc -17.58 0.20
## schoolV:classc -13.38 0.20
## schoolVI:classc -5.58 0.20
## schoolI:classd -30.00 0.20
## schoolII:classd -24.57 0.20
## schoolIII:classd -20.64 0.20
## schoolIV:classd -16.60 0.20
## schoolV:classd -12.04 0.20
## ---
## n = 1200, k = 27
## residual deviance = 1135.9, null deviance = 104432.7 (difference = 103296.8)
## overdispersion parameter = 1.0
## residual sd is sqrt(overdispersion) = 0.98
AIC(MLexamp.4)
## [1] 3396
這非常有用，但如果我們想了解學校的影響和課堂的影響，以及學校和班級的影響，該怎么辦？
MLexamp.5 <- glm(extro ~ open + agree + social + school * class - 1, data = lmm.data)
display(MLexamp.5)
## glm(formula = extro ~ open + agree + social + school * class -
## 1, data = lmm.data)
## coef.est coef.se
## open 0.01 0.00
## agree -0.01 0.01
## social 0.00 0.00
## schoolI 39.96 0.36
## schoolII 52.21 0.36
## schoolIII 56.78 0.36
## schoolIV 60.60 0.36
## schoolV 64.86 0.36
## schoolVI 69.90 0.36
## classb 5.79 0.20
## classc 8.53 0.20
## classd 10.39 0.20
## schoolII:classb -4.40 0.28
## schoolIII:classb -4.80 0.28
## schoolIV:classb -4.74 0.28
## schoolV:classb -4.60 0.28
## schoolVI:classb -3.87 0.28
## schoolII:classc -6.02 0.28
## schoolIII:classc -6.54 0.28
## schoolIV:classc -6.36 0.28
## schoolV:classc -6.41 0.28
## schoolVI:classc -3.65 0.28
## schoolII:classd -6.81 0.28
## schoolIII:classd -7.45 0.28
## schoolIV:classd -7.24 0.28
## schoolV:classd -6.93 0.28
## schoolVI:classd 0.06 0.28
## ---
## n = 1200, k = 27
## residual deviance = 1135.9, null deviance = 4463029.9 (difference = 4461894.0)
## overdispersion parameter = 1.0
## residual sd is sqrt(overdispersion) = 0.98
AIC(MLexamp.5)
## [1] 3396
探索隨機斜率
另一種選擇是為每個學校和班級組合安裝一個單獨的模型。如果我們認為我們的變量之間的關系可能高度依賴於學校和班級組合，我們可以簡單地擬合一系列模型並探索它們之間的參數變化：
require(plyr)
display(modellist[[1]])
## glm(formula = extro ~ open + agree + social, data = x)
## coef.est coef.se
## (Intercept) 35.87 5.90
## open 0.05 0.09
## agree 0.02 0.10
## social 0.01 0.03
## ---
## n = 50, k = 4
## residual deviance = 500.2, null deviance = 506.2 (difference = 5.9)
## overdispersion parameter = 10.9
## residual sd is sqrt(overdispersion) = 3.30
display(modellist[[2]])
## glm(formula = extro ~ open + agree + social, data = x)
## coef.est coef.se
## (Intercept) 47.96 2.16
## open -0.01 0.03
## agree -0.03 0.03
## social -0.01 0.01
## ---
## n = 50, k = 4
## residual deviance = 47.9, null deviance = 49.1 (difference = 1.2)
## overdispersion parameter = 1.0
## residual sd is sqrt(overdispersion) = 1.02
我們將在未來的教程中更深入地討論此策略，包括如何在此命令中生成的模型列表中進行性能推斷。
使用lmer安裝不同的斜率模型
輸入lme4。雖然上述所有技術都是解決這一問題的有效方法，但當我們明確感興趣的是群體之間的變化時，它們並不一定是最好的方法。這是混合效果建模框架有用的地方。現在我們使用lmer具有熟悉的公式接口的函數，但現在使用特殊語法指定組級變量：(1|school) tell 使用變量lmer擬合具有變量截距組效果的線性模型school。
display(MLexamp.6)
## lmer(formula = extro ~ open + agree + social + (1 | school),
## data = lmm.data)
## coef.est coef.se
## (Intercept) 59.12 4.10
## open 0.01 0.01
## agree 0.03 0.02
## social 0.00 0.00
##
## Error terms:
## Groups Name Std.Dev.
## school (Intercept) 9.79
## Residual 2.67
## ---
## number of obs: 1200, groups: school, 6
## AIC = 5836.1, DIC = 5789
## deviance = 5806.5
我們可以使用多個組效果術語來適應多個組效果。
display(MLexamp.7)
## lmer(formula = extro ~ open + agree + social + (1 | school) +
## (1 | class), data = lmm.data)
## coef.est coef.se
## (Intercept) 60.20 4.21
## open 0.01 0.01
## agree -0.01 0.01
## social 0.00 0.00
##
## Error terms:
## Groups Name Std.Dev.
## school (Intercept) 9.79
## class (Intercept) 2.41
## Residual 1.67
## ---
## number of obs: 1200, groups: school, 6; class, 4
## AIC = 4737.9, DIC = 4683
## deviance = 4703.6
最后，我們可以通過以下語法擬合嵌套的組效果術語：
display(MLexamp.8)
## lmer(formula = extro ~ open + agree + social + (1 | school/class),
## data = lmm.data)
## coef.est coef.se
## (Intercept) 60.24 4.01
## open 0.01 0.00
## agree -0.01 0.01
## social 0.00 0.00
##
## Error terms:
## Groups Name Std.Dev.
## class:school (Intercept) 2.86
## school (Intercept) 9.69
## Residual 0.98
## ---
## number of obs: 1200, groups: class:school, 24; school, 6
## AIC = 3568.6, DIC = 3508
## deviance = 3531.1
在這里(1|school/class)，我們想要為1|學校和學校內嵌的課程設置不同截取的混合效應術語。
用lmer擬合變化的斜率模型
但是，如果我們想要探索不同學生水平指標的影響，因為它們因教室而異。我們可以適應不同的坡度模型，而不是按學校（或學校/班級）擬合獨特的模型。在這里，我們修改我們的隨機效應項，以在分組術語之前包含變量：(1 + open|school/class)告訴R適合變化的斜率和不同的學校和學校類別的攔截模型，並允許open變量的斜率因學校而異。
display(MLexamp.9)
## lmer(formula = extro ~ open + agree + social + (1 + open | school/class),
## data = lmm.data)
## coef.est coef.se
## (Intercept) 60.26 3.93
## open 0.01 0.01
## agree -0.01 0.01
## social 0.00 0.00
##
## Error terms:
## Groups Name Std.Dev. Corr
## class:school (Intercept) 2.62
## open 0.01 1.00
## school (Intercept) 9.51
## open 0.00 1.00
## Residual 0.98
## ---
## number of obs: 1200, groups: class:school, 24; school, 6
## AIC = 3574.7, DIC = 3506
## deviance = 3529.3
結論
在R語言和生態系統中，擬合混合效應模型和探索組級變異非常容易。在以后的教程中，我們將探索跨模型的比較，使用混合效果模型進行推理，以及創建混合效果模型的圖形表示以了解它們的效果。
附錄
## R version 3.0.1 (2013-05-16)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] plyr_1.8 arm_1.6-10 MASS_7.3-29 lme4_1.0-5
## [5] Matrix_1.1-0 lattice_0.20-24 knitr_1.5
##
## loaded via a namespace (and not attached):
## [1] abind_1.4-0 coda_0.16-1 evaluate_0.5.1 formatR_0.10
## [5] grid_3.0.1 minqa_1.2.1 nlme_3.1-113 splines_3.0.1
## [9] stringr_0.6.2 tools_3.0.1