How to report negative binomial regression results from R

I have fit a negative binomial model in R, and would like to report the findings, but I'm unsure how (or if) I should convert the estimates to reportable coefficients. Here is my output:

> summary(negbinom)

Call:
glm.nb(formula = Frequency ~ Groups, data = Data, init.theta = 4.577879741,

Deviance Residuals:
Min       1Q   Median       3Q      Max
-2.3978  -0.8346   0.0000   0.4622   2.4414

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  1.32176    0.14217   9.297   <2e-16 ***
Group1        -0.47446    0.20440  -2.321   0.0203 *
Group2        -0.01117    0.20137  -0.055   0.9558
Group3         0.06454    0.19221   0.336   0.7370
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(4.5779) family taken to be 1)

Null deviance: 139.05  on 105  degrees of freedom
Residual deviance: 129.82  on 102  degrees of freedom
AIC: 478.1

Number of Fisher Scoring iterations: 1

Theta:  4.58
Std. Err.:  1.70

2 x log-likelihood:  -468.104

My questions are:

1) is it okay to simply report that Group 1 reported stuff I'm looking at with a lower frequency than Group 0 (the reference group), B = - .475, z = -2.32, p = .020? Or do I need to do something (exp() maybe?) with the coefficients to be able to report? Even then is this the right notation to use (B, z, etc.)?

2) can I use Anova() or anova() to get the omnibus effect of Groups? e.g., can I do this:

> Anova(negbinom)
Analysis of Deviance Table (Type II tests)

Response: Frequency
LR Chisq Df Pr(>Chisq)
Groups   9.2233  3    0.02646 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

...and then report it as: ChiSq(3) = 9.22, p = .026?

• Who is your audience for the report? – James Sep 6 '18 at 10:13
• Cognitive psychologists, mostly. – MGy Sep 6 '18 at 13:59

It is typical to report the back-transformed rate ratios (i.e. exp(coefficient) is a ratio of the rate in group X versus group 0). I typical set of information to report would be rate ratio, confidence interval and p-value (with any intended adjustments for multiple testing, if needed - what you see is not taking into account that you are doing multiple comparisons).

In principle you can do this kind of test for the overall effect of groups, but I do not know whether the particular R commands do this correctly.

Also, be aware that the standard errors, p-values etc. reported by glm.nb are generally too liberal, see this question on how to fix that.

• So I can report exp(coeff) - e.g., exp(-.47446) in the above example -, and call it a "rate ratio"? The z-value and all that other stuff is not typically necessary? Thank you for your final point too - once I use that function to get a vcov matrix, how do I get the more conservative p-values? – MGy Sep 6 '18 at 13:59
• You get the standard error as the square root of the appropriate diagonal entry of the matrix this produces. Then you get the standard normal percentile corresponding to estimate / SE to get the one sided o-value. – Björn Sep 6 '18 at 17:06
• Interestingly, when I look at the matrix with vcov(negbinom) and with glm.nb.cov(negbinom) - the suggested function - the two return the exact same values (except for the additional row and column in glm.nb.cov). I might be misunderstanding something, but this is surprising to me. – MGy Sep 8 '18 at 17:34

It is typical to report a percentage reduction:

round((negbinom$family$linkinv(-.47466)-1)*100)
 -38

So you would say that the estimate for group 1 is 38% smaller (p=0.02) than the reference group. Also mention that the other groups were not significantly different from the reference group (assuming you are happy with the p-values).

Note that I used the linkinv function that comes with the model fit. This is a safer way to back transform from glm fits as you can't accidentally use the wrong function if you change your link function.

As an aside, the effects package has some useful functions for presenting GLMs, and takes away much of the pain in back-transforming to get estimates.

# some synthetic data
x <- data.frame(Frequency=c(rnbinom(30,100,0.1),rnbinom(30,100,0.7),rnbinom(30,100,0.1),rnbinom(30,100,0.5)),Groups=rep(paste0("Group",0:3),each=30))

# model

library(effects)
summary(allEffects(model1))
model: Frequency ~ Groups

Groups effect
Groups
Group0    Group1    Group2    Group3
895.00000  42.20000 908.33333  98.46667

Lower 95 Percent Confidence Limits
Groups
Group0    Group1    Group2    Group3
861.69018  39.54580 874.55058  93.61605

Upper 95 Percent Confidence Limits
Groups
Group0    Group1    Group2    Group3
929.59746  45.03234 943.42108 103.56861

The exp(coef) is conventionally referred to as the Incidence Rate Ratio (IRR). It is easier to interpret because the size of the effect is expressed as a percentage. A good reference is Hilbe (2011).

The Anova() function in the car package produces p-values from Type II LR tests, but there is a caveat: theta is assumed to be fixed.

Here is an example:

library(MASS)
quine.mod1 <- glm.nb(Days ~ Eth + Age, data = quine)
car::Anova(quine.mod1)
## Analysis of Deviance Table (Type II tests)
##
## Response: Days
##     LR Chisq Df Pr(>Chisq)
## Eth   12.496  1  0.0004077 ***
## Age   11.220  3  0.0105926 *
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

First, theta is extracted from quine.mod1; then an NB model dropping one predictor is fitted with the same theta; and finally Anova() carries out a LR test to compare the two models.

I'll calculate the p-value for Eth as a demonstration:

quine.mod3 <- glm(Days ~ Age, data = quine,
# theta is extracted form quine.mod1
family=negative.binomial(quine.mod1$theta)) 1 - pchisq(quine.mod1$twologlik - 2 * as.numeric(logLik(quine.mod3)), 1)
##  0.0004077261

The variance of a negative binomial distribution is $\mu + \mu^2/\theta$, and theta accommodates the Poisson overdisperison. Dropping a predictor from the full model changes the MLE of theta. This is why a p-value produced by car::Anova() is different to that from the LR test of two individually fitted models.

quine.mod2 <- glm.nb(Days ~ Age, data = quine)
anova(quine.mod1, quine.mod2)
## Likelihood ratio tests of Negative Binomial Models
##
## Response: Days
##       Model    theta Resid. df    2 x log-lik.   Test    df LR stat.      Pr(Chi)
## 1       Age 1.147298       142       -1107.816
## 2 Eth + Age 1.249143       141       -1095.801 1 vs 2     1 12.01442 0.0005279048

As shown in this example: theta is 1.249 in quine.mod1 and 1.147 in quine.mod2. The p-value of LR test for Eth is 0.0005.

Reference:

Hilbe, J.M., 2011. Negative Binomial Regression, 2nd ed, Cambridge University Press.

• Thank you for a very useful response - I've always wondered why Anova(modelX), and anova(modelY, modelX) don't give the same p-value in situations like the one you described. Is there a similar issue with using just anova(modelX), so not car::Anova() but simple, lower case anova()? – MGy Sep 8 '18 at 17:46
• Yes, it has a similar issue. For glm.nb, anova() is the generic function for anova.negbin() in the MASS package. It first extracts theta from the full model, then fits a null model and add the predictors sequentially (theta is specified to be the same as the full model), and conducts the LR test based on the sequence. – HanZA Sep 9 '18 at 9:23