Reporting glmer.nb Results

Question

I'm running a mixed negative binomial GLM that looks like this:

Niche2 <- glmer.nb(log_density ~ height * factor(Year) + (1 | Grouping), data = NicheData2)

To see if the way sward height determines the density (log transformed for normality) of a herbivorous insect has changed over time (so I'm particularly interested in the interactions between years and height). The random effect of grouping is of the different sites sampled in different years (to account for temporal pseudoreplication).

> summary(Niche2)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(225251.6)  ( log )
Formula: log_density ~ height * factor(Year) + (1 | Grouping)
 Data: NicheData2

     AIC      BIC   logLik deviance df.resid 
   341.4    364.3   -162.7    325.4      122 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.2436 -0.4351 -0.1003  0.3726  2.1061 

Random effects:
 Groups   Name        Variance  Std.Dev. 
 Grouping (Intercept) 1.591e-10 1.261e-05
Number of obs: 130, groups:  Grouping, 45

Fixed effects:
                        Estimate Std. Error z value Pr(>|z|)  
(Intercept)              0.65834    0.27708   2.376   0.0175 *
height                  -0.08548    0.04371  -1.956   0.0505 .
factor(Year)2010         0.17000    0.43243   0.393   0.6942  
factor(Year)2018        -0.40936    0.46493  -0.880   0.3786  
height:factor(Year)2010  0.03534    0.05402   0.654   0.5130  
height:factor(Year)2018  0.08860    0.05360   1.653   0.0983 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) height f(Y)2010 f(Y)2018 h:(Y)2010
height      -0.873                                   
fctr(Y)2010 -0.650  0.567                            
fctr(Y)2018 -0.555  0.485  0.361                     
hgh:(Y)2010  0.713 -0.814 -0.871   -0.396            
hgh:(Y)2018  0.684 -0.791 -0.444   -0.853    0.644

I'd like to know the correct way to report these results, (presumably the coefficients, the standard errors, and the p values, while also explaining how much variation is fuelled by the random effect), but I understand that the coefficients need transforming.

Could someone please advise me on the right way to transform these, and how to report the variation due to the random effect?

I've had a look at other questions particularly this one (How to report negative binomial regression results from R) but haven't managed to apply their answers to my model.

EDIT: Following advice, re-ran the model with length as an offset and the response variable as adjusted count data, am keen to know the best way to report the co-efficients (and does the offset change the way I should do this).

Generalized linear mixed model fit by maximum likelihood (Laplace 
Approximation) ['glmerMod']
 Family: Negative Binomial(1.7478)  ( log )
Formula: adjusted ~ height * factor(Year) + (1 | Grouping) + 
offset(log(length))
   Data: NicheData2

     AIC      BIC   logLik deviance df.resid 
  1105.8   1128.2   -544.9   1089.8      114 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.2672 -0.6135 -0.2009  0.5125  3.6734 

Random effects:
 Groups   Name        Variance Std.Dev.
 Grouping (Intercept) 0.553    0.7436  
Number of obs: 122, groups:  Grouping, 39

Fixed effects:
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)             -2.77889    0.33069  -8.403   <2e-16 ***
height                  -0.09984    0.04466  -2.235   0.0254 *  
factor(Year)2010         0.11267    0.45323   0.249   0.8037    
factor(Year)2018        -0.47509    0.51415  -0.924   0.3555    
height:factor(Year)2010  0.03472    0.05173   0.671   0.5021    
height:factor(Year)2018  0.08253    0.05486   1.504   0.1325    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Correlation of Fixed Effects:
             (Intr) height f(Y)2010 f(Y)2018 h:(Y)2010
 height      -0.827                                   
 fctr(Y)2010 -0.520  0.458                            
 fctr(Y)2018 -0.446  0.397  0.388                     
 hgh:(Y)2010  0.625 -0.749 -0.862   -0.400            
 hgh:(Y)2018  0.589 -0.715 -0.413   -0.866    0.613

Answer 1

A couple of comments:

• The estimated variance of the random effect is extremely low. This could indicate that you do not need to include the random effect (i.e., there are correlations in the measurements within each group) or could be an artifact of the estimation procedure. In particular, glmer.nb() fits the negative binomial mixed model using the Laplace approximation, which is known not to be optimal. As an alternative, you can try fitting the same model using the GLMMadaptive package, which uses the adaptive Gaussian quadrature rule; for example, check here.
• In mixed models with nonlinear link functions, as in your case with the negative binomial mixed model, you have to be aware that the estimated coefficients you obtain have an interpretation conditional on the random effects. Most often this is not the interpretation you are looking for, but instead, the interest most often lies in a marginal / population interpretation. For more on this, you may check the discussion in this question. When you fit the model with GLMMadaptive you can obtain the coefficients with a marginal interpretation using the function marginal_coefs(); for example, check here.