# Null GLMM Poisson underestimate the mean of the response variable. Is it indicative of poor fitting?

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?

summary(banco2$caes) 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 Coefficients: 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 exp(mod1$coefficients[1])
(Intercept)
0.7803772

(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 ***

exp(fixef(mod2))
(Intercept)
0.6595108


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.