How to find transformed glm.nb bivariate mean and CI for each factor in R?

I want to run a negative binomial regression on various treatments in tmt. For each treatment, I want to find a binomial regression mean and 95% CI and generate a table as shown below.

library(MASS)
df <- data.frame(x = rnorm(12, 4,1), y = rnorm(12, 6,4), tmt = rep(c("A","B","C"), each = 4))

library(ggplot2)
library(dplyr)

mod <- glm.nb(x ~ y + tmt,df)

#Need to store the mean and 95% CI for each treatment in a dataframe df.mean
#Dummy data
x.mean.data y.mean.data yaxis.CI.low yaxis.CI.up xaxis.CI.low xaxis.CI.high  tmt
     1            2         1          3               0          3           A
     2            1         1          1               1          3           B
     1            2         1          0              -1          3           C

I am not 100% sure if regression is the correct choice, but in the end I want to generate a plot that looks like this with means and bivariate CI

To get the 95% CI, you could use confint.default function:

Computes confidence intervals for one or more parameters in a fitted model. There is a default and a method for objects inheriting from class "lm".

Code:

library(MASS)
df <- data.frame(x = rnorm(12, 4,1), y = rnorm(12, 6,4), tmt = rep(c("A","B","C"), each = 4))
mod <- glm.nb(x ~ y + tmt,df)
exp(confint.default(mod))

Output:

               2.5 %   97.5 %
(Intercept) 1.377593 5.405914
y           0.925069 1.104674
tmtB        0.667262 3.067027
tmtC        0.507602 2.849177