I am trying to graph some results from a multinomial regression. Unfortunately, I am having issues including the confidence intervals in them. I have tried the following codes with the following results:

library(sjPlot)
library(nnet)

# Example Data
set.seed(123)

dat <- expand.grid(
  VetoType = c("Total Veto", "No Veto", "Partial Veto"),
  "CoalescenceAmorim" = c(0, 1)
)
dat <- replicate(10, dat, simplify = FALSE)
dat <- do.call(rbind.data.frame, dat)
dat$CoalitionPercentageAmorim[dat$CoalescenceAmorim == 0] <- 
  rnorm(3 * 10, 60, 20)

dat$CoalitionPercentageAmorim[dat$CoalescenceAmorim == 1] <- 
  rnorm(3 * 10, 40, 10)

model9 <- multinom(
  VetoType ~ CoalescenceAmorim + CoalitionPercentageAmorim,
  data = dat
)

plot_model(model9,
  type = "pred",
  terms = c("CoalitionPercentageAmorim", "CoalescenceAmorim [0,1]"),
  mdrt.values = "meansd",
  ci.lvl = 0.95,
  title = "Figure: Predicted Probability of Partial Vetoes of PLs",
  axis.title = c("Coalition Size", "Predicted Probability"),
  legend.title = "COALITION COALESCENCE"
)

Similarly, I have tried the following for what is my main graph of interest:

model_data <- get_model_data(model9,
  type = "pred",
  terms = c("CoalitionPercentageAmorim", "CoalescenceAmorim [0,1]"),
  mdrt.values = "meansd",
  ci.lvl = 0.95
)


library(ggplot2)

model_data |> 
  subset(response.level == "Total Veto") |> 
  ggplot(aes(x , predicted, color = factor(group))) +
  geom_line() +
  scale_color_brewer(type = "qual", palette = 6) +
  scale_y_continuous(
    labels = scales::percent,
    limits = c(0, NA)
  ) +
  facet_wrap(~response.level)

But once again, these do not include the confidence intervals despite the specification ci.lvl being included.