Plotting confidence intervals with ggeffects for multinom

[janix]

Hello R-people,

I have a question regarding the ggeffects package and its use with multinom functions (from nnet package): I am trying to plot marginal effects for a multinomial regression model. Though ggeffects() should be compatible with multinom, the plot does not display confidence intervals. If I plot the same data with effects(), I do get the CIs.

I couldn't find any example for the use of ggeffects with multinom, so I'd be grateful for any suggestion that works for numerical and factor predictor variables. (I also tried the MNLpred package and its mnl_pred_ova function as an alternative, but this one doesn't seem to be able to deal with factors...)

Here's a minimal example:

library(nnet)
library(ggeffects)
library(effects)

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 = "count")  # Numerical 
plot(pred1)
pred2 <- ggpredict(fit, terms = "case")
plot(pred2)

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

Thanks!

[joels]

I'm not sure why ggpredict isn't adding confidence intervals. However, with some additional coding you can roll your own ggplot effects plot using the output of ggpredict and effects. (Ideally, you would generate the predictions and confidence bands directly from the model object, but it was quicker to just get them from ggpredict and effects.)

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)
}

Now, generate the effects plots:

# 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

To create the ci.fnc function, I looked at the structure of the objects returned by effects and ggpredict. For example, ggpredict returns a data frame (run class(pred2) or str(pred2)), and we print it using as.data.frame to override the default print method for ggpredict objects.

as.data.frame(pred2)

  x predicted response.level group
1 N 0.3614441       1 Better     1
2 N 0.2951495       2 Medium     1
3 N 0.3434064        3 Worse     1
4 Y 0.3448997       1 Better     1
5 Y 0.3290675       2 Medium     1
6 Y 0.3260327        3 Worse     1

effects returns a list. Below is the output of str(pred2e). You can see that the x element contains the values of term variable (the x-axis value in each plot. Then if you scroll to the bottom, you can see that the lower.prob and upper.prob elements contain the confidence limits. So, ci.fnc is set up to extract these. You can also manually extract these values. For example:

pred2e$lower.prob
pred2e$upper.prob

pred2e$upper.prob
pred2e$x

with(pred2e, cbind(x, upper.prob, lower.prob))

str(pred2e)

List of 19
 $ term            : chr "case"
 $ formula         :Class 'formula'  language rating ~ count + case
  .. ..- attr(*, "variables")= language list(rating, count, case)
  .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:3] "rating" "count" "case"
  .. .. .. ..$ : chr [1:2] "count" "case"
  .. ..- attr(*, "term.labels")= chr [1:2] "count" "case"
  .. ..- attr(*, "order")= int [1:2] 1 1
  .. ..- attr(*, "intercept")= int 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(rating, count, case)
  .. ..- attr(*, "dataClasses")= Named chr [1:3] "factor" "numeric" "factor"
  .. .. ..- attr(*, "names")= chr [1:3] "rating" "count" "case"
 $ response        : chr "rating"
 $ y.levels        : chr [1:3] "1 Better" "2 Medium" "3 Worse"
 $ variables       :List of 1
  ..$ case:List of 3
  .. ..$ name     : chr "case"
  .. ..$ is.factor: logi TRUE
  .. ..$ levels   : chr [1:2] "N" "Y"
 $ x               :'data.frame':	2 obs. of  1 variable:
  ..$ case: Factor w/ 2 levels "N","Y": 1 2
 $ model.matrix    : num [1:2, 1:3] 1 1 3.56 3.56 0 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:2] "1" "2"
  .. ..$ : chr [1:3] "(Intercept)" "count" "caseY"
  ..- attr(*, "assign")= int [1:3] 0 1 2
  ..- attr(*, "contrasts")=List of 1
  .. ..$ case: chr "contr.treatment"
 $ data            :'data.frame':	18 obs. of  3 variables:
  ..$ rating: Factor w/ 3 levels "1 Better","2 Medium",..: 1 1 1 2 2 2 3 3 3 1 ...
  ..$ count : num [1:18] 2 0 5 8 10 3 2 1 0 0 ...
  ..$ case  : Factor w/ 2 levels "N","Y": 2 1 2 2 1 2 1 2 1 1 ...
  ..- attr(*, "terms")=Classes 'terms', 'formula'  language rating ~ count + case + (count + case)
  .. .. ..- attr(*, "variables")= language list(rating, count, case)
  .. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:3] "rating" "count" "case"
  .. .. .. .. ..$ : chr [1:2] "count" "case"
  .. .. ..- attr(*, "term.labels")= chr [1:2] "count" "case"
  .. .. ..- attr(*, "order")= int [1:2] 1 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: 0x7fddc4d92df8> 
  .. .. ..- attr(*, "predvars")= language list(rating, count, case)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:3] "factor" "numeric" "factor"
  .. .. .. ..- attr(*, "names")= chr [1:3] "rating" "count" "case"
 $ discrepancy     : num 0
 $ model           : chr "multinom"
 $ prob            : num [1:2, 1:3] 0.361 0.345 0.295 0.329 0.343 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:3] "prob.X1.Better" "prob.X2.Medium" "prob.X3.Worse"
 $ logit           : num [1:2, 1:3] -0.569 -0.642 -0.871 -0.712 -0.648 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:3] "logit.X1.Better" "logit.X2.Medium" "logit.X3.Worse"
 $ se.prob         : num [1:2, 1:3] 0.17 0.165 0.167 0.168 0.17 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:3] "se.prob.X1.Better" "se.prob.X2.Medium" "se.prob.X3.Worse"
 $ se.logit        : num [1:2, 1:3] 0.735 0.729 0.803 0.76 0.754 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:3] "se.logit.X1.Better" "se.logit.X2.Medium" "se.logit.X3.Worse"
 $ lower.logit     : num [1:2, 1:3] -2.01 -2.07 -2.44 -2.2 -2.13 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:3] "L.logit.X1.Better" "L.logit.X2.Medium" "L.logit.X3.Worse"
 $ upper.logit     : num [1:2, 1:3] 0.872 0.787 0.704 0.778 0.829 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:3] "U.logit.X1.Better" "U.logit.X2.Medium" "U.logit.X3.Worse"
 $ lower.prob      : num [1:2, 1:3] 0.1181 0.112 0.0798 0.0995 0.1067 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:3] "L.prob.X1.Better" "L.prob.X2.Medium" "L.prob.X3.Worse"
 $ upper.prob      : num [1:2, 1:3] 0.705 0.687 0.669 0.685 0.696 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:3] "U.prob.X1.Better" "U.prob.X2.Medium" "U.prob.X3.Worse"
 $ confidence.level: num 0.95
 - attr(*, "class")= chr "effpoly"

[janix]

Thank you @joels, this is already very helpful!

I am wondering, however, whether it would not be more adequate in this case to take all the data from effect given that effect and ggpredict seem to use slightly different methods of prediction. That starts with picking different default values of X for prediction (resulting in NAs in the data created by ci.fnc, which has to be manually fixed by defining the xlevels in effect).

More importantly, the two methods predict different values for P(Y) - which results in the effect CIs being somewhat "off" on the graph. The differences are small in the above example, but in my larger dataset, I get a predicted probability for rating=better at case=0 of 0.24 with ggpredict and of 0.20 with effect.

[PoWiManu]

Hey @janix,
yes, you are right: the MNLpred package (so far) only works with numeric independent variables. If the use of factors is important, feel free to open a new issue.

[janix]

In case this is of further interest to anyone: I made some minor modifications to @joels' code to fetch the predicted values (plus CIs) from effects only, which can then be fed into a ggplot. As long as ggeffects does not display confidence intervals, this allows to use the internally consistent predictions from effects while getting a nice and customizable ggplot.
Thank you @joels!

ci.fnc = function(effects.obj) {
  col = sym(names(effects.obj$x)) # Capture the name of the effects term
  cis = with(effects.obj, cbind(x, prob, upper.prob, lower.prob)) %>% # Generate a tidy data frame with predictions, lower and upper confidence intervals
    pivot_longer(cols=-col) %>% 
    separate(name, into=c('type', 'response.level'), sep="prob.X") %>% 
    mutate(response.level=gsub("\\.", " ", response.level)) %>% 
    mutate(type=gsub("\\.", "", type)) %>% 
    spread(type, value)
  return(cis)
}

matches = ci.fnc(pred1e) %>% 
  ggplot(aes(matches, V1)) +
  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:14) +
  labs(y="probability", x="matches", title="Matches effect plot")