1

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

在处理数据之前，我检查了它的正态性（它是零膨胀的）、方差同质性、相关性和异常值。数据有两个异常值，我从上面的数据集 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) 中的错误：optim 中的梯度评估为长度 1 而不是 15 另外：有 50 个或更多警告（使用 warnings() 查看前 50 个）

warnings()给出：

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

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

我怎样才能解决这个问题？我究竟做错了什么？

将不胜感激详细的解释，以便我将来可以学习如何更好地解决此问题。或者，我对MCMCglmm解决方案持开放态度，因为该功能也可以处理此类数据（尽管运行时间更长）。

4

1 回答 1

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: In fitTMB(TMBStruc) : 模型收敛问题；非正定 Hessian 矩阵。见 vignette('troubleshooting')
2: In fitTMB(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)（即使用带有标准日志链接的泊松）有效！using 也是如此family=nbinom2，这要好得多。

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

最后的几点说明：

  • 与 quadrat 相关的方差实际上为零
  • 我认为您不一定需要零通货膨胀；低均值和过度分散（即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。

于 2020-05-19T01:12:46.190 回答