关于数据挖掘:R语言用线性混合效应多水平层次嵌套模型分析声调高低与礼貌态度的关系附代码数据

全文下载链接：http://tecdat.cn/?p=23681

最近咱们被客户要求撰写对于线性混合效应的钻研报告，包含一些图形和统计输入。

线性混合效应模型与咱们曾经晓得的线性模型有什么不同？

线性混合模型（有时被称为 “ 多层次模型 “ 或 “ 层次模型 ”，取决于上下文）是一种回归模型，它同时思考了（1）被感兴趣的自变量（如 lm()）所解释的变动 – 固定效应，以及（2）不被感兴趣的自变量解释的变动 – 随机效应。因为该模型包含固定效应和随机效应的混合，所以被称为混合模型。这些随机效应实质上赋予误差项 ϵ 构造。

固定效应和随机效应的定义可能会有所不同，所以要留神你在文献中的解释；然而，对于大多数目标来说，如果从所有感兴趣的层面收集了数据，你能够把一个变量视为固定效应因素（例如。性别：男 / 女，条件：易 / 中 / 难，剂量：低 / 高），如果变量有一堆可能的程度，但你只对一个随机的汇合（如受试者、刺激物、教室）进行采样，只管这些样本会有一些特异性，但你个别不会关怀它们，目标是对更宽泛的人群进行概括（如所有的人、所有的场景、所有的教室）。

例子

比方说，你对语言感兴趣，更确切地说，是对声音的高下与礼貌态度的关系感兴趣。你要求你的受试者对假如的场景（IV，受试者外部）做出反馈，这些场景要么是须要礼貌态度的正式场合（例如，给传授一个早退的借口），要么是比拟非正式的场合（例如，向敌人解释你为什么早退），并测量他们的音调（DV）。每个受试者都会失去一份所有场景的清单，因而每个受试者都会给出多个礼貌态度的或非正式的答复。你还留神到每个受试者的性别（IV，受试者之间），因为这是对腔调的另一个重要影响。

在迄今为止咱们所看到的线性模型中，咱们将建设这样的模型。

腔调 = 礼貌态度 + 性别 +ϵ

其中最初一项是咱们的误差项。这个误差项代表了因为咱们无奈在试验中管制的 “ 随机 “ 因素而导致的与咱们预测的偏差。

对于这种数据，因为每个受试者都给出了多个反馈（” 反复测量 “ 设计），咱们能够看到，这将违反线性建模中重要的独立性假如：同一受试者的多个反馈不能被视为彼此独立。在咱们的计划中，每个人的腔调都略有不同，这将成为影响同一受试者所有反馈的特异性因素，从而使这些不同的反馈相互依赖（相干）而非独立。

随机效应

咱们要解决这种状况的办法是为主体增加一个随机效应。这使咱们可能通过为每个受试者假如不同的 “ 基准 “ 音高值来解决这种非独立性。因而，受试者 1 在不同的话语中可能有 233 赫兹的均匀腔调，而受试者 2 可能有 210 赫兹的均匀腔调。在咱们的模型中，咱们通过对受试者的随机效应来解释这些腔调的个体差异。

咱们将一些数据为例进行剖析。

table(subject)

把数据可视化。

qplot(condition, pitch, facets = . ~ subject)

受试者 “F#” 为女性受试者。对象 “M#” 是男性对象。你马上就会发现，男性的声音比女性低（这是能够预期的）。但除此之外，在男性和女性群体中，你会看到很多个体差异，一些人的性别值绝对较高，而另一些人的性别值绝对较低。

来自同一主体的样本的相关性

另一种说法是，在受试者外部，不同条件下的音高存在着相关性。让咱们把它形象化。

用随机截距对个体平均值进行建模

咱们能够通过为每个参与者假如不同的随机截距来建设这些个体差异的模型；每个参与者都被调配了不同的截距值（即不同的均匀腔调），而混合模型基本上是为你预计这些截距。

回过头来看咱们的模型，咱们以前的公式是。

腔调 = 截距 + 礼貌 + 性别 +ϵ

咱们更新后的公式是这样的。

腔调 = 截距 + 礼貌 + 性别 +(1| 个体)+ϵ

“(1|subject) “ 是随机截距的 R 语法。这句话的意思是 “ 假如每个主体的截距都不同 ”…… 而 “1 “ 代表这里的截距。你能够认为这个公式是通知你的模型，它应该冀望每个受试者会有多个反馈，而这些反馈将取决于每个受试者的基准程度。这就无效地解决了因同一受试者有多个反馈而产生的非独立性问题。

请留神，该公式依然蕴含一个个别误差项 ϵ。这是必要的，因为即便咱们思考到了每个主体的变动，同一主体的不同音高之间依然会存在 “ 随机 “ 差别。

点击题目查阅往期内容

R 语言 LME4 混合效应模型钻研老师的受欢迎水平

左右滑动查看更多

01

02

03

04

对不同条件下的不同参与者的平均值有一个概念。

 aggregate(pitch ~ subject, FUN = "mean")

当初用 lmer()，咱们能够预计每个参与者的平均值。为了做到这一点，咱们将为每个受试者蕴含一个随机截距，而后看一下预计的截距。

coef(lmer(pitch ~ (1 | subject))

# 固定效应 + 随机效应的主体
['（截距）'] + subject

请留神，估计值与理论均匀音高相当靠近，咱们能够看到，各受试者的理论均匀音高是估计值（Intercept），而各受试者均匀音高的标准差是随机效应的标准差（Std.Dev）。

# 应用原始数据
mean

## [1] 193

sd

## [1] 63.47

# 应用每个子项目的预计截距
mean(subject[1][,'(Intercept)'])

## [1] 193

sd

## [1] 62.4

# 这也是模型输入中的总结
summary(res1)

包含固定效应

因为咱们预测假如状态的条件（” 非正式 “ 与 “ 礼貌态度 ”）会影响音调（兴许在非正式状态下音调会更高），此外还有受试者的性别（女性的音调可能会更高），让咱们把这些条件纳入模型，同时也思考到每个受试者的随机截距（让截距因受试者而异）。

lmer(音调~ 礼貌 + 性别 +(1| 个体))

音调 j = 截距 + 截距 j + 状态 + 性别

音调 A= 截距 + A 截距的变动〔_Intercept_ _Shift_〕+ 状态 + 性别

ggplot(x=condition, y=mean)

在这里，咱们还将比照状态和性别，这样咱们就能够看到条件在女性和男性之间的 “ 平均值 “ 的影响，以及性别在 “ 非正式 “ 和 “ 不礼貌 “ 之间的平均值的影响。

能够看到，咱们的均匀音调是 192.88，非正式状态的音调比礼貌状态的音调高，b=9.71，t=3.03，女性的音调比男性高，b=54.10，t=5.14。咱们能够应用一个教训法令，如果 t 大于 2，就可能是显著的。

更多模型信息

咱们能够用许多不同类型的信息来评估模型。在比拟模型的时候，这些信息可能很有用 一个有用的衡量标准是 AIC，即偏差 +2∗(p+1)，其中 p 是模型中的参数数量（这里，咱们将参数合成，所以 1 是预计的残差，p 是所有其余参数，例如，固定效应系数 + 预计的随机效应的方差等）。较低的 AIC 比拟好，因为较高的偏差意味着模型不能很好地拟合数据。因为 AIC 随着 p 的减少而减少，所以 AIC 会因为更多的参数而受到惩办。

deviance = -2*logLikelihood; deviance

## [1] 789.5

# 用手计算 AIC
p = 4 # 参数 = 3（固定效应）+1（随机效应
deviance + 2*(p+1) # 总参数 = 4 + 1 用于预计残差

## [1] 799.5

AIC(res2)

## [1] 799.5

提取所有系数

ggplot(aes(x=factor(subject), y=mean_pitch))

coef(res2)

在这里，咱们能够看到，这个模型为每个受试者产生了一个独自的截距，此外还有一个参数估计值 / 斜率的条件和性别，在各受试者中是恒定的。从这里，咱们能够尝试依据这些系数来预计一个给定样本的均匀音高。例如，让咱们尝试用他们预计的截距和作为女性的影响来预计受试者 F1 的平均数（x¯=232.0357）。

179.3003 + 0*(9.7) + 1*(54.10244)

## [1] 233.4

相当靠近

当初是 M3（x¯=168.9786）:

220.3196 + 0*(9.7) + -1*(54.10244)

## [1] 166.2

还不错

随机斜率

在下面的模型中，咱们假如礼貌态度的影响对所有被试都是一样的，因而，礼貌态度的系数是一个。然而，礼貌态度的影响对于不同的受试者主体可能是不同的；也就是说，可能存在着受试者主体礼貌态度的相互作用。例如，咱们能够预期，有些人在有须要礼仪的场景下更有礼貌，有些人则不那么有礼貌。因而，咱们须要的是一个随机斜率模型，在这个模型中，不仅容许主体有不同的截距，而且还容许它们对礼貌的影响有不同的斜率（即状态对音调的不同影响）。

让咱们开始将数据可视化。

依据这幅图，看起来各受试者的斜率是否很不稳固？

当初退出随机斜率。

coef(res3)

anova(res2, res3, refit=FALSE)

减少每个受试者的随机斜率又占用了 2 个自由度，并没有明显改善模型拟合，χ2(2)=0.024,P=0.99。留神 df=2，因为咱们同时退出了斜率方差和截距与斜率之间的相干关系。看一下 AIC 值，更简单的模型的 AIC 值更高，所以咱们想用不太简单（更扼要）的模型。如果咱们看一下预计的斜率，咱们能够看到它们在不同的样本中是十分类似的。所以，咱们不须要在模型中退出条件的随机斜率。

然而，其他人会认为，咱们应该放弃咱们的模型最全面，并且应该包含随机斜率，即便它们没有改善模型!

测试显著性

尽管对是否应该取得 lmer()模型的 p 值有一些争执（例如，这个；大多数争执围绕着如何计算 dfs），但你能够应用 {lmerTest} 包取得 df 的近似值（以及因而取得 p 值）。

获取 P 值

summary(res3b)

将模型输入与 SS/Kenward-Roger appox 进行比拟

anova

anova(res2b)

模型比拟

另一方面，有些人认为，用似然比测验进行模型比拟是测验一个参数是否显著的更好办法。也就是说，如果在你的模型中退出该参数能显著进步模型的拟合度，那么该参数就应该被纳入模型中。

似然比测验实质上通知咱们，数据在更简单模型下的可能性比在简略模型下的可能性大多少（这些模型须要嵌套！）：

D 的散布大概是 χ2，自由度为 df2-df1。咱们要么 “ 手动 “ 做这个计算，要么就间接应用 anova()函数!

anova(res4, res4b)

# 手动计算
dev0 <- -2*logLik(res4) # 较简略的偏差模型
devdiff <- (dev0-dev1) # 偏差差值
dfdiff  # 参数的差别（应用 dfs）。

在这里，咱们比拟了两个嵌套模型，一个没有条件，另一个有条件。通过模型比拟，咱们得出结论，在咱们的模型中退出条件是有必要的，因为它明显改善了模型的拟合，χ2(1)=8.79,P<0.01。

REML 与 ML

让咱们从一个统计模型开始，指定（i）固定效应和（ii）各种随机效应的正态分布的变异和协方差。

• 在 ML（最大似然）预计中，咱们计算上述（i）和（ii）组中任意抉择的参数值的数据的对数（似然）（LL）。而后，咱们寻找能使 L 最大化（或最小化 -L）的参数值。这些最佳参数值被称为 ML 参数估计值。L 的最大值（- 2 倍）被称为模型的偏差。对于某些目标，如形容数据，咱们关注 ML 参数估计值；对于其余目标，如模型比拟，咱们关注偏差。在比拟固定效应不同的模型时，你应该应用 ML，而且你必须包含 lmer(, REML=FALSE)。此外，如果你要比拟一个 lm()和 lmer()模型（即测试是否有必要应用任何随机效应），你也应该应用 ML 预计。
• 在 REML（REstricted ML）预计中，咱们的次要趣味是预计随机效应，而不是固定效应。设想一下，咱们通过将下面第 (i) 组中的固定效应参数设置为某些正当的值来限度咱们的参数空间。在这个限度的空间里，咱们寻找集 (ii) 中的随机效应参数值，使数据的 LL 值最大化；同时留神 LL 值的最大值。而后多次重复这个过程。而后对固定效应参数值、随机效应参数的估计值和 LL 的最大值进行均匀。这种平均法能够失去 REML 参数估计值和 REML 偏差值。因为这个过程对固定效应参数的关注度很低，所以它不应该被用来比拟固定效应构造不同的模型。你应该在比拟随机效应不同的模型时应用这个办法。要做到这一点，你应该养成在运行模型比拟时包含 lmer(, REML=TRUE)，并应用 anova(, refit=FALSE)的习惯。

例子

测试随机效应是否有必要

dev1 <- -2*logLik(res5)
dev0 <- -2*logLik(res5b)
devdiff <- as.numeric(dev0-dev1); devdiff

# 比拟 AICs
AIC(res5b, res5)

AICtab(res5b, res5) #  AIC

用 anova 间接比拟 lm 和 nlme 模型

看起来，是的，纳入受试者的随机截距是有情理的，χ2(1) = 19.51, P < 0.001.

在这里，χ2 散布并不总是一个很好的有效散布的近似值（在测试一些随机效应时过于激进，而在测试一些固定效应时不够激进）。

场景效应

不同的场景可能会引起不同的 “ 音调 “ 值；因而，在不同的受试者之间，甚至在一个受试者外部，在礼貌和非正式的状态下，特定场景的音调可能是相干的。咱们能够把这作为一个随机效应的模型。

ggplot(d, aes(x=scenario, y=pitch)

anova

coef(res4)

与受试者的随机截距类似，当初咱们有了每种场景下的均匀音高程度。

那么，每个场景的斜率变动呢？

  ggplot(byscenario, aes(x=condition, y=pitch))

anova

没有改善模型。这阐明咱们的计划在非正式和礼貌状态下的差别方面可能是类似的。

混合模型阐明

这里有几件重要的事件要讲。你可能会问 “ 我应该指定哪些随机斜率？” 甚至是 “ 随机斜率到底有没有必要？”

从概念上讲，将随机斜率和随机截距一起包含进来是十分有意义的。毕竟，你能够认为人们对试验操纵的反馈不同！同样，你能够认为人们对试验操纵的反馈不同。同样，你总是能够认为，试验操纵的成果对试验中的所有我的项目都不一样。

在上述模型中，咱们的整个钻研关键在于阐明无关礼貌的货色。咱们对性别差异不感兴趣，但它们是很值得管制的。这就是为什么咱们对礼貌态度的影响有随机斜率（按被试和我的项目），而不是性别。换句话说，在礼貌态度对音调的影响方面，咱们只模仿了按主体和按我的项目的变动。

在线性模型背景下探讨的所有都间接实用于混合模型。因而，你还必须放心共线性和异样值。你还得放心同方差（方差相等）和潜在的正态性缺失问题。

独立性，作为最重要的假如，须要一个非凡的词。咱们转向混合模型而不是仅仅应用线性模型的次要起因之一是为了解决咱们数据中的非依赖性。然而，混合模型依然能够违反独立性。如果你短少重要的固定或随机效应。因而，例如，如果咱们用一个不包含随机效应 “ 主体 “ 的模型来剖析咱们的数据，那么咱们的模型就不会 “ 晓得 “ 每个主体有多个反馈。这就相当于违反了独立假如。因而，要认真抉择你的固定效应和随机效应，解决非独立性问题。

其余一些阐明。

如果你的因变量是 …

• 间断：应用混合效应的线性回归模型
• 二元：应用混合效应的 Logistic 回归模型

函数 lmer 用于拟合线性混合模型，函数 glmer 用于拟合狭义（非高斯）线性混合模型。

点击文末 “浏览原文”

获取全文残缺材料。

本文选自《R 语言用线性混合效应（多程度 / 档次 / 嵌套）模型剖析腔调高下与礼貌态度的关系》。

点击题目查阅往期内容

R 语言建设和可视化混合效应模型 mixed effect model
R 语言 线性混合效应模型实战案例
R 语言用潜类别混合效应模型 (Latent Class Mixed Model ,LCMM) 剖析老年痴呆年龄数据
R 语言贝叶斯狭义线性混合（多层次 / 程度 / 嵌套）模型 GLMM、逻辑回归剖析教育留级影响因素数据 R 语言预计多元标记的潜过程混合效应模型（lcmm）剖析心理测试的认知过程
R 语言因子实验设计 nlme 拟合非线性混合模型剖析有机农业施氮程度
R 语言非线性混合效应 NLME 模型 (固定效应 & 随机效应) 反抗哮喘药物茶碱动力学钻研
R 语言用线性混合效应（多程度 / 档次 / 嵌套）模型剖析腔调高下与礼貌态度的关系
R 语言 LME4 混合效应模型钻研老师的受欢迎水平 R 语言 nlme、nlmer、lme4 用（非）线性混合模型 non-linear mixed model 剖析藻类数据实例
R 语言混合线性模型、多层次模型、回归模型剖析学生均匀问题 GPA 和可视化
R 语言线性混合效应模型（固定效应 & 随机效应）和交互可视化 3 案例
R 语言用 lme4 多层次（混合效应）狭义线性模型（GLM），逻辑回归剖析教育留级考察数据 R 语言 线性混合效应模型实战案例
R 语言混合效应逻辑回归（mixed effects logistic）模型剖析肺癌数据
R 语言如何用潜类别混合效应模型（LCMM）剖析抑郁症状
R 语言基于 copula 的贝叶斯分层混合模型的诊断准确性钻研
R 语言建设和可视化混合效应模型 mixed effect model
R 语言 LME4 混合效应模型钻研老师的受欢迎水平
R 语言 线性混合效应模型实战案例
R 语言用 Rshiny 摸索 lme4 狭义线性混合模型（GLMM）和线性混合模型（LMM）
R 语言基于 copula 的贝叶斯分层混合模型的诊断准确性钻研
R 语言如何解决线性混合模型中畸形拟合 (Singular fit) 的问题
基于 R 语言的 lmer 混合线性回归模型
R 语言用 WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型
R 语言分层线性模型案例
R 语言用 WinBUGS 软件对学术能力测验（SAT）建设分层模型
应用 SAS，Stata，HLM，R，SPSS 和 Mplus 的分层线性模型 HLM
R 语言用 WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型
SPSS 中的多层（等级）线性模型 Multilevel linear models 钻研整容手术数据
用 SPSS 预计 HLM 多层（档次）线性模型模型