如何在二项式 GLMM（lme4）中确定正确的混合效应结构？

Question

有人可以帮我在 lme4 中的二项式 GLMM 中确定正确的随机变量结构吗？

我将首先尽力解释我的数据。我有被吃掉（1）或没有被吃掉（0）的幼苗的二项式数据，以及植被覆盖的数据。我试图弄清楚植被覆盖率与树木被吃掉的概率之间是否存在关系，因为其他植被是一种食物来源，可以将食草动物吸引到某个森林斑块。

数据收集在分散在国家公园内的约 90 个地块中，历时 9 年。有些测量了所有年份，有些只测量了几年（销毁/新添加的地块）。原始数据集分为 2 个（落叶树与针叶树），均包含约 55.000 个条目。每个地块每次测量大约 100 棵树苗，因此两个单独的数据集可能包含每个地块大约 50 棵树（尽管情况并非总是如此，因为 decid:conif 比率并不总是相等）。每个图由 4 个子图组成。 我知道由于绘图位置可能存在空间自相关，但我们不会对此进行纠正。

每年在同一时期对植被进行调查。植被覆盖在地块级别估计，单个树木（二元）在子地块级别测量。 所有树木都经过测量，因此随着森林自然再生，每个子样地的响应量在子样地和年份之间会有所不同。

很遗憾，我无法分享我的原始数据，但我尝试创建一个包含要点的示例：

#set seed for whole procedure
addTaskCallback(function(...) {set.seed(453);TRUE})
# Generate vector containing individual vegetation covers (in %) 
cover1vec <- c(sample(0:100,10, replace = TRUE)) #the ',number' is amount of covers generated

# Create dataset
DT <- data.frame(
                 eaten = sample(c(0,1), 80, replace = TRUE),
                 plot = as.factor(rep(c(1:5), each = 16)), 
                 subplot = as.factor(rep(c(1:4), each = 2)),
                 year = as.factor(rep(c(2012,2013), each = 8)),
                 cover1 = rep(cover1vec, each = 8)
)

这将生成这个数据集：

>DT
   eaten plot subplot year cover1
1      0    1       1 2012      4
2      0    1       1 2012      4
3      1    1       2 2012      4
4      1    1       2 2012      4
5      0    1       3 2012      4
6      1    1       3 2012      4
7      0    1       4 2012      4
8      1    1       4 2012      4
9      1    1       1 2013     77
10     0    1       1 2013     77
11     0    1       2 2013     77
12     1    1       2 2013     77
13     1    1       3 2013     77
14     0    1       3 2013     77
15     1    1       4 2013     77
16     0    1       4 2013     77
17     0    2       1 2012     46
18     0    2       1 2012     46
19     0    2       2 2012     46
20     1    2       2 2012     46
....etc....
80     0    5       4 2013     82    

注意 1： 再次澄清一下，在这个例子中，每个 subplot:year 组合的响应数量是相同的，从而使数据平衡，而原始数据集中的情况并非如此。 注意 2： 这个例子不能在 GLMM 中运行，因为我收到了一个奇点警告并且我所有的随机效应测量值都为零。显然我的示例不适合实际使用（因为使用 sample() 导致 0 和 1 的数量太均匀而无法产生足够大的影响？）。

从示例中可以看出，每个 plot:year 组合的封面数据都是相同的。 图是多年测量的（示例中只有 2012 年和 2013 年），因此存在重复测量。 此外，考虑到我们有例如 年份效应，干燥/潮湿年份。

首先我想到了以下模型结构：

library(lme4)
mod1 <- glmer(eaten ~ cover1 + (1 | year) + (1 | plot), data = DT, family = binomial)
summary(mod1)

其中 (1 | year) 应校正年份之间的差异，而 (1 | plot) 应校正重复测量。

但后来我开始思考：与 2013 年在地块 1 中测量的所有树木相比，与 2013 年地块 1 中的（部分相同的）树木相比，它们之间的相似度会更高。 所以，我怀疑这种随机模型结构能否纠正这种图内时间效应。

所以我最好的猜测是添加另一个随机变量，其中考虑了这种“交互”。 我知道有两种方法可以实现这一点：

方法一。 添加随机变量“+ (1 | year:plot)”

方法二。 添加随机变量“+ (1 | year/plot)”

从其他人告诉我的情况来看，我仍然不知道两者之间的区别。 我看到方法 2 相比 方法 1 添加了一个额外的随机变量（year.1），但我不知道如何解释这个额外的随机变量。

作为示例，我使用 方法 2 添加了随机效应摘要（由于我的示例数据存在奇异性问题，因此为零）：

Random effects:
 Groups    Name        Variance Std.Dev.
 plot.year (Intercept) 0        0       
 plot      (Intercept) 0        0       
 year      (Intercept) 0        0       
 year.1    (Intercept) 0        0       
Number of obs: 80, groups:  plot:year, 10; plot, 5; year, 2

谁能解释一下方法 1 和 方法 2 之间的实际区别？ 我试图了解正在发生的事情，但无法掌握。

我已经尝试从一位同事那里获得建议，他提到使用 cbind(success, failure) per plot:year 组合可能更合适。 通过这个站点，我发现当 Ntrails > 1 时在二项式模型中使用 cbind，我认为考虑到我们的抽样程序，情况确实如此。

我想知道，如果 cbind 已经用于 plot:year 组合，是否需要添加 plot:year 随机变量？ 使用 cbind 时，示例数据如下所示：

>DT3
   plot year cover1 Eaten_suc Eaten_fail
8     1 2012      4         4          4
16    1 2013     77         4          4
24    2 2012     46         2          6
32    2 2013     26         6          2
40    3 2012     91         2          6
48    3 2013     40         3          5
56    4 2012     61         5          3
64    4 2013     19         2          6
72    5 2012     19         5          3
80    5 2013     82         2          6

正确的随机模型结构是什么？为什么？ 我在想：

可能性A

mod4 <- glmer(cbind(Eaten_suc, Eaten_fail) ~ cover1 + (1 | year) + (1 | plot), 
              data = DT3, family = binomial)

可能性 B

mod5 <- glmer(cbind(Eaten_suc, Eaten_fail) ~ cover1 + (1 | year) + (1 | plot) + (1 | year:plot), 
              data = DT3, family = binomial)

但是 cbind(success, failure) 不是已经正确的年份：情节依赖？

可能性C

mod6 <- glmer(cbind(Eaten_suc, Eaten_fail) ~ cover1 + (1 | year) + (1 | plot) + (1 | year/plot), 
              data = DT3, family = binomial)

因为我还不明白 year:plot 和 year/plot 之间的区别

因此：使用 cbind 方法确实比原始二进制数据更合适吗？为了防止伪复制和其他依赖关系，需要什么随机模型结构？

提前感谢您的时间和投入！

编辑 7/12/20：我添加了一些关于原始数据的额外信息

Answer 1

您在问题中提出了很多问题。我将尝试涵盖所有内容，但我建议阅读来自lme4 和glmmFAQ 页面的文档和小插图以获取更多信息。此外，我强烈建议您在 google 学者上搜索这些主题，因为它们的内容相当丰富。

我会从一个简单的地方开始

注 2（为什么我的模型是单数的？）

您的模型是高度单一的，因为您模拟数据的方式并不表明数据本身之间存在任何依赖关系。如果您想模拟二项式模型，您可以使用g(eta) = X %*% beta 来模拟您的线性预测器，从而模拟成功的概率。然后可以使用此概率来模拟您的二元结果。因此，这将是一个两步过程，首先使用一些已知的X 或随机模拟的X，给定我们选择的一些先验分布。在第二步中，我们将使用rbinom 来模拟二进制结果，同时使其依赖于我们的预测器X。

在您的示例中，您正在模拟独立的X 和y，其中概率也独立于X。因此，当我们查看结果y 时，对于某个常数c，所有子组的成功概率等于p=c。

谁能解释一下方法 1 和方法 2 之间的实际区别？ （(1| year:plot) 与 (1|year/plot)）

这在第 7 页表格中的包小插图fitting linear mixed effects models with lme4 中进行了说明。

1. (1|year/plot) 表示我们有 2 个混合拦截效果，year 和 plot 和 plot 嵌套在 year 中。
2. (1|year:plot) 表示单个混合拦截效果，plot 嵌套在 year 中。例如。我们不包括year 的主要效果。这有点类似于没有截距的模型（虽然不那么激烈，并且解释不会被破坏）。

第一个比第二个更常见，但我们可以将第一个写成第二个(1|year) + (1|year:plot)的函数。

因此：使用 cbind 方法确实比原始二进制数据更合适吗？

公式中的cbind 用于二项式数据（或多变量分析），而对于二进制数据，我们使用原始向量或0/1 表示成功/失败，例如。聚合二进制数据（类似于我们使用 glm 的方式）。如果您对子图的随机/固定效应不感兴趣，您也许可以跨图聚合数据，然后这可能是有意义的。否则留在你身边0/1 指示成功或失败的结果向量。

正确的随机模型结构是什么？为什么？

这是一个很难给出明确答案的话题，并且仍在积极研究中。根据您的统计范式，意见差异很大。

方法一：经典方法

经典混合建模基于对您正在处理的数据的了解。一般来说，选择这些参数有几个“经验法则”。我在my answer here 中经历了一些。一般来说，如果您对系统效应“不感兴趣”并且可以将其视为某些人群的随机样本，那么它可能是随机效应。如果是人口，例如。如果重复该过程，样本不会改变，那么它可能不应该改变。

这种方法通常会为那些不熟悉混合效果模型的人提供“不错的”选择，但受到倾向于使用类似于我们在非混合模型中使用的方法的作者的强烈批评（例如，可视化以基于我们的选择和检验显着性）。

方法二：使用可视化

如果您能够将数据分成独立的子组并保持固定效应结构，则检查潜在随机效应的合理方法是跨这些子组估计边际模型（例如使用glm）并查看固定效应是否在这些观察值之间“正态分布”。函数lmList（在lme4 中）就是为这种特定方法设计的。在线性模型中，我们确实希望这些是正态分布的，因此我们可以得到一个指示，一个特定的分组“可能”是一个有效的随机效应结构。我相信在广义线性模型的情况下也是如此，但我缺乏参考资料。我知道 Ben Bolker 在我在论文中使用的他之前的一篇文章（下面的第一个参考资料）中提倡这种方法。然而，这只是对严格可分数据的有效方法，并且在因子水平未在所有组之间共享的情况下，实施并不稳健。

简而言之：如果您拥有正确的数据，这种方法简单、快速且看似高度可靠。

方法 3：基于 AIC 或 AICc（或 p 值检验或替代指标）拟合最大/最小模型和减小/扩展模型

最后是使用类似“逐步”程序的替代方法。有人提倡从最大模型和最小模型开始（我确信下面的参考文献中至少有一篇谈到了两者的问题，否则请查看 glmmFAQ），然后测试您的随机效应的有效性。就像经典回归一样，这在某种程度上是一把双刃剑。原因既非常容易理解，又非常复杂。

要使此方法成功，您必须执行交叉验证或样本外验证以避免选择偏差，就像标准模型一样，但与标准模型不同，抽样变得复杂，因为：

1. 固定效应取决于随机结构。
2. 您需要独立的训练和测试样本
3. 由于这取决于您的随机结构，并且这是逐步选择的，因此很难避免某些模型中的信息泄漏。
4. 此处避免出现问题的唯一确定方法是定义空间 您将根据最 限制性模型定义。

接下来，我们在选择评估指标时也遇到了问题。如果对随机效应感兴趣，则使用 AICc（条件模型的 AIC 估计）是有意义的，而对于固定效应，优化 AIC（边际模型的 AIC 估计）可能更有意义。我建议在 glmmFAQ 上检查对 AIC 和 AICc 的引用，并且要小心，因为在一组非常严格的混合模型（即“足够的独立样本超过随机效应”）之外，这些的大样本结果可能是不确定的。

这里的另一种方法是使用 p 值而不是过程的某些度量。但人们可能对随机效应的测试更加谨慎。即使使用贝叶斯方法或以难以置信的大量重采样进行引导，有时这些也不是很好。同样，我们需要“足够多的独立样本而不是随机效应”来确保准确性。

DHARMA 为可能更适合的混合效果提供了一些非常有趣的测试方法。当我在该领域工作时，作者仍在（似乎）开发一篇文章，记录他们选择的方法的有效性。即使不使用它进行初始选择，我也只能建议检查一下并决定是否相信他们的方法。这是迄今为止最简单的视觉测试方法，具有简单的解释（例如，解释图几乎不需要先验知识）。

因此，关于此方法的最后一点说明是：这确实是一种方法，但我个人不推荐。它需要非常小心或作者接受对模型假设的无知。

结论

混合效果参数选择是困难。我的经验告诉我，主要使用方法 1 和 2 的组合，而方法 3 似乎主要由较新的作者使用，这些往往忽略样本外错误（根据用于训练的数据测量模型指标） )，在拟合随机效应时忽略样本的独立性问题或限制自己仅使用此方法测试固定效应参数。然而，所有 3 个都具有一定的有效性。我自己往往属于第一类，我的决定基于我在该领域的“经验”、经验法则和我的数据限制。

您的具体问题。

鉴于您的具体问题，我假设(1|year/plot/subplot) 的混合效果结构将是正确的结构。如果您添加自回归（时间-空间）效果，year 可能会消失。这种结构的原因是，在地块的地理分析和分析中，经典的方法是包括每个地块的影响。如果每个图可以进一步索引到子图，那么很自然地认为“子图”嵌套在“图”中。假设您不对自回归效应建模，我会认为 time 是随机的，原因您已经陈述过。有些年份我们的天气会比其他年份更干燥、更热。由于测量的地块必须存在于给定年份，因此这些地块将嵌套在年份中。

这就是我所说的maximal 模型，根据您的数据量，它可能不可行。在这种情况下，我会尝试使用(1|time) + (1|plot/subplot)。如果两者都可行，我会比较这些模型，使用自举方法或近似 LRT 测试。

注意：(1|time/plot/subplot) 似乎不太可能导致“个体水平效应”。例如，数据中每行 1 个随机效应。由于我早已忘记（但一旦阅读）的原因，在二元混合模型中具有个体（也称为主题级别）效果是不合理的。在这种情况下，使用替代方法或测试在从随机效应中保留 subplot 时是否保留模型假设也可能有意义。

下面我添加了一些有用的参考资料，其中一些与问题直接相关。此外，请查看 Ben Bolker 的 glmmFAQ 网站等。

参考文献

1. Bolker, B. 等人。 （2009 年）。 “广义线性混合模型：生态学和进化的实用指南”。在：生态与进化趋势 24.3，p. 127–135。
2. Bolker, B. 等人。 （2011）。 “GLMMs in action：基因与环境的相互作用在拟南芥野生种群果实总产量中的作用”。在：修订版，第 1 部分 1，p。 127–135。
3. 渴望，C. og J. Roy (2017)。 “混合效果模型有时很糟糕”。在：arXiv 预印本 arXiv:1701.04858。网址：https://arxiv.org/abs/1701.04858（上次访问时间为 19.09.2019）。
4. 冯，辛迪等人。 （2017）。 “随机分位数残差：具有统一参考分布的综合模型诊断工具”。在：arXiv 预印本 arXiv:1708.08527。 （最后一次出现在 2019 年 9 月 19 日）。
5. Gelman, A. og Jennifer Hill (2007)。使用回归和多级/分层模型的数据分析。剑桥大学出版社。
6. Hartig, F. (2019)。 DHARMa：分层（多级/混合）回归模型的残差诊断。 R 包版本 0.2.4。网址：http://florianhartig.github.io/DHARMa/（最后一次出现时间为 19.09.2019）。
7. Lee, Y. og J. A. Nelder (2004)。 “条件和边际模型：另一种观点”。在：统计科学 19.2，p。 219–238.
   doi: 10.1214/088342304000000305。网址：https://doi.org/10.1214/088342304000000305
8. Lin, D. Y. 等人。 （2002 年）。 “基于累积残差的模型检查技术”。在：生物识别 58.1，p。 1-12。 （最后一次看到 19.09.2019）。 林 X. (1997)。 “具有随机效应的广义线性模型中的方差分量测试”。在：Biometrika 84.2，p。 309–326。 issn：00063444。网址：http://www.jstor.org/stable/2337459 （最后一次出现在 2019 年 9 月 19 日）。
9. Stiratelli, R. 等人。 （1984 年）。 “具有二元响应的串行观察的随机效应模型”。在：
   生物识别，p。 961–971。