R中的零膨胀过度分散计数数据glmmTMB错误

Question

我正在使用零膨胀和过度分散且具有随机效应的计数数据（可用 here）。最适合处理此类数据的软件包是 glmmTMB（详情 here 和故障排除 here）。

在处理数据之前，我检查了它的正态性（它是零膨胀的）、方差同质性、相关性和异常值。数据有两个异常值，我从上面的数据集 linekd 中删除了它们。来自 18 个位置的 351 个观测值 (prop_id)。

数据如下：

euc0 ea_grass ep_grass np_grass np_other_grass month year precip season   prop_id quad
 3      5.7      0.0     16.7            4.0     7 2006    526 Winter    Barlow    1
 0      6.7      0.0     28.3            0.0     7 2006    525 Winter    Barlow    2
 0      2.3      0.0      3.3            0.0     7 2006    524 Winter    Barlow    3
 0      1.7      0.0     13.3            0.0     7 2006    845 Winter    Blaber    4
 0      5.7      0.0     45.0            0.0     7 2006    817 Winter    Blaber    5
 0     11.7      1.7     46.7            0.0     7 2006    607 Winter    DClark    3

响应变量为euc0，随机效应为prop_id 和quad_id。其余变量是固定效应（均代表不同植物物种的覆盖百分比）。

我要运行的模型：

library(glmmTMB)
seed0<-glmmTMB(euc0 ~ ea_grass + ep_grass + np_grass + np_other_grass + month + year*precip + season*precip + (1|prop_id)  + (1|quad), data = euc, family=poisson(link=identity))

fit_zinbinom <- update(seed0,family=nbinom2) #allow variance increases quadratically

运行seed0代码后我得到的错误是：

optimHess(par.fixed, obj$fn, obj$gr) 中的错误：优化中的梯度 评估为长度 1 而不是 15 另外：有 50 或更多 警告（使用 warnings() 查看前 50 个）

warnings() 给出：

1. In (function (start, objective, gradient = NULL, hessian = NULL,  ... :
  NA/NaN function evaluation

我通常也指中心化和标准化我的数值变量，但这只会删除第一个错误并保留NA/NaN 错误。我尝试添加 glmmTMBControl 语句，例如 this OP，但它打开了一个全新的错误世界。

我该如何解决这个问题？我做错了什么？

非常感谢您提供详细的解释，以便我将来可以自己学习如何更好地解决此问题。 或者，我对MCMCglmm 解决方案持开放态度，因为该函数也可以处理此类数据（尽管运行时间更长）。

Answer 1

一个不完整的答案...

• 有限域响应分布（例如 Gamma 或 Poisson，其中不可能出现负值）的身份链接模型在计算上存在问题；在我看来，它们通常在概念上也存在问题，尽管有一些合理的论据支持它们。您有充分的理由这样做吗？
• 对于您要拟合的模型，这是一个非常小的数据集：13 个固定效应预测变量和 2 个随机效应预测变量。经验法则是，您需要大约 10 到 20 倍的观察结果：似乎符合您的 345 个左右观察结果，但是……您的观察结果中只有 40 个是非零！这意味着您的“有效”观察次数/信息量会少得多（有关这一点的更多讨论，请参阅 Frank Harrell 的回归建模策略）。

也就是说，让我回顾一下我尝试过的一些事情以及我最终的结果。

• GGally::ggpairs(euc, columns=2:10) 没有检测到任何明显可怕的数据（我确实用 euc0==78 丢弃了数据点）

为了尝试使身份链接模型工作，我在 glmmTMB 中添加了一些代码。你应该可以通过remotes::install_github("glmmTMB/glmmTMB/glmmTMB@clamp") 安装（注意你需要安装编译器等来安装它）。此版本采用负预测值并强制它们为非负，同时对负对数似然增加相应的惩罚。

使用新版本的 glmmTMB 我没有收到错误，但确实收到了以下警告：

警告信息： 1：适合TMB（TMBStruc）： 模型收敛问题；非正定 Hessian 矩阵。见小插图（'疑难解答'）
2：适合TMB（TMBStruc）： 模型收敛问题；错误收敛 (8)。见小插图（'疑难解答'）

Hessian（二阶导数）矩阵是非正定的，这意味着存在一些（仍然难以解决）问题。 heatmap(vcov(f2)$cond,Rowv=NA,Colv=NA) 让我看看协方差矩阵。 （我也喜欢corrplot::corrplot.mixed(cov2cor(vcov(f2)$cond),"ellipse","number")，但当vcov(.)$cond 是非正定时，这不起作用。在紧要关头，您可以使用sfsmisc::posdefify() 强制它为正定......）

尝试缩放：

eucsc <- dplyr::mutate_at(euc1,dplyr::vars(c(ea_grass:precip)), ~c(scale(.)))

这将有所帮助 - 现在我们仍在做一些愚蠢的事情，例如将年份视为数字变量而不将其居中（因此模型的“截距”位于公历的第 0 年......）

但这仍然不能解决问题。

更仔细地观察ggpairs 图，看起来season 和year 混淆了：with(eucsc,table(season,year)) 表明观察发生在一年的春季和冬季，另一年的秋季。 season 和 month 也很困惑：如果我们知道月份，那么我们就会自动知道季节。

此时我决定放弃身份链接，看看发生了什么。 update(<previous_model>, family=poisson)（即使用带有标准日志链接的泊松）有效！使用family=nbinom2 也是如此，效果要好得多。

我查看了结果，发现 precip X 季节系数的 CI 很疯狂，因此删除了交互项 (update(f2S_noyr_logNB, . ~ . - precip:season))，此时结果看起来很合理。

最后几点说明：

• 与样方相关的方差实际上为零
• 我认为你不一定需要零通货膨胀；低均值和过度分散（即family=nbinom2）可能就足够了。
• 残差的分布看起来不错，但似乎仍然存在一些模型不适合 (library(DHARMa); plot(simulateResiduals(f2S_noyr_logNB2)))。我会花一些时间针对各种预测变量组合绘制残差和预测值，看看您是否可以定位问题。

PS 一种更快的方法来查看固定效应是否存在问题（多重共线性）：

X <- model.matrix(~ ea_grass + ep_grass +
                   np_grass + np_other_grass + month +
                   year*precip + season*precip,
                  data=euc)
ncol(X)  ## 13
Matrix::rankMatrix(X) ## 11

lme4 有这样的测试，以及自动删除别名列的机制，但它们目前在glmmTMB 中没有实现。