【问题标题】:formula or matrix underlying count predictions using predprob in pscl在 pscl 中使用 predprob 的公式或矩阵基础计数预测
【发布时间】:2015-12-29 08:21:00
【问题描述】:

R 中的 pscl 包用于拟合零膨胀泊松模型等模型。

我有一个关于每周食物消费量的数据框,它是零膨胀的。我在拟合零膨胀泊松 (ZIP) 模型时没有任何问题,该模型成功地预测了我的数据中零计数的百分比。数据包含大约 30% 的零计数,只有 4600 多个观测值。

ZIP 模型是:

 zeroinfl(LightC.home~SEX+AGEDET+Inhabitants, data=FoodAnalysis)

其中LightC 是每周消费计数(范围 0 - 74),SEX 是受访者性别(2 个级别的因子),AGEDET 是受访者年龄范围(12 个级别的因子,例如“5-9 岁"),Inhabitants 是家庭人数(整数,范围 1 到 8)。

ZIP 模型总结为:

Call:
zeroinfl(formula = LightC.home ~ SEX + AGEDET + Inhabitants, data = FoodAnalysis)

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-3.3479 -1.2781 -0.5170  0.6342 12.2912 

Count model coefficients (poisson with log link):
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)             2.133213   0.029229  72.983  < 2e-16 ***
SEXFemale              -0.079156   0.010284  -7.697 1.39e-14 ***
AGEDET5-9 years        -0.038153   0.032745  -1.165 0.243963    
AGEDET10-14 years       0.014199   0.032411   0.438 0.661317    
AGEDET15-17 years       0.253942   0.035419   7.170 7.52e-13 ***
AGEDET18-24 years       0.148395   0.029089   5.101 3.37e-07 ***
AGEDET25-34 years       0.158506   0.026291   6.029 1.65e-09 ***
AGEDET35-44 years       0.157821   0.026043   6.060 1.36e-09 ***
AGEDET45-54 years       0.307299   0.026153  11.750  < 2e-16 ***
AGEDET55-64 years       0.340590   0.026913  12.655  < 2e-16 ***
AGEDET65-74 years       0.361976   0.027260  13.278  < 2e-16 ***
AGEDET75-84 years       0.251614   0.039866   6.311 2.76e-10 ***
AGEDET85 years or more  0.606829   0.083345   7.281 3.32e-13 ***
Inhabitants             0.014697   0.004302   3.416 0.000635 ***

Zero-inflation model coefficients (binomial with logit link):
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)             -0.82861    0.17389  -4.765 1.89e-06 ***
SEXFemale               -0.08333    0.07196  -1.158 0.246893    
AGEDET5-9 years         -0.27280    0.18082  -1.509 0.131390    
AGEDET10-14 years       -0.32820    0.18390  -1.785 0.074308 .  
AGEDET15-17 years       -0.06472    0.20589  -0.314 0.753276    
AGEDET18-24 years       -0.06973    0.16047  -0.435 0.663896    
AGEDET25-34 years       -0.27054    0.14791  -1.829 0.067385 .  
AGEDET35-44 years       -0.71412    0.15648  -4.564 5.03e-06 ***
AGEDET45-54 years       -0.50510    0.15585  -3.241 0.001191 ** 
AGEDET55-64 years       -0.65281    0.16541  -3.947 7.92e-05 ***
AGEDET65-74 years       -1.09276    0.18313  -5.967 2.41e-09 ***
AGEDET75-84 years       -1.28092    0.34508  -3.712 0.000206 ***
AGEDET85 years or more -12.75039  279.21031  -0.046 0.963577    
Inhabitants              0.01733    0.02872   0.603 0.546296    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Number of iterations in BFGS optimization: 35 
Log-likelihood: -2.136e+04 on 28 Df

我使用pscl 包中的predprob 命令来构造食物的消费概率向量。这构建了一个矩阵,每个观察一行,每个观察给出的每个计数的消耗概率在 0(最小值)和 74(最大值)之间。此矩阵中有 165 个唯一行,这是由于零膨胀泊松模型中的 165 个唯一预测变量组合。

如何为每种预测变量组合(例如女性、25-34 岁、4 个居民)构建支持这 75 个概率的公式?从数学的角度,我不明白这 75 个概率之间的关系。

例如,预测变量组合 SEX=="Female"AGEDET=="25-34 years"Inhabitants == 4 为数据中具有该预测变量组合的每个人提供了这些概率(仅显示了 75 个中的前 6 个):

        X0            X1           X2           X3          X4          X5
 0.2473285  0.0004504804  0.002183136  0.007053337  0.01709108  0.03313101

考虑到原始zeroinfl 公式看起来只是在估计预测变量组合而不是所有 75 个概率的消耗平均值,这些概率是如何估计的?

【问题讨论】:

    标签: r poisson


    【解决方案1】:

    零膨胀泊松模型是零点质量和泊松分布的混合。具体公式请参见vignette("countreg", package = "pscl") 的第 2.3 节。

    在您的模型中,有两个线性预测变量:一个用于 Poisson 分量的平均值(带有对数链接),一个用于零通胀概率(带有 logit 链接)。使用您问题中指出的特定回归向量,我们可以计算两个线性预测变量,然后应用各自的反向链接函数来获得平均值 mu 和零通胀概率 pr

    mu <- exp(2.133213 -0.079156 + 0.158506 + 0.014697 * 4)
    mu
    ## [1] 9.692487
    pr <- plogis(-0.82861 -0.08333 -0.27054 + 0.01733 * 4)
    pr
    ## [1] 0.2472822
    

    然后我们可以将这些输入到混合密度中(参见小插图中的公式 7:

    y <- 0:5
    pr * I(y == 0) + (1-pr) * dpois(y, mu)
    ## [1] 0.2473287134 0.0004504784 0.0021831278 0.0070533124 0.0170910338
    ## [6] 0.0331309230
    

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2018-01-03
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多