原文链接：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