R语言广义线性混合模型GLMMs在生态学中应用可视化2实例合集|附数据代码1：https://developer.aliyun.com/article/1501226

3. 参数自助法似然比检验：对新的固定效应模型进行了参数自助法似然比检验。
   

# 模拟过程  
  for(i in 1:n.sim){  
    # 从零模型中模拟新的观测值  
    dattemp$ysim <- unlist(simulate(mod3_lmer))  
      
    
  # 计算p值  
  (sum(lrt.sim >= lrt.obs) + 1) / (n.sim + 1)

解释：

• 绘制直方图：通过直方图，我们可以直观地看到模拟的似然比检验统计量的分布情况，并标出观察到的统计量，从而评估其显著性。
  
• 固定效应部分的模型选择：这部分代码通过拟合不同的混合效应模型来比较固定效应部分的影响。lme函数用于拟合线性混合效应模型，而lmer函数用于拟合线性混合效应模型，但使用的是lme4包。两种模型都考虑了随机效应，但使用了不同的估计方法（lme使用最大似然法，lmer使用REML或ML）。
  
• 参数自助法似然比检验：与之前的自助法类似，但这次是针对固定效应部分的模型进行比较。代码从mod3_lmer（只包含NAP作为固定效应的模型）中模拟新的观测值，然后拟合零模型和替代模型，并计算似然比检验统计量。最后，基于模拟的统计量计算p值，以评估固定效应Exposure是否显著。
  

注意：在代码中，simulate函数用于从模型生成模拟数据，而anova函数用于比较模型的差异。此外，unlist函数用于将列表转换为向量，因为simulate函数返回的可能是一个列表

# 使用蓝色绘制直方图，直方图的x轴范围为0到lrt.sim和lrt.obs中的最大值，并设置x轴和y轴的标签大小  
    
  # 在直方图上添加一条垂直于x轴的橙色线，表示lrt.obs的值，线宽为3  
  abline(v=lrt.obs, col="orange", lwd=3)  
    
  # 下一步将是先删除NAP，然后看似然比检验是否显著，以及是否总是先删除Exposure导致更高的LRT统计量  
  # 其他方法，如AIC...  
  # GLMM的R平方计算，参见Nakagawa 2013 MEE的补充材料  
  # 计算固定效应方差  
  # VarCorr()函数用于提取方差分量  
  # attr(VarCorr(lmer.model),’sc’)^2提取残差方差，VarCorr()$plot提取plot效应的方差  
  # 计算条件R平方  
  #conditionnal R-square 的计算公式

解释：

• hist函数用于绘制直方图，展示lrt.sim（可能是似然比检验的统计量模拟值）的分布情况。直方图的x轴范围设置为从0到lrt.sim和lrt.obs（观察到的似然比检验统计量）中的最大值。同时设置了直方图的颜色、x轴和y轴的标签以及标签的大小。
  
• abline函数在直方图上添加了一条垂直于x轴的线，线的位置为lrt.obs的值，线的颜色为橙色，线宽为3。这通常用于在直方图上标识某个特定的观察值或阈值。
  
• 注释部分说明了下一步的分析计划，即先删除NAP变量，然后检查似然比检验是否显著，以及比较先删除NAP还是先删除Exposure变量对LRT统计量的影响。同时提到了其他分析方法，如AIC（赤池信息准则）。
  
• 接下来的代码计算了线性混合效应模型mod1_lmer的条件R平方。这包括计算固定效应的方差（VarF），提取模型的方差分量（VarCorr），以及计算条件R平方的值。条件R平方是解释模型中固定效应和随机效应共同解释的变异比例。
  

R复制代码
  # 从模型中推断  
  # lme 和 glmer 可以获取 p 值，但 lmer 不行  
    
  # 使用 glmer 拟合模型  
  summary(mod1_glmer)  
    
  # 使用 arm 包中的 sim 函数  
  n.sim <- 1000  
  head(simu@fixef)  
  # 95% 可信区间  
  # 绘制 NAP 对丰富度的影响  
  predmat <- matrix(ncol = nsim, nrow = nrow(newdat))  
    
  lines(newdat$NAP, newdat$lower, col = "red", lty = 2, lwd = 1.5)  
    
  # 比较不同项的系数，标准化变量  
  data$stdNAP <- scale(data$NAP)  
  # 95% 可信区间  
  
    
  # 绘图  
  plot(1:3, coeff[, 2], ylim = c(-0.8, 2), xaxt = "n", xlab = "", ylab = "Estimated values")

解释

这段代码的主要目的是从广义线性混合效应模型（GLMM）中推断结果，并绘制相关效应的可视化图形。

1. 模型推断：
   

• summary(mod1) 和 summary(mod1_lmer) 用于展示模型 mod1 和 mod1_lmer 的统计摘要，包括系数估计、标准误差、z值或t值以及对应的p值。
  

2. 拟合GLMM模型：
   

• glmer 函数用于拟合广义线性混合效应模型，这里以物种丰富度（Richness）作为响应变量，NAP和Exposure作为固定效应，同时考虑到Beach作为随机效应。
  

3. 模拟后验分布
   

2.生态学模拟对广义线性混合模型GLMM进行功率（功效、效能、效力）分析power analysis环境监测数据|附代码数据

假设检验的功效定义为假设原假设为假，检验拒绝原假设的概率。换句话说，如果一个效应是真实的，那么分析判断该效应具有统计显着性的概率是多少？

概括

• r 语言允许用户计算 lme 4 包中广义线性混合模型的功效。功率计算基于蒙特卡罗模拟。
  
• 它包括用于 (i) 对给定模型和设计进行功效分析的工具；(ii) 计算功效曲线以评估功效和样本量之间的权衡。
  
• 本文提供了一个教程，使用具有混合效果的计数数据的简单示例（具有代表环境监测数据的结构）。
  

介绍

如果一项研究的功效不足，资源可能被浪费，真正的效果可能被遗漏。另一方面，一项大型研究的花费可能过大，因此其费用也会超过必要的范围。因此，在收集数据之前进行功效分析是一个很好的做法，以确保样本具有适当的规模来回答正在考虑的任何研究问题。

Kaizong Ye

分析师

广义线性混合模型 (GLMM) 在生态学中很重要，它允许分析计数和比例以及连续数据，并控制空间非独立性.

蒙特卡罗模拟是一种灵活且准确的方法，适用于现实的生态研究设计。在某些情况下，我们可以使用解析公式来计算功效，但这些通常是近似值或需要特殊形式的设计。仿真是一种适用于各种模型和方法的单一方法。即使公式可用于特定模型和设计，定位和应用适当的公式也可能非常困难，因此首选仿真。

对于对 r 不够熟悉的研究人员，设置模拟实验可能太复杂了。在本文中，我们介绍了一个工具来自动化这个过程。

r 包

有一系列的 r 包目前可用于混合模型的功效分析 。然而，没有一个可以同时处理非正态因变量和广泛的固定和随机效应规范。

图1

r 旨在与任何可以与 lme 4 中的 lmer 或 glmer 配合的线性混合模型 (LMM) 或 GLMM 一起使用。这允许具有不同固定和随机效应规范的各种模型。还支持在 r 中使用 lm 和 glm 的线性模型和广义线性模型，以允许没有随机效应的模型。

r 中的功效分析从适合 lme 4 的模型开始。

在 r 中，通过重复以下三个步骤来计算功效：(i) 使用提供的模型模拟因变量的新值；(ii) 将模型重新拟合为模拟因变量；(iii) 对模拟拟合应用统计检验。在此设置中，已知存在测试效果，因此每个阳性测试都是真正的阳性，每个阴性测试都是 II 类错误。可以根据步骤 3 的成功和失败次数计算测试的功效。

教程

本教程使用包含的数据集。该数据集代表环境监测数据，在连续固定效应变量_x _（例如研究年份）的10 个水平上测量三个组 _g _（例如研究地点）的因变量 _z _（例如鸟类丰度 ）。还有一个连续因变量 _y _，在本教程中没有使用。

拟合模型

我们首先将 lme 4 中的一个非常简单的泊松混合效应模型拟合到数据集。在这种情况下，我们有一个随机截距模型，其中每个组 ( g ) 都有自己的截距，但这些组共享一个共同的趋势。

glm
summary

本教程重点介绍关于_x _趋势的推断 。在这种情况下，_x _的估计效应大小为 -0.11，使用默认_z_检验在 0.01 水平上显着 。

请注意，我们特意使用了一个非常简单的模型来使本文易于理解。例如，适当的分析会包含更多的组，并会考虑过度分散等问题。

简单的功率分析

假设我们想重复这项研究。如果效果是真实的，我们是否有足够的功效来期待积极的结果？

指定效应量

在开始功效分析之前，重要的是要考虑您感兴趣的效果大小类型。功效通常随效果大小而增加，较大的效果更容易检测。回顾性“观察功效”计算，其中目标效应大小来自数据，给出误导性结果.

对于此示例，我们将考虑检测 -0.05 斜率的功效。可以使用 lme 4 函数拟合 glmer 模型中的固定效应。然后可以更改固定效应的大小。变量_x _的固定效应的大小 可以从 -0.11 更改为 -0.05，如下所示：

fixe<‐ ‐0.05

在本教程中，我们只更改变量_x _的固定斜率 。但是，我们也可以更改随机效应参数或残差方差（适用于合适的模型）。

运行功效分析

一旦指定了模型和效应大小，在 r 中进行功效分析就非常容易了。由于这些计算基于蒙特卡罗模拟，因此您的结果可能略有不同。如果你想得到和本文一样的结果，你可以使用 set.seed(123)。

power

于此特定设置，拒绝_x _中零趋势的零假设的 能力约为 33%。这几乎总是被认为是不够的；传统上，80% 的功率被认为是足够的.

在实践中，  z_检验可能不适合这样一个小例子。参数引导测试 可能是最终分析的首选。但是，更快的 _z -test 更适合学习使用该包以及在功效分析期间进行初始探索性工作。

增加样本量

在第一个示例中，估计功率很低。小型试点研究通常没有足够的功效来检测微小的影响，但更大的研究可能会。

试点研究对_x 的 _10 个值进行了观察， 例如代表研究第 1 年到第 10 年。在此步骤中，我们将计算将其增加到 20 年的影响。

modl2 <‐ extend
power(modl2)

沿参数指定要扩展的变量，n 指定要替换它的级别。扩展模型 2 现在将具有 从 1 到 20 的_x _值，与以前一样分为三组，总共 60 行（与模型 1 中的 30 行相比）。

通过观察_x 的 _20 个值 ，我们将有足够的能力来检测大小为 -0.05 的效应。

各种样本量的功效分析

当数据收集成本高昂时，用户可能只想收集达到一定统计能力所需的数据量。功效曲线 函数可用于探索样本大小和功效之间的权衡。

确定所需的最小样本量

在前面的示例中，当对变量_x 的_20 个值进行观察时，我们发现了非常高的 _功效 _。我们能否减少这个数字，同时保持我们的功效高于通常的 80% 阈值？

poerCure
print
plot

请注意，我们已将此结果保存到变量 pc2 以匹配模型 2 中的编号。由于模型 1 没有足够的功率，我们没有通过 powerCurve 运行它。绘制的输出如图所示。 我们可以看到，检测_x _趋势的 能力随着采样大小的增加而增加。这里的结果基于将模型拟合到 10 个不同的自动选择的子集。最小的子集仅使用前 3 年（即 9 个观测值），最大的子集使用所有 20 个假设研究年份（即 60 行数据）。该分析表明，该研究必须运行 16 年才能有≥80% 的功效来检测指定大小的影响。

图2

检测大小为 -0.05 的固定效应的功效 (±95% CI)，使用 powerCurve 函数在一系列样本大小上计算。变量_x 的不同值的数量 从 3 ( _n  = 9) 到 20 ( n  = 60) 不等。

改变组的数量和大小

增加观察到的_x _值的数量可能不可行 。例如，如果 _x _是研究年份，我们可能不愿意等待更长时间的结果。在这种情况下，增加研究地点的数量或每个地点的测量数量可能是更好的选择。这两项分析从我们的原始模型 1 开始，该模型已有 10 年的研究时间。

添加更多组

我们可以像为_x _添加额外值一样 为_g _添加额外级别 。例如，如果变量 _g _代表我们的研究站点，我们可以将站点数量从 3 增加到 15。

extend(n=15)
plot(pc3)

与上一个示例的主要变化是我们将变量_g _传递 给了沿参数。该分析的输出如图 1 所示。要达到 80% 的功率，我们至少需要 11 个站点。

图 3

检测大小为 -0.05 的固定效应的功效 (±95% CI)，使用 powerCurve 在一系列样本大小上计算。因子_g 的级别数 从 3 ( _n  = 30) 到 15 ( n  = 150) 不等。

增加组内的大小

我们可以用内参数替换扩展和 powerCurve 的沿参数以增加组内的样本大小。每个组在_x _和 _g 的 _每个水平上只有一个观察值 。我们可以将其扩展到每个站点每年 5 次观测，如下所示：

extend( n=5)
plot(p4)

请注意 powerCurve 的breaks 参数。为_x _和 _g 的 _每个组合提供一到五个观察结果 。图表明每年每个站点 4 次观测会给我们 80% 的效力。

图 4

检测大小为 -0.05 的固定效应的功效 (±95% CI)，使用 powerCurve 函数在一系列样本大小上计算。_x 和 _g 的 _每个组合的观察数 从 1 ( _n  = 30) 到 5 ( n  = 150) 不等。