31 Multinomial Logit (`mlogit`)

This vignette shows how to estimate and interpret multinomial logit (mlogit) models with the marginaleffects package. Multinomial logit models are useful when we wish to model choices among three or more unordered alternatives—for example, commuters selecting a mode of transport (car, bus, bicycle, walk). In these models the data are arranged in “long” format, with one row for every unit–option combination. Because this structure differs from the single‑row‑per‑unit layout used by most other models, current support in marginaleffects is limited to the predictions() and avg_predictions() functions. The examples below show that, by using the hypothesis argument, we can still obtain marginal effects and treatment effects that answer substantive questions about policy interventions.

31.1 Note on implementation

This tutorial demonstrates how to work with multinomial logit models using predictions() together with custom hypothesis functions, rather than comparisons(). While comparisons() would be more complex to implement for multinomial models, the predictions() approach with custom hypotheses provides a flexible and powerful alternative. With minimal additional work, this approach can be extended to support more advanced analyses while maintaining the core functionality needed for most use‑cases.

31.2 Fit the model

library(dplyr)
library(mlogit)
library(marginaleffects)

# Fishing mode‑choice data from mlogit
 data(Fishing)
 Fish <- dfidx(Fishing, varying = 2:9, shape = "wide", choice = "mode")

 m <- mlogit(mode ~ price + catch | income, data = Fish)

31.3 Contrast analysis

Assume we want to know how a 100‑dollar increase in beach‐trip prices affects mode choice. We create two replicas of the original data: one with the observed prices (lo) and one with the higher beach price (hi). The term variable tags each scenario; later we will contrast them.

# Prepare data for the two scenarios
lo <- transform(as.data.frame(Fish), term = "lo")
lo$id1 <- Fish$idx$id1
lo$id2 <- Fish$idx$id2
lo$idx <- NULL

hi <- transform(lo,
    price = ifelse(id2 == "beach", price + 100, price),
    term = "hi")

dat <- rbind(lo, hi)
dat$term <- factor(dat$term, levels = c("lo", "hi"))

31.4 Average predicted probabilities

The call below obtains unit‑level fitted values and then averages them by choice group (the alternative) and term (scenario):

predictions(m,
    newdata = dat,
    by = c("group", "term"))


 Term   Group Estimate Std. Error     z Pr(>|z|)     S   2.5 % 97.5 %
   lo beach     0.1134    0.00843 13.45   <0.001 134.5 0.09684 0.1299
   lo boat      0.3536    0.01304 27.13   <0.001 535.9 0.32809 0.3792
   lo charter   0.3824    0.01338 28.58   <0.001 594.2 0.35617 0.4086
   lo pier      0.1506    0.00928 16.22   <0.001 194.1 0.13239 0.1688
   hi beach     0.0125    0.00233  5.38   <0.001  23.7 0.00796 0.0171
   hi boat      0.3781    0.01382 27.36   <0.001 545.3 0.35098 0.4051
   hi charter   0.4104    0.01411 29.08   <0.001 615.3 0.38274 0.4381
   hi pier      0.1990    0.00996 19.97   <0.001 292.4 0.17949 0.2186

Type: response

31.5 Alternative: custom hypothesis function

The same result can be produced with a custom summary function supplied to the hypothesis argument. Custom functions offer considerable flexibility:

h <- \(x) aggregate(estimate ~ term + group, data = x, FUN = mean)
predictions(m, newdata = dat, hypothesis = h)


 Term   Group Estimate Std. Error     z Pr(>|z|)     S   2.5 % 97.5 %
   lo beach     0.1134    0.00843 13.45   <0.001 134.5 0.09684 0.1299
   hi beach     0.0125    0.00233  5.38   <0.001  23.7 0.00796 0.0171
   lo boat      0.3536    0.01304 27.13   <0.001 535.9 0.32809 0.3792
   hi boat      0.3781    0.01382 27.36   <0.001 545.3 0.35098 0.4051
   lo charter   0.3824    0.01338 28.58   <0.001 594.2 0.35617 0.4086
   hi charter   0.4104    0.01411 29.08   <0.001 615.3 0.38274 0.4381
   lo pier      0.1506    0.00928 16.22   <0.001 194.1 0.13239 0.1688
   hi pier      0.1990    0.00996 19.97   <0.001 292.4 0.17949 0.2186

Type: response

Instead of using the aggregate() function from base R, we can use the dplyr package to summarize the data.

library(dplyr) 
h <- function(x) {
  x |> summarize(estimate = mean(estimate), .by = c("term", "group"))
}
predictions(m, newdata = dat, hypothesis = h)

31.6 Average marginal effects (risk differences)

Finally, we estimate the average effect of the price increase on the expected probability of each mode. We proceed in three steps inside a custom hypothesis function:

Average the unit‑level fitted values within each term and group.
Compute the difference between the hi and lo averages.
Rename the grouping variable to term, as expected by marginaleffects.

h <- function(x) {
    x |>
        aggregate(estimate ~ term + group, data = _, FUN = mean) |>
        aggregate(estimate ~ group, data = _, FUN = diff) |>
        setNames(c("term", "estimate"))
}
predictions(m, newdata = dat, hypothesis = h)


    Term Estimate Std. Error     z Pr(>|z|)     S   2.5 %  97.5 %
 beach    -0.1008    0.00749 -13.5   <0.001 134.9 -0.1155 -0.0862
 boat      0.0244    0.00226  10.8   <0.001  88.2  0.0200  0.0288
 charter   0.0280    0.00257  10.9   <0.001  89.5  0.0230  0.0330
 pier      0.0484    0.00426  11.4   <0.001  97.1  0.0401  0.0568

Type: response

Interpretation: a 100‑dollar increase in beach‑trip prices lowers the probability of choosing the beach by roughly ten percentage points while modestly raising the probabilities of the other modes.

Or the same results can be obtained using dplyr:

h <- function(x) {
  x |>
    summarize(estimate = mean(estimate), .by = c("term", "group")) |>
    summarize(estimate = diff(estimate), .by = "group") |>
    rename(term = group)
}
predictions(m, newdata = dat, hypothesis = h)