Solved – Why is Poisson regression different with glmer and gamlss

count-datagamlssglmmpoisson distributionrandom-effects-model

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! 😉

Best Answer

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.

Related Solutions

Solved – Help with zero-inflated generalized linear mixed models with random factor in R

No, zeroinfl() currently does not support random effects. So the formula you specified actually means something different: You use a fixed treatment effect in the count part and a fixed site effect in the zero-inflation part. See vignette("countreg", package = "pscl") for more details.
If you want random effects, then no. If you use fixed interaction effects instead, you could still try to find a suitable model with zeroinfl(). But with your number of observations this is probably not the best solution.
As the model is not the one you would want to fit, this is not relevant here.
For zeroinfl() there would be and I suppose that for glmmADMB there are as well. But I'm not an expert on that.
You could employ effect plots for the covariate effects or rootograms for the goodness of fit. It depends on what you really want to show.

Solved – Binomial glmer with data between 0-1, not count data, not normal proportion

In my opinion, you are modeling the establishment limitation, you just made a mistake with the weights argument. From R help:

For the binomial and quasibinomial families the response can be specified in one of three ways:

As a factor: ‘success’ is interpreted as the factor not having the first level (and hence usually of having the second level).

As a numerical vector with values between 0 and 1, interpreted as the proportion of successful cases (with the total number of cases given by the weights).

As a two-column integer matrix: the first column gives the number of successes and the second the number of failures.

Your data clearly fits in the second case, but you have to supply the total number of cases to weights.

The following code works for me:

model1 <- glmer(Est.Limit ~ Treatment + log(Size) + (1|Species),
                data    = Limitdr.Est,
                family  = "binomial",
                weights = Seed.sites)

the result:

> summary(model1)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: Est.Limit ~ Treatment + log(Size) + (1 | Species)
   Data: x
Weights: Seed.sites

     AIC      BIC   logLik deviance df.resid 
   110.3    117.8    -50.2    100.3       28 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7705 -0.4858 -0.1659  0.6247  2.0272 

Random effects:
 Groups  Name        Variance Std.Dev.
 Species (Intercept) 4.519    2.126   
Number of obs: 33, groups:  Species, 11

Fixed effects:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)           0.9418     0.8370   1.125 0.260523    
TreatmentIsland       0.5172     0.3926   1.317 0.187700    
TreatmentPlantation   1.6060     0.4670   3.439 0.000583 ***
log(Size)            -1.0333     0.3793  -2.724 0.006442 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) TrtmnI TrtmnP
TrtmntIslnd -0.230              
TrtmntPlntt -0.189  0.466       
log(Size)   -0.479 -0.041 -0.100

Hope this helps, Cheers!

Best Answer

Related Solutions

Solved – Help with zero-inflated generalized linear mixed models with random factor in R

Solved – Binomial glmer with data between 0-1, not count data, not normal proportion

Related Question