处理分组数据和复杂层次结构的分析师，从嵌入在参与者中的测量，嵌套在州内的县或嵌套在教室内的学生，经常发现他们需要建模工具来反映他们数据的这种结构。
在R中，有两种主要的方法来拟合多级模型，这些模型考虑了数据中的这种结构。
可下载资源
介绍
这些教程将向用户展示如何使用lme4R中的包来拟合线性和非线性混合效果模型，以及如何使用rstan以完全适合贝叶斯多级模型。这里的重点是如何使模型适合R而不是模型背后的理论。有关多级建模的背景知识，请参阅参考资料。
本教程将介绍如何lme4 设置和运行一些基本模型，其中包括：
- 在R中构造变化的截距，变化的斜率以及变化的斜率和截距模型
- 从混合效应模型中生成预测和解释参数
- 广义和非线性多层次模型
- 完全贝叶斯多级模型适合
rstan或其他MCMC方法
设置 环境
在R中开始多级建模很简单。lme4是在R中实现多级模型的规范包，尽管有许多包依赖并增强其功能集，包括贝叶斯扩展。lme4 最近已被重写以提高速度并整合C ++代码库，因此封装的功能有些不断变化。
要安装lme4，我们只需运行：
# 主要版本
install.packages("lme4")
# 或安装开发版本
library(devtools)
install_github("lme4", user = "lme4")
×
如果对不满足正态性、独立性、方差齐性三个适用条件的资料采用传统的分析方法，对所有样本一视同仁，建立回归方程，就会带来如下问题：
(1)参数估计值不再具有最小方差线性无偏性。
(2)会严重低估回归系数的标准误差。
(3)容易导致估计值过高，使常用的检验失效，从而增加统计检验I型错误发生的概率。
如果我们对不同的群体分别建立各自的回归模型，当群体数较少，群体内样本容量较大，传统的分析方法可能是有效的。或者，我们的兴趣仅在于对这些群体分别做一些统计推断时，也适合用这种方法。但是如果我们把这些群体看成是从总体中抽样来的一个样本(例如多阶段抽样和重复测度数据)，并想分析不同群体之间的总体差异，那么简单地使用传统的统计方法是不够的。同样，如果一些群体包含的样本容量较少，对这些群体做出的推断也不可靠。因此，我们需要把这些群体看成是从总体抽样来的样本，并使用样本总体的信息来进行推断。
本文所讨论的混合线性模型既保留了传统线性模型中的正态性假定条件，又对独立性和方差齐性不作要求，从而扩大了传统线性模型的适用范围 。
读入数据
多级模型适用于特定类型的数据结构，其中单元嵌套在组内（通常为5个以上组），并且我们希望对数据的组结构进行建模。对于我们的介绍性示例，我们将从lme4文档中的一个简单示例开始，并解释模型正在执行的操作。
library(lme4) # 加载库
library(arm) # 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具有熟悉的公式接口的函数， 使用特殊语法指定组级变量：(1|school) ，使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语言拟合线性混合效应模型、固定效应、随机效应参数估计可视化生物生长、发育、繁殖影响因素
【视频】因子分析简介及R语言应用实例：对地区经济研究分析重庆市经济指标
R语言宏观经济学：IS-LM曲线可视化货币市场均衡
R语言代做编程辅导和解答GLM Coursework
