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

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

Solution

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