2
$\begingroup$

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!

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$
    – whuber ♦
    Commented Oct 21, 2020 at 15:51
  • $\begingroup$ 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? $\endgroup$
    – Frith Nanth
    Commented Oct 23, 2020 at 7:05
  • $\begingroup$ 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? $\endgroup$
    – Frith Nanth
    Commented Oct 23, 2020 at 7:05

1 Answer 1

Reset to default
0
$\begingroup$

If anybody ever encounters this issue again, please check out the solution here: R glm Coefficient Slightly Off Weighted Effect Coding Binomial Logistic Regression

This has to do with nonlinear averaging.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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