乐趣区

关于算法:R语言线性混合效应模型固定效应随机效应和交互可视化3案例

jiezi

原文链接：http://tecdat.cn/?p=23050

原文出处：拓端数据部落公众号

在本文中，咱们将用 R 语言对数据进行线性混合效应模型的拟合，而后可视化你的后果。

线性混合效应模型是在有随机效应时应用的，随机效应产生在对随机抽样的单位进行屡次测量时。来自同一天然组的测量后果自身并不是独立的随机样本。因而，这些单位或群体被假设为从一个群体的 “ 人口 “ 中随机抽取的。示例状况包含

  • 当你划分并对各局部进行独自试验时（随机组）。
  • 当你的抽样设计是嵌套的，如横断面内的四分仪；林地内的横断面；地区内的林地（横断面、林地和地区都是随机组）。
  • 当你对相干个体进行测量时（家庭是随机组）。
  • 当你反复测量受试者时（受试者是随机组）。

混合效应的线性模型在 R 命令 lme4 和 lmerTest 包中实现。另一个抉择是应用 nmle 包中的 lme 办法。lme4 中用于计算近似自由度的办法比 nmle 包中的办法更精确一些，特地是在样本量不大的时候。

测量斑块长度

这第一个数据集是从 Griffith 和 Sheldon（2001 年，《动物行为学》61：987-993）的一篇论文中提取的，他们在两年内对瑞典哥特兰岛上的 30 只雄性领头鶲的红色额斑进行了测量。该斑块在吸引配偶方面很重要，但其大小每年都有变动。咱们在这里的指标是预计斑块长度（毫米）。

读取和检查数据

  1. 从文件中读取数据。
  2. 查看数据的前几行，看是否正确读取。
  3. 创立一个显示两年钻研中每只飞鸟的测量对图。能够尝试制作点阵图。是否有证据表明不同年份之间存在着测量变异性？

构建线性混合效应模型

  1. 对数据进行线性混合效应模型，将单个鸟类视为随机组。注：对每只鸟的两次测量是在钻研的间断年份进行的。为了简略起见，在模型中不包含年份。在 R 中把它转换成一个字符或因子，这样它就不会被当作一个数字变量。依照上面步骤（2）和（3）所述，用这个模型从新计算可重复性。重复性的解释如何扭转？
  2. 从保留的 lmer 对象中提取参数估计值（系数）。查看随机效应的输入。随机变异的两个起源是什么？固定效应指的是什么？
  3. 在输入中，查看随机效应的标准差。应该有两个标准差：一个是 ”（截距）”，一个是 “ 残差 ”。这是因为混合效应模型有两个随机变异的起源：鸟类外部反复测量的差别，以及鸟类之间额斑长度的实在差别。这两个起源中的哪一个对应于 ”（截距）”，哪一个对应于 “ 残差 ”？
  4. 同时查看固定效应后果的输入。模型公式中惟一的固定效应是所有长度测量的平均值。它被称为 ”（截距）”，但不要与随机效应的截距相混同。固定效应输入给了你平均值的估计值和该估计值的标准误差。留神固定效应输入是如何提供均值估计值的，而随机效应输入则提供方差（或标准差）的估计值。
  5. 从拟合模型中提取方差重量，预计各年斑块长度的可重复性 *。
  6. 解释上一步中取得的重复性测量后果。如果你失去的重复性小于 1.0，那么个体内测量后果之间的变动起源是什么。仅是测量误差吗？
  7. 产生一个残差与拟合值的图。留神到有什么问题？仿佛有一个轻微的正向趋势。这不是一个谬误，而是最佳线性无偏预测器（BLUPs）” 膨胀 “ 的后果。

剖析步骤

读取并检查数据。

head(fly)

# 点阵图
chart(patch ~ bird)

# 但显示成对数据的更好办法是用成对的交互图来显示
plot(res=patch, x = year)

# 优化版本
plot(y = patch, x = factor(year), theme_classic)

拟合一个线性混合效应模型。summary()的输入将显示两个随机变异的起源：单个鸟类之间的变异（鸟类截距），以及对同一鸟类进行的反复测量之间的变异（残差）。每个起源都有一个预计的方差和标准差。固定效应只是所有鸟类的平均值 – 另一个 “ 截距 ”。

# 1. 混合效应模型
# 2. 参数估计
summary(z)

# 5. 方差重量
VarCorr(z)

# 可重复性
1.11504^2/(1.11504^2 + 0.59833^2)
## \[1\] 0.7764342
# 7. 残差与拟合值的关系图
plot(z)

金鱼视觉

Cronly-Dillon 和 Muntz(1965; J. Exp. Biol 42: 481-493)用视静止反馈来测量金鱼的色觉。在这里，咱们将对数据进行拟合，包含测试的全副波长。5 条鱼中的每一条都以随机的程序在所有的波长下被测试。敏感度的值大表明鱼能够检测到低的光强度。视静止反馈的一个重要特点是，鱼不习惯，在一个波长下的视觉敏感度的测量不太可能对起初在另一个波长下的测量产生影响。

读取和检查数据

  1. 读取文件中的数据，并查看前几行以确保读取正确。
  2. 应用交互图来比拟不同光波长试验下的个体鱼的反馈。
  3. 应用什么类型的实验设计？* 这将决定在拟合数据时应用的线性混合模型。

构建线性混合效应模型

  1. 对数据拟合一个线性混合效应模型。能够用 lmer()来实现。发现“畸形拟合”，“boundary (singular) fit: see ?isSingular
    ”
  2. 绘制拟合（预测）值 **。每条鱼的预测值和察看值之间的差别代表残差。
  3. 你在（1）中做了什么假如？创立一个残差与拟合值的图，以查看这些假如之一。
  4. 从保留的 lmer 对象中提取参数估计值。查看固定效应的后果。给出的系数与应用 lm 剖析的分类变量的解释雷同。
  5. 查看随机效应的输入。咱们的混合效应模型中再次出现了两个随机误差的起源。它们是什么？其中哪个对应于输入中的 ”（截距）”，哪个对应于 “ 残差 ”？留神，在这个数据集中，其中一个变动源的预计标准差十分小。这就是畸形拟合信息背地的起因。鱼类之间的方差不太可能真的为零，然而这个数据集十分小，因为抽样误差，可能会呈现低方差预计。
  6. 生成基于模型的每个波长的均匀敏感度的预计。
  7. 各个波长之间的差别是否显著？生成 lmer 对象的方差分析表。这里测试的是什么效应，随机效应还是固定效应？解释方差分析后果。

* 这是一个 “ 按试验对象 “ 的反复测量设计，因为每条鱼在每个试验下被测量一次。它实质上与随机齐全区块设计雷同（把每条鱼看作是 “ 区块 ”）。

* 可视化是首选，因为数据和拟合值都被绘制进去。请留神鱼与鱼之间的预测值是如许的类似。这表明在这项钻研中，个体鱼之间的预计差别十分小。

*一般来说，在方差分析表中只测试固定效应。应用测试随机效应中没有方差的无效假设是可能的。

剖析步骤

读取并检查数据。

x <- read.csv("fish.csv", 
        stringsAsFactors = FALSE)
head(x)

拟合一个线性混合效应模型。

该模型假如所有拟合值的残差为正态分布，方差相等。该办法还假如个体鱼之间的随机截距为正态分布。该办法还假如组（鱼）的随机抽样，对同一鱼的测量之间没有影响。

# # 1. 拟合混合效应模型。
## boundary (singular) fit: see ?isSingular

# 2. 这就为每条鱼别离绘制了拟合值。vis(z)

 

# 3. 测试假如
plot(z)

# 4. 提取参数估计值
summary(z)

# 6.  基于模型的均匀敏感度预计 
means(z)

# 7. ANOVA 方差分析

蓍草酚类物质的浓度

我的项目实验性地考察了国家公园的南方森林生态系统中施肥和食草的影响（Krebs, C.J., Boutin, S. & Boonstra, R., eds (2001a) Ecosystem dynamics of the Boreal Forest.Kluane 我的项目. 牛津大学出版社，纽约）)，目前的数据来自于一项对于动物资源和食草动物对底层动物物种防御性化学的影响的钻研。

16 个 5 ×5 米的小区中的每一个都被随机调配到四个试验之一。1）用栅栏围起来排除食草动物；2）用 N -P- K 肥料施肥；3）用栅栏和施肥；4）未试验的对照。而后，16 块地中的每一块被分成两块。每块地的一侧（随机抉择）在 20 年的钻研中继续承受试验。每块地的另一半在头十年承受试验，之后让它复原到未试验的状态。这里要剖析的数据记录了欧蓍草（Achillea millefolium）中酚类物质的浓度（对动物进攻化合物的粗略测量），欧蓍草是地块中常见的草本植物。测量单位是每克干重毫克丹宁酸当量。

可视化数据

  1. 从文件中读取数据。
  2. 查看前几行的数据。试验是作为一个有四个档次的繁多变量给出的（而不是作为两个变量，围墙和肥料，用 2 ×2 因子设计的模型）。持续时间示意半块土地是否承受了整整 20 年的试验，或者是否在 10 年后进行试验。变量 “ch “ 是蓍草中酚类物质的浓度。
  3. 画一张图来阐明不同试验和持续时间类别中蓍草中的酚类物质的浓度。在每个试验和持续时间程度的组合中没有很多数据点，所以按组画条形图可能比按组画箱形图更好。
  4. 增加线段来连接成对的点。

拟合一个线性混合效应模型

  1. 应用的是什么类型的实验设计？* 这将决定对数据的线性混合模型的拟合。
  2. 在没有试验和持续时间之间的交互作用的状况下，对数据进行线性混合模型拟合。应用酚类物质的对数作为因变量，因为对数转换改善了数据与线性模型假如的拟合。
  3. 可视化模型对数据的拟合。按持续时间（如果 xvar 是试验）或试验（如果 xvar 是持续时间）离开面板。visreg()不会保留配对，但会容许你查看残差。
  4. 当初反复模型拟合，但这次包含试验和持续时间之间的相互作用。将模型与数据的拟合状况可视化。两个模型拟合之间最显著的区别是什么，一个有交互作用，另一个没有？形容包含交互项的模型 “ 容许 “ 什么，而没有交互项的模型则不容许。判断，哪个模型最适宜数据？
  5. 应用诊断图查看包含交互项的模型的线性混合模型的一个要害假如。
  6. 应用拟合模型对象预计线性模型的参数（包含交互作用）。请留神，当初固定效应表中有许多系数。
  7. 在上一步的输入中，你会看到 “ 随机效应 “ 标签下的 “Std.Dev “ 的两个数量。解释一下这些数量指的是什么。
  8. 来预计所有固定效应组合的模型拟合平均值。
  9. 生成固定效应的方差分析表。哪些项在统计学上是显著的？
  10. 默认状况下，lmerTest 将应用 Type 3 的平方和来测试模型项，而不是按程序（Type 1）。用类型 1 来反复方差分析表。后果有什么不同吗？**

* 试验采纳了分块设计，即整个块被随机调配到不同的试验，而后将第二种试验（持续时间）的不同程度调配到块的一半。

* 应该没有差异，因为设计是齐全均衡的。

剖析步骤

浏览并检查数据。

一个好的策略是对试验类别进行排序，把对照组放在后面。这将使线性模型的输入更加有用。

# 1. 读取数据
# 2. 查看
head(x)

# 3. 分组带状图
# 首先，重新排列试验类别
factor(treat,levels=c("cont","exc","fer","bo"))
plot(data = x, y = log(phe), x = treat, fill = dura, color = dura)

# 4. 在多个面板上别离绘制成对的数据
plot(data = x,y = log(ach, x = dur, fill = dur, col = dur)

拟合一个线性混合效应模型。固定效应是 “ 试验 “ 和 “ 持续时间 ”，而 “ 块 ” 是随机效应。拟合交互作用时，试验程度之间的差别大小在持续时间程度之间会有所不同。

因为随机效应也存在（块），系数表将显示两个随机变动起源的方差预计。一个是拟合模型的残差的方差。第二个是（随机）块截距之间的方差。

# 2. 拟合混合效应模型 - 无交互作用
# 3. 可视化
vis(z)

# 4. 包含交互项和再次视觉化
z.int <- lmer(log(phen.ach) ~ treatment * duration + (1|plot), data=x)

vis(z.int, overlay = TRUE)

# 5. 绘制图表以测验方差齐性（以及正态性）plot(z)

# 6. 系数
summary(z)

# 8. 模型拟合平均值
means(z, data = x)

# 9. 方差分析表
anova(z) # lmerTest 中默认为 3 类平方和。

# 10.  改为 1 类
anova(z, type = 1)

最受欢迎的见解

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

退出移动版