原文链接： http://tecdat.cn/?p=3015
介绍
首先，请注意，围绕多级模型的术语非常不一致。例如，多级模型本身可以称为分级线性模型，随机效应模型，多级模型，随机截距模型，随机斜率模型或汇集模型。根据学科，使用的软件和学术文献，许多这些术语可能指的是相同的一般建模策略。
读入数据
多级模型适用于特定类型的数据结构，其中单元嵌套在组内（通常为5个以上组），并且我们希望对数据的组结构进行建模。
## 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
在这里，我们有关于嵌套在课堂内和学校内的科目外向性的数据。
在开始之前让我们先了解一下数据的结构：
## 'data.frame': 1200 obs. of 7 variables:
## $ id : int 1 2 3 4 5 6 7 8 9 10 ...
## $ extro : num 63.7 69.5 79.7 63 64.2 ...
## $ open : num 43.4 46.9 32.3 44.4 36.9 ...
## $ agree : num 38 31.5 40.2 30.5 37.4 ...
## $ social: num 75.1 98.1 116.3 90.5 98.5 ...
## $ class : Factor w/ 4 levels "a","b","c","d": 4 1 4 3 4 4 4 4 1 2 ...
## $ school: Factor w/ 6 levels "I","II","III",..: 4 6 6 4 4 1 3 4 3 1 ...
在这里，我们看到我们有两个可能的分组变量 - class和school。让我们进一步探索它们：
##
## a b c d
## 300 300 300 300
##
## I II III IV V VI
## 200 200 200 200 200 200
##
## I II III IV V VI
## a 50 50 50 50 50 50
## b 50 50 50 50 50 50
## c 50 50 50 50 50 50
## d 50 50 50 50 50 50
这是一个完美平衡的数据集。您很可能没有使用完美平衡的数据集，但我们将在未来探讨其中的含义。现在，让我们稍微绘制一下数据。使用包中的优秀xyplot功能lattice，我们可以跨变量探索学校和班级之间的关系。
在这里，我们看到，班内有明显的分层，我们也看到，social变量是从强不同open和agree变量。我们还看到，班级a和班级d在最低和最高频段分别有更多的传播。让我们接下来绘制数据school。
通过学校我们看到学生紧密分组，但学校I和学校的VI分散程度远远高于其他学校。我们的预测因子中的相同模式在学校之间就像在课堂之间一样。
在这里我们可以看到，学校和阶级似乎在密切区分我们的预测者和外向性之间的关系。
探索merMod对象的内部
在上一个教程中，我们为嵌套数据拟合了一系列随机拦截模型。我们lmerMod将更深入地研究在拟合此模型时生成的对象，以便了解如何使用R中的混合效果模型。我们首先拟合下面按类分组的基本示例：
## [1] "lmerMod"
## attr(,"package")
## [1] "lme4"
首先，我们看到它MLexamp1现在是该类的R对象lmerMod。这个lmerMod对象是一个S4类，为了探索它的结构，我们使用slotNames：
## [1] "resp" "Gp" "call" "frame" "flist" "cnms" "lower"
## [8] "theta" "beta" "u" "devcomp" "pp" "optinfo"
在lmerMod对象中，我们看到了许多我们可能希望探索的对象。要查看其中的任何一个，我们只需键入MLexamp1@然后键入插槽名称本身即可。例如：
## lmer(formula = extro ~ open + agree + social + (1 | school),
## data = lmm.data)
## [1] 59.116514 0.009751 0.027788 -0.002151
## [1] "data.frame"
## extro open agree social school
## 1 63.69 43.43 38.03 75.06 IV
## 2 69.48 46.87 31.49 98.13 VI
## 3 79.74 32.27 40.21 116.34 VI
## 4 62.97 44.41 30.51 90.47 IV
## 5 64.25 36.86 37.44 98.52 IV
## 6 50.97 46.26 38.83 75.22 I
该merMod对象有许多可用的方法 - 在这里枚举太多。但是，我们将在下面的列表中介绍一些更常见的内容：
## [1] anova.merMod* as.function.merMod* coef.merMod*
## [4] confint.merMod deviance.merMod* df.residual.merMod*
## [7] drop1.merMod* extractAIC.merMod* extractDIC.merMod*
## [10] family.merMod* fitted.merMod* fixef.merMod*
## [13] formula.merMod* isGLMM.merMod* isLMM.merMod*
## [16] isNLMM.merMod* isREML.merMod* logLik.merMod*
## [19] model.frame.merMod* model.matrix.merMod* ngrps.merMod*
## [22] nobs.merMod* plot.merMod* predict.merMod*
## [25] print.merMod* profile.merMod* qqmath.merMod*
## [28] ranef.merMod* refit.merMod* refitML.merMod*
## [31] residuals.merMod* sigma.merMod* simulate.merMod*
## [34] summary.merMod* terms.merMod* update.merMod*
## [37] VarCorr.merMod* vcov.merMod weights.merMod*
##
## Non-visible functions are asterisked
常见的需求是从merMod对象中提取固定效果。fixef提取固定效果的命名数字向量，这很方便。
## (Intercept) open agree social
## 59.116514 0.009751 0.027788 -0.002151
如果您想了解这些参数的p值或统计显着性，请首先查看lme4帮助，?mcmcsamp了解各种方法的执行情况。内置于包中的一种便捷方式是：
## 0.5 % 99.5 %
## .sig01 4.91840 23.88758
## .sigma 2.53287 2.81456
## (Intercept) 46.27751 71.95610
## open -0.02465 0.04415
## agree -0.01164 0.06721
## social -0.01493 0.01063
从这里我们可以首先看到我们的固定效果参数重叠0表示没有效果的证据。我们还可以看到.sig01，这是我们对随机效应变化的估计，是非常大且非常广泛的定义。这表明我们的团队之间可能缺乏精确性 - 要么是因为群体之间的群体影响很小，要么得到更精确的估计的群体太少，我们每个群体中的单位太少，或者所有群体的组合都是以上。
另一个常见的需求是提取残余标准误差，这是计算效果大小所必需的。要获得残差标准误的命名向量：
## [1] 2.671
例如，教育研究中的常见做法是通过将固定效应参数除以残差标准误差来将固定效果标准化为“效果大小”，这可以lme4很容易地实现：
## (Intercept) open agree social
## 22.1336705 0.0036508 0.0104042 -0.0008055
从这一点，我们可以看到，我们对开放性，宜人性和社交性的预测因素在预测外向性方面几乎毫无用处 - 正如我们的情节所显示的那样。让我们把注意力转向下一个随机效应。
探索组变化和随机效果
您很可能适合混合效果模型，因为您直接对模型中的组级变化感兴趣。目前还不清楚如何从结果中探索这种群体水平的变化summary.merMod。我们从这个输出得到的是组效应的方差和标准偏差，但我们不会对个别组产生影响。这是ranef功能派上用场的地方。
## $school
## (Intercept)
## I -14.091
## II -6.183
## III -1.971
## IV 1.966
## V 6.331
## VI 13.948
运行该ranef功能为我们提供了每个学校的拦截，但没有太多额外的信息 - 例如这些估计的精确度。为此，我们需要一些额外的命令：
## [1] "ranef.mer"
## , , 1
##
## [,1]
## [1,] 0.03565
##
## , , 2
##
## [,1]
## [1,] 0.03565
##
## , , 3
##
## [,1]
## [1,] 0.03565
##
## , , 4
##
## [,1]
## [1,] 0.03565
##
## , , 5
##
## [,1]
## [1,] 0.03565
##
## , , 6
##
## [,1]
## [1,] 0.03565
该ranef.mer对象是一个列表，其中包含每个组级别的data.frame。数据框包含每个组的随机效果（这里我们只对每个学校进行拦截）。当我们要求lme4随机效应的条件方差时，它被存储在attribute那些数据帧的一个中作为方差 - 协方差矩阵的列表。
这种结构确实_很复杂_，但它很强大，因为它允许嵌套，分组和跨级随机效果。此外，创建者lme4已经为用户提供了一些简单的快捷方式，以便从ranef.mer对象中获得他们真正感兴趣的内容。
## $school
## (Intercept)
## I -14.091
## II -6.183
## III -1.971
## IV 1.966
## V 6.331
## VI 13.948
##
## with conditional variances for "school"
## $school
此图显示了dotplot随机效应。
使用模拟和图来探索随机效应
一种常见的计量经济学方法是创建所谓的集团级术语的经验贝叶斯估计。不幸的是，关于什么构成随机效应项的适当标准误差甚至如何一致地定义经验贝叶斯估计，没有太多的一致意见。
## X1 X2 mean median sd
## 1 I (Intercept) -14.053 -14.093 3.990
## 2 II (Intercept) -6.141 -6.122 3.988
## 3 III (Intercept) -1.934 -1.987 3.981
## 4 IV (Intercept) 2.016 2.051 3.993
## 5 V (Intercept) 6.378 6.326 3.981
## 6 VI (Intercept) 13.992 14.011 3.987
该REsim函数为每个学校返回级别名称X1，估计名称，X2估计值的平均值，中位数和估计的标准差。
另一个便利功能可以帮助我们绘制这些结果，看看他们如何与以下结果进行比较dotplot：
这提供了对随机效应分量之间的变化的更保守的观点。根据您的数据收集方式和研究问题，可以采用其他方法来估算这些影响大小。但是，请谨慎行事。
作者推荐的另一种方法lme4涉及RLRsim包。使用该软件包，我们可以测试随机效应的包含是否改善了模型拟合，我们可以使用基于模拟的似然比检验来评估其他随机效应项的p值。
##
## simulated finite sample distribution of LRT. (p-value based on
## 10000 simulated values)
##
## data:
## LRT = 2958, p-value < 2.2e-16
这里exactLRT发出警告，因为我们最初使用REML而不是完全最大可能性来拟合模型。幸运的是，该refitML功能lme4允许我们使用完全最大可能性轻松地重新调整我们的模型，以便轻松地进行精确测试。
##
## simulated finite sample distribution of LRT. (p-value based on
## 10000 simulated values)
##
## data:
## LRT = 2958, p-value < 2.2e-16
在这里，我们可以看到包含我们的分组变量是重要的，即使每个单独的组的影响可能实质上很小和/或不精确地测量。这对于理解模型的正确规范很重要。我们的下一个教程将更详细地介绍这样的规范测试。
随机效应至关重要？
如何解释我们随机效应的_实质性_影响？这在尝试使用多级结构来理解分组可能对个体观察产生的影响的观察工作中通常是至关重要的。为此，我们选择了12个随机病例，然后我们模拟了extro它们在6所学校中的每一所学校的预测值。注意，这是一个非常简单的模拟，仅使用固定效应的平均值和随机效应的条件模式，而不是复制或采样以获得可变性的感觉。这将留给读者和/或未来的教程练习！
现在我们已经建立了一个模拟数据帧，
这个图表告诉我们，在每个情节中，代表一个案例，学校有很大的变化。因此，将每个学生转移到另一所学校对外向学分数有很大影响。但是，每个学校的每个案例都有所不同吗？
在这里我们可以清楚地看到，在每个学校中，案例相对相同，表明群体效应大于个体效应。
这些图可用于以实质方式证明群体和个体效果的相对重要性。可以做更多的事情来使图表更具信息性，例如放置对结果的总可变性的参考，并且还观察距离，移动组将每个观察值从其真实值移开。
结论
lme4提供了一个非常强大的面向对象的工具集，用于处理R中的混合效果模型。理解lme4对象的模型拟合和置信区间需要一些勤奋的研究和使用各种函数和扩展lme4本身。在下一个教程中，我们将探索如何lme4为难以指定的模型确定随机效应模型的适当规范和框架的贝叶斯扩展。我们还将探讨广义线性模型框架和glmer多级广义线性建模的功能。
如果您有任何疑问，请在下面发表评论。
大数据部落 -中国专业的第三方数据服务提供商，提供定制化的一站式数据挖掘和统计分析咨询服务
统计分析和数据挖掘咨询服务：y0.cn/teradat（咨询服务请联系官网客服）
QQ：3025393450
QQ交流群：186388004 
【服务场景】
科研项目; 公司项目外包;线上线下一对一培训;数据爬虫采集;学术研究;报告撰写;市场调查。
【大数据部落】提供定制化的一站式数据挖掘和统计分析咨询
欢迎选修我们的R语言数据分析挖掘必知必会课程！
欢迎关注 微信公众号，了解更多数据干货资讯！
