Plot predicted probabilities (with confidence intervals) for multinomial models with interaction terms

[glittergirl]

I am looking for a solution to plot predicted probabilities with confidence intervals for multinomial logic regression models with interaction terms. I've tried several attempts and the example below seems the furthest I've gone, but I am open to all suggestions even if it does not follow the steps below.

In addition, I would be curious as to how to remove one of the facets. (For example, if I am not interested in showing the graph representing the "2 Medium" rating.)

I based my example off of the following post: Plotting confidence intervals with ggeffects for multinom

If count and case were an interaction term instead, it appears that I can change ggpredict() to account for it...

df <- data.frame(rating = c("1 Better", "1 Better", "1 Better", "2 Medium", "2 Medium", "2 Medium", "3 Worse", "3 Worse", "3 Worse","1 Better", "1 Better", "1 Better", "2 Medium", "2 Medium", "2 Medium", "3 Worse", "3 Worse", "3 Worse"),
                 count = c(2,0,5,8,10,3,2,1,0,0,9,1,0,5,7,2,9,0),
                 case = c("Y","N","Y","Y","N","Y","N","Y","N","N","Y","N","Y","N","N","Y","N","Y"))
fit <- multinom(rating ~ count*case,
                data = df)
summary(fit)

pred1 <- ggpredict(fit, terms = c("count", "case"))  # Numerical 
plot(pred1)

I changed the effect() as well:

pred1e <- effect(c("count", "case"), fit)
plot(pred1e)

And ran the rest of the code:

library(tidyverse)
library(patchwork)

# Function to combine output of effects and ggpredict into a single 
#  tidy data frame
ci.fnc = function(effects.obj, ggpredict.obj) {
  
  # Capture the name of the effects term
  col = sym(names(effects.obj$x))
  
  # Generate a tidy data frame with lower and upper confidence intervals
  cis = with(effects.obj, cbind(x, upper.prob, lower.prob)) %>% 
    pivot_longer(cols=-col) %>% 
    separate(name, into=c('type', 'response.level'), sep=".prob.X") %>% 
    mutate(response.level=gsub("\\.", " ", response.level)) %>% 
    spread(type, value)
  
  # Rename the term column from "x" to the actual term name
  ggpredict.obj = ggpredict.obj %>% rename(!!col := x)
  
  # joint the predictions and confidence intervals into a single data frame
  cis = left_join(cis, ggpredict.obj) %>% as.data.frame()
  
  return(cis)
}

# Numeric x-axis variable
count = ci.fnc(pred1e, pred1) %>% 
  ggplot(aes(count, predicted)) +
  geom_ribbon(aes(ymin=L, ymax=U), fill="blue", alpha=0.15) +
  geom_line(aes(group=response.level), colour="blue", size=1) +
  facet_grid(. ~ response.level) +
  theme_bw() +
  theme(plot.title=element_text(hjust=0.5),
        panel.grid.minor.x=element_blank()) +
  scale_y_continuous(limits=c(0,1), expand=c(0,0)) +
  scale_x_continuous(breaks=0:10) +
  labs(y="probability", x="count", title="Count effect plot")

# Categorical x-axis variable
case = ci.fnc(pred2e, pred2) %>% 
  ggplot(aes(case, predicted)) +
  geom_line(aes(group=response.level), colour="grey70") +
  geom_pointrange(aes(ymin=L, ymax=U), fatten=7, shape=21, 
                  fill="blue", colour="red", stroke=0) +
  facet_grid(. ~ response.level) +
  theme_bw() +
  theme(plot.title=element_text(hjust=0.5)) +
  scale_y_continuous(limits=c(0,1), expand=c(0,0)) +
  labs(y="probability", x="case", title="Case effect plot")

# Lay out both plots together
count / case

I got the following error:

Error in sym():
! Can't convert a character vector to a symbol.
Run rlang::last_trace() to see where the error occurred.

[joesho112358]

Hi,

Looks like the root of the problem is in col = sym(names(effects.obj$x)) where

print(names(effects.obj$x))
[1] "count" "case"   
Error in `sym()`:
! Can't convert a character vector to a symbol.
Run `rlang::last_trace()` to see where the error occurred.

A similar example:

> z <- list(a = 1, b = 2, c = 3)
> names(z)
[1] "a" "b" "c"
> sym(names(z))
Error in `sym()`:
! Can't convert a character vector to a symbol.
Run `rlang::last_trace()` to see where the error occurred.

I think this might solve that particular problem:

    col = c()
    for (ex in names(effects.obj$x)) {
        append(col, sym(ex))
    }

But it might just push the problem down the road.