零膨胀泊松模型无法拟合
Zero inflated poisson model fails to fit
这是数据和设置:
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)
第一次调用成功。决赛说它未能估计,我不确定为什么。谢谢!
编辑:
假设我正确地编写了代码(不确定),那么我唯一能想到的是没有足够个零来适合模型?
您的数据并非真正零膨胀,因此拟合模型不会带来改进。我没有使用 fitdistr
方法,而是使用 glm
和下面的扩展回归模型。不过,所有回归(子)模型都只使用一个常数(或截距)而没有任何真正的回归变量。对于可视化,我通过 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
因此,计数均值与之前基本相同,零-inflation 的概率基本为零(plogis(-9.151)
约为 0.01%)。
障碍模型的工作原理类似,但可以使用零删失泊松模型来表示 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
这再次具有几乎相同的均值和足够接近无穷大(= 泊松)的巨大 theta 参数。因此,总体而言,泊松模型就足够了,不需要考虑任何扩展。可能性几乎没有变化,附加参数(零-inflation、零障碍、theta-分散)没有产生任何改进:
AIC(mp, mzip, mhp, mnb)
## df AIC
## mp 1 418.2993
## mzip 2 420.2996
## mhp 2 420.2631
## mnb 2 420.2959
这是数据和设置:
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)
第一次调用成功。决赛说它未能估计,我不确定为什么。谢谢!
编辑:
假设我正确地编写了代码(不确定),那么我唯一能想到的是没有足够个零来适合模型?
您的数据并非真正零膨胀,因此拟合模型不会带来改进。我没有使用 fitdistr
方法,而是使用 glm
和下面的扩展回归模型。不过,所有回归(子)模型都只使用一个常数(或截距)而没有任何真正的回归变量。对于可视化,我通过 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
因此,计数均值与之前基本相同,零-inflation 的概率基本为零(plogis(-9.151)
约为 0.01%)。
障碍模型的工作原理类似,但可以使用零删失泊松模型来表示 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
这再次具有几乎相同的均值和足够接近无穷大(= 泊松)的巨大 theta 参数。因此,总体而言,泊松模型就足够了,不需要考虑任何扩展。可能性几乎没有变化,附加参数(零-inflation、零障碍、theta-分散)没有产生任何改进:
AIC(mp, mzip, mhp, mnb)
## df AIC
## mp 1 418.2993
## mzip 2 420.2996
## mhp 2 420.2631
## mnb 2 420.2959