# Modeling continuous abundance data with a GLM in R: how to select the correct distribution family?

I have abundance data (counts) that I have standardized by area sampled, making them continuous. I would like to explain them with my two independent variables using a GLM but I am having trouble specifying a model distribution. The data are derived from raw counts of salamanders at 40 sites standardized by area of pond sampled. The raw counts are as follows:

0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1
1  1  1  1  2  2  2  2  3  3  3  4  5  8 10 11 12 21


And after standardizing the counts by area sampled:

0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
0.0000000  0.0000000  0.1754590  0.4491828  0.8517423  1.4341965
2.2777698  4.0467065  0.3889454  0.5935273  1.9376223  2.0924642
0.5110034  0.5418544  6.9358962  1.5491324  2.2689315 14.2278592
1.4483645  0.3947695  6.2244910  8.7240609


It seems that poisson and negative binomial are now out of the question because my data are not integers. My data has zeros so I don't think I can use gamma distribution without transforming it.

I used the fitdist function in R (package fitdistrplus) to generate parameters for continuous distributions (exponential and normal). I then randomly sampled true exponential and normal distributions with the generated parameters [rexp(n,rate) and rnorm(n,mean,sd)], respectively. Using a two-sided KS test [ks.test(x,y)], I compared my data with the generated data and the distributions were significantly different (p<<<0.05), ruling out my data being normal or exponential. I transformed the data with an ln+1 transformation (to keep zero values) and they still deviate significantly from the hypothesized distributions.

Because my data has no negative values, contains zeros, and is not right and/or left bound, I'm not sure what distribution to specify for the glm.

My questions are:

1. How can I determine the distribution of my response variable, and if I need to transform, can I use a transformation that turns my zeros into non-zeros or even negative values?
2. If the model family depends on the distribution of errors, how can I know their distribution without first removing the potentially explainable variation? It seems like I would need to explain some variation with my predictor variables, and then model the distribution of the remaining error. Is there a better approach to this?

Any suggestions are helpful, specifically those for R or excel. My knowledge of statistical theory is limited.

• Is there a specific reason for why you use fitdistrplus package instead of MASS package (fitdistr())? You don't seem to have censored data. Dec 19, 2014 at 14:03
• Aleksandr, I was unaware of the difference between the two functions but you're right that I don't have censored data. I just used the MASS package and the parameter estimates with fitdistr are equal to those of fitdist for the two distributions I tried.
– kyle
Dec 19, 2014 at 14:57
• That's expected - for non-censored data the results should be the same. If you're not getting good fit, consider performing mixture analysis (of original or transformed data) to see, if you have a mixture distribution (however, I doubt that it could be useful, as your data set is really small). Should you decide to do that, check my different answer. Dec 19, 2014 at 15:08
• 1. It's not strictly necessary for your data to be integer for a poisson-like glm to apply, if you are OK with not having a generative model for you data. See quasi poisson models. 2. If your normalization is just dividing your counts by something akin to an area or an time-interval, you could include an offset in your model which will have the same effect as normalization by division and will keep the counts preserved as counts. Dec 19, 2014 at 15:43
• Thank you Aleksandr and Andrew. I had not thought of these options (like I said, my stats knowledge=weak) and I will explore them asap. I'm also looking into an inverse hyperbolic sine transformation. I will update this post when I find something.
– kyle
Dec 19, 2014 at 15:53

The problem arose from standardizing the raw count values first. The best way to handle this situation is to model the raw counts directly and use an offset in the model to correct for area sampled. With a log link in your GLM model, your offset would be a term offset(log(area)). The coefficient of the offset is forced to be 1. Then you can try Poisson or negative binomial models.