# Count Data - Gaussian, poisson or quasipoisson?

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 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? Please help. 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()


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)

• 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?? – The learning R Kid Aug 26 '16 at 13:01
• quasipoisson also works when the data is underdispersed. – JDL Aug 26 '16 at 14:29

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.

• 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. – The learning R Kid Aug 26 '16 at 17:02
• 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. – Ben Bolker Aug 26 '16 at 17:06
• 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. – The learning R Kid Aug 26 '16 at 20:40
• 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 ... – Ben Bolker Aug 26 '16 at 21:11
• 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. – The learning R Kid Aug 28 '16 at 14:44