Plotting multinomial choice regression model in R I generated myself some data to try some data with 4 levels (health status based on body-mass-index, please don't judge me if the values don't make sense):
library(tidyverse)
library(ZeligChoice)

set.seed(100)
health_excellent <- rnorm(20, mean = 20, sd = 3)
health_good <- rnorm(20, mean = 25, sd = 3)
health_fair <- rnorm(20, mean = 30, sd = 3)
health_poor <- rnorm(20, mean = 35, sd = 3)

health_condition_data <- data.frame(health_excellent, health_good, health_fair, health_poor) %>%
  pivot_longer(cols = health_excellent:health_poor, names_to = "health_condition", values_to = "bmi") %>%
  mutate(health_condition_nr = case_when(health_condition == "health_excellent" ~ 1,
                                         health_condition == "health_good" ~ 2,
                                         health_condition == "health_fair" ~ 3,
                                         health_condition == "health_poor" ~ 4)) %>%
  arrange(health_condition_nr)

health_condition_data %>%
  ggplot(aes(x = bmi, y = health_condition_nr))+
  geom_point()

The plot looks like this:

Now I run an ordered choice logistic regression
health_condition_regression <- zelig(formula = as.factor(health_condition_nr) ~ bmi,
           data = health_condition_data,
           model = "ologit")

summary(health_condition_regression)

And get the following output:
Coefficients:
     Value Std. Error t value
bmi 0.6146     0.0893   6.883

Intercepts:
    Value   Std. Error t value
1|2 13.7532  2.0962     6.5610
2|3 16.7425  2.4167     6.9277
3|4 20.2020  2.9289     6.8975

So far, so good.
I now want to plot this, but I don't know how and also can't find a good example of how such a plot would look. I tried
coef_health <- coefficients(health_condition_regression)
# coef_heatl is 0.6146474

health_condition_data %>%
  ggplot(aes(x = bmi, y = health_condition_nr))+
  geom_point()+
  stat_function(fun = function(x) {exp(coef_health*x)/(1+exp(coef_health*x))})

But that does not give me the wanted result:

I also tried adding in an intercept from the output above like this
health_condition_data %>%
  ggplot(aes(x = bmi, y = health_condition_nr))+
  geom_point()+
  stat_function(fun = function(x) {exp(15.4+coef_health*x)/(1+exp(15.4+coef_health*x))})

but without success.
Maybe I have a wrong understanding of the matter. Help (regarding the plotting and my perhaps wrong understanding) would be very much appreciated.
Thank you
 A: (Note that you're fitting an ordinal model ("ologit"), not a multinomial.) You could plot the predicted probability of falling into a certain category with respect to BMI (see this tutorial). Here's the plot:

Here's the R code I used (I used MASSs polr function to fit the model):
library(tidyverse)
library(MASS)

set.seed(100)
health_excellent <- rnorm(20, mean = 20, sd = 3)
health_good <- rnorm(20, mean = 25, sd = 3)
health_fair <- rnorm(20, mean = 30, sd = 3)
health_poor <- rnorm(20, mean = 35, sd = 3)

health_condition_data <- data.frame(health_excellent, health_good, health_fair, health_poor) %>%
  pivot_longer(cols = health_excellent:health_poor, names_to = "health_condition", values_to = "bmi") %>%
  mutate(health_condition_nr = case_when(health_condition == "health_excellent" ~ 1,
                                         health_condition == "health_good" ~ 2,
                                         health_condition == "health_fair" ~ 3,
                                         health_condition == "health_poor" ~ 4)) %>%
  arrange(health_condition_nr)


health_condition_regression <- polr(formula = as.factor(health_condition_nr) ~ bmi,
                                     data = health_condition_data)
# Set up new data frame for prediction
newdat <- data.frame(
  bmi = seq(17, 43, length = 1000)
)
# Predict the probabilities
newdat <- cbind(newdat, predict(health_condition_regression, newdat, type = "probs"))
# Convert wide dataset to long format
lnewdat <- pivot_longer(newdat, cols = c("1", "2", "3", "4"), names_to = "Level", values_to = "Probability")
# Plot using ggplot2
theme_set(theme_bw())
ggplot(lnewdat, aes(x = bmi, y = Probability, colour = Level)) +
  geom_line(linewidth = 1) +
  labs(
    x = "BMI"
    , y = "Probability"
  ) +
  theme(
    axis.title.y=element_text(colour = "black", size = 17, hjust = 0.5, margin=margin(0,12,0,0)),
    axis.title.x=element_text(colour = "black", size = 17, margin=margin(10,0,0,0)),
    axis.text.x=element_text(colour = "black", size=15),
    axis.text.y=element_text(colour = "black", size=15),
    legend.position="right",
    legend.text=element_text(size=15),
    legend.title = element_text(size = 15),
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    legend.key=element_blank(),
    plot.title = element_text(face = "bold"),
    strip.text.x=element_text(size=15),
    legend.key.size = unit(2,"line")
  )

If you wanted to use the coefficients directly and use functions to plot the predicted probabilities, here's a way to do it (with anonymous functions):
beta <- coef(health_condition_regression) # Coefficient beta
ints <- health_condition_regression$zeta # Intercepts

theme_set(theme_bw())
ggplot(health_condition_data, aes(x = bmi, y = NULL)) +
  geom_function(fun = {\(x) plogis(ints[1] - beta*x)}, linewidth = 1, aes(col = "Excellent")) +
  geom_function(fun = {\(x) (plogis(ints[2] - beta*x) - plogis(ints[1] - beta*x))}, linewidth = 1, aes(col = "Good")) +
  geom_function(fun = {\(x) (plogis(ints[3] - beta*x) - plogis(ints[2] - beta*x))}, linewidth = 1, aes(col = "Fair")) +
  geom_function(fun = {\(x) (1 - plogis(ints[3] - beta*x))}, linewidth = 1, aes(col = "Poor")) +
  scale_colour_manual(name = "Health condition", values = c("#f8766d", "#7cae00", "#00bfc4", "#c77cff"), breaks = c("Excellent", "Good", "Fair", "Poor")) +
  labs(
    x = "BMI"
    , y = "Probability"
  ) +
  theme(
    axis.title.y=element_text(colour = "black", size = 17, hjust = 0.5, margin=margin(0,12,0,0)),
    axis.title.x=element_text(colour = "black", size = 17, margin=margin(10,0,0,0)),
    axis.text.x=element_text(colour = "black", size=15),
    axis.text.y=element_text(colour = "black", size=15),
    legend.position="right",
    legend.text=element_text(size=15),
    legend.title = element_text(size = 15),
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    legend.key=element_blank(),
    plot.title = element_text(face = "bold"),
    strip.text.x=element_text(size=15),
    legend.key.size = unit(2,"line")
  )

