嵌套设置中的混合效应模型或多重回归比较

Posted 2023-02-14

【中文标题】嵌套设置中的混合效应模型或多重回归比较

【英文标题】：Mixed effect model or multiple regressions comparison in nested setup

【发布时间】：2022-01-23 05:39:29

【问题描述】：

我有一个回复Y，这是一个介于 0-1 之间的百分比。我的数据按分类或进化关系嵌套，例如 phylum/genus/family/species，我有一个连续协变量 temp 和一个分类协变量 fac，级别为 fac1 和 fac2。

我有兴趣估算：

Y 在fac1 和fac2 之间是否存在差异（截距）以及由此解释的差异程度 每个级别的 fac 对 temp 的响应是否不同（线性斜率） Y 对于我的分类的每个级别是否存在差异，以及它们解释了多少差异（请参阅 varcomp） 我的分类的每个级别对temp 的响应是否不同（线性斜率）

一个蛮力的想法是将我的数据拆分为这里分类最低的物种，为每个物种 i 做一个线性 beta 回归 betareg(Y(i)~temp) 。然后提取每个物种的斜率和截距，并将它们分组到每个 fac 更高的分类水平，并比较斜率（截距）的分布，例如，通过 Kullback-Leibler 散度与我在引导我的 Y 值时得到的分布。或者分别比较分类水平或我的因子 fac 之间的斜率（或截距）分布。或者只比较分类水平或我的因子水平之间的平均斜率和截距。 不确定这是不是个好主意。而且也不确定如何回答我的分类级别解释了多少方差的问题，例如在嵌套随机混合效应模型中。

另一种选择可能只是那些混合模型，但我怎样才能在一个模型中包含我想要测试的所有方面

说我可以使用“gamlss”包来做：

library(gamlss)

model<-gamlss(Y~temp*fac+re(random=~1|phylum/genus/family/species),family=BE)

但在这里我看不到任何方法可以合并随机斜率，或者我可以这样做：

model<-gamlss(Y~re(random=~temp*fac|phylum/genus/family/species),family=BE)

但是对 lme 的内部调用对此有一些问题，并且猜测这无论如何都不是正确的表示法。 有什么方法可以实现我想要测试的东西，不一定是 gamlss 而是任何其他包含嵌套结构和 beta 回归的包？ 谢谢！

【问题讨论】：

您的数据中是否有精确的 0/1 值（尤其是精确的 1 值）？ glmmTMB 应该能够处理具有随机斜率的 Beta 分布响应......请注意，在简单的方差分解不一定有效的混合模型世界中，“解释了多少方差”可能是一个棘手的问题...

brms 也可以做到这一点，我认为，虽然你必须弄清楚整个贝叶斯/MCMC 的事情（值得注意的是，它也可以处理零一膨胀的 Beta）

非常感谢您的帮助！你知道带有 glmmTMB 包的模型是什么样子吗？我不确定如何将因子 fac 放在那里同时回答上面的 2) 和 4)。你会说蛮力方法通常也可以接受吗？再次感谢

【参考方案1】：

在glmmTMB 中，如果您的响应中没有确切的 0 或 1 值，这样的事情应该可以工作：

library(glmmTMB)
glmmTMB(Y ~ temp*fac + (1 + temp | phylum/genus/family/species),
           data = ...,
           family = beta_family)

如果您的值为零，您将需要做一些事情。例如，您可以在glmmTMB 中添加一个零通胀术语； brms 可以处理零一膨胀的 Beta 响应；您可以稍微“挤压” 0/1 值（参见 Smithson 和 Verkuilen 关于 Beta 回归的论文的附录）。如果你只有几个 0/1 值，那么你做什么并不重要。如果你有很多，你需要花一些时间认真思考它们的含义，这将影响你处理它们的方式。它们是否代表审查（即不完全是 0/1 但太接近边界而无法测量差异的值）？它们是质量不同的反应吗？等等...） 正如我在评论中所说，计算 GLMM 的方差分量非常棘手 - 不一定有简单的分解，例如见here。但是，您可以计算每个分类级别的截距和斜率的方差并进行比较（您可以使用标准差与固定效应的大小进行比较...） 这里给出的模型可能要求很高，具体取决于系统发育的大小 - 例如，您可能没有足够的门级复制（在这种情况下，您可以拟合模型 ~ temp*(fac + phylum) + (1 + temp | phylum:(genus/family/species))，即拔出门效应作为固定效应）。 这是假设您愿意假设fac 的影响以及它与temp 的相互作用在整个系统发育过程中不会有所不同...

【讨论】：

非常感谢这个详细的回答！这是超级有帮助的。

对于方差分解，我找到了这个包：rdrr.io/cran/MuMIn/man/r.squaredGLMM.html 这是什么意思，在这里？但也会听从你的建议。您介意简要说明如何将随机效应的 sd 与答案中固定效应的大小进行比较吗？我只是没有得到每个 fac 级别的斜率/截距吗？ - 还有一个小问题。为什么我需要 temp 作为“全局”固定效果？非常感谢，从中学到了很多关于混合模型的知识。

- 我想知道是否（在我没有数据的先验信息的情况下）对每个 fac 级别 j 中的每个物种 i 使用 beta 回归所以 Yij ~tempij，然后分别提取斜率然后使用说条件推理树，例如ctree();这样 ctree(slopes~phylum+genus+class+family+fac, 可以给我任何关于在哪里期望不同的分组响应 temp over taxonomy 和 fac ....