# Intercept of negative binomial Regression with weighted effect coded variables does not equal sample mean

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!

• In general, it is never the case that the intercept of any GLM (mixed or not) must equal the grand mean of the responses (especially when using a non-linear link). That result generally holds only for ordinary least squares regression.
– whuber
Commented Oct 21, 2020 at 15:51
• Thanks for your fast reply, @whuber! I am a bit confused by your answer because Agresti (2013, p.552 ff.) states that the grand mean of nbinom-regression always equals the sample mean (his formulation sounds like a theorem). Although the example he is providing has only one dummy coded predictor. So I thought, if the coefficient for the dummy variable there expresses the difference between the two categories of the predictor so that the intercept is the mean of one of the dummy category this must work for weighted effect coding with the sample grand mean, too. Or is this naive? Commented Oct 23, 2020 at 7:05
• Further I found out that if I remove values with residual outliers the mean of the sample and the intercept are closer to each other. But why is this not the case for ordinary least squares regression? Commented Oct 23, 2020 at 7:05