Confidence interval for the difference of two negative binomial rates I have a model with a negative binomial distribution using the glm.nb function from R:
Call:
MASS::glm.nb(formula = Counts ~ Gender + offset(log(Offset_Days)), 
    data = r_dataset, init.theta = 1.023811633, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.5973  -0.5973  -0.5594  -0.5594   3.7897  

Coefficients:
            Estimate Std. Error  z value Pr(>|z|)    
(Intercept) -7.67759    0.02613 -293.787  < 2e-16 ***
GenderMale   0.14226    0.03450    4.124 3.73e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(1.0238) family taken to be 1)

    Null deviance: 12494  on 22136  degrees of freedom
Residual deviance: 12477  on 22135  degrees of freedom
AIC: 22558

Number of Fisher Scoring iterations: 1


              Theta:  1.0238 
          Std. Err.:  0.0863 

 2 x log-likelihood:  -22551.7100

From here I get that the ratio of the two gender groups shows the males having rates about 14% greater than females.
My question is if there is a way to derive a confidence interval for the difference of the rates? From the model I get that the rates for the male and female groups are 0.1948688 0.1690280 respectively. From this data I can get a confidence interval of the ratio, but I'm wondering if there is a way to get a confidence interval for the difference, 0.1948688 - 0.1690280 = 0.0258408?
Also, I only have one covariate here, but in the final model there will be several, both continuous and discrete.

Two ways I see.
Marginal Effects
You're basically asking about a contrast.  The {marginaleffects} library handles this elegantly.
library(tidyverse)
library(rsample)
library(marginaleffects)

y <- c(7, 5, 4, 7, 5, 2, 11, 5, 5, 4, 2, 3, 4, 3, 5, 9, 6, 7, 10, 6, 12,
       6, 3, 5, 3, 9, 13, 0, 6, 1, 2, 0, 1, 0, 0, 4, 5, 1, 5, 3, 3, 4)

sex <- sample(c("M","F"), size = length(y), replace = T)

d <- tibble(y, sex)


# Marginal Effects

fit <- MASS::glm.nb(y ~ ., data = d)

marginaleffects::marginaleffects(fit, variables = 'sex') %>% 
  summary


  Term Contrast  Effect Std. Error z value Pr(>|z|)  2.5 % 97.5 %
1  sex    M - F -0.6364      1.024 -0.6216  0.53419 -2.643   1.37

Model type:  negbin 
Prediction type:  response 

Here, the marginaleffects function call computes the Male/Female contrast.
Bootstrapping
# bootstrap

compute_differences <- function(data, ...){
  d <- analysis(data)
  fit <- MASS::glm.nb(y ~ ., data = d)
  
  # Predict for males
  newdata <-  tibble(sex = c('F','M'))
  preds <- predict(fit, newdata = newdata, type='response')
  
  tibble(term = 'difference', estimate = diff(preds))
}


bootstraps(d, times = 250, apparent=TRUE) %>% 
  mutate(
    difference = map(splits, compute_differences)
  ) %>% 
  int_bca(difference, .fn = compute_differences)

 term       .lower .estimate .upper .alpha .method
  <chr>       <dbl>     <dbl>  <dbl>  <dbl> <chr>  
1 difference  -2.85     -1.03   1.26   0.05 BCa   


Note that both of these methods only really work conditional on all the other variables being the same.  If Sex is associated with a 14% change in the rate, then 14% of 100 is a lot more than 14% of 10.
