Displaying one response level when using plot_model function of sjplot

Question

I have successfully used the plot_model function of sjplot to plot a multinomial logistic regression model. The regression contains an outcome (Veto Type, with 3 levels) and 2 continuous predictors (Coalition Size and coalition Coalescence). I have also changed the values of coalition Coalescence in the plot_model so as to plot predicted effect outcomes based on the coalition Coalescence mean value and SDs

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

However, I want the graph to only display one Veto Type level at a time. I have tried using rm.terms, and subsetting in multiple ways but so far I have not been successful. I would appreciate any help or advice regarding how to display only one level at a time.

Thanks!

addition

I tried using gtable as suggested here but this produced a 1/3 graph rather than one that cover the whole graph. I would still prefer a full graph for each of the levels.

Answer 1

Instead of trying to fiddle around with the plot object via the gtableone option would be to extract the model data using get_model_data and build your desired plot from scratch using ggplot2.

Using some fake random example data.

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
)
#> # weights:  12 (6 variable)
#> initial  value 65.916737 
#> iter  10 value 65.390983
#> final  value 65.390969 
#> converged

plot_model(model9,
  type = "pred",
  terms = c("CoalitionPercentageAmorim", "CoalescenceAmorim [0,1]"),
  mdrt.values = "meansd",
  ci.lvl = 0.00,
  title = "Figure: Predicted Probability of Partial Vetoes of PLs",
  axis.title = c("Coalition Size", "Predicted Probability"),
  legend.title = "COALITION COALESCENCE"
)
#> Data were 'prettified'. Consider using `terms="CoalitionPercentageAmorim
#>   [all]"` to get smooth plots.

model_data <- get_model_data(model9,
  type = "pred",
  terms = c("CoalitionPercentageAmorim", "CoalescenceAmorim [0,1]"),
  mdrt.values = "meansd",
  ci.lvl = 0.00
)
#> Data were 'prettified'. Consider using `terms="CoalitionPercentageAmorim
#>   [all]"` to get smooth plots.

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)