# Plotting probabilities from multinomial regression output

I generated some data to visualize a multinomial logistic regression, where individuals choose a mode of transportation based on their income. I then set up a regression and predicted the probabilities to then plot them. Here's my code:

library(tidyverse)
library(ggpubr)
library(nnet)

# Generating the data --------------------------
set.seed(100)
helicopter <- rnorm(20, mean = 35, sd = 3)
car <- rnorm(20, mean = 30, sd = 3)
bus <- rnorm(20, mean = 25, sd = 3)
bike <- rnorm(20, mean = 20, sd = 3)

transportation_data <- data.frame(helicopter, car, bus, bike) %>%
pivot_longer(cols = 1:4, values_to = "income", names_to = "mode")

# Plotting the data ---------------------------
transportation_plot <- transportation_data %>%
ggplot(aes(x = income, y = mode, color = mode))+
geom_point()+
coord_cartesian(xlim = c(0,50))

# Setting up the regression -------------------
transportation_regression <- multinom(mode~income, data = transportation_data)
summary(transportation_regression)

# Predicting the probabilities ---------------
new_data <- data.frame(income = seq(0,50,0.1))
prediction <- as.data.frame(predict(transportation_regression, new_data, type = "probs"))
new_data <- cbind(new_data, prediction)

# Plotting the probabilities -----------------
prob_plot <- new_data %>%
pivot_longer(2:5, names_to = "mode", values_to = "prob") %>%
ggplot(aes(x = income, y = prob, color = mode))+
geom_line()

# Merging the two plots -----------------------
ggarrange(transportation_plot, prob_plot, nrow = 2)


The regression output is:

Coefficients:
(Intercept)    income
bus          -8.458078 0.3821646
car         -26.317817 1.0150949
helicopter  -69.080279 2.3148401


And the plot looks like this:

What I would like to do now is to plot the same probabilities, but not using the predict() function. I want to use stat_function() and the coefficients of the regression output. My uni-script says that the probability of the choice of alternative j is $$Pr(y_i = j | \bf{x_i}) = \frac{exp(\bf{x_i\prime \beta_j})}{\sum_{h=1}^J exp(\bf{x_i\prime \beta_h})}$$ so I guess I need this function. But I have trouble understanding and implementing this.

EDIT

I tried the following, but it does not yield reasonable results:

ins <- coef(transportation_regression)[1:3]
betas <- coef(transportation_regression)[4:6]

transportation_data %>%
ggplot(aes(x = income))+
stat_function(fun = function(x) { exp(1) / (1 + sum(exp(ins + betas * x))) }, aes(color = "bike"))+
stat_function(fun = function(x) { exp(ins[1] + betas[1] * x) / (1 + sum(exp(ins + betas * x))) }, aes(color = "bus"))+
stat_function(fun = function(x) { exp(ins[2] + betas[2] * x) / (1 + sum(exp(ins + betas * x))) }, aes(color = "car"))+
stat_function(fun = function(x) { exp(ins[3] + betas[3] * x) / (1 + sum(exp(ins + betas * x))) }, aes(color = "helicopter"))


However, if I take one of the function, let's say for car and plug in 30 for x in the console, I get a sensible result (compare to the plot above):

> exp(ins[2] + betas[2] * 30) / (1 + sum(exp(ins + betas * 30)))
[1] 0.73385


So why won't it work as a function of x in ggplot?

• $\exp z = e^z$, $x_i$ is income for the ith person, and $\beta_j$ is the coefficient on income for mode $j$. So take the intercept, add income multiplied by the coefficient, and wrap an exponential around it. The $\Sigma$ in the denominator means do it for each mode and add together. Commented Feb 4, 2023 at 16:38
• This will give you four probabilities for each person that add up to one. Commented Feb 4, 2023 at 16:41
• Try replacing 1 + sum(exp(ins + betas * x)) with 1 + sapply(x, \(x){sum(exp(ins + betas * x))}) in your formulas. It should work then. Commented Feb 4, 2023 at 20:46
• @COOLSerdash thank you, this worked! Could you briefly explain why this is the case? Or point me somewhere where I can learn more about why my approach didn't work? Commented Feb 4, 2023 at 20:51
• Test your code with a small x. For example: x <- c(20, 30); 1 + sum(exp(ins + betas * x)) gives and error. You have to do the sum for each value of x, leading to the use of sapply, which applies the function for each x. Another possibility would be to just explicitly list all terms: (1 + exp(ins[1] + betas[1]*x) + exp(ins[2] + betas[2]*x) + exp(ins[3] + betas[3]*x)). Commented Feb 4, 2023 at 20:57

As COOLSerdash suggested, swapping 1 + sum(exp(ins + betas * x)) with
1 + sapply(x, \(x){sum(exp(ins + betas * x))}) worked:

transportation_data %>%
ggplot(aes(x = income))+
stat_function(fun = function(x) { 1 / (1 + sapply(x, \(x){sum(exp(ins + betas * x))})) }, aes(color = "bike"))+
stat_function(fun = function(x) { exp(ins[1] + betas[1] * x) / (1 + sapply(x, \(x){sum(exp(ins + betas * x))})) }, aes(color = "bus"))+
stat_function(fun = function(x) { exp(ins[2] + betas[2] * x) / (1 + sapply(x, \(x){sum(exp(ins + betas * x))})) }, aes(color = "car"))+
stat_function(fun = function(x) { exp(ins[3] + betas[3] * x) / (1 + sapply(x, \(x){sum(exp(ins + betas * x))})) }, aes(color = "helicopter"))