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

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

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

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

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

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

测量斑块长度

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

读取和检查数据

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

构建线性混合效应模型

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

剖析步骤

读取并检查数据。

head(fly)

# 点阵图chart(patch ~ bird)

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

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

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

点击题目查阅往期内容

R语言用lme4多层次（混合效应）狭义线性模型（GLM），逻辑回归剖析教育留级考察数据

左右滑动查看更多

01

02

03

04

# 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条鱼中的每一条都以随机的程序在所有的波长下被测试。敏感度的值大表明鱼能够检测到低的光强度。视静止反馈的一个重要特点是，鱼不习惯，在一个波长下的视觉敏感度的测量不太可能对起初在另一个波长下的测量产生影响。

读取和检查数据

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

构建线性混合效应模型

  • 对数据拟合一个线性混合效应模型。能够用lmer()来实现。发现“畸形拟合”，“boundary (singular) fit: see ?isSingular
    ”
  • 绘制拟合（预测）值**。每条鱼的预测值和察看值之间的差别代表残差。
  • 你在（1）中做了什么假如？创立一个残差与拟合值的图，以查看这些假如之一。
  • 从保留的lmer对象中提取参数估计值。查看固定效应的后果。给出的系数与应用lm剖析的分类变量的解释雷同。
  • 查看随机效应的输入。咱们的混合效应模型中再次出现了两个随机误差的起源。它们是什么？其中哪个对应于输入中的"（截距）"，哪个对应于 "残差"？留神，在这个数据集中，其中一个变动源的预计标准差十分小。这就是畸形拟合信息背地的起因。鱼类之间的方差不太可能真的为零，然而这个数据集十分小，因为抽样误差，可能会呈现低方差预计。
  • 生成基于模型的每个波长的均匀敏感度的预计。
  • 各个波长之间的差别是否显著？生成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个5x5米的小区中的每一个都被随机调配到四个试验之一。1）用栅栏围起来排除食草动物；2）用N-P-K肥料施肥；3）用栅栏和施肥；4）未试验的对照。而后，16块地中的每一块被分成两块。每块地的一侧（随机抉择）在20年的钻研中继续承受试验。每块地的另一半在头十年承受试验，之后让它复原到未试验的状态。这里要剖析的数据记录了欧蓍草（Achillea millefolium）中酚类物质的浓度（对动物进攻化合物的粗略测量），欧蓍草是地块中常见的草本植物。测量单位是每克干重毫克丹宁酸当量。

可视化数据

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

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

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

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

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

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

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

# 4. 包含交互项和再次视觉化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)

本文摘选 《 R语言线性混合效应模型（固定效应&随机效应）和交互可视化3案例 》 ，点击“浏览原文”获取全文残缺材料。

点击题目查阅往期内容

非线性混合效应 NLME模型反抗哮喘药物茶碱动力学钻研
生态学模仿对狭义线性混合模型GLMM进行功率（效用、效力、效劳）剖析power analysis环境监测数据
无限混合模型聚类FMM、狭义线性回归模型GLM混合利用剖析威士忌市场和钻研专利申请数据
如何用潜类别混合效应模型(Latent Class Mixed Model ,LCMM)剖析老年痴呆年龄数据
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多层（档次）线性模型模型