将概率分布拟合到 R 中的零膨胀数据

机器算法验证 r 分布 可能性 零通胀
2022-04-12 18:41:28

我正在尝试学习如何使用程序 R 将概率分布拟合到数据向量，但是有很多潜在的概率分布可供使用！所以我的问题是，如何为我的数据找到最佳分布，以及如何证明我选择了正确的分布？我可以获取一整套不同分布的 AIC 值吗？

这些数据是蜜蜂访问花朵的观察计数数据。每个物种都有一定数量的访问，因此频率不同。目标是找到描述蜜蜂访问的最佳分布，表明我选择了正确的分布，然后使用该分布随机抽样进行一组模拟。

这是数据的样子，它是计数观察的向量。它是零膨胀的，具有长尾分布（可能是零膨胀的负二项式？）。

i.vec=c(0,63,1,4,1,44,2,2,1,0,1,0,0,0,0,1,0,0,3,0,0,2,0,0,0,0,0,2,0,0,0,0,
0,0,0,0,0,0,0,0,6,1,11,1,1,0,0,0,2)

这是我计算的一些基本参数。我使用 sigma 的标准差，phi 是数据中零的比例。

m=mean(i.vec)
#[1] 3.040816
sig=sd(i.vec)
#[1] 10.86078
tab<-table(i.vec)
zero.prop<-as.numeric(tab[1])/sum(as.numeric(tab))
#[1] 0.6122449

如您所见，标准差远大于平均值，而且我的零比例非常高。

3个回答

您可以使用pscl 包中的Vuong测试来比较非嵌套模型。这是一个例子

> m1 <- zeroinfl(i.vec ~ 1 | 1, dist = "negbin")
> summary(m1)

Call:
zeroinfl(formula = i.vec ~ 1 | 1, dist = "negbin")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-0.3730 -0.3730 -0.3730 -0.2503  7.3544 

Count model coefficients (negbin with log link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.1122     0.3831   2.903  0.00369 ** 
Log(theta)   -1.9256     0.2839  -6.784 1.17e-11 ***

Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   -9.815     96.462  -0.102    0.919
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta = 0.1458 
Number of iterations in BFGS optimization: 579 
Log-likelihood: -80.51 on 3 Df


> m2 <- zeroinfl(i.vec ~ 1 | 1, dist = "poisson")
> summary(m2)

Call:
zeroinfl(formula = i.vec ~ 1 | 1, dist = "poisson")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-0.7242 -0.7242 -0.7242 -0.4860 14.2795 

Count model coefficients (poisson with log link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.05911    0.08205    25.1   <2e-16 ***

Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.4561     0.2933   1.555     0.12
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Number of iterations in BFGS optimization: 11 
Log-likelihood: -233.7 on 2 Df


> vuong(m1, m2)
Vuong Non-Nested Hypothesis Test-Statistic: 1.946095 
(test-statistic is asymptotically distributed N(0,1) under the
 null that the models are indistinguishible)
in this case:
model1 > model2, with p-value 0.02582165

Vuong 检验还表明，与普通的负二项式相比，零膨胀负二项式更适合您的数据（此处未显示，但您可以拟合两个模型并进行比较）。

我认为您不一定需要夸大零。

您的数据似乎与负二项式非常一致：

> library(MASS)
> table(rnegbin(49,mu=3.1,theta=0.075))

 0  1  2  3  5 18 20 21 31 61 
36  4  2  1  1  1  1  1  1  1 
> table(i.vec)
i.vec
 0  1  2  3  4  6 11 44 63 
30  8  5  1  1  1  1  1  1

所以让我们得到一些估计：

> mean(i.vec)
[1] 3.040816
> theta.ml(i.vec,3.041)
[1] 0.145777
attr(,"SE")
[1] 0.04136887

所以让我们看一下 mu = 3.041 和 theta 在 0.14 （这里的 theta 有很多不确定性）：

这是来自该分布的三个随机样本：

> table(rnegbin(49,mu=3.041,theta=0.14))

 0  1  2  3  4  5  6  7  8 18 49 
30  5  1  3  2  2  1  2  1  1  1 
> table(rnegbin(49,mu=3.041,theta=0.14))

 0  1  2  3  4  7  9 15 29 31 47 56 
33  2  4  1  1  1  1  2  1  1  1  1 
> table(rnegbin(49,mu=3.041,theta=0.14))

 0  1  2  3  4  7  9 12 16 48 66 
36  4  1  1  1  1  1  1  1  1  1

它们看起来与您的数据足够相似。负二项式似乎至少是合理的。

您可能喜欢玩 MASS 中的功能

我不确定在这种情况下，如果没有关于您的数据的更多信息（尤其是因为您的观察结果很少），那么在这种情况下，您是否可以做得比仅插入经验测量更好。在这种情况下，你的误差方差应该是你观察次数的倒数（通过 Efron-Stein）。

也许您可以使用一些基于卷积的估计器（例如 R 中的密度函数，但使用支持整数的内核）来稍微平滑一些事情。或将事物视为混合物。但是，如果您不知道数据的来源，则没有理由这样做。

其它你可能感兴趣的问题
上一篇是一次输入所有预测变量还是执行分层回归？ 下一篇lme4中的不同组系数