Statistical test in R for count data I am struggling with using the correct test for count data in R. The dataset I have is the number of nymphs produced by three aphid species on wheat.
So this is count data with one response variable (number of nymphs) and one explanatory variable with three levels (i.e. three different species).
I want to test whether there are any differences between the three species with the number of nymphs they produce.
I think that I should not use an ANOVA with count data. I have tried a GLM with QuasiPoisson, but I am struggling to interpret the output and be sure that I have used the correct test. Essentially I want to equivalent for a one-way anova but for count data.
 A: I suggest two consideration before modelling.


*

*Are you or not comparing fecundities of species through their lifetime? 
Poisson distribution might not be the case here if you do compare in this way, because species might have different lifespan on the same host plant, Poisson however refer to occurrence of events at the same time interval. Non-parametric comparison can be used here but please try to have enough replicates because generally non-parametric methods are more conservative. 

*If hypothesis on Poisson is met here, it will be better to check whether, overall, the average value of response variable equal to their dispersion. You might consider to use "Negative binomial regression" if your data is over-dispersed. 
A: I would suggest a Kruskal-Wallis test followed by Dunn Test for multiple comparisons:
library(dunn.test) #available on CRAN

## count data
df <- data.frame(group=LETTERS[1:3], count = c(1,3,5,4,4,6,1,3,5,2,3,5,1,5,3,4,2,5))

kruskal.test(count~group, data=df)

#Kruskal-Wallis rank sum test
#
#data:  count by group
#Kruskal-Wallis chi-squared = 8.7041, df = 2, p-value = 0.01288






# non-formula usage (default "holm" method)
dunnTest(df$count,df$group)

#Dunn (1964) Kruskal-Wallis multiple comparison
#  p-values adjusted with the Holm method.
#
#  Comparison         Z     P.unadj      P.adj
#1      A - B -1.131517 0.257837400 0.25783740
#2      A - C -2.925386 0.003440288 0.01032086
#3      B - C -1.793869 0.072834082 0.14566816

A: I like your original approach. There is a worked example here. I have never done Poisson regression, but that seems to me what you have here, where number of nymphs is your Poisson distributed y and species is your factor variable x.
glm(y~x, family="poisson")

You would interpret the output similarly to how you would interpret an ordinary ANOVA output from R.
