处理分组数据和复杂层次结构的分析师，从嵌入在参与者中的测量，嵌套在州内的县或嵌套在教室内的学生，经常发现他们需要建模工具来反映他们数据的这种结构。

由Kaizong Ye，Coin Ge撰写

在R中，有两种主要的方法来拟合多级模型，这些模型考虑了数据中的这种结构。

可下载资源

完整程序、数据和文档（word）

介绍

这些教程将向用户展示如何使用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

自适应网页宽度的 Bilibili 视频

视频

线性混合效应模型(LMM,Linear Mixed Models)和R语言实现案例

探索见解 ➜

去bilibili观看 ➜

探索更多视频 ➜

模型

让我们首先拟合一个简单的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

最受欢迎的见解

1.基于R语言的lmer混合线性回归模型

2.R语言用Rshiny探索lme4广义线性混合模型（GLMM）和线性混合模型（LMM）

3.R语言线性混合效应模型实战案例

4.R语言线性混合效应模型实战案例2

5.R语言线性混合效应模型实战案例

6.线性混合效应模型Linear Mixed-Effects Models的部分折叠Gibbs采样

7.R语言LME4混合效应模型研究教师的受欢迎程度

8.R语言中基于混合数据抽样(MIDAS)回归的HAR-RV模型预测GDP增长

9.使用SAS，Stata，HLM，R，SPSS和Mplus的分层线性模型HLM

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