Why is Poisson regression different with glmer and gamlss? I have a set of count data that seems to fit "Poisson" = not overdispersed, alpha = 0.
The problem is, I get different results using gamlss vs glmer. Any help explaining the difference would be appreciated:
Treatment is a factor with 6 levels. Trial is a random intercept factor with 3 levels.
Glmer.POI <- glmer( Y ~ Treatment + offset(log(Total)) + (1|Trial),      family=poisson))
Gamlss.PO <- gamlss(Y ~ Treatment + offset(log(Total)), random=~1|Trial, family=PO)

# Glmer result:  AIC = 284. deviance = 270
# Gamlss result: AIC = 388. deviance = 376

I get different coefficients and different p-values with the 2 models. When I look at the diagnostic plots, I see residuals that increase with fitted values in the gamlss model "heteroscedasticity", but not with the glmer.
My intuition tells me this is related to how random effects are handled, but how?
I think I had similar problems with another dataset using negative binomial. I pulled the question from this site to concentrate on this one: If gamlss and glmer behaviour with random effects explains both problems, Great! But I'd like to know how a biologist (non statistician) "poor soul" like me can understand how to use these tools and choose the right one! ;-)
 A: It may be that it has changed since this question was written, but it looks as though the random effect is not coded correctly for gamlss.  You have it written as "random=~1|Trial," but when I try to run that through gamlss it states that the "|" is not valid for factors.  More details on how to code random effects for gamlss is in the manual: http://www.gamlss.org/wp-content/uploads/2013/01/gamlss-manual.pdf
The benefit of gamlss is that you can model different aspects (I use a zero-inflated beta distribution and can model both the distribution of non-zero values and the probability of zeroes).  It could be that, as coded, the random effect is not influencing the same part of the model as it is in glmer.  When I use a ~ for the random effect in my gamlss models, it doesn't contribute to the distribution of non-zero values anymore; it is transferred to the probability of zero.
Unfortunately I have not figured out how to properly code mixed models for gamlss, but hopefully this at least clears up why you're getting different results.
