How to plot confidence intervals for multinominal logit predicted values in R? I am going to use UCLA example on multinominal logit as a running example. I want to analyze the predicted probability of each response based on social-economic status (ses) and the writing grad (write)---
ml <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")

ml$prog2 <- relevel(ml$prog, ref = "academic")
test <- multinom(prog2 ~ ses + write, data = ml)

dwrite <- data.frame(ses = rep(c("low", "middle", "high"), each = 41), write = rep(c(30:70),
    3))

## store the predicted probabilities for each value of ses and write
pp.write <- cbind(dwrite, predict(test, newdata = dwrite, type = "probs", se = TRUE))

lpp <- melt(pp.write, id.vars = c("ses", "write"), value.name = "probability")

ggplot(lpp, aes(x = write, y = probability, colour = ses)) + geom_line() + facet_grid(variable ~
    ., scales = "free")


For that, my objective is to replicate the following plot with confidence intervals.

I have already tried the following without success: R How to get confidence interval for multinominal logit?
 A: You can use emmeans to get predicted probabilities and standard errors for each response level, and ggeffects to plot them.
library("tidyverse")

data <- foreign::read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")

# reference categories:
data$prog <- relevel(data$prog, ref = "academic")
data$ses <- relevel(data$ses, ref = "high")

model <- nnet::multinom(prog ~ write + ses, data = data)
#> # weights:  15 (8 variable)
#> initial  value 219.722458 
#> iter  10 value 179.985215
#> final  value 179.981726 
#> converged

predictions <- ggeffects::ggemmeans(model, terms = c("write [all]", "ses"))
#> Loading required namespace: emmeans

# plot(predictions)

predictions %>% ggplot(aes(x = x, y = predicted, colour = group)) +   
  geom_line(aes(colour = group), size = 1) +
  xlab("write") +
  facet_grid(response.level ~ ., scales = "free") + 
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), 
              alpha = 0.2) + 
  theme_bw()



Created on 2022-06-08 by the reprex package (v2.0.1)
