Multinomial Choice Models and the Indepdence of Irrelevant Alternatives#

import arviz as az
import matplotlib.pyplot as plt
import pandas as pd

from pymc_marketing.customer_choice.mnl_logit import MNLogit
from pymc_marketing.paths import data_dir

az.style.use("arviz-darkgrid")
plt.rcParams["figure.figsize"] = [12, 7]
plt.rcParams["figure.dpi"] = 100
%config InlineBackend.figure_format = "retina"

Discrete choice models are a class of statistical models used to analyze and predict choices made by individuals among a finite set of alternatives. The set of alternatives ussually represent a choice between “products” broadly construed.

These models are grounded in utility maximization theory, where each option provides a certain level of utility to the decision-maker, and the option with the highest perceived utility is chosen.

Discrete choice models are widely applied in fields like transportation, marketing, and health economics to understand behavior and inform policy or design. Common variants include the multinomial logit, nested logit, and mixed logit models, each capturing different aspects of choice behavior, such as similarity among alternatives or individual heterogeneity. In this notebook we will demonstrate how to specify the multinomial logit model and highlight a property of this model known as the Independence of Irrelevant Alternatives.

The Data#

We will examine a case of choices between heating systems. This data is drawn from an example in the R package mlogit.

Note the data has been formatted in a “wide” fashion where each row represents a choice scenario characterised by (a) the chosen outcome depvar, (b) the product attributes: installation costs (ic_x) and operating costs (oc_x) and (c) fixed attributs of the agent making the choice e.g. the consumer’s income and rooms.

data_path = data_dir / "choice_wide_heating.csv"
df = pd.read_csv(data_path)
df
idcase depvar ic_gc ic_gr ic_ec ic_er ic_hp oc_gc oc_gr oc_ec oc_er oc_hp income agehed rooms region
0 1 gc 866.00 962.64 859.90 995.76 1135.50 199.69 151.72 553.34 505.60 237.88 7 25 6 ncostl
1 2 gc 727.93 758.89 796.82 894.69 968.90 168.66 168.66 520.24 486.49 199.19 5 60 5 scostl
2 3 gc 599.48 783.05 719.86 900.11 1048.30 165.58 137.80 439.06 404.74 171.47 4 65 2 ncostl
3 4 er 835.17 793.06 761.25 831.04 1048.70 180.88 147.14 483.00 425.22 222.95 2 50 4 scostl
4 5 er 755.59 846.29 858.86 985.64 883.05 174.91 138.90 404.41 389.52 178.49 2 25 6 valley
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
895 896 gc 766.39 877.71 751.59 869.78 942.70 142.61 136.21 474.48 420.65 203.00 6 20 4 mountn
896 897 gc 1128.50 1167.80 1047.60 1292.60 1297.10 207.40 213.77 705.36 551.61 243.76 7 45 7 scostl
897 898 gc 787.10 1055.20 842.79 1041.30 1064.80 175.05 141.63 478.86 448.61 254.51 5 60 7 scostl
898 899 gc 860.56 1081.30 799.76 1123.20 1218.20 211.04 151.31 495.20 401.56 246.48 5 50 6 scostl
899 900 gc 893.94 1119.90 967.88 1091.70 1387.50 175.80 180.11 518.68 458.53 245.13 2 65 4 ncostl

900 rows × 16 columns

Seeing that we have 5 alternatives to choose from:

df["depvar"].value_counts()

depvar
gc    573
gr    129
er     84
ec     64
hp     50
Name: count, dtype: int64

The idea of discrete choice models is to determine how propensity to choose these products is driven by their observable attributes.

Utility Theory and Maximisation#

The general idea of these models is that an individuals’ choices among observed alternatives reflect the maximization of an underlying utility function. In this framework, each alternative in a choice set is associated with a latent utility composed of observable components (e.g., cost, travel time) and unobservable factors captured as random errors. Daniel McFadden formalized this approach using random utility models, including the multinomial logit model, which assumes that the unobserved components of utility follow an extreme value distribution. His models provided a rigorous econometric foundation for analyzing choice behavior from observed decisions, enabling estimation of how changes in attributes influence choice probabilities. This work earned McFadden the Nobel Prize in Economics in 2000 and has had lasting influence in transportation, marketing, and labor economics.

To specify our multinomial logit model we therefore specify the form of the utility equations that might drive consumers to purchase particular goods amongst a choice set. Note how we distinguish each of the individual alternatives and the alternative-specific-covariates and the individual-specific-covariates by using the | notation. This is a variant of the standard Wilkinson stlye formula notation for specifying regression models. Under the hood we fit N regression models for each of the N goods. These regressions are used to generate a Utility score which is fed through a softmax transform to give us probabilities of each particular choice. The most probable choice should be the one with the greatest predicted utility.

\[\begin{split} \begin{split} \begin{pmatrix} u_{gc} \\ u_{gr} \\ u_{ec} \\ u_{er} \\ u_{hp} \\ \end{pmatrix} = \begin{pmatrix} gc_{ic} & gc_{oc} \\ gr_{ic} & gr_{oc} \\ ec_{ic} & ec_{oc} \\ er_{ic} & er_{oc} \\ hp_{ic} & hp_{oc} \\ \end{pmatrix} \begin{pmatrix} \beta_{ic} \\ \beta_{oc} \\ \end{pmatrix} \end{split} \end{split}\]
\[\text{softmax}(u)_{j} = \frac{\exp(u_{j})}{\sum_{q=1}^{J}\exp(u_{q})}\]

For more details on the background theory and varieties of utility specification see the PyMC example here

utility_formulas = [
    "gc ~ ic_gc + oc_gc | income + rooms + agehed",
    "gr ~ ic_gr + oc_gr | income + rooms + agehed",
    "ec ~ ic_ec + oc_ec | income + rooms + agehed",
    "er ~ ic_er + oc_er | income + rooms + agehed",
    "hp ~ ic_hp + oc_hp | income + rooms + agehed",
]

mnl = MNLogit(df, utility_formulas, "depvar", covariates=["ic", "oc"])
mnl

<pymc_marketing.customer_choice.mnl_logit.MNLogit at 0x12b8cea50>

Once we have initialised our model class we can estimate the model with the sample command

mnl.sample(
    fit_kwargs={
        "target_accept": 0.97,
        "tune": 2000,
        "idata_kwargs": {"log_likelihood": True},
        "nuts_sampler": "nutpie",
        "progressbar": True,
    }
)

Sampling: [alphas_, betas, betas_fixed_, likelihood]

NUTS[nutpie]: [alphas_, betas, betas_fixed_]



/Users/carlostrujillo/Documents/GitHub/pymc-nb-review/mnl_logit/.venv/lib/python3.13/site-packages/pymc/sampling/mcmc.py:995: FutureWarning: Passing `log_likelihood` via `idata_kwargs` is deprecated and will be removed in future versions. Call `pm.compute_log_likelihood(idata)` instead.
  return _sample_external_nuts(





Sampling: [likelihood]







<pymc_marketing.customer_choice.mnl_logit.MNLogit at 0x12b8cea50>









We can then access the fitted inference data object to inspect the posterior estimates for the model parameters of interest.






az.summary(mnl.idata, var_names=["alphas", "betas", "betas_fixed"])









/Users/carlostrujillo/Documents/GitHub/pymc-nb-review/mnl_logit/.venv/lib/python3.13/site-packages/arviz_stats/base/diagnostics.py:90: RuntimeWarning: invalid value encountered in scalar divide
  (between_chain_variance / within_chain_variance + num_samples - 1) / (num_samples)
/Users/carlostrujillo/Documents/GitHub/pymc-nb-review/mnl_logit/.venv/lib/python3.13/site-packages/arviz_stats/base/diagnostics.py:313: RuntimeWarning: invalid value encountered in scalar divide
  varsd = varvar / evar / 4
/Users/carlostrujillo/Documents/GitHub/pymc-nb-review/mnl_logit/.venv/lib/python3.13/site-packages/arviz_stats/base/diagnostics.py:313: RuntimeWarning: invalid value encountered in scalar divide
  varsd = varvar / evar / 4









  
    

      
      mean
      sd
      eti89_lb
      eti89_ub
      ess_bulk
      ess_tail
      r_hat
      mcse_mean
      mcse_sd
    

  
  
    

      alphas[gc]
      1.44
      0.71
      0.31
      2.6
      739
      1051
      1.01
      0.026
      0.019
    

    

      alphas[gr]
      0.42
      0.8
      -0.84
      1.7
      821
      1218
      1.00
      0.028
      0.021
    

    

      alphas[ec]
      0.86
      0.99
      -0.71
      2.5
      962
      1530
      1.01
      0.032
      0.023
    

    

      alphas[er]
      2.38
      0.87
      1
      3.8
      889
      1346
      1.00
      0.029
      0.021
    

    

      alphas[hp]
      0
      0
      0
      0
      4000
      4000
      NaN
      0
      NaN
    

    

      betas[ic]
      -0.00154
      0.00063
      -0.0025
      -0.00053
      5873
      3175
      1.00
      8.2e-06
      5.8e-06
    

    

      betas[oc]
      -0.00691
      0.00154
      -0.0094
      -0.0044
      2327
      2484
      1.00
      3.2e-05
      2.3e-05
    

    

      betas_fixed[gc, income]
      -0.061
      0.085
      -0.2
      0.075
      1028
      1512
      1.00
      0.0027
      0.0018
    

    

      betas_fixed[gc, rooms]
      -0.008
      0.082
      -0.14
      0.12
      1096
      1496
      1.00
      0.0025
      0.0017
    

    

      betas_fixed[gc, agehed]
      0.0151
      0.0104
      -0.0011
      0.032
      970
      1622
      1.00
      0.00033
      0.00024
    

    

      betas_fixed[gr, income]
      -0.168
      0.096
      -0.32
      -0.016
      1145
      1684
      1.00
      0.0028
      0.002
    

    

      betas_fixed[gr, rooms]
      -0.023
      0.093
      -0.17
      0.13
      1240
      1709
      1.00
      0.0026
      0.0018
    

    

      betas_fixed[gr, agehed]
      0.0186
      0.0118
      -0.00027
      0.038
      1040
      1597
      1.00
      0.00037
      0.00026
    

    

      betas_fixed[ec, income]
      -0.054
      0.111
      -0.23
      0.13
      1301
      1978
      1.00
      0.0031
      0.0021
    

    

      betas_fixed[ec, rooms]
      0.046
      0.108
      -0.13
      0.22
      1458
      1973
      1.00
      0.0028
      0.002
    

    

      betas_fixed[ec, agehed]
      0.0196
      0.0133
      -0.0015
      0.041
      1154
      1637
      1.00
      0.00039
      0.00028
    

    

      betas_fixed[er, income]
      -0.088
      0.105
      -0.26
      0.083
      1259
      1691
      1.00
      0.003
      0.0021
    

    

      betas_fixed[er, rooms]
      0.023
      0.1
      -0.14
      0.19
      1342
      1816
      1.00
      0.0027
      0.002
    

    

      betas_fixed[er, agehed]
      -0.0057
      0.0128
      -0.026
      0.015
      1185
      1724
      1.00
      0.00037
      0.00027
    

    

      betas_fixed[hp, income]
      0
      0
      0
      0
      4000
      4000
      NaN
      0
      NaN
    

    

      betas_fixed[hp, rooms]
      0
      0
      0
      0
      4000
      4000
      NaN
      0
      NaN
    

    

      betas_fixed[hp, agehed]
      0
      0
      0
      0
      4000
      4000
      NaN
      0
      NaN
    

  







Note here how the beta coefficients attaching to the installation and operating costs are negative implying that unit increases in cost drive a decrease in latent subjective utility across the consumer base. It’s this kind of behaviourial insight about the attributes of each product that these models are intended to uncover.


However, some care must be taken in the interpretation of the parameters in these models. Fundamentally discrete choice models are species of causal-inference model where we aim to not merely make inferences about relations in data but make claims about the causal effects of changes in product attributes e.g. how does the market respond to price changes? The criteria of adequacy of a discrete choice model is keyed to how plausible they are as guides to future action and action under counterfactual settings.






New Pricing Intervention#


Consider the following intervention? How might the consumer popoulation respond to a price intervention which targets a particular market segment? We can fit the multinomial logit model as before under a pricing intervention.






new_policy_df = df.copy()
new_policy_df[["ic_ec", "ic_er"]] = new_policy_df[["ic_ec", "ic_er"]] * 1.5

idata_new_policy = mnl.apply_intervention(new_choice_df=new_policy_df)









Sampling: [likelihood]











Then we can inspect how the market share allocated to each of the products is updated after this pricing intervention. On the one hand it makes sense to see a drop in demand for the electrical heating systems which just saw a price increase. However, we can also note a curious fact about the share allocation for the three remaining goods. Each of the goods absorbs an equal 0.08 lift in their market share. This is not a fluke.


The observation stems from IIA property of the multinomial logit model. Under counterfactual interventions market changes lead to somewhat counter intuitive patterns of market substitution. The remaining goods on the market exhibit proportional substitution, where all remaining choices “compete” equally. This is often a quite implausible assumption about market behaviour because it fails to respect properties of market structure or preference. In our case we may wonder if consumers inclined towards electrical heating systems would all equally flock to gas systems of if their is an inherent market bias among that consumer base towards the more eco-friendly heat pump systems?






change_df = mnl.calculate_share_change(mnl.idata, mnl.intervention_idata)
change_df













  
    

      
      policy_share
      new_policy_share
      relative_change
    

    

      product
      
      
      
    

  
  
    

      gc
      0.636687
      0.689231
      0.082527
    

    

      gr
      0.143201
      0.155084
      0.082979
    

    

      ec
      0.071197
      0.043524
      -0.388686
    

    

      er
      0.093125
      0.051522
      -0.446737
    

    

      hp
      0.055791
      0.060640
      0.086910
    

  











fig = mnl.plot_change(change_df, figsize=(20, 8))








../../_images/986b213440e4ecd369f9ebc191f6629c98ca8f187dbae83754b0e17acdd98dd6.png






The idea of structural preference or non-proportional substitution patterns simply cannot be accounted for with the multinomial logit model. This is occasionally called the Red-Bus, Blue Bus paradox because the model property implies that even in a transport market with a 50% between a car and bus. If we introduce another bus similar in all respects to the original, then the new market is split equally between the three goods yielding a 33% share for each. This seems wrong because a new bus should only cannibalise market share from the bus consumers.


Let’s see how this pattern plays out in our example. We can remove a product as follows and re-estimate our multinomial logit.






Market Product Intervention#






new_policy_df = df.copy()
new_policy_df = new_policy_df[new_policy_df["depvar"] != "hp"]

new_utility_formulas = [
    "gc ~ ic_gc + oc_gc | income + rooms + agehed",
    "gr ~ ic_gr + oc_gr | income + rooms + agehed",
    "ec ~ ic_ec + oc_ec | income + rooms + agehed",
    "er ~ ic_er + oc_er | income + rooms + agehed",
]

idata_new_policy = mnl.apply_intervention(
    new_choice_df=new_policy_df, new_utility_equations=new_utility_formulas
)









Sampling: [alphas_, betas, betas_fixed_, likelihood]





NUTS[nutpie]: [alphas_, betas, betas_fixed_]







/Users/carlostrujillo/Documents/GitHub/pymc-nb-review/mnl_logit/.venv/lib/python3.13/site-packages/pymc/sampling/mcmc.py:995: FutureWarning: Passing `log_likelihood` via `idata_kwargs` is deprecated and will be removed in future versions. Call `pm.compute_log_likelihood(idata)` instead.
  return _sample_external_nuts(





Sampling: [likelihood]











Here agains we see a basically proportional re-allocation where the market demand is re-allocated equally across the remaining goods.






change_df = mnl.calculate_share_change(mnl.idata, mnl.intervention_idata)
change_df













  
    

      
      policy_share
      new_policy_share
      relative_change
    

    

      product
      
      
      
    

  
  
    

      gc
      0.636687
      0.673671
      0.058089
    

    

      gr
      0.143201
      0.152152
      0.062507
    

    

      ec
      0.071197
      0.075685
      0.063040
    

    

      er
      0.093125
      0.098492
      0.057635
    

    

      hp
      0.055791
      0.000000
      -1.000000
    

  











fig = mnl.plot_change(change_df, title="Product Removal", figsize=(20, 8))








../../_images/3254bccb23fed0ae1df935a2ac285a1683ceacbbb117299e49e456703a1cd19f.png










Counterfactual Plausibility as a Criteria of Adequacy#


Kenneth Train observes that




“[The IIA] property can be seen either as a restriction imposed by the model or as the natural outcome of a wellspecified model that captures all sources of correlation over alternatives into representative utility, so that only white noise remains. Often the researcher is unable to capture all sources of correlation explicitly, so
that the unobserved portions of utility are correlated and IIA does not hold.” - pg 76 in “Discrete Choice Methods with Simulation”




which suggests that the Multinomial Logit model can have a compelling role in markets where the structure of the preferences aren’t governed by anything we haven’t included in the utility equations. But given the difficulty of such a comprehensive model specification we might want to seek alternative model specifications that can support more plausible counterfactual inference about the patterns of market substitution under intervention. One such model specification is the nested logit model.






%load_ext watermark
%watermark -n -u -v -iv -w -p pymc_marketing,pytensor









Last updated: Wed, 08 Jul 2026

Python implementation: CPython
Python version       : 3.13.9
IPython version      : 9.15.0

pymc_marketing: 1.0.0.dev0
pytensor      : 3.0.7

arviz         : 1.2.0
matplotlib    : 3.10.9
pandas        : 2.3.3
pymc_marketing: 1.0.0.dev0

Watermark: 2.6.0