pscl与countreg包中zeroinfl函数结果差异的原因探究

我以教学为目的，使用R语言COUNT包中的mdvis数据集，分别调用pscl和countreg包的zeroinfl函数拟合零膨胀负二项模型，但两者输出结果差异显著。其中Stata的拟合结果与pscl包的zeroinfl结果高度相似，尤其是计数部分。我已通过zeroinfl.control()设置相同参数，特询问产生该差异的原因。

pscl包zeroinfl拟合结果 #

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

countreg包zeroinfl拟合结果 #

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

差异原因分析 #

• 负二项分布参数化差异：两个包对负二项分布的过度离散参数theta定义不同。pscl采用方差 = 均值 + 均值²/theta的参数化方式，countreg则使用方差 = 均值*(1 + theta*均值)，这种参数化差异会直接导致计数模型系数、theta估计值及似然计算结果出现偏差。

• 优化初始化策略不同：尽管都设置了BFGS优化方法且关闭EM算法，两个包的默认起始值选择逻辑不同。pscl的起始值更贴近Stata的拟合策略，易收敛到与Stata一致的局部最优解；而countreg的起始值可能引导优化过程收敛到另一个似然值略高的局部最优（从log-likelihood值-4554 vs -4557可看出）。

• 零膨胀模型的数值稳定性差异：当前数据中零膨胀部分的系数估计极不稳定（两个结果中零膨胀截距都极大，标准误也异常大），说明零膨胀部分的识别性较弱。此时，两个包在似然函数梯度计算、数值优化精度控制（如收敛阈值、步长调整）上的细微差异，会被放大为显著的结果差异。

• 底层实现逻辑差异：pscl的zeroinfl是较早实现的零膨胀模型工具，基于独立的优化逻辑；而countreg的zeroinfl基于广义线性模型框架扩展，两者在似然函数计算细节、优化算法底层实现上存在差异，即使设置了相同控制参数，默认的数值计算细节也会影响最终结果。

内容的提问来源于stack exchange，提问作者Ayalew A.