How to make a plot of a negative binomial regression (fitted values and 95% confidence bands) when using offset variables?

Question

I am trying to make a plot of a negative binomial model in R but I cannot even extract the confidence limits for the fitted values when using offset variables. Is there any way to make a plot without trying to get the 95% CIs manually?

y <- c(602,38,616,256,21,723,245,176,89,1614,31,27,313,251,345)
x <- c(31,35,21,30,37,26,45,21,74,27,37,37,31,37,25)
offset_1 <-c(72,50,31,30,16,25,75,16,78,40,68,25,71,52,17) 
newdata <- data.frame(y,x,offset_1)
nb.model <- MASS::glm.nb(y ~ x + offset(log(offset_1)), data=newdata)
summary(nb.model)
predict(nb.model, type="response")/newdata$offset_1

Answer 1

Here's how to do it with emmeans + ggplot (ggeffects/sjPlot might automate the process a bit further, but I like the control of doing it myself). Handling offsets in emmeans is discussed here.

library(emmeans)
## generated predicted values / CIs
ee <- emmeans(nb.model, ~x,
       at = list(x = seq(20,75, length.out = 51)), 
       type = "response", offset = log(1))
library(ggplot2)
ggplot(as.data.frame(ee), aes(x, response)) +
    geom_line() +
    geom_ribbon(aes(ymin = asymp.LCL, ymax = asymp.UCL), 
                colour = NA, alpha = 0.2) +
    ## add observed data
    geom_point(data = newdata, aes(y=y/offset_1))

Answer 2

I'm not aware of any direct method (apart from maybe the ggeffects package) of plotting the confidence intervals from a negative binomial model with an offset though we can use broom to return the standard error of the fit and calculate the confidence intervals using a normal approximation (Wald confidence interval).

library(broom)
library(dplyr)
library(ggplot2)
conf_level <- 0.95
cval <- qnorm(1 - ( (1 - conf_level) / 2))
augment(nb.model, se_fit = TRUE) %>%
  mutate(lcl = .fitted - cval * .se.fit,
         ucl = .fitted + cval * .se.fit) %>%
  mutate(across(c(.fitted, lcl, ucl), exp)) %>%
  mutate(across(c(y, .fitted, lcl, ucl), ~ .x / exp(offset(log(offset_1))))) %>%
  ggplot(aes(x = x, y = y)) + 
  geom_point() +
  geom_line(aes(y = lcl), lty = 2) + 
  geom_line(aes(y = ucl), lty = 2)