# Use of Gamma Distribution for count data

I am working on my data including the insect abundance in dependence of landscape variables with a nested random effect. Since i collected the individuals in the field i have count data and thought a poisson distribution would fit best (using lme4). Currently my model looks like this example glmer(insects~landscape1*landscape2 +(1|region/location), family="poisson", data=..). (no overdisperion)

Following "A practical guide to mixed models in R" http://ase.tufts.edu/gsc/gradresources/guidetomixedmodelsinr/mixed%20model%20guide.html i "should use" a lognormal or gamma distribution since they fit best. I tried lme with log-transformed response and glmer with gamma even if I have no continous data and both show similiar results in contrast to the glmer with poisson distribution. I also tried using ord_plot() and distplot() (pos. slope, negative intercept) which showed that log.series would be the best choice again (according to Friendly http://www.datavis.ca/courses/grcat/grcat.pdf chapter 2.3.).

I don't want to use the log-transformed approach but was wondering if I could also use gamma for discrete count data or only for continuous? Or do you can suggest any alternatives not using poisson for count data? Hope to get some new insights. Thank you

Here are the two outputs of the methods i was following to decide for the right distribution. There is one count with 0. I am not that familiar with statistics yet but i read a lot that log-transformed data are difficult to interpret and that it is better to go for another method if the data are not normally distributed. So, i thought there should be another way.

• It might help those who would like to try to answer if you could show results that document why certain models "fit best" and the differences between other results and the glmer/poisson model. Also, why don't you want to use the log-transformed approach; are there cases where there are 0 insects? – EdM Feb 3 '17 at 14:51

The negative binomial distribution can be considered a discrete analog to the gamma distribution. Thus, if a gamma gives a good fit but you're concerned about the continuous approximation, then fitting for a negative binomial response would fix that problem and maintain a good fit. Unlike the Poisson distribution, it allows for overdispersion.

I think that you have been misled by "A practical guide to mixed models in R". The tests for normal, Poisson, gamma, etc. distributions proposed on that site look at the distribution of the response variable without reference to the values of the predictors. That's not helpful, and if that's what you did the results are not pointing you in the correct direction.

For example, in a simple linear regression with a true linear relation between predictor and response, the response variable will tend to follow the distribution of the predictor variable. If the predictor variable doesn't have a normal distribution there is no reason to suspect that the response variable has a normal distribution. The errors between predicted and observed response-variable values ideally are normal, but little is to be gained by looking at the distribution of response-variable values while ignoring the corresponding predictors.

This becomes even more of an issue when you have multiple predictors and interactions, as in your case. There is no reason to suspect that your response variable values overall show any single distribution; the tests you performed based on that website's suggestions thus may be leading you astray. Your responses may actually be Poisson-distributed, but with the mean value of the Poisson distribution depending on the values of the predictors. As I read that website and look at your plots, it seems that the values of the predictors are nowhere considered. The points with highest quantiles on your Poisson plot may well be those whose predictor variables predict higher insect numbers.

Ignore those plots, and instead fit a model that shows well-behaved residuals. A generalized linear model based on the Poisson family would be a good place to start; that will allow for Poisson-based responses whose mean values are related to the predictors. As @Björn suggested in another answer, related approaches might be useful if Poisson in its simplest form fails.

Finally, I don't at all agree that log-transformed variables are hard to interpret. Changes in log-scale units are simply fractional or percentage changes in the original units. Sometimes percentage changes are even easier to think about.

• Thank you very much! Poisson was my first try fitting the model but i was really confused by gamma, lognormal,.. So i will just go back to poisson and probably try Björns suggested approaches! – svenja Feb 6 '17 at 8:37

Technically the gamma distribution (and log-normal distribution to name another similar one) is just suitable for continuous data. If counts are really high though (let's say all counts >50 or so), this will in practice not matter so much. However, in that case the Poisson distribution should also look very similar. I am not sure why you would want to use some sensible discrete distribution for count data, there are plenty of ways of extending the Poisson distribution, if you need more flexibility. E.g. negative binomial (if there is more variability between units than the Poisson would suggest, in case of the negative binomial distribution this is assumed to vary according to a gamma distribution across units) or zero inflated version of these two.

In fact, one thing that may or may not make sense in your case is to fit a hierarchical model with e.g. a Poisson likelihood and random effects (e.g. normal distributed on the log-mean rate scale or gamma distributed on the mean rate scale) for region and location nested within region. This would assume that locations within a region are sort of similar in terms of the expected mean counts once you have taken model covariates (such as landscape) into account, and that regions are sort of similar to each other (again, once you have taken model covariates into account). "Sort of similar" here means that distribution parameters are exchangeable, i.e. after accounting for your model covariates, you would not know which one you would expect to have higher or lower counts. The main effect of doing this is that if your regions and locations are indeed somewhat similar, you end up getting some borrowing of strength across them. I.e. estimation becomes more efficient, because you exploit that you do not have estimate completely separate parameters for all locations (but also relies more on the assumption of exachangeability).

• Thanks for your answer. I collected the insects in ten different regions which means they differ in their landscape variables on a certain scale (e.g. % of agricultural area). All insects where collected on the same field type at the edge and interior (location). I thought i have to account for this by using the random effect anyway. Isn't that right? – svenja Feb 6 '17 at 8:45
• Accounting for region could be via a fixed (=completely separate parameter for each region) or a random effect (=parameters are pulled towards each other according to the random effects distribution). The second option technically requires that you cannot predict, which of the regions would have a higher insect count - at least not beyond the variables in the model (this assumption is called "exchangeability"). E.g. if a higher % agricultural area would likely result in fewer (or more) insects, then it needs to be in the model (and anything else that would predict the count). – Björn Feb 6 '17 at 11:23
• Now i am bit confused whether i need to include region as random. Region includes different % of landuse types (agriculture, forest,..). But in the models I include only one type and the interaction with local flower cover respectively: e.g. response~agricultural area* flower cover crop. I was always including region as random but is this really necessary if agric. area is already part of it and has the same value for the whole region? At first i also wanted to include the effect of different landuse types in one model but i have always error warnings (in glmer, model failed to convergence) – svenja Feb 6 '17 at 12:59
• I compared the model with and without the random (region) and Pr(>Chisq) is significant (glm.nb and glmer.nb). That means that region as random has an influence and should stay in the model right? – svenja Feb 6 '17 at 13:41