How to report negative binomial regression results from R

Question

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, 
    link = log)

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?

Answer 1

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.

Answer 2

It is typical to report a percentage reduction:

round((negbinom$family$linkinv(-.47466)-1)*100)
[1] -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

model1 <- glm.nb(Frequency~Groups,data=x,link=log)

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 

Answer 3

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)
## [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.