时间戳

首页 > 解决方案 > zeroinfl 模型 - 警告消息：在 sqrt(diag(object$vcov)) 中：产生了 NaN

r - zeroinfl 模型 - 警告消息：在 sqrt(diag(object$vcov)) 中：产生了 NaN

问题描述

我正在尝试运行零膨胀负二项式，但在检查阻止我看到结果的模型时，我遇到了“产生的 NaN”警告。这是一些模拟数据，它们是我的实际数据的简化版本——我的真实数据对每个物种有更多的观察结果 + 更多的物种：

df1 <- data.frame(species = c("Rufl","Rufl","Rufl","Rufl","Assp","Assp","Assp","Assp","Elre", "Elre","Elre", "Elre","Soca","Soca","Soca","Soca"),
                  state = c("warmed","ambient","warmed","ambient","warmed","ambient","warmed","ambient","warmed","ambient","warmed","ambient","warmed","ambient","warmed","ambient"),
                  p_eaten = c(0, 0, 3, 0, 0, 1, 15, 0, 20, 0, 0, 2, 0, 3, 87, 0))

这是我试图运行的模型，状态和物种之间存在交互：

library(pscl)
mod1 <- zeroinfl(p_eaten ~ state * species,
                     dist = 'negbin',
                     data = df1)
summary(mod1)

这是我得到的时候Warning message: In sqrt(diag(object$vcov)) : NaNs produced。如何修复此警告消息，以便能够查看模型结果？谢谢！

使用 R 版本 4.0.2，Mac OS X 10.13.6

标签： r

解决方案

对于如此复杂的模型来说，这是非常薄的数据，但是如果您xtabs对您的数据框版本执行操作，您会发现您的参考类别之一的计数为零。如果您交换state变量的级别，则 NA 会消失，尽管一些较大的标准错误仍然存​​在。

xtabs(p_eaten~ state + species, data=df1)
         species
state     Assp Elre Rufl Soca
  ambient    1    2    0    3
  warmed    15   20    3   87

未清理的控制台输出如下：

df1 <- data.frame(species = c("Rufl","Rufl","Rufl","Rufl","Assp","Assp","Assp","Assp","Elre", "Elre","Elre", "Elre","Soca","Soca","Soca","Soca"),
+                   state = factor(c("warmed","ambient","warmed","ambient","warmed","ambient","warmed","ambient","warmed","ambient","warmed","ambient","warmed","ambient","warmed","ambient"), levels=c("warmed","ambient")),
+                   p_eaten = c(0, 0, 3, 0, 0, 1, 15, 0, 20, 0, 0, 2, 0, 3, 87, 0))
> xtabs(p_eaten~ state + species, data=df1)
         species
state     Assp Elre Rufl Soca
  warmed    15   20    3   87
  ambient    1    2    0    3

试图：

> library(pscl)
> mod1 <- zeroinfl(p_eaten ~ state * species,
+                      dist = 'negbin',
+                      data = df1)
> summary(mod1)

Call:
zeroinfl(formula = p_eaten ~ state * species, data = df1, dist = "negbin")

Pearson residuals:
     Min       1Q   Median       3Q      Max 
-0.98868 -0.80384 -0.00342  0.80387  0.98872 

Count model coefficients (negbin with log link):
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               2.708e+00  2.582e-01  10.488  < 2e-16 ***
stateambient             -3.401e+00  1.033e+00  -3.292 0.000994 ***
speciesElre               2.877e-01  3.416e-01   0.842 0.399623    
speciesRufl              -1.671e+00  6.874e-01  -2.431 0.015068 *  
speciesSoca               1.758e+00  2.796e-01   6.288 3.22e-10 ***
stateambient:speciesElre  8.714e-01  1.400e+00   0.622 0.533627    
stateambient:speciesRufl -3.972e-04         NA      NA       NA    
stateambient:speciesSoca -2.763e-02  1.218e+00  -0.023 0.981906    
Log(theta)                1.501e+01  1.848e+02   0.081 0.935234    

Zero-inflation model coefficients (binomial with logit link):
                           Estimate Std. Error z value Pr(>|z|)
(Intercept)              -1.327e-04  1.414e+00   0.000    1.000
stateambient             -8.538e+00  1.206e+02  -0.071    0.944
speciesElre               1.379e-04  2.000e+00   0.000    1.000
speciesRufl              -1.267e-01  2.083e+00  -0.061    0.952
speciesSoca               1.757e-04  2.000e+00   0.000    1.000
stateambient:speciesElre  8.016e+00  1.206e+02   0.066    0.947
stateambient:speciesRufl  1.757e+01         NA      NA       NA
stateambient:speciesSoca  8.411e+00  1.206e+02   0.070    0.944
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta = 3316847.5216 
Number of iterations in BFGS optimization: 41 
Log-likelihood: -21.87 on 17 Df
Warning message:
In sqrt(diag(object$vcov)) : NaNs produced

方差-协方差矩阵应该是正定的，负值可能会妨碍反演工作。当在 CrossValidated.com 上提出类似问题时，Ben Bolker 建议使用 brglm2 版本的 glm：

> library(brglm2)
> summary(m1 <- glm(p_eaten~ state * species, data=df1,
+             family=poisson,
+             method="brglmFit"))

Call:
glm(formula = p_eaten ~ state * species, family = poisson, data = df1, 
    method = "brglmFit")

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-9.3541  -1.8708  -0.7071   0.8567   5.7542  

Coefficients:
                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)                2.0477     0.2540   8.062 7.52e-16 ***
stateambient              -2.3354     0.8551  -2.731  0.00631 ** 
speciesElre                0.2796     0.3366   0.831  0.40619    
speciesRufl               -1.4881     0.5918  -2.514  0.01192 *  
speciesSoca                1.7308     0.2756   6.281 3.37e-10 ***
stateambient:speciesElre   0.2312     1.0863   0.213  0.83142    
stateambient:speciesRufl   0.3895     1.7369   0.224  0.82258    
stateambient:speciesSoca  -0.8835     1.0141  -0.871  0.38362    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 443.21  on 15  degrees of freedom
Residual deviance: 183.08  on  8  degrees of freedom
AIC: 225.38

Number of Fisher Scoring iterations: 1

事实上，互动的必要性似乎值得怀疑。删除交互项后偏差的变化是微乎其微的。很难知道您的完整数据集是否也会出现这种情况。

推荐阅读