Nested Logit and Non-Proportional Patterns of Substitution#

import arviz as az
import arviz_plots as azp
import matplotlib.pyplot as plt
import pandas as pd
from pymc_extras.prior import Prior

from pymc_marketing.customer_choice.nested_logit import NestedLogit
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"

We’ve seen in the other example notebook for consumer choice how the Multinomial Logit model suffers from the IIA limitation that leads to implausible counterfactual inferences regarding market behaviour. If you haven’t read that one, we advise you read it before continuing. We will now show how the nested logit model specification avoids this property by adding more explicit structure to the choice scenarios that are modelled.

In this notebook we will re-use the same heating choice data set seen in the Multinomial Logit example and apply a few different specifications of a nested logit model.

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

Single Layer Nesting#

The important addition gained through the nested logit specification is the ability to specify “nests” of products in this way we can partition the market into “natural” groups of competing products ensuring that there is an inherent bias in the model towards a selective pattern of preference. As before we specify the models using formulas, but now we also add a nesting structure.

## No Fixed Covariates
utility_formulas = [
    "gc ~ ic_gc + oc_gc |",
    "ec ~ ic_ec + oc_ec | ",
    "gr ~ ic_gr + oc_gr | ",
    "er ~ ic_er + oc_er | ",
    "hp ~ ic_hp + oc_hp | ",
]


nesting_structure = {"central": ["gc", "ec", "hp"], "room": ["gr", "er"]}


nstL_1 = NestedLogit(
    df,
    utility_formulas,
    "depvar",
    covariates=["ic", "oc"],
    nesting_structure=nesting_structure,
    alphas_nests=True,
)
nstL_1

<pymc_marketing.customer_choice.nested_logit.NestedLogit at 0x1296967b0>

We will dwell a bit on the manner in which these nests are specified and why. The nested logit model partitions the choice set into nests of alternatives that share common unobserved attributes (i.e., are more similar to each other). It computes the overall probability of choosing an alternative as the product of (1) The probability of choosing a nest (marginal probability), and (2) the probability of choosing an alternative within that nest (conditional probability, given that nest). In our case we want to isolate the probability of choosing a central heatings system and a room based heating system.

Each of the alternatives alt is indexed to a nest. So that we can determine (§) the marginal probability of choosing room or central and (2) conditional probability of choosing ec given that you have chosen central. Our utilities are decomposed into contributions from fixed_covariates and alternative specific covariates:

\(U = Y + W + \epsilon\)

and the probabilities are derived from these decomposed utilitiies in the following manner.

\[ P(i) = P( \text{choose nest B}) \cdot P( \text{choose i} | \text{ i } \in \text{B}) \]

where

\[ P(\text{choose nest B}) = \dfrac{e^{W + \lambda_{k}I_{k}}}{\sum_{l=1}^{K} e^{W + \lambda_{l}I_{l}}} \]

and

\[ P(\text{choose i} | \text{ i} \in \text{B}) = \dfrac{e^{Y_{i} / \lambda_{k}}}{\sum_{j \in B_{k}} e^{Y_{j} / \lambda_{k}}} \]

while the inclusive term \(I_{k}\) is:

\[I_{k} = ln \sum_{j \in B_{k}} e^{Y_{j} / \lambda_{k}} \text{ and } \lambda_{k} \sim Beta(1, 1)\]

such that \(I_{k}\) is used to “aggregate” utilities within a nest a “bubble up” their contribution to decision through the product of the marginal and conditional probabilities. More extensive details of this mathematical formulation can be found in Kenneth Train’s “Discrete Choice Methods with Simulation”.

nstL_1.sample(
    fit_kwargs={
        "idata_kwargs": {"log_likelihood": True},
        "nuts_sampler": "nutpie",
        "progressbar": True,
    }
)

Sampling: [alphas_, alphas_nests, betas, lambdas_nests, likelihood]

NUTS[nutpie]: [alphas_, betas, lambdas_nests, alphas_nests]



/Users/carlostrujillo/Documents/GitHub/pymc-nb-review/nested_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.nested_logit.NestedLogit at 0x1296967b0>









The model structure is quite a bit more complicated now than the simpler multinomial logit as we need to calculate the marginal and conditional probabilities within each of the nests seperately and then “bubble up” the probabilies as a product over the branching nests. These probabilities are deterministic functions of the summed utilities and are then fed into our categorical likelihood to calibrate our parameters against the observed data.






nstL_1.graphviz()








../../_images/0323972c5c6c6b0b795c6ce6f08e28076825d35144b2d9050e5274611fb7542c.svg





But again we are able to derive parameter estimates for the drivers of consumer choice and consult the model implications is in a standard Bayesian model.






az.summary(
    nstL_1.idata,
    var_names=["betas", "alphas", "lambdas_nests", "alphas_nests"],
    round_to=5,
)









/Users/carlostrujillo/Documents/GitHub/pymc-nb-review/nested_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/nested_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
    

  
  
    

      betas[ic]
      -0.00246
      0.00070
      -0.00358
      -0.00138
      2756.01000
      2463.28407
      1.00108
      0.00001
      0.00001
    

    

      betas[oc]
      -0.00623
      0.00166
      -0.00900
      -0.00361
      1745.41915
      1778.70984
      1.00255
      0.00004
      0.00003
    

    

      alphas[gc]
      1.12047
      0.29626
      0.64565
      1.59000
      2308.55102
      2379.51187
      1.00123
      0.00618
      0.00412
    

    

      alphas[ec]
      1.22247
      0.43125
      0.54930
      1.92093
      1851.99416
      1817.08597
      1.00314
      0.01002
      0.00721
    

    

      alphas[gr]
      -0.78266
      3.54242
      -6.69734
      4.89425
      526.98814
      694.94045
      1.00557
      0.15468
      0.11135
    

    

      alphas[er]
      0.59676
      3.53590
      -5.27067
      6.20324
      527.98508
      706.86694
      1.00522
      0.15411
      0.11165
    

    

      alphas[hp]
      0.00000
      0.00000
      0.00000
      0.00000
      4000.00000
      4000.00000
      NaN
      0.00000
      NaN
    

    

      lambdas_nests[central]
      0.83023
      0.12616
      0.59841
      0.98679
      1679.37477
      1696.26872
      1.00211
      0.00295
      0.00225
    

    

      lambdas_nests[room]
      0.79502
      0.14439
      0.53370
      0.98101
      2219.73062
      1735.07855
      1.00174
      0.00281
      0.00192
    

    

      alphas_nests[central]
      0.66517
      0.70919
      -0.46359
      1.79575
      1582.93866
      1733.05377
      1.00379
      0.01796
      0.01215
    

    

      alphas_nests[room]
      -0.69281
      0.70803
      -1.82051
      0.43299
      1406.64895
      1830.68341
      1.00421
      0.01917
      0.01276
    

  







Here we see additionally lambda parameters for each of the nests. These terms measure how strongly correlated the unobserved utility components are for alternatives within the same nest. Closer to 0 indicates a high correlation, substitution happens mostly within the nest. Whereas a value closer to 1 implies lower within nest correlation suggesting IIA approximately holds within the nest. The options available for structuring a market can be quite extensive. As we might have “nests within nests” where the conditional probabilities flow through successive choices within segments of the market.






azp.plot_trace(nstL_1.idata, var_names=["betas"]);








../../_images/b4d78fac06200dbe243e7144eb915f44b1a6d4445b84645761504dcf149e7740.png






Again the parameter estimates seem to be recovered sensibly on the beta coefficients, but the question remains as to whether this additional nesting structure will help support plausible counterfactual reasoning. Note how the model struggles to identify the the intercept terms and places weight on






Making Interventions in Structured Markets#






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

idata_new_policy_1 = nstL_1.apply_intervention(new_choice_df=new_policy_df)









Sampling: [likelihood]











Here we can see that both nesting structures recover non-proportional patterns of product substitution. We have elided the IIA feature of the multinomial logit and can continue to assess whether or not the behaviour implications of these utility theory makes sense.






change_df_1 = nstL_1.calculate_share_change(nstL_1.idata, nstL_1.intervention_idata)
change_df_1













  
    

      
      policy_share
      new_policy_share
      relative_change
    

    

      product
      
      
      
    

  
  
    

      gc
      0.635503
      0.703662
      0.107252
    

    

      ec
      0.071248
      0.026335
      -0.630377
    

    

      gr
      0.150178
      0.180035
      0.198811
    

    

      er
      0.087495
      0.028262
      -0.676986
    

    

      hp
      0.055576
      0.061707
      0.110304
    

  











Visualising the Substitution Patterns#






fig = nstL_1.plot_change(change_df=change_df_1, figsize=(10, 6))

# Nudge overlapping share annotations apart (several products have near-identical shares)
ax = fig.axes[0]
labels = [t for t in ax.texts if t.get_fontsize() == 8]
for x_side in sorted({round(t.get_position()[0], 2) for t in labels}):
    side_labels = sorted(
        (t for t in labels if round(t.get_position()[0], 2) == x_side),
        key=lambda t: t.get_position()[1],
    )
    min_gap = 0.022
    y_prev = None
    for label in side_labels:
        x, y = label.get_position()
        if y_prev is not None and y - y_prev < min_gap:
            y = y_prev + min_gap
            label.set_position((x, y))
        y_prev = y








../../_images/c1f61b7e91adf8c97365f57f53e1e7e9e04c2b1226f65c1405fe1341ee1921bd.png






We can additionally test counterfactuals about the removal of products from the market.






fit_kwargs = {
    "target_accept": 0.97,
    "tune": 2000,
    "nuts_sampler": "nutpie",
    "idata_kwargs": {"log_likelihood": True},
    "progressbar": False,
}

new_policy_df = nstL_1.choice_df.copy()
new_policy_df = new_policy_df[
    (new_policy_df["depvar"] != "hp") & (new_policy_df["depvar"] != "gr")
]

new_utility_formulas = [
    "gc ~ ic_gc + oc_gc | income",
    "ec ~ ic_ec + oc_ec | income ",
    "er ~ ic_er + oc_er | income  ",
]
nstL_1.nesting_structure = {"central": ["gc", "ec"], "room": ["er"]}
nstL_1.alternatives = ["gc", "ec", "er"]
idata_new_policy_3 = nstL_1.apply_intervention(
    new_choice_df=new_policy_df,
    new_utility_equations=new_utility_formulas,
    fit_kwargs=fit_kwargs,
)

change_df_3 = nstL_1.calculate_share_change(nstL_1.idata, nstL_1.intervention_idata)
change_df_3









Sampling: [alphas_, alphas_nests, betas, betas_fixed_, lambdas_nests, likelihood]





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





/Users/carlostrujillo/Documents/GitHub/pymc-nb-review/nested_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]











  
    

      
      policy_share
      new_policy_share
      relative_change
    

    

      product
      
      
      
    

  
  
    

      gc
      0.635503
      0.792799
      0.247514
    

    

      ec
      0.071248
      0.088894
      0.247671
    

    

      gr
      0.150178
      0.000000
      -1.000000
    

    

      er
      0.087495
      0.118307
      0.352158
    

    

      hp
      0.055576
      0.000000
      -1.000000
    

  











A Different Market#


Let’s now briefly look at a different market that highlights a limitation of the nested logit model.






data_path = data_dir / "choice_crackers.csv"
df_new = pd.read_csv(data_path)
last_chosen = pd.get_dummies(df_new["lastChoice"]).drop("private", axis=1).astype(int)
last_chosen.columns = [col + "_last_chosen" for col in last_chosen.columns]
df_new[last_chosen.columns] = last_chosen
df_new













  
    

      
      personId
      disp_sunshine
      disp_keebler
      disp_nabisco
      disp_private
      feat_sunshine
      feat_keebler
      feat_nabisco
      feat_private
      price_sunshine
      price_keebler
      price_nabisco
      price_private
      choice
      lastChoice
      personChoiceId
      choiceId
      keebler_last_chosen
      nabisco_last_chosen
      sunshine_last_chosen
    

  
  
    

      0
      1
      0
      0
      0
      0
      0
      0
      0
      0
      0.99
      1.09
      0.99
      0.71
      nabisco
      nabisco
      1
      1
      0
      1
      0
    

    

      1
      1
      1
      0
      0
      0
      0
      0
      0
      0
      0.49
      1.09
      1.09
      0.78
      sunshine
      nabisco
      2
      2
      0
      1
      0
    

    

      2
      1
      0
      0
      0
      0
      0
      0
      0
      0
      1.03
      1.09
      0.89
      0.78
      nabisco
      sunshine
      3
      3
      0
      0
      1
    

    

      3
      1
      0
      0
      0
      0
      0
      0
      0
      0
      1.09
      1.09
      1.19
      0.64
      nabisco
      nabisco
      4
      4
      0
      1
      0
    

    

      4
      1
      0
      0
      0
      0
      0
      0
      0
      0
      0.89
      1.09
      1.19
      0.84
      nabisco
      nabisco
      5
      5
      0
      1
      0
    

    

      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
    

    

      3151
      136
      0
      0
      0
      0
      0
      0
      0
      0
      1.09
      1.19
      0.99
      0.55
      private
      private
      9
      3152
      0
      0
      0
    

    

      3152
      136
      0
      0
      0
      1
      0
      0
      0
      0
      0.78
      1.35
      1.04
      0.65
      private
      private
      10
      3153
      0
      0
      0
    

    

      3153
      136
      0
      0
      0
      0
      0
      0
      0
      0
      1.09
      1.17
      1.29
      0.59
      private
      private
      11
      3154
      0
      0
      0
    

    

      3154
      136
      0
      0
      0
      0
      0
      0
      0
      0
      1.09
      1.22
      1.29
      0.59
      private
      private
      12
      3155
      0
      0
      0
    

    

      3155
      136
      0
      0
      0
      0
      0
      0
      0
      0
      1.29
      1.04
      1.23
      0.59
      private
      private
      13
      3156
      0
      0
      0
    

  



3156 rows × 20 columns










utility_formulas = [
    "sunshine ~ disp_sunshine + feat_sunshine + price_sunshine ",
    "keebler ~ disp_keebler + feat_keebler + price_keebler  ",
    "nabisco ~ disp_nabisco + feat_nabisco + price_nabisco  ",
    "private ~ disp_private + feat_private + price_private  ",
]


nesting_structure = {
    "private": ["private"],
    "brand": ["keebler", "sunshine", "nabisco"],
}


nstL_3 = NestedLogit(
    choice_df=df_new,
    utility_equations=utility_formulas,
    depvar="choice",
    covariates=["disp", "feat", "price"],
    nesting_structure=nesting_structure,
    model_config={
        "alphas_": Prior("Normal", mu=0, sigma=1, dims="alts"),
        "betas": Prior("Normal", mu=0, sigma=1, dims="alt_covariates"),
        "lambdas_nests": Prior("Beta", alpha=1, beta=1, dims="nests"),
    },
    alphas_nests=False,
)
nstL_3









<pymc_marketing.customer_choice.nested_logit.NestedLogit at 0x12e1c0190>













nstL_3.sample(
    fit_kwargs={
        "target_accept": 0.95,
        "tune": 2000,
        "nuts_sampler": "nutpie",
        "progressbar": True,
    }
)









Sampling: [alphas_, betas, lambdas_nests, likelihood]





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







Sampling: [likelihood]







<pymc_marketing.customer_choice.nested_logit.NestedLogit at 0x12e1c0190>













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









/Users/carlostrujillo/Documents/GitHub/pymc-nb-review/nested_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/nested_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[sunshine]
      -0.98
      0.6
      -1.9
      -0.0099
      436
      459
      1.01
      0.029
      0.021
    

    

      alphas[keebler]
      -0.41
      0.6
      -1.4
      0.56
      437
      505
      1.01
      0.029
      0.021
    

    

      alphas[nabisco]
      1.35
      0.59
      0.4
      2.3
      427
      473
      1.01
      0.029
      0.022
    

    

      alphas[private]
      0
      0
      0
      0
      4000
      4000
      NaN
      0
      NaN
    

    

      betas[disp]
      0.098
      0.09
      -0.047
      0.24
      4116
      2905
      1.00
      0.0014
      0.00097
    

    

      betas[feat]
      0.336
      0.146
      0.1
      0.57
      4185
      2796
      1.00
      0.0023
      0.0016
    

    

      betas[price]
      -3.6
      0.327
      -4.1
      -3.1
      3307
      2601
      1.00
      0.0057
      0.004
    

  











pc = azp.plot_forest(
    nstL_3.idata,
    var_names=["alphas", "betas"],
    combined=True,
)
ax = pc.viz["/"]["figure"].values.item().axes[0]
ax.axvline(0, color="k", linestyle="--")
ax.set_title("Parameter Estimates from the Nested Logit Model");








../../_images/f19032701e9fc005f89ecdd3e64944429829b261ce965a03d357a46dd7fe33d7.png






This model suggests plausibly that price increases have a strongly negative effect of the purchase probability of crackers. However, we should note that we’re ignoring information when we use the nested logit in this context. The data set records multiple purchasing decisions for each individual and by failing to model this fact we lose insight into individual heterogeneity in their responses to price shifts. To gain more insight into this aspect of the purchasing decisions of individuals we might augment our model with hierarchical components or switch to alternative Mixed logit model.






Choosing a Market Structure#


Causal inference is hard and predicting the counterfactual actions of agents in a competitive market is very hard. There is no guarantee that a nested logit model will accurately represent the choice of any particular agent, but you can be hopeful that it highlights expected patterns of choice when the nesting structure reflects the natural segmentation of a market. Nests should group alternatives that share unobserved similarities, ideally driven by a transparent theory of the market structure. Well-specified nests should show stronger substitution within nests than across nests, and you can inspect the substitution patterns as above. Ideally you can always try and hold out some test data to evaluate the implications of your fitted model. Discrete choice models are causal inference models and their structural specification should support generalisable inference across future and counterfactual situations.






%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
arviz_plots   : 1.2.0
matplotlib    : 3.10.9
pandas        : 2.3.3
pymc_extras   : 0.12.2.dev1+gee8cc37df
pymc_marketing: 1.0.0.dev0

Watermark: 2.6.0