I am facing a really puzzling issue. I am analysing data of employees' absentism because of illness during four months with a negative binomial GLMM with the employees' department as a random effect. My predictors are sex, age, income group (factor with four levels), status group (civil servants vs. employees) and the position (superior vs. non-superior). I want to present the results to people who are absolutely unfamiliar with multivariate statistics and therefore estimate regression coefficients for all levels of the factor variables so I decided to use weighted effect coding rather than dummy coding.

The grand mean of the sample is

> mean(df$WorkingDaysIll)
[1] 9.567867 

The estimates of the GLMM negative binomial regression with weighted effect coded variables are

> Model1.6 <- glmer.nb(WorkingDaysIll ~ 1 + sex + age_scaled + incomegroup + statusgroup + position + 
(1 | department), df)      
> summary(Model1.6)$coefficients
                              Estimate 
(Intercept)                  2.1072362 
sexF                         0.1342055 
age_scaled                   0.2605714 
incomegroup0101              0.6820038 
incomegroup0102              0.1950740 
incomegroup0201             -0.1159360 
statusgroup                  0.3243721 
superiorNsup                 0.1376719 

So the intercept of the model is 2.1072, but

> exp(summary(Modell1.6)$coefficients[1,1])
[1] 8.225477

8.2255 days is relatively far away from the sample grand mean of 9.5679 days.

The mean of the predicted values is relatively close to the grand mean of the sample:

> mean(fitted(Model1.6))
[1] 9.511731

If I only include sex as a predictor the intercept is closer to the grand mean:

> summary(Modell1.1)$coefficients[,1]
   (Intercept)           SexF
     2.2333826      0.1790173 
> exp(summary(Modell1.1)$coefficients[,1])
   (Intercept)           SexF
      9.331377       1.196041 

The more predictors I include, the smaller the intercept and the bigger the difference to the observed grand mean becomes. However if I run Model1.6 as an ordinary regression I get

> Model1.6lmer <- lmer(WorkingDaysIll ~ 1 + sex + age_scaled + incomegroup + statusgroup + position + 
(1 | department), df) 
> summary(Modell1.6lmer)$coefficients
                            Estimate 
(Intercept)                   9.5722    
sexF                          1.2784     
age_scaled                    2.1945     
incomegroup0101               7.6714    
incomegroup0102               1.8612     
incomegroup0201              -1.5037     
statusgroup                   3.2919     
superiorNsup                  0.9916    

The intercept is very close to the sample grand mean of 9.5679 what seems ok to me and could differ slightly because the relationship between age and WorkingDaysIll could be a bit too far away from being linear. It is interesting that the ordinary regression gets a more accurate estimate of the intercept than the negative binomial although the latter describes the data better (the dependent variable is a negative binomial distributed count variable and the deviance is smaller 2268.3 against 2911.9 if I run Model1.6lmer with log-likelihood optimitisation). The results are equal even if I use other functions for negative binomial such as glmmTMB or mixed_model from package GLMMadaptive.

How do I have to interpret the difference between regression intercept and sample grand mean? Is weighted effect coding not reliable for negative binomial regressions (although I do not know any reason why this could be the case), is that a normal behaviour, is the model against my supposition not well specified or do I oversee anything?

Thanks for help!