I want to test the fixed and random effects of some covariates on a discrete variable with non negative values. In exploratory analysis I fitted a null Poisson GLM and an null Poisson GLMM. However, the GLMM underestimated the mean value of the response variable even after inclusion of fixed and/or random covariates. I also tried Bayesian approaches, zero-inflated models and negative binomial distributions but the "problem" remains.

Response variable mean: 0.7804
GLM intercept: 0.7803772
GLMM intercept: 0.6595108

Is the estimated intercept of the GLMM an indicative of poor fitting of the model?

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.   
 0.0000  0.0000  0.0000  0.7804  1.0000 12.0000  

mod1 <- glm(caes ~ 1, poisson, banco2)  
Deviance Residuals:  
    Min       1Q   Median       3Q      Max    
-1.2493  -1.2493  -1.2493   0.2381   6.5689  
            Estimate Std. Error z value Pr(>|z|)      
(Intercept) -0.24798    0.01078  -23.01   <2e-16 ***  
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  
(Dispersion parameter for poisson family taken to be 1)  
    Null deviance: 15304  on 11027  degrees of freedom  
Residual deviance: 15304  on 11027  degrees of freedom  
AIC: 27654  
Number of Fisher Scoring iterations: 5  


(mod2 <- lmer(caes ~ 1 + (1 | setor), poisson, data = banco2))  
Generalized linear mixed model fit by the Laplace approximation  
Formula: caes ~ 1 + (1 | setor)  
   Data: banco2  
   AIC   BIC logLik deviance  
 13575 13590  -6785    13571  
Random effects:  
 Groups Name        Variance Std.Dev.  
 setor  (Intercept) 0.39817  0.63101  
Number of obs: 11028, groups: setor, 559  
Fixed effects:  
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -0.41626    0.02937  -14.18   <2e-16 ***  


First I'll outline the logic behind why the exponentiation of the intercept is not sufficient, then I'll present a simple method of calculating the model estimate of the mean of the series.

When calculating the model-based estimate of the overall mean of the series in the random effects case, you have to not only exponentiate the log of the intercept, but take the variability of the observation-specific means into account. The individual observations do not all have mean equal to $\exp\{\text{intercept}\}$, due to that random effect term. See for a not-altogether-dissimilar example the mean of the lognormal distribution, which is not just $\exp\{\text{the mean of the log of the lognormal variate}\}$. The reason for this is that the mean of the transform is not the transform of the mean, unless the transform is linear or you got lucky with the transform, which you haven't with exponentiation.

If the categories are roughly the same size, the variance of the random effects can be used, using the same formula as in the link:

$$\mathbb{E}[\text{caes}] = \exp\{\mu + \sigma^2/2)$$

which translates to $\exp\{-0.41626 + 0.39817/2\} = 0.8048$.

In practice, you can use the predict function to predict the individual values on the response scale, then calculate the mean of the predictions and use that as the basis for your comparison. That's what I would do.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.