`pscl`和`countreg`包中的`zeroinfl`给予了非常不同的结果,为什么?

8gsdolmq  于 2023-04-27  发布在  其他
关注(0)|答案(1)|浏览(113)

为了教学目的,我正在使用RCOUNT包中的mdvis数据集进行实验。我使用psclcountreg包中的zeroinfl函数拟合了零膨胀负二项模型。然而,pscl包和countreg包中的zeroinfl的结果有很大不同。
模型和输出如下所示。

zeroinflpscl

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`

zeroinflcountreg

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()指定相同的选项来完全相似地拟合模型。我还试图在网上搜索,看看其他人以前是否报告过类似的问题,以及是否已经有答案。我也在stackoverflowCV上搜索答案。但我无法得到任何答案。

62lalag4

62lalag41#

问题来源:

问题是,没有零通货膨胀,但实际上比负二项模型预期的零更少-特别是在一些关于reformagegrp的子群中。

症状:

由于没有零膨胀,二元零膨胀部分尽可能地产生预测的零膨胀概率,这些概率对于某些子群非常接近于零,特别注意模型中这部分的截距非常小,标准误差很大。
总的来说,这看起来类似于二进制回归模型中的准完全分离。可能性变得非常平坦,并且取决于优化器的设置,它确切地停止在哪里(这些设置在psclcountreg之间有所不同)。
由于零通货膨胀模型的两个组成部分(计数与二进制部分)不能单独估计,一个组成部分估计中的问题也可能导致另一个组成部分估计中的问题/差异。

备选:

障碍模型在这里有几个优点:(1)它不仅可以处理零通货膨胀,而且可以处理缺零的情况。(2)模型的两个组成部分可以分别估计,因此问题不会溢出。

图示:

让我们将基本负二项模型与零通胀和障碍模型进行比较:

data("mdvis", package = "COUNT")
mdvis <- transform(mdvis,
  reform = factor(reform),
  badh = factor(badh),
  agegrp = factor(agegrp),
  educ = factor(educ)
)
library("countreg")
f <- numvisit ~ reform + badh + agegrp + educ + loginc
m <- glm.nb(f, data = mdvis)
m2 <- zeroinfl(f, dist = "negbin", data = mdvis)
m3 <- hurdle(f, dist = "negbin", data = mdvis)

就对数似然性而言,零通胀模型在基本模型的基础上只得到了非常轻微的改进,尽管需要几乎两倍的参数。障碍模型的改进稍微多一点,但也不多:

c(logLik(m), logLik(m2), logLik(m3))
## [1] -4560.910 -4554.146 -4550.520

就AIC和BIC而言,零通胀模型是这三个模型中最差的。AIC稍微倾向于障碍模型,而BIC更倾向于基本模型。

AIC(m, m2, m3)
##    df      AIC
## m   9 9139.819
## m2 17 9142.293
## m3 17 9135.040
BIC(m, m2, m3)
##    df      BIC
## m   9 9191.195
## m2 17 9239.336
## m3 17 9232.083

将根图作为基本模型的诊断显示,我们可以看到,总体而言,实际上观察到的零比负二项模型的预期要少,但基本模型基本上已经相当好地拟合了。

rootogram(m)

  • (注意:具有置信限的Rootogram版本实际上是由另一个正在开发的软件包生成的。countreg中的版本没有这些置信限。)*

相关问题