零膨胀泊松模型无法拟合该数据

6

以下是数据和设置:

library(fitdistrplus)
library(gamlss)

finalVector <- dput(finalVector)
c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 2, 1, 1, 1, 1, 1, 
1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 2, 1, 3, 2, 1, 1, 1, 1, 1, 1, 2, 
2, 1, 4, 2, 3, 1, 2, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
2, 1, 2, 2, 1, 1, 4, 1, 2, 2, 1, 1, 1, 1, 1, 2, 1, 1, 2, 2, 1, 
2, 1, 1, 4, 2, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1)

countFitPoisson <- fitdist(finalVector,"pois",  method = "mle", lower = 0)
countFitZeroPoisson <- fitdist(finalVector, 'ZIP', start = list( ##mu  = mean of poisson, sigma = prob(x = 0))
                                                           mu = as.numeric(countFitPoisson$estimate), 
                                                            sigma = as.numeric(as.numeric(countFitPoisson$estimate))
                                                          ), method = "mle", lower= 0) 

第一次调用成功。最后一次调用失败,我不确定原因。谢谢!
编辑:
假设我编写的代码是正确的(不确定),那么我唯一能想到的就是模型拟合时零值不够?

1
当您填充几个0(5)并将起始参数更改为更接近mle值(0.7,0.05)时,它似乎可以工作。 - Vlo
1个回答

8
您的数据不真正存在零膨胀,因此拟合模型并不能带来改进。我将使用glm和扩展回归模型,而不是使用fitdistr方法。所有回归(子)模型仅使用常数(或截距),没有实际的回归变量。为了可视化,我使用根图通过R-Forge提供的countreg包(其中包含pscl计数数据回归的后继函数)。
首先,让我们看看泊松拟合:
(mp <- glm(finalVector ~ 1, family = poisson))
## Call:  glm(formula = finalVector ~ 1, family = poisson)
## 
## Coefficients:
## (Intercept)  
##      -0.284  
## 
## Degrees of Freedom: 181 Total (i.e. Null);  181 Residual
## Null Deviance:      200.2 
## Residual Deviance: 200.2        AIC: 418.3

这对应于一个拟合的平均值为exp(-0.284),即约为0.753。如果您比较观察到的和拟合的频率,这非常适合数据:

library("countreg")
rootogram(mp)

这表明0、1、2的计数适应性基本完美,3、4、5只有轻微偏差,但这些频率非常低。因此,从这个结果来看,似乎不需要扩展模型。
但为了正式比较模型与其他模型,可以考虑使用零膨胀泊松分布(如您尝试过的)、障碍泊松分布或负二项分布。零膨胀模型如下:
(mzip <- zeroinfl(finalVector ~ 1 | 1, dist = "poisson"))
## Call:
## zeroinfl(formula = finalVector ~ 1 | 1, dist = "poisson")
## 
## Count model coefficients (poisson with log link):
## (Intercept)  
##     -0.2839  
## 
## Zero-inflation model coefficients (binomial with logit link):
## (Intercept)  
##      -9.151  

因此,计数均值基本上与之前相同,零膨胀的概率基本上为零(plogis(-9.151)约为0.01%)。
障碍模型的工作方式类似,但可以使用零截断泊松模型进行0与大于0的区分,并使用截断泊松模型进行正计数。然后,这也嵌套在泊松模型中,因此可以轻松进行Wald检验:
mhp <- hurdle(finalVector ~ 1 | 1, dist = "poisson", zero.dist = "poisson")
hurdletest(mhp)
## Wald test for hurdle models
## 
## Restrictions:
## count_((Intercept) - zero_(Intercept) = 0
## 
## Model 1: restricted model
## Model 2: finalVector ~ 1 | 1
## 
##   Res.Df Df Chisq Pr(>Chisq)
## 1    181                    
## 2    180  1 0.036     0.8495

这也清楚地表明,没有多余的零,并且简单的泊松模型就足够了。

最后也可以考虑负二项式模型:

(mnb <- glm.nb(finalVector ~ 1))
## Call:  glm.nb(formula = finalVector ~ 1, init.theta = 125.8922776, link = log)
## 
## Coefficients:
## (Intercept)  
##      -0.284  
## 
## Degrees of Freedom: 181 Total (i.e. Null);  181 Residual
## Null Deviance:      199.1 
## Residual Deviance: 199.1        AIC: 420.3

这个模型的平均值与泊松分布非常相似,而且大的θ参数足够接近无穷大(=泊松分布)。因此,总体来说,泊松模型已经足够简单,不需要考虑任何扩展。似然函数几乎没有改变,额外的参数(零膨胀、零障碍、θ离散度)也没有带来任何改进:

AIC(mp, mzip, mhp, mnb)
##      df      AIC
## mp    1 418.2993
## mzip  2 420.2996
## mhp   2 420.2631
## mnb   2 420.2959

我完全同意你的分析 -- 但是我想知道,包是否应该返回一个常规的泊松模型,而不是崩溃? - user1357015
2
它不会“崩溃”,但会以足够信息的错误终止:sigma 变为 0,但必须在 0 和 1 之间。这表明存在退化模型,因此需要出现错误。zeroinfl 函数通过使用 logit 链接来避免这种情况,以便概率必须严格在 0 和 1 之间。但您仍然可以看到退化,因为在 logit 比例上的 -9.2 的截距已经非常接近于 -Infinity。使用 reltol = 1e-12 将导致 -11.7,更低的 reltol 值也会导致无法收敛... - Achim Zeileis

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接