# Reporting glmer.nb Results

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


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.
• Maybe I am missing something, but doesn't a mixed effects negative binomial regression model need a response variable which is a count variable (possibly with inclusion of an offset)? The response variable in this model is log density, which is a continuous rather than a count response. Nov 22 '18 at 14:19
• @IsabellaGhement yes, you're right. Nonetheless, AFAIK it is also used with non-integer response variables. An example I know of is RNAseq data, which are normalized counts (i.e., originally you have integer RNA counts, but then they normalize them to account for what is called library size). In standard software for RNAseq (e.g., bioconductor.org/packages/release/bioc/html/edgeR.html) they still use the negative binomial distribution. I thought of this case to be something similar. Nov 22 '18 at 14:24
• Thanks, Dimitris! As always, I'm learning a lot from you! ❤️ Nov 22 '18 at 14:33
• Indeed, also in my GLMMadaptive package, this version is programmed, and hence it also works for the "continualized" negative binomial distribution. Nov 22 '18 at 14:49
• I think the motivation is that is has a specific mean-variance relationship that you will not have if your fit it as a normal distribution. Nov 22 '18 at 15:36