为了教学目的,我正在使用R
的COUNT
包中的mdvis
数据集进行实验。我使用pscl
和countreg
包中的zeroinfl
函数拟合了零膨胀负二项模型。然而,pscl
包和countreg
包中的zeroinfl
的结果有很大不同。
模型和输出如下所示。
zeroinfl
从pscl
:
mdvisit_zeroinf<-pscl::zeroinfl(numvisit ~ reform + badh + agegrp + educ + loginc | reform + badh + agegrp + educ + loginc, dist = "negbin", data = mdvis, control = zeroinfl.control(method = "BFGS", EM=F)
summary(mdvisit_zeroinf)
Call:
pscl::zeroinfl(formula = numvisit ~ reform + badh + agegrp + educ + loginc |
reform + badh + agegrp + educ + loginc, data = mdvis, dist = "negbin")
Pearson residuals:
Min 1Q Median 3Q Max
-0.9667 -0.8037 -0.3597 0.3154 9.4878
Count model coefficients (negbin with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.30316 0.56565 -0.536 0.5920
reform -0.09916 0.05543 -1.789 0.0736 .
badh 1.11555 0.07517 14.841 <2e-16 ***
agegrp2 0.11376 0.06990 1.628 0.1036
agegrp3 0.21126 0.08620 2.451 0.0143 *
educ2 0.07605 0.07209 1.055 0.2914
educ3 -0.03979 0.07295 -0.545 0.5854
loginc 0.13209 0.07499 1.761 0.0782 .
Log(theta) 0.04747 0.05815 0.816 0.4144
Zero-inflation model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -30.7110 772.2930 -0.040 0.968
reform 14.3274 763.7283 0.019 0.985
badh -1.9503 3.5392 -0.551 0.582
agegrp2 10.9236 113.6292 0.096 0.923
agegrp3 9.7723 113.6558 0.086 0.931
educ2 -0.6632 1.9878 -0.334 0.739
educ3 -0.8743 1.3501 -0.648 0.517
loginc 0.5437 1.9718 0.276 0.783
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Theta = 1.0486
Number of iterations in BFGS optimization: 106
Log-likelihood: -4557 on 17 Df`
zeroinfl
从countreg
:
mdvisit_zeroinf2<-countreg::zeroinfl(numvisit ~ reform + badh + agegrp + educ + loginc | reform + badh + agegrp + educ + loginc, dist = "negbin", data = mdvis, control = zeroinfl.control(method = "BFGS", EM=F))
summary(mdvisit_zeroinf2)
Call:
countreg::zeroinfl(formula = numvisit ~ reform + badh + agegrp + educ + loginc |
reform + badh + agegrp + educ + loginc, data = mdvis, dist = "negbin",
control = zeroinfl.control(method = "BFGS", EM = F))
Pearson residuals:
Min 1Q Median 3Q Max
-0.9523 -0.8057 -0.3615 0.3106 9.6300
Count model coefficients (negbin with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.44473 0.53752 -0.827 0.40803
reform -0.12962 0.05106 -2.538 0.01113 *
badh 1.13558 0.07413 15.319 < 2e-16 ***
agegrp2 0.05851 0.06359 0.920 0.35756
agegrp3 0.18678 0.07176 2.603 0.00925 **
educ2 0.06398 0.06621 0.966 0.33385
educ3 -0.04825 0.07087 -0.681 0.49599
loginc 0.15338 0.07056 2.174 0.02973 *
Log(theta) 0.01232 0.04778 0.258 0.79652
Zero-inflation model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -119.137 95.259 -1.251 0.2111
reform 4.821 2.646 1.822 0.0684 .
badh -21.461 4614.989 -0.005 0.9963
agegrp2 9.261 28.341 0.327 0.7438
agegrp3 3.502 27.984 0.125 0.9004
educ2 -19.790 203.802 -0.097 0.9226
educ3 -7.894 4.758 -1.659 0.0971 .
loginc 13.010 10.463 1.243 0.2137
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Theta = 1.0124
Number of iterations in BFGS optimization: 197
Log-likelihood: -4554 on 17 Df
Stata
的结果与pscl
包中zeroinfl
的结果非常相似(特别是模型的计数部分的结果)。
我通过使用zeroinfl.control()
指定相同的选项来完全相似地拟合模型。我还试图在网上搜索,看看其他人以前是否报告过类似的问题,以及是否已经有答案。我也在stackoverflow
和CV
上搜索答案。但我无法得到任何答案。
1条答案
按热度按时间62lalag41#
问题来源:
问题是,没有零通货膨胀,但实际上比负二项模型预期的零更少-特别是在一些关于
reform
和agegrp
的子群中。症状:
由于没有零膨胀,二元零膨胀部分尽可能地产生预测的零膨胀概率,这些概率对于某些子群非常接近于零,特别注意模型中这部分的截距非常小,标准误差很大。
总的来说,这看起来类似于二进制回归模型中的准完全分离。可能性变得非常平坦,并且取决于优化器的设置,它确切地停止在哪里(这些设置在
pscl
和countreg
之间有所不同)。由于零通货膨胀模型的两个组成部分(计数与二进制部分)不能单独估计,一个组成部分估计中的问题也可能导致另一个组成部分估计中的问题/差异。
备选:
障碍模型在这里有几个优点:(1)它不仅可以处理零通货膨胀,而且可以处理缺零的情况。(2)模型的两个组成部分可以分别估计,因此问题不会溢出。
图示:
让我们将基本负二项模型与零通胀和障碍模型进行比较:
就对数似然性而言,零通胀模型在基本模型的基础上只得到了非常轻微的改进,尽管需要几乎两倍的参数。障碍模型的改进稍微多一点,但也不多:
就AIC和BIC而言,零通胀模型是这三个模型中最差的。AIC稍微倾向于障碍模型,而BIC更倾向于基本模型。
将根图作为基本模型的诊断显示,我们可以看到,总体而言,实际上观察到的零比负二项模型的预期要少,但基本模型基本上已经相当好地拟合了。