Several packages in the R ecosystem allow users to estimate multinomial logit model and discrete choice models. This case study illustrates the use of marginaleffects with the nnet and mlogit packages.

We begin by loading two libraries:

nnet package

The multinom function of the nnet package allows users to fit log-lienar models via neural networks. The data used for this function is a data frame with one observation per row, and the response variable is coded a factor. All the marginaleffects package function work seamlessly with this model. For example, we can estimate a model and compute average marginal effects as follows:

library(nnet)

head(mtcars)
#>                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
#> Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
#> Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
#> Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
#> Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
#> Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
#> Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

mod <- multinom(factor(gear) ~ hp + mpg, data = mtcars, trace = FALSE)

marginaleffects(mod, type = "probs") %>% summary() 
#> Average marginal effects 
#>   Group Term     Effect Std. Error  z value   Pr(>|z|)     2.5 %     97.5 %
#> 1     3   hp -3.438e-05   0.002253 -0.01526 0.98782440 -0.004450  0.0043815
#> 2     3  mpg -7.131e-02   0.026441 -2.69703 0.00699606 -0.123134 -0.0194886
#> 3     4   hp -4.667e-03   0.002199 -2.12294 0.03375920 -0.008976 -0.0003583
#> 4     4  mpg  1.591e-02   0.020003  0.79525 0.42646997 -0.023298  0.0551125
#> 5     5   hp  4.702e-03   0.001304  3.60439 0.00031289  0.002145  0.0072584
#> 6     5  mpg  5.540e-02   0.016469  3.36416 0.00076776  0.023126  0.0876827
#> 
#> Model type:  multinom 
#> Prediction type:  probs

Notice that in such models, we get one marginal effect for each term, for each level of the response variable.

mlogit package

The mlogit package uses data in a slightly different structure, with one row per observation-choice combination. For example, this data on choice of travel mode includes 4 rows per individual, one for each mode of transportation:

library("AER")
library("mlogit")
library(tidyverse)
data("TravelMode", package = "AER")

head(TravelMode)
#>   individual  mode choice wait vcost travel gcost income size
#> 1          1   air     no   69    59    100    70     35    1
#> 2          1 train     no   34    31    372    71     35    1
#> 3          1   bus     no   35    25    417    70     35    1
#> 4          1   car    yes    0    10    180    30     35    1
#> 5          2   air     no   64    58     68    68     30    2
#> 6          2 train     no   44    31    354    84     30    2

mod <- mlogit(choice ~ wait + gcost | income + size, TravelMode)

marginaleffects(mod, variables = c("income", "size")) %>% summary()
#> Average marginal effects 
#>   Group   Term     Effect Std. Error z value Pr(>|z|)      2.5 %     97.5 %
#> 1   air income  2.534e-04  0.0001096   2.312 0.020766  0.0000386  0.0004682
#> 2   bus income -2.104e-05  0.0001394  -0.151 0.880006 -0.0002943  0.0002522
#> 3   car income -2.534e-04  0.0001096  -2.312 0.020766 -0.0004682 -0.0000386
#> 4 train income  2.104e-05  0.0001394   0.151 0.880006 -0.0002522  0.0002943
#> 5   air   size -6.048e-03  0.0030190  -2.003 0.045156 -0.0119646 -0.0001305
#> 6   bus   size  1.142e-03  0.0029362   0.389 0.697262 -0.0046125  0.0068970
#> 7   car   size  6.048e-03  0.0030190   2.003 0.045156  0.0001305  0.0119646
#> 8 train   size -1.142e-03  0.0029362  -0.389 0.697262 -0.0068970  0.0046125
#> 
#> Model type:  mlogit 
#> Prediction type:  response

Note that the marginaleffects function will always return estimates of zero for regressors before the vertical bar in the formula. This is because the predict() function supplied by the mlogit package does not produce different predictions for different values of these variables.

To compute different kinds of marginal effects, we can construct customized data frames and feed them to the newdata argument of the marginaleffects function.

Important: The newdata argument for mlogit models must be a “balanced” data frame, that is, it must have a number of rows that is a multiple of the number of choices.

If we want to compute the slope of the response function (marginal effecs) when each of the predictors is fixed to its global mean, we can do:

nd <- TravelMode %>%
    summarize(across(c("wait", "gcost", "income", "size"),
              function(x) rep(mean(x), 4)))
nd
#>       wait    gcost   income     size
#> 1 34.58929 110.8798 34.54762 1.742857
#> 2 34.58929 110.8798 34.54762 1.742857
#> 3 34.58929 110.8798 34.54762 1.742857
#> 4 34.58929 110.8798 34.54762 1.742857

marginaleffects(mod, newdata = nd, variables = c("income", "size")) |> summary()
#> Average marginal effects 
#>   Group   Term     Effect Std. Error z value  Pr(>|z|)      2.5 %     97.5 %
#> 1   air income  6.656e-03  2.426e-03   2.743 0.0060816  1.901e-03  0.0114108
#> 2   bus income -1.141e-03  9.454e-04  -1.207 0.2273591 -2.994e-03  0.0007117
#> 3   car income  6.480e-06  2.032e-05   0.319 0.7497346 -3.334e-05  0.0000463
#> 4 train income -5.521e-03  1.910e-03  -2.890 0.0038527 -9.265e-03 -0.0017767
#> 5   air   size -1.694e-01  5.877e-02  -2.883 0.0039386 -2.846e-01 -0.0542485
#> 6   bus   size  4.672e-02  2.723e-02   1.716 0.0862434 -6.656e-03  0.1000991
#> 7   car   size  1.358e-03  8.808e-04   1.542 0.1230361 -3.680e-04  0.0030845
#> 8 train   size  1.214e-01  4.447e-02   2.729 0.0063528  3.420e-02  0.2085119
#> 
#> Model type:  mlogit 
#> Prediction type:  response

If we want to compute marginal effects with the gcost and wait fixed at their mean value, conditional on the choice of transportation mode:

nd <- TravelMode %>%
    group_by(mode) %>%
    summarize(across(c("wait", "gcost", "income", "size"), mean))
nd
#> # A tibble: 4 × 5
#>   mode   wait gcost income  size
#>   <fct> <dbl> <dbl>  <dbl> <dbl>
#> 1 air    61.0 103.    34.5  1.74
#> 2 train  35.7 130.    34.5  1.74
#> 3 bus    41.7 115.    34.5  1.74
#> 4 car     0    95.4   34.5  1.74

marginaleffects(mod, newdata = nd, variables = c("income", "size")) |> summary()
#> Average marginal effects 
#>   Group   Term     Effect Std. Error z value   Pr(>|z|)      2.5 %    97.5 %
#> 1   air income  0.0060149   0.002332  2.5793  0.0098996  0.0014443  0.010586
#> 2   bus income -0.0007128   0.001461 -0.4878  0.6256800 -0.0035768  0.002151
#> 3   car income  0.0054450   0.002288  2.3800  0.0173140  0.0009609  0.009929
#> 4 train income -0.0107471   0.002563 -4.1925 2.7592e-05 -0.0157713 -0.005723
#> 5   air   size -0.2329273   0.056622 -4.1137 3.8936e-05 -0.3439050 -0.121950
#> 6   bus   size  0.0204397   0.034364  0.5948  0.5519797 -0.0469129  0.087792
#> 7   car   size  0.0678200   0.041226  1.6451  0.0999517 -0.0129810  0.148621
#> 8 train   size  0.1446676   0.047752  3.0295  0.0024493  0.0510747  0.238261
#> 
#> Model type:  mlogit 
#> Prediction type:  response

We can also explore more complex alternatives. Here, for example, only one alternative is affected by cost reduction:

nd <- datagrid(mode = TravelMode$mode, newdata = TravelMode)
nd <- lapply(1:4, function(i) mutate(nd, gcost = ifelse(1:4 == i, 30, gcost)))
nd <- bind_rows(nd)
nd
#>     individual choice     wait    vcost   travel    gcost   income     size
#>  1:          1     no 34.58929 47.76071 486.1655  30.0000 34.54762 1.742857
#>  2:          1     no 34.58929 47.76071 486.1655 110.8798 34.54762 1.742857
#>  3:          1     no 34.58929 47.76071 486.1655 110.8798 34.54762 1.742857
#>  4:          1     no 34.58929 47.76071 486.1655 110.8798 34.54762 1.742857
#>  5:          1     no 34.58929 47.76071 486.1655 110.8798 34.54762 1.742857
#>  6:          1     no 34.58929 47.76071 486.1655  30.0000 34.54762 1.742857
#>  7:          1     no 34.58929 47.76071 486.1655 110.8798 34.54762 1.742857
#>  8:          1     no 34.58929 47.76071 486.1655 110.8798 34.54762 1.742857
#>  9:          1     no 34.58929 47.76071 486.1655 110.8798 34.54762 1.742857
#> 10:          1     no 34.58929 47.76071 486.1655 110.8798 34.54762 1.742857
#> 11:          1     no 34.58929 47.76071 486.1655  30.0000 34.54762 1.742857
#> 12:          1     no 34.58929 47.76071 486.1655 110.8798 34.54762 1.742857
#> 13:          1     no 34.58929 47.76071 486.1655 110.8798 34.54762 1.742857
#> 14:          1     no 34.58929 47.76071 486.1655 110.8798 34.54762 1.742857
#> 15:          1     no 34.58929 47.76071 486.1655 110.8798 34.54762 1.742857
#> 16:          1     no 34.58929 47.76071 486.1655  30.0000 34.54762 1.742857
#>      mode
#>  1:   air
#>  2: train
#>  3:   bus
#>  4:   car
#>  5:   air
#>  6: train
#>  7:   bus
#>  8:   car
#>  9:   air
#> 10: train
#> 11:   bus
#> 12:   car
#> 13:   air
#> 14: train
#> 15:   bus
#> 16:   car

marginaleffects(mod, newdata = nd, variables = c("income", "size")) |> summary()
#> Average marginal effects 
#>   Group   Term     Effect Std. Error    z value Pr(>|z|)      2.5 %    97.5 %
#> 1   air income  2.223e-15  3.618e-11  6.143e-05  0.99995 -7.091e-11 7.092e-11
#> 2   bus income  2.176e-15  7.163e-11  3.038e-05  0.99998 -1.404e-10 1.404e-10
#> 3   car income  9.781e-15  7.024e-11  1.393e-04  0.99989 -1.377e-10 1.377e-10
#> 4 train income  6.010e-15  6.306e-11  9.530e-05  0.99992 -1.236e-10 1.236e-10
#> 5   air   size -2.196e-15  6.238e-10 -3.520e-06  1.00000 -1.223e-09 1.223e-09
#> 6   bus   size  3.578e-14  9.634e-10  3.714e-05  0.99997 -1.888e-09 1.888e-09
#> 7   car   size  1.952e-14  6.045e-10  3.228e-05  0.99997 -1.185e-09 1.185e-09
#> 8 train   size -2.602e-15  8.319e-10 -3.128e-06  1.00000 -1.631e-09 1.631e-09
#> 
#> Model type:  mlogit 
#> Prediction type:  response