# Negative binomial pairwise comparison in R

I have fitted a negative binomial regression model to my data, and the summary of this compares latency of 3 resources to that of burrows:

NegativeBinomalLatencyModel <- glm.nb(Latency_s ~ Resource, data = Cricket)

summary(NegativeBinomalLatencyModel)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)      4.9416     0.3055  16.178   <2e-16 ***
ResourceFemale   0.3292     0.4895   0.672    0.501
ResourceFood     0.1878     0.4228   0.444    0.657
ResourceNone    -0.2179     0.4231  -0.515    0.606


I was wondering how to produce a pairwise comparison of this, comparing each resource to all the other resources.

You could do it using the emmeans package.

Then simply do:

m_means <- emmeans(NegativeBinomalLatencyModel, ~ Resource)
#TO GET PAIRWISE COMPARISONS WITH DIFFERENCES INDICATED AS LETTERS
cld(m_means, Letters = letters)


The emmeans package has a very good documentation (see link above).

If you want to plot the data, you can do it simply via the emmip() function (from the emmeans package). Have look at ?emmip for details. Using your specific example a basic plot could be generated like this:

#BASIC PLOT
emmip(m_means, ~ Resource)
#BASIC PLOT WITH CONFIDENCE LIMITS
emmip(m_means, ~ Resource, CIs=T)
#BASIC PLOT WITH CONFIDENCE LIMITS ON THE RESPONSE SCALE
emmip(m_means, ~ Resource, CIs=T, type="response")


Another way of plotting can be achieved by simply using the plot() function. For that have a look at ?plot.emmGrid.

If you want more control, you can store the output of cld() in an object such as this:

m_means_table <- cld(m_means, Letters = letters)


This can then be used in ggplot2 for example:

require(ggplot2)
ggplot(m_means_table, aes(x=Resource, y=emmean)) + geom_point() +
geom_errorbar(aes(ymin=emmean-SE, ymax=emmean+SE))


If you want upper and lower confidence limits, you can simply replace emmean-SE with asymp.LCL and emmean+SE with asymp.UCL, respectively (from the m_means_table object).

Also have a look at my answer on poisson and glm.nb models here: poisson glm to observe whether effects of artificial light on the number of bat passes in each location were significant

• Thanks, I'm having some trouble interpreting the output. It looks like this Resource emmean SE df asymp.LCL asymp.UCL .group 4 None 4.723694 0.2927040 Inf 4.150004 5.297383 a 1 Burrow 4.941642 0.3054607 Inf 4.342951 5.540334 a 3 Food 5.129405 0.2922820 Inf 4.556543 5.702268 a 2 Female 5.270799 0.3825413 Inf 4.521032 6.020567 a Is there a way to plot this? – Harry Mar 9 '18 at 11:40
• @Harry I added some more detail. Have a look and see if this answers you questions. – Stefan Mar 9 '18 at 13:12
• that's exactly what I needed. Thank you very much for your help. – Harry Mar 9 '18 at 13:37