R语言生态学模拟对广义线性混合模型GLMM进行功率（功效、效能、效力）分析power analysis环境监测数据

最近我们被客户要求撰写关于生态学的研究报告，包括一些图形和统计输出。

相关视频：线性混合效应模型(LMM,Linear Mixed Models)和R语言实现

线性混合效应模型(LMM,Linear Mixed Models)和R语言实现案例

时长12:13

相关视频：马尔可夫链蒙特卡罗方法MCMC原理与R语言实现

马尔可夫链蒙特卡罗方法MCMC原理与R语言实现

，时长08:47

概括

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

介绍

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

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

广义线性混合模型 (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) 不等。

---