# Understanding emmeans outputs for poisson and negative binomial GLM fitted on count data with or without offset

I fitted a poisson and negative binomial GLM on count data (=larva) and try to explain it as a function of a factor (=modality). However, as the traps used to trap larvae from inflorescences were exposed to different number of inflorescences, I than added an offset (=inflo) to the models for count normalization (in this case I am modeling the number of larvae per inflorescence).

I used the following data:

mydata <- data.frame(
modality = c("ER","ER","ER","ER","ER","ER","ER","ER","ER","ER","ER","ER","ER","ER","EH","EH","EH","EH","EH","EH","EH","EH","EH","EH","EH","EH","EH","EH","EH","EH","EH","EH","TS","TS","TS","TS","TS","TS","TS","TS","TS","TS","TS","TS","TS","TS"),
larva = c(149,184,51,35,10,6,102,29,151,37,57,95,44,38,2,245,29,22,30,32,124,42,17,49,39,36,60,14,73,22,16,21,21,54,53,39,41,58,47,42,12,1,4,2,3,1),
inflo = c(61,48,68,28,33,15,49,31,87,40,21,27,25,20,9,10,28,23,20,32,102,34,18,31,33,32,86,17,77,58,20,22,50,53,38,31,21,35,78,48,15,22,3,10,13,5)
) %>%
mutate(rate = larva / inflo)


Mean values per modality for count and rate data are :

mydata_stat <- mydata %>%
group_by(modality) %>%
summarise(mean_larva = mean(larva),
mean_inflo = mean(inflo),
mean_rate = mean(rate),
n = length(larva)
)
mydata_stat

# A tibble: 3 x 5
modality mean_larva mean_inflo mean_rate     n
<chr>         <dbl>      <dbl>     <dbl> <int>
1 EH             48.5       36.2     2.28     18
2 ER             70.6       39.5     1.75     14
3 TS             27         30.1     0.856    14


From what I understand, a count GLM models the conditional mean of counts. If μ is the mean response conditional on the level of modality and the log-link function is η = log(μ), then the model is:
η = β0 + β1 modalityER + β2 modalityTS, that is:
μ = exp(β0 + β1 modalityER + β2 modalityTS)
And for count GLM with an offset models the conditional mean of normalized counts. In this case, the model is:
η = β0 + β1 modalityER + β2 modalityTS + log(inflo), that is:
μ / inflo = exp(β0 + β1 modalityER + β2 modalityTS)

1/ Here the models and emmeans outputs for the Poisson and negative binomial GLMs without offset:

fit1 <- glm(larva ~ modality, family = "poisson", data = mydata)
emmeans(fit1, ~ modality, type="response")

modality rate   SE  df asymp.LCL asymp.UCL
EH       48.5 1.64 Inf      45.4      51.8
ER       70.6 2.25 Inf      66.3      75.1
TS       27.0 1.39 Inf      24.4      29.9

Confidence level used: 0.95
Intervals are back-transformed from the log scale

fit3 <- glm.nb(larva ~ modality, data = mydata)
emmeans(fit3, ~ modality, type="response")

modality response    SE  df asymp.LCL asymp.UCL
EH           48.5 10.59 Inf      31.6      74.4
ER           70.6 17.40 Inf      43.5     114.4
TS           27.0  6.75 Inf      16.5      44.1

Confidence level used: 0.95
Intervals are back-transformed from the log scale


I was expecting that the response means provided by emmeans was equal to observed means (=mean_larva in mydata_stat), which is the case for both models.

2/ Here the models and emmeans outputs for the Poisson and negative binomial GLMs with offset:

fit1_offset <- glm(larva ~ modality + offset(log(inflo)), family = "poisson", data = mydata)
emmeans(fit1_offset, ~ modality, type="response", offset = log(1))

modality  rate     SE  df asymp.LCL asymp.UCL
EH       1.339 0.0453 Inf      1.25     1.431
ER       1.787 0.0568 Inf      1.68     1.902
TS       0.896 0.0461 Inf      0.81     0.991

fit3_offset <- glm.nb(larva ~ modality + offset(log(inflo)), data = mydata)
emmeans(fit3_offset, ~ modality, type="response", offset = log(1))

modality response    SE  df asymp.LCL asymp.UCL
EH          2.245 0.478 Inf     1.478      3.41
ER          1.755 0.424 Inf     1.093      2.82
TS          0.861 0.215 Inf     0.528      1.41

Confidence level used: 0.95
Intervals are back-transformed from the log scale


In the case of count GLMs with offset, I was also expecting that the response rates provided by emmeans was equal to observed rates (=mean_rate in mydata_stat), which is not the case.
In the case of the negative binomial GLM, it was close to mean_rate but not exactly the same, particularly for EH modality. And in the case of the Poisson GLM, it was instead exactly equal to mean_rate_2 = mean_larva / mean_inflo.

mydata_stat %>% mutate(mean_rate_2 = mean_larva / mean_inflo)

# A tibble: 3 x 6
modality mean_larva mean_inflo mean_rate     n mean_rate_2
<chr>         <dbl>      <dbl>     <dbl> <int>       <dbl>
1 EH             48.5       36.2     2.28     18       1.34
2 ER             70.6       39.5     1.75     14       1.79
3 TS             27         30.1     0.856    14       0.896


I can't understand why...

I realized that there is a possible oulier in modality EH:

boxplot(rate ~ modality, data=mydata) When removing the outlier, response rates provided by emmeans were much more closer between the two models and also closer the observed rates (=mean_rate in mydata_stat_cor - see below) :

mydata_cor <- mydata %>% filter(rate < 20)
mydata_stat_cor <- mydata_cor %>%
group_by(modality) %>%
summarise(mean_larva = mean(larva),
mean_inflo = mean(inflo),
mean_rate = mean(rate),
n = length(larva)
)
mydata_stat_cor %>% mutate(mean_rate_2 = mean_larva / mean_inflo)

# A tibble: 3 x 6
modality mean_larva mean_inflo mean_rate     n mean_rate_2
<chr>         <dbl>      <dbl>     <dbl> <int>       <dbl>
1 EH             36.9       37.8     0.976    17       0.978
2 ER             70.6       39.5     1.75     14       1.79
3 TS             27         30.1     0.856    14       0.896

fit1_offset <- glm(larva ~ modality + offset(log(inflo)), family = "poisson", data = mydata_cor)
fit3_offset <- glm.nb(larva ~ modality + offset(log(inflo)), data = mydata_cor)

as.data.frame(emmeans(fit1_offset, ~ modality, type="response", offset = log(1)))$rate  0.9781931 1.7866184 0.8957346 as.data.frame(emmeans(fit3_offset, ~ modality, type="response", offset = log(1)))$response
 0.9833245 1.7576010 0.8702802


Does it mans that mean_rates are managed differently in Poisson and Negative binomila GLMs? How can these differences in results be explained ? Thanks in advance for your helpful explanations.

• Note also that the SEs of the estimates for the NB models are 5--10 times greater than those for the Poisson models. Seems like the NB models are numerically poorly conditioned. Jul 28 at 2:39

First of all, the one question I actually can answer is the one about dividing by the mean inflo. This is indeed the case; look at:

> emmeans(fit3_offset, "modality")@grid
modality .offset. .wgt.
1       EH 3.565852    18
2       ER 3.565852    14
3       TS 3.565852    14

> log(mean(mydata\$inflo))
 3.565852


This happens because inflo is considered a covariate in the model, and emmeans() by default sets covariates to their mean, making the offset equal to the log of that mean.. If we do not specify offset in the emmeans() call, then we are adding the .offset. to the linear estimates. If on the other hand we specify offset = 0, then all the estimates will be smaller by that amount; so on the response scale, we'll get the estimates divided by the mean inflo.

I think what is happening here is that, while we indeed have a model for the mean rates (inasmuch that the log offsets are subtracted from the log responses, log(larva) - log(inflo) = log(larva/inflo)), the estimates are just not the same mean rates you computed; it's more complicated than that. If we go back to the the model predictions and compute the rates from those, we can reproduce the estimated rates obtained from emmeans():

> tapply(exp(predict(fit1_offset)) / mydata$$inflo, mydata$$modality, mean)
EH        ER        TS
1.3389571 1.7866184 0.8957346

> tapply(exp(predict(fit3_offset)) / mydata$$inflo, mydata$$modality, mean)
EH        ER        TS
2.2445143 1.7552025 0.8613856


Both of these are legitimate estimates of mean rates; they are just different models so they yield different results.

• Thank you for your explanation. I have a good understanding of how mean rates are calculated from parameter estimates. However, I was expecting that estimates would be such that both models predict the same mean rates as the observed one, but that only their standard errors would be different (which is indeed the case: due to overdispersion, the SE is underestimated for Poisson model). However, when the outlier from EH is removed, my expectation seems to be correct. It looks as if the estimated parameters are differently sensitive to the presence of outliers between the two models? Jul 30 at 12:32
• It's reasonable to expect that removing outliers leads to more stable results. But even with the outlier removed, the two models don't produce exactly the same mean rates as you calculated "by hand." It is more complex than that. Jul 31 at 16:35

TLDR: I don't address the question how emmeans handles an offset in a count model, which is answered by @RussLenth (+1). Instead I plot the data which suggests that the counts are overdispersed, with a different dispersion in each modality.

Let's start by plotting the data. We can do better than a boxplot by showing larva against inflo by modality. We learn lots from these three scatterplots:

• There is a linear relationship between inflo and larva, so the model should include inflo as an offset. (You have already done this.)
• The counts are overdispersed, esp. for the ER modality. A Poisson model won't be appropriate for this data; instead consider a quasipoisson model or a negative binomial model.
• The dispersion (ie. the variance) is different for the different modalities. For example, even by eye, ER has about twice as much variability than EH.
• There is an outlier in the EH subset. You should check whether this might be a data entry issue.
• But even more intriguingly, there seems to be a shift in the rate above inflo = 50 or so. If we ignore the outlier, notice how a straight line that fits the points below inflo = 50 will pass above the observations for inflo = 50. Something might have happened during the experiment or else there might be something that happenes to larva when the inflo is too high? There are few observations collected at high inflo, so the data can't tell what's going on.

So to summarize: The data suggests to consider a quasipoisson or a negative binomial model and model the dispersion as well as the mean rate as a function of modality. Since modality is a factor variable, I make a first stab at this analysis by fitting three separate models, one for each modality. I fit both the quasipoisson and the negative binomial regression; the results are similar.

Note: I remove the outlier before I do the analysis.

Here are the estimated dispersion parameters from the quasipoisson models of each modality. These are substantially different and indicate that it may not be reasonable to assume the same dispersion across modalities (as fit1 and fit3 do implicitly).

#>   modality dispersion
#> 1 EH             3.88
#> 2 ER            24.4
#> 3 TS             9.08


Estimates of the rate $$\mu$$ with quasipoisson regression:

#>   modality estimate conf.low conf.high
#> 1 EH          0.978    0.835      1.14
#> 2 ER          1.79     1.29       2.39
#> 3 TS          0.896    0.650      1.20


as well as with negative binomial regression:

#>   modality estimate conf.low conf.high
#> 1 EH          0.989    0.842      1.17
#> 2 ER          1.76     1.28       2.51
#> 3 TS          0.864    0.589      1.32


Notice that — should we assume to trust the confidence intervals produced by the two analyses — the CI for the ER rate doesn't overlap with the CIs for the EH and TS rates.

And finally, the data (with one "unusual" data point removed) and with the mean rate estimated by the quasipoisson regression. The model for the EH modality is not a great fit: notice how the line does its best to fit the observations for low inflo (< 50) and high inflo (> 50) with the end result that it fits neither particularly convincingly. R code to reproduce the analysis and the figures:


library("MASS")
library("broom")
library("tidyverse")

calculate_dispersion <- function(object) {
# NOTE: This snippet is copied from summary.glm.
df.r <- object$$df.residual fam <- object$$family
if (!is.null(fam$$dispersion) && !is.na(fam$$dispersion)) {
fam$$dispersion } else if (fam$$family %in% c("poisson", "binomial")) {
1
} else if (df.r > 0) {
est.disp <- TRUE
if (any(object$$weights == 0)) { warning("observations with zero weight not used for calculating dispersion") } sum((object$$weights * object$$residuals^2)[object$$weights > 0]) / df.r
} else {
est.disp <- TRUE
NaN
}
}

mydata <-
tibble(
modality = c("ER", "ER", "ER", "ER", "ER", "ER", "ER", "ER", "ER", "ER", "ER", "ER", "ER", "ER", "EH", "EH", "EH", "EH", "EH", "EH", "EH", "EH", "EH", "EH", "EH", "EH", "EH", "EH", "EH", "EH", "EH", "EH", "TS", "TS", "TS", "TS", "TS", "TS", "TS", "TS", "TS", "TS", "TS", "TS", "TS", "TS"),
larva = c(149, 184, 51, 35, 10, 6, 102, 29, 151, 37, 57, 95, 44, 38, 2, 245, 29, 22, 30, 32, 124, 42, 17, 49, 39, 36, 60, 14, 73, 22, 16, 21, 21, 54, 53, 39, 41, 58, 47, 42, 12, 1, 4, 2, 3, 1),
inflo = c(61, 48, 68, 28, 33, 15, 49, 31, 87, 40, 21, 27, 25, 20, 9, 10, 28, 23, 20, 32, 102, 34, 18, 31, 33, 32, 86, 17, 77, 58, 20, 22, 50, 53, 38, 31, 21, 35, 78, 48, 15, 22, 3, 10, 13, 5)
) %>%
arrange(
modality
)

mydata %>%
ggplot(
aes(inflo, larva)
) +
geom_point() +
facet_grid(
~modality
)

# Remove 1 outlier.
mydata <- mydata %>%
filter(
larva < 200
)

fits <- mydata %>%
group_by(
modality
) %>%
nest(
data = c(larva, inflo)
) %>%
mutate(
fit.qp = map(data, ~ glm(larva ~ 1 + offset(log(inflo)), family = quasipoisson, data = .)),
fit.nb = map(data, ~ glm.nb(larva ~ 1 + offset(log(inflo)), data = .))
)

fits %>%
mutate(
dispersion = map_dbl(fit.qp, ~ calculate_dispersion(.))
) %>%
select(
modality, dispersion
)

mus.qp <- fits %>%
mutate(
tidy = map(fit.qp, tidy, conf.int = TRUE, exponentiate = TRUE)
) %>%
unnest(
tidy
) %>%
select(
modality, estimate, conf.low, conf.high
)
mus.qp

mus.nb <- fits %>%
mutate(
tidy = map(fit.nb, tidy, conf.int = TRUE, exponentiate = TRUE)
) %>%
unnest(
tidy
) %>%
select(
modality, estimate, conf.low, conf.high
)
mus.nb

mydata %>%
ggplot(
aes(inflo, larva)
) +
geom_point() +
facet_grid(
~modality
) +
geom_abline(
aes(intercept = 0, slope = estimate, group = modality),
data = mus.qp
)

• Thank you for your feedback. I was well aware of the overdispersion (I was thinking of using the NB model) and the presence of the outlier for the EH modality. However, you raised two interesting points that I hadn't noticed before: 1) the different response pattern before and after inflos = 50, and 2) the difference in the dispersion parameter between the three modalities. Jul 30 at 12:34
• My concern is how can I deal with this? To model the dispersion as a function of modality, you fit 3 models. But in this case, how do I get the p-value that gives the significance of the modality effect? Can I get it from the 3 fitted models or can I fit a single model where the dispersion parameter depends on modality? With the fit3_offset model I would do this: library(car) Anova(fit3_offset, test="LR") fit_emm <- emmeans(fit3_offset, specs=c("modality"), type="response", offset = log(1)) contrast(fit_emm, method = "pairwise", adjust="tukey") %>% summary(infer = TRUE) Jul 30 at 12:35
• (Caveat: Personally I think significance tests and p-values are not very informative, so take this with a grain of salt.) (a) I would look at the confidence intervals; they tell us something about each modality. Isn't confidence intervals what you are getting with emmeans as well? (b) I got interested as well about how to do the "different-dispersions-in-one-model" analysis. This CV thread suggests VGAM for standard analysis and brms for Bayesian analysis. I haven't tried it though. Jul 30 at 12:47
• However, I expect that since modality is a categorical variable, you would get very similar results with the "three separate models" approach and the "estimate the mean rate and the dispersion for each modality in one big model" approach. On the other hand, if you had other variables (covariates), I would expect a bigger model to be a better choice because the big model would learn about the effects of the covariates from the entire dataset. Jul 30 at 12:51
• For point 2, the only solution would be to fit one model for inflos <50 and a second for inflos >50 and test the modality effect for both situations? Jul 30 at 13:07