Negative binomial pairwise comparison in R

Question

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.

Answer 1

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

Edit to address OPs comments:

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