9
$\begingroup$

I have a fairly simple dataset looking at the relationship between the first nesting date of a bird in a given year (Date) and the birds overall fledgling production from that year (Fledge; count data from 0-3 fledglings). I want to determine the optimal laying date for this species (i.e. the date where fledgling production is highest); however, I have been struggling to work out which statistical analysis is most appropriate for my data.

Most birds produce no fledglings in a year, so there are many zero values in the data. With this in mind, I thought that a zero inflated poisson regression might be most appropriate. To test this in R, I fitted a regular glm with poisson distribution (model1 below) and a zero inflated poisson model using zeroinfl() from the pscl library (model2 below). I then compared the two using vuong test statistic (output below).

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

According to the vuong output, a regular non-zero inflated poisson regression is most appropriate. This wasn't that surprising as, from my understanding, the zeros in my data are 'true zeroes' i.e. they are legitimate data points and not due to sampling technique or design. Next I tested for overdispersion by fitting a glm with a quasipoisson distribution to check the dispersion parameter (model3).

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

As you can see, the output showed a dispersion parameter close to one (1.297), and a null deviance (455.36) very close to the residual deviance (443.40). So I concluded that overdispersion was not a big problem, and a poisson, rather than negative binomial, distribution was required.

After doing all this, I was fairly confident that my regular poisson regression (model1) was best for my data. HOWEVER, when I plotted the model outputs of the zero inflated and non-zero inflated data, the non-zero inflated model output didn't appear to fit my data at all while the zero inflated model fitted much better.

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)

enter image description here

As you can see, the non-zero inflated model (red line) doesn't seem to fit the data at all, as it is predicting fledgling production less than 0 on all dates (obviously not biologically possible)! Conversely, the zero inflated model (blue line) seems to fit the data very well. So I have two questions to ask;

  1. Were the methods I used to test zero inflation (and overdispersion) appropriate and interpreted correctly?

  2. If so, is my application of a poisson regression appropriate? Or are there other more appropriate options for my data outside of the two I've tried here?

I've included a dput() of my data below to allow for replication.

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))
Improve this question
$\endgroup$

1 Answer 1

Reset to default
6
$\begingroup$

I suspect that your problem may be that the default behavior of predict.glm isn't what you think it is.

Specifically, predict used on a glm object will by default gives a response on the scale of the linear predictors, not the response.

This is quite clearly stated in the help (?predict.glm) but seems to trip people up very often (suggesting the default ought to be changed, perhaps; you might like to raise it on the relevant mailing list).

To get the values you want, try predict(model1,type="response")

Improve this answer
$\endgroup$
5
  • $\begingroup$ that makes things look much more sensible, thank you! As for my selection process, is my use of a vuong test and calculation of dispersion parameters appropriate when deciding if I should use zero-inflation and/or negative binomial? Would you agree with the use of a poisson regression for this analysis? I don't have much experience analysing such data, so I just want to make sure I haven't made any glaring errors! $\endgroup$
    – Liam B.
    Commented Sep 19, 2014 at 3:56
  • $\begingroup$ Your general approach seems to be broadly-speaking reasonable. If you mean this test, that's probably a reasonable thing to try, but it's probably not what I'd have done. I'd fit the model I felt was most interpretable and check it for reasonableness. If the model proportion of zeros and the residual variance were both reasonably consistent with a Poisson, I'd probably have used that, but another criterion would be what similar studies have used (as long as their model choices seem adequate). You should be a bit cautious about...(ctd) $\endgroup$
    – Glen_b
    Commented Sep 19, 2014 at 6:06
  • $\begingroup$ (ctd)... the consequences of choosing between competing models and doing estimation or testing of parameters on the same data. $\endgroup$
    – Glen_b
    Commented Sep 19, 2014 at 6:19
  • 2
    $\begingroup$ Another possibility when it comes to choosing between NB or ZIP regression models is to compare predicted vs. actual fledgling counts (perhaps with RMSE) using 5- or 10-fold cross validated samples. $\endgroup$
    – RobertF
    Commented Sep 19, 2014 at 19:05
  • $\begingroup$ Thanks for your responses. Makes much more sense to me now. Cheers. $\endgroup$
    – Liam B.
    Commented Sep 21, 2014 at 23:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.