机器算法验证 - 何时使用零膨胀泊松回归和负二项分布 - 吾爱随笔录

我有一个相当简单的数据集，查看给定年份（日期）中鸟类的第一次筑巢日期与该年鸟类总体雏鸟产量（雏鸟；计数 0-3 雏鸟的数据）之间的关系。我想确定该物种的最佳产蛋日期（即雏鸟产量最高的日期）；然而，我一直在努力找出最适合我的数据的统计分析。

大多数鸟类在一年内不会产生雏鸟，因此数据中有许多零值。考虑到这一点，我认为零膨胀泊松回归可能是最合适的。为了在 R 中测试这一点，我使用 pscl 库中的 zeroinfl() 拟合了一个带有泊松分布的常规 glm（下面的模型 1）和一个零膨胀泊松模型（下面的模型 2）。然后我使用 vuong 测试统计数据（下面的输出）比较了两者。

model1<-glm(Fledge~Date,data=OPT,family="poisson")
model2<-zeroinfl(Fledge~Date,data=OPT,dist="poisson")
vuong(model1,model2)

Vuong Non-Nested Hypothesis Test-Statistic: 4.25169 
(test-statistic is asymptotically distributed N(0,1) under the
null that the models are indistinguishible)
in this case:
model1 > model2, with p-value 1.0608e-05

根据 vuong 输出，常规的非零膨胀泊松回归是最合适的。这并不奇怪，因为据我了解，我的数据中的零是“真零”，即它们是合法的数据点，而不是由于采样技术或设计。接下来，我通过拟合具有准泊松分布的 glm 来测试过度分散，以检查分散参数（模型 3）。

model3<-glm(Fledge~Date,data=OPT,family="quasipoisson")
summary(model3)

Call:
glm(formula = Fledge ~ Date, family = "quasipoisson", data = OPT)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.7877  -0.6258  -0.5578  -0.4504   3.3648  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -0.15109    0.54874  -0.275  0.78315   
Date        -0.03288    0.01133  -2.901  0.00385 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 1.296628)

Null deviance: 455.36  on 642  degrees of freedom
Residual deviance: 443.40  on 641  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 6

如您所见，输出显示的离差参数接近 1 (1.297)，零偏差 (455.36) 非常接近剩余偏差 (443.40)。所以我得出结论，过度分散不是一个大问题，需要泊松分布，而不是负二项分布。

完成所有这些之后，我相当有信心我的常规泊松回归（模型 1）最适合我的数据。但是，当我绘制零膨胀和非零膨胀数据的模型输出时，非零膨胀模型输出似乎根本不适合我的数据，而零膨胀模型拟合得更好。

ggplot(OPT,aes(x=Date,y=Fledge))+
geom_point()+
theme_bw()+
geom_line(data=cbind(OPT,optpred=predict(model1)),aes(y=optpred),size=1,colour="red")+
geom_line(data=cbind(OPT,optpred2=predict(model2)),aes(y=optpred2),size=1,colour="blue)

在此处输入图像描述

如您所见，非零膨胀模型（红线）似乎根本不适合数据，因为它预测所有日期的新兴产量都小于 0（显然在生物学上不可能）！相反，零膨胀模型（蓝线）似乎非常适合数据。所以我有两个问题要问；

我用来测试零通胀（和过度分散）的方法是否恰当且解释正确？
如果是这样，我应用泊松回归是否合适？或者除了我在这里尝试的两个之外，还有其他更适合我的数据的选项吗？

我在下面包含了我的数据的 dput() 以允许复制。

structure(list(Date = c(45L, 40L, 42L, 41L, 39L, 34L, 40L, 44L, 
36L, 32L, 33L, 89L, 58L, 50L, 46L, 56L, 48L, 69L, 64L, 56L, 61L, 
58L, 66L, 63L, 57L, 58L, 60L, 65L, 31L, 48L, 42L, 41L, 46L, 38L, 
59L, 41L, 65L, 41L, 34L, 41L, 36L, 60L, 42L, 39L, 43L, 46L, 47L, 
38L, 38L, 71L, 65L, 51L, 42L, 37L, 51L, 41L, 65L, 59L, 44L, 50L, 
51L, 47L, 40L, 53L, 56L, 62L, 50L, 46L, 51L, 55L, 50L, 46L, 45L, 
39L, 36L, 52L, 50L, 73L, 42L, 38L, 51L, 49L, 43L, 45L, 44L, 76L, 
68L, 65L, 70L, 56L, 40L, 45L, 49L, 52L, 66L, 80L, 45L, 42L, 44L, 
37L, 48L, 43L, 53L, 31L, 47L, 49L, 44L, 46L, 54L, 55L, 48L, 53L, 
55L, 72L, 54L, 45L, 83L, 59L, 48L, 47L, 52L, 72L, 51L, 70L, 48L, 
44L, 42L, 38L, 48L, 43L, 45L, 39L, 45L, 42L, 64L, 46L, 56L, 34L, 
50L, 48L, 47L, 47L, 60L, 50L, 61L, 40L, 72L, 63L, 55L, 66L, 69L, 
66L, 61L, 60L, 60L, 40L, 70L, 45L, 40L, 41L, 42L, 71L, 54L, 45L, 
52L, 48L, 40L, 39L, 49L, 42L, 43L, 53L, 38L, 53L, 52L, 68L, 61L, 
62L, 87L, 41L, 45L, 37L, 44L, 45L, 43L, 72L, 39L, 56L, 34L, 74L, 
62L, 46L, 43L, 47L, 35L, 54L, 61L, 44L, 49L, 54L, 61L, 37L, 51L, 
48L, 52L, 48L, 48L, 44L, 45L, 44L, 45L, 68L, 61L, 87L, 51L, 52L, 
50L, 50L, 56L, 55L, 56L, 57L, 65L, 41L, 63L, 76L, 52L, 62L, 50L, 
50L, 54L, 63L, 48L, 54L, 46L, 57L, 54L, 52L, 45L, 41L, 54L, 74L, 
69L, 68L, 51L, 60L, 54L, 44L, 67L, 52L, 49L, 43L, 41L, 44L, 49L, 
46L, 43L, 46L, 49L, 46L, 47L, 54L, 55L, 67L, 52L, 55L, 52L, 49L, 
50L, 51L, 57L, 48L, 34L, 54L, 49L, 47L, 71L, 62L, 43L, 45L, 45L, 
49L, 58L, 57L, 55L, 54L, 52L, 51L, 41L, 54L, 70L, 52L, 53L, 53L, 
50L, 71L, 56L, 48L, 33L, 43L, 41L, 68L, 42L, 38L, 39L, 46L, 55L, 
64L, 62L, 56L, 69L, 44L, 49L, 54L, 86L, 46L, 46L, 50L, 44L, 45L, 
55L, 55L, 52L, 49L, 49L, 56L, 41L, 34L, 50L, 62L, 39L, 41L, 56L, 
42L, 40L, 43L, 44L, 45L, 43L, 48L, 41L, 45L, 62L, 49L, 47L, 49L, 
63L, 69L, 46L, 53L, 49L, 59L, 54L, 33L, 46L, 44L, 49L, 36L, 41L, 
33L, 66L, 56L, 67L, 43L, 66L, 31L, 51L, 59L, 57L, 51L, 39L, 44L, 
31L, 40L, 39L, 42L, 27L, 43L, 42L, 78L, 60L, 70L, 64L, 67L, 66L, 
67L, 66L, 62L, 58L, 51L, 50L, 60L, 38L, 45L, 34L, 69L, 38L, 45L, 
39L, 44L, 39L, 44L, 43L, 46L, 37L, 59L, 74L, 59L, 39L, 43L, 40L, 
38L, 45L, 45L, 42L, 36L, 33L, 51L, 64L, 52L, 40L, 89L, 49L, 37L, 
51L, 70L, 65L, 71L, 62L, 61L, 68L, 59L, 54L, 75L, 57L, 55L, 58L, 
52L, 58L, 45L, 50L, 41L, 64L, 49L, 50L, 67L, 54L, 43L, 49L, 54L, 
55L, 53L, 53L, 59L, 47L, 47L, 48L, 45L, 50L, 39L, 48L, 51L, 54L, 
44L, 43L, 56L, 51L, 38L, 71L, 62L, 56L, 65L, 69L, 68L, 52L, 47L, 
47L, 47L, 80L, 51L, 48L, 36L, 32L, 39L, 45L, 31L, 43L, 57L, 65L, 
60L, 62L, 36L, 53L, 64L, 57L, 43L, 71L, 66L, 63L, 49L, 39L, 49L, 
43L, 32L, 47L, 44L, 35L, 35L, 41L, 54L, 50L, 44L, 44L, 48L, 50L, 
41L, 40L, 46L, 48L, 38L, 43L, 54L, 52L, 36L, 62L, 72L, 47L, 66L, 
50L, 51L, 50L, 56L, 47L, 67L, 50L, 35L, 40L, 43L, 42L, 31L, 35L, 
43L, 46L, 45L, 46L, 39L, 40L, 40L, 39L, 36L, 45L, 43L, 44L, 44L, 
44L, 38L, 49L, 52L, 49L, 43L, 42L, 47L, 56L, 51L, 51L, 51L, 59L, 
64L, 46L, 40L, 75L, 65L, 51L, 91L, 56L, 83L, 56L, 57L, 58L, 51L, 
50L, 56L, 40L, 69L, 54L, 45L, 35L, 41L, 48L, 60L, 54L, 39L, 39L, 
31L, 92L, 39L, 66L, 56L, 48L, 44L, 40L, 42L, 47L, 51L, 47L, 45L, 
49L, 69L, 48L, 42L, 58L, 56L, 58L, 61L, 42L, 36L, 47L, 52L, 45L, 
54L, 55L, 62L, 48L, 44L, 54L, 51L, 46L, 44L, 50L, 37L, 33L, 40L, 
57L, 54L, 64L, 59L, 69L, 46L, 40L, 51L, 53L, 79L, 60L), Fledge = c(1L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 3L, 0L, 1L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
0L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 2L, 0L, 0L, 0L, 0L, 
0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 1L, 2L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
0L, 1L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 3L, 0L, 0L, 1L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
2L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 1L, 2L, 0L, 0L, 0L, 1L, 1L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 
1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 3L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 2L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 
0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 2L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
2L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 
1L, 0L, 2L, 0L, 0L, 0L, 1L, 0L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 1L, 1L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L)), .Names = c("Date", "Fledge"), class = "data.frame", row.names = c(NA, 
    -643L))