0
$\begingroup$

So my dataset has columns Male, Female, Month and Region (East or West).

And as this is a count data (number of recordings of either male or female cannot be less than 0 for any entry) I am supposed to use poisson.

First with gaussian distribution I get significant p-values about these birds being spotted in certain months.

eb1$bird_cnt <- eb1$Male + eb1$Female

summary(glm(bird_cnt ~ Region, data = eb1[eb1$DATE >= "2005-01-01",]))

Call:
glm(formula = bird_cnt ~ Region + Month, data = eb1[eb1$DATE >= "2005-01-01", ])

Deviance Residuals: 
 Min        1Q    Median        3Q       Max  
-0.40541  -0.24594  -0.16190  -0.01577   1.76026  

Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
(Intercept)     1.01577    0.08501  11.948  < 2e-16 ***
RegionWEST     -0.12613    0.06070  -2.078  0.03855 *  
MonthAugust     0.16911    0.12107   1.397  0.16345    
MonthDecember   0.13641    0.10872   1.255  0.21055    
MonthFebruary   0.25450    0.11875   2.143  0.03288 *  
MonthJanuary    0.22397    0.10932   2.049  0.04132 *  
MonthJuly       0.04157    0.15236   0.273  0.78518    
MonthJune       0.15991    0.15877   1.007  0.31464    
MonthMarch      0.14613    0.11059   1.321  0.18737    
MonthMay        0.01226    0.16225   0.076  0.93980    
MonthNovember   0.27872    0.10377   2.686  0.00763 ** 
MonthOctober    0.35631    0.10715   3.325  0.00099 ***
MonthSeptember  0.38964    0.11940   3.263  0.00122 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.1720771)

    Null deviance: 57.879  on 321  degrees of freedom
Residual deviance: 53.172  on 309  degrees of freedom
AIC: 361.87

Number of Fisher Scoring iterations: 2

I see that residual deviance by degree of freedom shows its underdispersed.

So I try poisson and this is what I get -

Call:
glm(formula = bird_cnt ~ Region + Month, family = "poisson", data = eb1[eb1$DATE >= "2005-01-01", ])

Deviance Residuals: 
 Min        1Q    Median        3Q       Max  
-0.36585  -0.22290  -0.15306  -0.01284   1.33480  

Coefficients:
           Estimate Std. Error z value Pr(>|z|)
(Intercept)     0.01281    0.20472   0.063    0.950
RegionWEST     -0.10741    0.13711  -0.783    0.433
MonthAugust     0.15643    0.28058   0.558    0.577
MonthDecember   0.12780    0.25589   0.499    0.617
MonthFebruary   0.22705    0.27229   0.834    0.404
MonthJanuary    0.20210    0.25377   0.796    0.426
MonthJuly       0.03459    0.36665   0.094    0.925
MonthJune       0.14555    0.36954   0.394    0.694
MonthMarch      0.13644    0.25944   0.526    0.599
MonthMay        0.01008    0.39106   0.026    0.979
MonthNovember   0.24614    0.24109   1.021    0.307
MonthOctober    0.30528    0.24541   1.244    0.214
MonthSeptember  0.33203    0.26834   1.237    0.216

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 41.040  on 321  degrees of freedom
Residual deviance: 37.071  on 309  degrees of freedom
AIC: 746.26

Number of Fisher Scoring iterations: 4

And again I try it with quasipoisson and get significant p-values for certain months but underdispersed.

Call:
glm(formula = bird_cnt ~ Region + Month, family = "quasipoisson", data = eb1[eb1$DATE >= "2005-01-01", ])

Deviance Residuals: 
 Min        1Q    Median        3Q       Max  
-0.36585  -0.22290  -0.15306  -0.01284   1.33480  

Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
(Intercept)     0.01281    0.07582   0.169 0.865926    
RegionWEST     -0.10741    0.05078  -2.115 0.035214 *  
MonthAugust     0.15643    0.10392   1.505 0.133284    
MonthDecember   0.12780    0.09477   1.348 0.178502    
MonthFebruary   0.22705    0.10085   2.251 0.025066 *  
MonthJanuary    0.20210    0.09399   2.150 0.032306 *  
MonthJuly       0.03459    0.13580   0.255 0.799133    
MonthJune       0.14555    0.13687   1.063 0.288410    
MonthMarch      0.13644    0.09609   1.420 0.156630    
MonthMay        0.01008    0.14484   0.070 0.944563    
MonthNovember   0.24614    0.08929   2.757 0.006189 ** 
MonthOctober    0.30528    0.09089   3.359 0.000881 ***
MonthSeptember  0.33203    0.09939   3.341 0.000938 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 0.1371779)

Null deviance: 41.040  on 321  degrees of freedom
Residual deviance: 37.071  on 309  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4

Sorry for the long post.

I am not able to understand the reason behind this. Also I am not able to understand which distribution should I stick to? Gaussian, poisson or quasipoisson?

I am extremely sorry for the long post but I didnt know how to explain things.

Please help.

Thanks.

Additional details -

The codes used here

model2<-(glm(bird_cnt ~ Region + Month , family='quasipoisson', data = eb1[eb1$DATE >= "2005-01-01",]))
summary(model2)

png("myplot.png"); par(mfrow=c(2,2)); plot(model2); dev.off()

PNG of the plot

$\endgroup$

migrated from stackoverflow.com Aug 26 '16 at 11:53

This question came from our site for professional and enthusiast programmers.

3
$\begingroup$

The answer is likely to be quasipoisson.

This will depend a bit on how much data you have. Is it only slightly more than the number of parameters (12)? Assuming you have at least, say, 24 counts:

When you model data with a poisson distribution, you are saying that the variance of that data is equal to its mean. In other words, if you predict a count of 10000, then the variance of that count is 10000 (std.dev 100).

In real life, that isn't always true. Some data have more variance than this, and some less. It looks like your data has less (if we predict a count of 10000, then the variance appears to be more like 1371 rather than 10000).

Your (non-quasi-)poisson model ignores that fact. It is taking the predictive variance to always be equal to the predictive mean even when the data indicates otherwise. This is why it thinks the parameters are insignificant, because it is highly overstating the predictive variance.

If you only have 13-15 rows of data then it might just be that the poisson glm happens to fit very well and the residuals were smaller than expected.

If the counts are reasonably large, the Gaussian distribution is a good approximation. If some counts are quite small (say, less than 25) then it works less well. Bear in mind also that if you use a Gaussian LM, the effects are additive (observing in November = +1000 birds against June, for example) rather than multiplicative (observing in November = x2 birds against June)

$\endgroup$
  • $\begingroup$ Hey there. Thanks for your reply. well I have 387 rows -.-. its really annoying coz I really dont know which would fit better. But quasi poisson is taken when the data is over dispersed correct?? $\endgroup$ – The learning R Kid Aug 26 '16 at 13:01
  • 3
    $\begingroup$ quasipoisson also works when the data is underdispersed. $\endgroup$ – JDL Aug 26 '16 at 14:29
2
$\begingroup$

To add to the other answer (i.e. I agree that family=quasipoisson is likely to be the right choice).

It does (all other things being equal) make sense to use a distribution that naturally applies to count data (Poisson, negative binomial; quasi-Poisson is not a distribution but fits in the same category). However, the most important step here is to check the diagnostic plots: plot(fitted_model,which=3) will show you a scale-location plot which will give you an idea of whether the model is correctly scaling the variance with the mean (this should be an approximately even cloud of points, with a relatively constant smoothed line going through it).

Based on your diagnostic plots above, I'm going to guess that most of your responses are 0 or 1 with only a few observations >1 (this pattern is what gives rise to the odd lines of points in your diagnostic plots). Under these circumstances, it might be best to collapse to the data to (0,>0) and do a binomial/Bernoulli model (e.g. eb1$bird_present <- as.numeric(eb1$bird_cnt>0)) and use family=binomial ...

It is not a good idea to try to pick the best distribution by looking at the parameter summaries, for two reasons: (1) the parameter estimates and p-values don't actually tell you which model fits best; (2) looking at parameter summaries for lots of models will tempt to you to pick the one with the largest number of stars (which may not actually be the appropriate model).

Also, it probably doesn't make too much sense to pick out the months that happen to be significant and say that there's something special about them; rather, test the overall effect of the month variable (e.g. via anova() or car::Anova()), then describe the overall seasonal pattern verbally/graphically.

$\endgroup$
  • $\begingroup$ Greetings Sir, Thank you for your reply. I tried to plot the model and I got a nearly straight line but I dont see random scattered points. I am just too confused. Do I check the distribution plot for all the models? My region is a Binomial variable (East or West) which is why the estimate for region is negative and shows only west. I have read about the dummy variable trap, but is that the only solution to get the value for East? Thank you for your help Sir. $\endgroup$ – The learning R Kid Aug 26 '16 at 17:02
  • $\begingroup$ please use png("myplot.png"); par(mfrow=c(2,2)); plot(fitted_quasipoisson_model); dev.off() and edit your question to post a link to the result. $\endgroup$ – Ben Bolker Aug 26 '16 at 17:06
  • $\begingroup$ Greetings Dr. Bolker, Thank you for you comments. I tried doing it and the results I am getting is bizarre. Which I feel isnt really the case with the data. Well if you dont mind and If you have some spare time to look into my data and R script, I can send it to you by email so that you can spot exactly where I am going wrong? Thanks. $\endgroup$ – The learning R Kid Aug 26 '16 at 20:40
  • $\begingroup$ Sorry, but you're just going to have to struggle to present your problem here; I can't answer questions by e-mail. Read stackoverflow.com/questions/5963269/… and [mwe] to get some ideas about how to post a reproducible example, and compose a new question ... $\endgroup$ – Ben Bolker Aug 26 '16 at 21:11
  • $\begingroup$ Greetings Dr. Bolker, Sorry I am a bit amateur so I didnt understand the part where you said collapse the data to (0,>0). I have been trying to understand this but couldnt. $\endgroup$ – The learning R Kid Aug 28 '16 at 14:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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