zeroinfl模型-警告消息:sqrt(diag(对象$vcov))：生成NaNs

Question

我试图运行零膨胀负二项环，但我遇到了一个"NaNs生产“警告时，检查模型，阻止我看到结果。下面是一些模拟数据，这是我实际数据的简化版本--我的真实数据有更多的每个物种的观察结果+更多的物种：

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

Answer 1

这很可能是完全分离的一种情况，尽管没有完整的数据集是不可能确定的。

这很可能发生在您的类别是全为零，或全部-非零。在上面给出的示例中：

with(df1,table(species,state,p_eaten==0))

在species=="Rufl“、state==”combination“和p==0是FALSE的情况下，这些都不是观测值；换句话说，对于这些因素组合，所有的观测结果都是零的。因此，与这种状态相比较的任何系数都有一个大幅度的参数值(即abs(β) >> 1)；理论上它应该是无限的，但通常在10到30之间(取决于数值方法放弃的地方)。这些系数要么具有非常大的标准误差和(Wald)置信区间，要么具有(在您的情况下) NaN值。

这种描述适用于speciesRufl和statewarmed:speciesRufl的零通货膨胀系数.计数模型系数不太大，但仍然存在NaN标准误差，我认为是因为它们的不确定性与零通货膨胀系数的不确定性有关。

Count model coefficients (negbin with log link):
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)             -0.69312    1.00016  -0.693 0.488307    
...
speciesRufl             -1.69427        NaN     NaN      NaN    

Zero-inflation model coefficients (binomial with logit link):
                        Estimate Std. Error z value Pr(>|z|)
...
speciesRufl               17.560        NaN     NaN      NaN

你能做些什么？

• 忽略这个问题。您将无法获得可靠的标准错误(因此p值、置信区间等)。对于受问题影响的系数，但是模型在原则上仍然是可以的。进行似然比检验(将模型的对数可能性与没有指定的预测器集合的模型进行比较)仍然可以，并可用于获得特定效果的p值，例如通过lmtest::lrtest()。
• 简化你的模型:把一些类别分类起来，决定你是否真的需要零通货膨胀，等等。
• 有许多其他方法涉及惩罚或在相关系数上强加贝叶斯优先级以保持它们的合理性(例如brglm2包)，但我不知道这些方法是否对zeroinfl模型实现/可用您可以这样做，例如通过brms包，但这需要大量的工作才能跟上包的基础。

Answer 2

对于这样一个复杂的模型来说，这是非常薄的数据，但是如果您在您的dataframe版本上做了一个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上提出类似的问题时，Bolker建议使用glm的brglm2版本：

> 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

事实上，是否需要进行互动似乎是有问题的。随着相互作用项的去除，偏差的变化很小。很难知道mich是否也会出现在您的完整数据集中。