Nested logit¶
Download the script: nested.py
The conditional logit’s substitution pattern is rigid: when an alternative
gets worse, its choosers go to the others in proportion to their shares (the
independence of irrelevant alternatives). The nested logit (experimental)
groups alternatives that are closer substitutes into nests, within which
their unobserved utilities are correlated, so that choosers who leave one
alternative go disproportionately to its nest-mates. nlogit_stream takes
the same long-format data as clogit_stream, a column naming the
alternatives, and the tree of nests.
import polars as pl
from logit_stream import clogit_stream, nlogit_stream
from logit_stream.simulate import simulate_choices
nests = {"rail": ["train", "metro"], "road": ["car"]}
data = (simulate_choices(50_000, nests=nests, lambdas=0.5, seed=1)
.with_columns(income_car=pl.col("income") * (pl.col("alt") == "car")))
FML = "chosen ~ cost + time + C(alt) + income_car"
fit = nlogit_stream(FML, data, case="case", alt="alt", nests=nests, vcov="hetero")
fit.summary()
### chosen ~ cost + time + C(alt) + income_car
Nested logit (logit_stream, Newton) vcov: hetero
choice sets: 50,000 alternatives (rows): 134,988
nests (alt): rail: train, metro; road: car
log-likelihood of the conditional logit: -34,936.848
log-likelihood: -34111.627 null: -48517.828 pseudo R2: 0.296926 AIC: 68235.255 BIC: 68288.173
┌─────────────────┬───────────┬────────────┬────────────┬─────────────┐
│ Coefficient ┆ Estimate ┆ Std. Error ┆ t value ┆ Pr(>|t|) │
╞═════════════════╪═══════════╪════════════╪════════════╪═════════════╡
│ cost ┆ -1.112621 ┆ 0.016045 ┆ -69.343914 ┆ 0.0 │
│ time ┆ -1.328957 ┆ 0.017765 ┆ -74.807139 ┆ 0.0 │
│ C(alt)[T.metro] ┆ 0.703156 ┆ 0.01266 ┆ 55.54271 ┆ 0.0 │
│ C(alt)[T.train] ┆ 0.293196 ┆ 0.013829 ┆ 21.201982 ┆ 9.1558e-100 │
│ income_car ┆ 0.292545 ┆ 0.011669 ┆ 25.070388 ┆ 1.0466e-138 │
│ lambda.rail ┆ 0.514243 ┆ 0.00864 ┆ 59.515987 ┆ 0.0 │
└─────────────────┴───────────┴────────────┴────────────┴─────────────┘
The data are drawn from a nested logit with train and metro in a nest with \(\lambda = 0.5\), and the car alone. The estimates recover it: -1.1 on cost, -1.3 on time, 0.7 and 0.3 for metro and train against the car, 0.3 on income for the car, and \(\lambda_\text{rail} = 0.514\).
The nested and mixed logit refuse collinear covariates rather than drop
them, as the conditional logit does, so a chooser’s characteristic enters
as the columns it should have: income_car here, instead of
income:C(alt) with one of its levels dropped.
The model¶
It uses the utility-consistent normalization (RU2, Stata’s nlogit default
and R mlogit’s unscaled = FALSE). With two levels, alternative \(r\) of nest
\(m\) has probability
with one \(\lambda\) per nest of two or more children (a nest of one has none). \(1 - \lambda_m^2\) is the correlation of the unobserved utilities within nest \(m\). With every \(\lambda = 1\) it is the conditional logit, which is where the fit starts; estimates outside \((0, 1]\), or above the parent nest’s, where the model is not consistent with utility maximization everywhere, get a warning.
The likelihood is not concave, so the fit is the mixed logit’s: each read evaluates the log-likelihood and its analytic gradient at several points, the information is their central differences, and between them the steps are quasi-Newton. The test suite checks it against R’s mlogit: the same estimates, and the BHHH standard errors mlogit reports to six digits (the iid variance here is the inverse of the observed information, as Stata’s), and a three-level fit against the likelihood written out recursively.
Against the conditional logit¶
The conditional logit is \(\lambda = 1\), which the likelihood-ratio test and the Wald test both reject decisively here:
cl = clogit_stream(FML, data, case="case")
lr = 2 * (fit.loglik - cl.loglik) # 1650.4, chi-squared(1)
fit.wald_test(keep="lambda.rail", q=1.0) # statistic 3160.6
Fitted to these data, the conditional logit’s coefficients are wrong too, the cost and time coefficients by about a third:
┌─────────────────┬─────────────┬───────────┐
│ coefficient ┆ conditional ┆ nested │
╞═════════════════╪═════════════╪═══════════╡
│ cost ┆ -1.443622 ┆ -1.112621 │
│ time ┆ -1.796145 ┆ -1.328957 │
│ C(alt)[T.metro] ┆ 0.680443 ┆ 0.703156 │
│ C(alt)[T.train] ┆ 0.087421 ┆ 0.293196 │
│ income_car ┆ 0.327711 ┆ 0.292545 │
└─────────────────┴─────────────┴───────────┘
Substitution¶
Raising metro fares by half, where do the metro’s lost riders go? Predicting the shares before and after with each model:
dearer = data.with_columns(
cost=pl.when(pl.col("alt") == "metro").then(1.5 * pl.col("cost")).otherwise("cost"))
fit.predict(dearer)
share of the metro's loss going to each alternative
┌─────────────┬──────────┬──────────┐
│ model ┆ train ┆ car │
╞═════════════╪══════════╪══════════╡
│ conditional ┆ 0.546189 ┆ 0.453811 │
│ nested ┆ 0.65264 ┆ 0.34736 │
└─────────────┴──────────┴──────────┘
The conditional logit splits them roughly by the train’s and car’s shares;
the nested logit sends more of them to the train, the metro’s nest-mate.
predict() and resid() work as for the conditional logit; the
probabilities need each whole choice set, so Polars evaluates them in
memory.
Deeper trees¶
A nest of nests is a dict within the dict. With four alternatives, public transport {rail {train, metro}, bus} against the car:
tree = {"public": {"rail": ["train", "metro"], "bus": ["bus"]}, "private": ["car"]}
four = simulate_choices(50_000, alternatives=("train", "metro", "bus", "car"), nests=tree,
lambdas={"public": 0.7, "rail": 0.4}, seed=2)
deep = nlogit_stream(FML, four, case="case", alt="alt", nests=tree)
┌─────────────────┬───────────┬────────────┐
│ Coefficient ┆ Estimate ┆ Std. Error │
╞═════════════════╪═══════════╪════════════╡
│ cost ┆ -1.098435 ┆ 0.016694 │
│ time ┆ -1.28752 ┆ 0.019098 │
│ C(alt)[T.car] ┆ -0.346578 ┆ 0.016252 │
│ C(alt)[T.metro] ┆ 0.37017 ┆ 0.011904 │
│ C(alt)[T.train] ┆ -0.044541 ┆ 0.011175 │
│ income_car ┆ 0.317571 ┆ 0.012138 │
│ lambda.public ┆ 0.698161 ┆ 0.012858 │
│ lambda.rail ┆ 0.386346 ┆ 0.007534 │
└─────────────────┴───────────┴────────────┘
log-likelihoods: two levels -44793.3, three -44247.6: LR 1091.4, p = 2.5e-239
Deeper trees repeat the two-level step: a nest’s value in its parent is \(\lambda_m I_m\), and an alternative’s probability is the product of the conditional probabilities along its path. The two-level tree with all public transport in one nest is the three-level one with \(\lambda_\text{rail} = \lambda_\text{public}\), rejected here. The bus is a nest of one, with no \(\lambda\) of its own. Every alternative must be in exactly one nest, a nest holds alternatives or nests but not both, and the top level needs at least two nests.
Covariates of a nest¶
Stata’s nlogit takes variables at each level of the tree; here every
covariate enters the alternatives’ utilities, and under RU2 that is the same
model. A term \(z_m \gamma\) that is the same for every alternative of nest
\(m\) leaves each conditional probability within the nest unchanged (it
cancels from \(V_r/\lambda_m - I_m\)) and adds exactly \(z_m \gamma\) to the
nest’s value \(\lambda_m I_m\) – which is what a nest-level variable with
coefficient \(\gamma\) does. So a nest-level variable is a column constant
within each set’s nest, entered as usual: a nest-specific constant is an
indicator of the nest, and a chooser’s characteristic with an effect on the
nest is its product with that indicator (from a column naming each row’s
nest; not <alt>_nest, the name the fit uses for its own). (Under RU1, the
other normalization, the two differ by the factor \(\lambda_m\).)
Inference, weights and offsets¶
vcov= is "iid" (the inverse of the observed information), "hetero"
(the sandwich with the choice sets’ scores) or clustered, {"CRV1": var}
(multi-way "a+b"), with the clusters constant within choice sets.
Weights must be constant within sets too. With five choice sets per person
and weights by person:
nlogit_stream(FML, panel, case="case", alt="alt", nests=nests,
vcov={"CRV1": "person"}, weights="w")
offset= adds a known term to each alternative’s utility, as for the
conditional logit. device="cuda" runs the passes on an NVIDIA GPU
(experimental).
The full script: examples/nested.py
"""The nested logit: nests of alternatives that are closer substitutes, the
dissimilarity parameters, a test against the conditional logit, deeper
trees, predictions and substitution.
python examples/nested.py
"""
import polars as pl
from scipy.stats import chi2
from logit_stream import clogit_stream, nlogit_stream
from logit_stream.simulate import simulate_choices
# Mode-choice data drawn from a nested logit: train and metro share a nest
# with dissimilarity parameter 0.5 (their unobserved utilities correlated,
# 1 - 0.5^2 = 0.75), the car is a nest of its own.
nests = {"rail": ["train", "metro"], "road": ["car"]}
data = (simulate_choices(50_000, nests=nests, lambdas=0.5, seed=1)
# the chooser's income affects the car only; a variable constant within
# a set has to be interacted with the alternatives, and the nested
# logit refuses collinear covariates, so make the one column
.with_columns(income_car=pl.col("income") * (pl.col("alt") == "car")))
FML = "chosen ~ cost + time + C(alt) + income_car"
fit = nlogit_stream(FML, data, case="case", alt="alt", nests=nests, vcov="hetero",
verbose=False)
fit.summary()
# ---- against the conditional logit, which is lambda = 1: the likelihood-
# ratio test, and the Wald test of lambda.rail = 1
cl = clogit_stream(FML, data, case="case", verbose=False)
lr = 2 * (fit.loglik - cl.loglik)
print(f"conditional logit {cl.loglik:.1f}, nested {fit.loglik:.1f}: LR {lr:.1f}, "
f"p = {chi2.sf(lr, 1):.2g}")
print(fit.wald_test(keep="lambda.rail", q=1.0))
print(pl.DataFrame({"coefficient": cl.coefnames, "conditional": cl.beta,
"nested": fit.beta[:len(cl.beta)]}))
# ---- substitution: raise metro fares and see where its riders go. The
# conditional logit sends them in proportion to the other shares (IIA); the
# nested logit sends more of them to the train, its nest-mate.
dearer = data.with_columns(
cost=pl.when(pl.col("alt") == "metro").then(1.5 * pl.col("cost")).otherwise("cost"))
def shares(model, frame):
return (model.predict(frame).group_by("alt").agg(pl.col("prob").sum())
.with_columns(pl.col("prob") / pl.col("prob").sum()).collect())
rows = []
for name, model in (("conditional", cl), ("nested", fit)):
now, then = shares(model, data), shares(model, dearer)
d = now.join(then, on="alt", suffix="_after").with_columns(
change=pl.col("prob_after") - pl.col("prob"))
lost = -d.filter(pl.col("alt") == "metro")["change"].item()
for alt in ("train", "car"):
gain = d.filter(pl.col("alt") == alt)["change"].item()
rows.append({"model": name, "to": alt, "share of metro's loss": gain / lost})
print(pl.DataFrame(rows).pivot(on="to", index="model", values="share of metro's loss"))
# ---- a deeper tree: nests of nests. Four alternatives, public transport
# {rail {train, metro}, bus} against the car, with lambda 0.4 in the rail
# nest and 0.7 in the public one.
tree = {"public": {"rail": ["train", "metro"], "bus": ["bus"]}, "private": ["car"]}
four = (simulate_choices(50_000, alternatives=("train", "metro", "bus", "car"), nests=tree,
lambdas={"public": 0.7, "rail": 0.4}, seed=2)
.with_columns(income_car=pl.col("income") * (pl.col("alt") == "car")))
deep = nlogit_stream(FML, four, case="case", alt="alt",
nests=tree, verbose=False)
print(deep.tidy().select("Coefficient", "Estimate", "Std. Error"))
# the two-level tree that puts all public transport in one nest is the
# deeper one with lambda.rail = lambda.public: rejected
flat = nlogit_stream(FML, four, case="case", alt="alt",
nests={"public": ["train", "metro", "bus"], "private": ["car"]},
verbose=False)
lr = 2 * (deep.loglik - flat.loglik)
print(f"log-likelihoods: two levels {flat.loglik:.1f}, three {deep.loglik:.1f}: "
f"LR {lr:.1f}, p = {chi2.sf(lr, 1):.2g}")
# ---- clustering beyond the choice set, and weights constant within sets
panel = (simulate_choices(50_000, nests=nests, lambdas=0.5, tasks=5, seed=3)
.with_columns(income_car=pl.col("income") * (pl.col("alt") == "car"),
w=1.0 + (pl.col("person") % 3)))
cl_fit = nlogit_stream(FML, panel, case="case", alt="alt", nests=nests,
vcov={"CRV1": "person"}, weights="w", verbose=False)
print(cl_fit.tidy().select("Coefficient", "Estimate", "Std. Error"))