R 中的零膨胀过度分散计数数据 glmmTMB 错误
zero-inflated overdispersed count data glmmTMB error in R
我正在处理计数数据（可用 here) that are zero-inflated and overdispersed and has random effects. The package best suited to work with this sort of data is the glmmTMB (details here and troubleshooting here）。
在处理数据之前，我检查了它的正态性（零膨胀）、方差同质性、相关性和异常值。数据有两个异常值，我从上面的数据集 linekd 中删除了它们。有来自 18 个位置 (prop_id) 的 351 个观测值。
数据如下所示：
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 代码后得到的错误是：
Error in optimHess(par.fixed, obj$fn, obj$gr) : gradient in optim
evaluated to length 1 not 15 In addition: There were 50 or more
warnings (use warnings() to see the first 50)
warnings() 给出：
1. In (function (start, objective, gradient = NULL, hessian = NULL, ... :
NA/NaN function evaluation
我通常也指居中和标准化我的数值变量，但这只会消除第一个错误并保留 NA/NaN 错误。我尝试添加 glmmTMBControl 语句，如 ，但它打开了一个全新的错误世界。
我该如何解决这个问题？我做错了什么？
如有详细说明，我将不胜感激，以便我以后可以更好地学习如何解决此问题。 或者，我对 MCMCglmm 解决方案持开放态度，因为该函数也可以处理此类数据（尽管 运行 需要更长的时间）。
答案不完整...
- identity-link 有限域响应分布模型（例如 Gamma 或 Poisson，其中不可能出现负值）在计算上存在问题；在我看来，它们在概念上也经常存在问题，尽管有一些合理的论据支持它们。你有充分的理由这样做吗？
- 对于您要拟合的模型，这是一个非常小的数据集：13 个固定效应预测变量和 2 个随机效应预测变量。经验法则是您需要大约 10-20 倍的观测值：似乎 适合您的 345 次左右的观测值，但是......只有 40 次观测值是非零的！这意味着您的 'effective' 数量的 observations/amount 信息将少得多（有关这一点的更多讨论，请参阅 Frank Harrell 的 回归建模策略 ）。
话虽如此，让我 运行 回顾一下我尝试过的一些事情以及最终的结果。
GGally::ggpairs(euc, columns=2:10) 没有检测到任何明显糟糕的数据（我确实用 euc0==78 丢弃了数据点）
为了尝试使 identity-link 模型工作，我在 glmmTMB 中添加了一些代码。您应该能够通过 remotes::install_github("glmmTMB/glmmTMB/glmmTMB@clamp") 安装（请注意，您需要安装编译器等才能安装它）。此版本采用负预测值并强制它们为非负，同时对负对数似然添加相应的惩罚。
使用新版本的 glmmTMB 我没有收到错误，但我确实收到了这些警告：
Warning messages:
1: In fitTMB(TMBStruc) :
Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
2: In fitTMB(TMBStruc) :
Model convergence problem; false convergence (8). See vignette('troubleshooting')
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(.)))
这会有所帮助 - 现在我们仍在做一些愚蠢的事情，例如将年份视为数字变量而不将其居中（因此模型的 'intercept' 是公历的第 0 年...)
但这仍然不能解决问题。
更仔细地观察 ggpairs 图，看起来 season 和 year 混淆了：with(eucsc,table(season,year)) 表明观察发生在 Spring 和一年是冬天，另一年是秋天。 season 和 month 也很困惑：如果我们知道月份，那么我们就自动知道季节。
此时我决定放弃身份 link 看看会发生什么。 update(<previous_model>, family=poisson)（即使用泊松和标准对数 link）有效！使用 family=nbinom2 也是如此。
我查看了结果，发现降水 X 季节系数的 CI 很疯狂，因此删除了交互项 (update(f2S_noyr_logNB, . ~ . - precip:season))，此时结果看起来很合理。
一些最后的说明：
- 与样方相关的方差实际上为零
- 我认为你不一定需要零-inflation；低均值和过度分散（即
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 中实现。
我正在处理计数数据（可用 here) that are zero-inflated and overdispersed and has random effects. The package best suited to work with this sort of data is the glmmTMB (details here and troubleshooting here）。
在处理数据之前，我检查了它的正态性（零膨胀）、方差同质性、相关性和异常值。数据有两个异常值，我从上面的数据集 linekd 中删除了它们。有来自 18 个位置 (prop_id) 的 351 个观测值。
数据如下所示：
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 代码后得到的错误是：
Error in optimHess(par.fixed, obj$fn, obj$gr) : gradient in optim evaluated to length 1 not 15 In addition: There were 50 or more warnings (use warnings() to see the first 50)
warnings() 给出：
1. In (function (start, objective, gradient = NULL, hessian = NULL, ... :
NA/NaN function evaluation
我通常也指居中和标准化我的数值变量，但这只会消除第一个错误并保留 NA/NaN 错误。我尝试添加 glmmTMBControl 语句，如
我该如何解决这个问题？我做错了什么？
如有详细说明，我将不胜感激，以便我以后可以更好地学习如何解决此问题。 或者，我对 MCMCglmm 解决方案持开放态度，因为该函数也可以处理此类数据（尽管 运行 需要更长的时间）。
答案不完整...
- identity-link 有限域响应分布模型（例如 Gamma 或 Poisson，其中不可能出现负值）在计算上存在问题；在我看来，它们在概念上也经常存在问题，尽管有一些合理的论据支持它们。你有充分的理由这样做吗？
- 对于您要拟合的模型，这是一个非常小的数据集：13 个固定效应预测变量和 2 个随机效应预测变量。经验法则是您需要大约 10-20 倍的观测值：似乎 适合您的 345 次左右的观测值，但是......只有 40 次观测值是非零的！这意味着您的 'effective' 数量的 observations/amount 信息将少得多（有关这一点的更多讨论，请参阅 Frank Harrell 的 回归建模策略 ）。
话虽如此，让我 运行 回顾一下我尝试过的一些事情以及最终的结果。
GGally::ggpairs(euc, columns=2:10)没有检测到任何明显糟糕的数据（我确实用euc0==78丢弃了数据点）
为了尝试使 identity-link 模型工作，我在 glmmTMB 中添加了一些代码。您应该能够通过 remotes::install_github("glmmTMB/glmmTMB/glmmTMB@clamp") 安装（请注意，您需要安装编译器等才能安装它）。此版本采用负预测值并强制它们为非负，同时对负对数似然添加相应的惩罚。
使用新版本的 glmmTMB 我没有收到错误，但我确实收到了这些警告：
Warning messages: 1: In fitTMB(TMBStruc) : Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
2: In fitTMB(TMBStruc) : Model convergence problem; false convergence (8). See vignette('troubleshooting')
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(.)))
这会有所帮助 - 现在我们仍在做一些愚蠢的事情，例如将年份视为数字变量而不将其居中（因此模型的 'intercept' 是公历的第 0 年...)
但这仍然不能解决问题。
更仔细地观察 ggpairs 图，看起来 season 和 year 混淆了：with(eucsc,table(season,year)) 表明观察发生在 Spring 和一年是冬天，另一年是秋天。 season 和 month 也很困惑：如果我们知道月份，那么我们就自动知道季节。
此时我决定放弃身份 link 看看会发生什么。 update(<previous_model>, family=poisson)（即使用泊松和标准对数 link）有效！使用 family=nbinom2 也是如此。
我查看了结果，发现降水 X 季节系数的 CI 很疯狂，因此删除了交互项 (update(f2S_noyr_logNB, . ~ . - precip:season))，此时结果看起来很合理。
一些最后的说明：
- 与样方相关的方差实际上为零
- 我认为你不一定需要零-inflation；低均值和过度分散（即
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 中实现。