sim_nested_data <- function(N, T, beta, lambda) {
J <- 3
n_rows <- N * T * J
df <- data.frame(
id = rep(seq_len(N), each = T * J),
task = rep(rep(seq_len(T), each = J), times = N),
alt = rep(c("Acer", "Dell", "Apple"), times = N * T),
ram = factor(sample(c("8","16","32"), n_rows, replace = TRUE),
levels = c("8","16","32")),
screen = factor(sample(c("13","15"), n_rows, replace = TRUE),
levels = c("13","15")),
price = round(runif(n_rows, min = 0.8, max = 2.4), 2)
)
df$alt <- factor(df$alt, levels = c("Acer", "Dell", "Apple"))
df$nest <- ifelse(df$alt == "Apple", "mac", "win")
X <- model.matrix(~ alt + ram + screen + price, data = df)[, -1]
V <- as.vector(X %*% beta)
task_id <- cumsum(!duplicated(df[, c("id", "task")]))
## nested logit probabilities, task by task (eq-nested-prob)
P <- numeric(n_rows)
for (i in unique(task_id)) {
rows <- which(task_id == i)
v <- V[rows]; nest <- df$nest[rows]
Iv <- tapply(v / lambda, nest, function(z) log(sum(exp(z))))
Pnest <- exp(lambda * Iv) / sum(exp(lambda * Iv))
Pwithin <- exp(v / lambda) / ave(exp(v / lambda), nest, FUN = sum)
P[rows] <- Pnest[nest] * Pwithin
}
## choices: one categorical draw per task
df$choice <- 0L
for (i in unique(task_id)) {
rows <- which(task_id == i)
df$choice[rows[sample.int(J, 1, prob = P[rows])]] <- 1L
}
df
}
set.seed(103)
beta_true <- c(dell = 0.5, apple = 1.0, ram16 = 0.6, ram32 = 0.9,
screen15 = 0.3, price = -1.2)
lambda_true <- 0.5
nested <- sim_nested_data(N = 500, T = 8, beta = beta_true,
lambda = lambda_true)10 Relaxing IIA: Nested Logit and GEV Models
The MNL bought its closed form with the IIA property (Equation 7.5), and Chapter 7 showed the bill can be absurd: paint a bus red and steal car riders. This chapter makes the first — and gentlest — repair. The nested logit lets alternatives share unobserved appeal within researcher-specified groups, breaking IIA across groups while preserving a closed-form probability, a tractable likelihood, and everything we know about estimating such things. It is the smallest possible step beyond the MNL, which makes it the ideal place to learn how model generalizations propagate through the estimation machinery: a new parameter, a new likelihood, the same workflow.
The chapter also marks a small narrative shift in the running example, explained below, and introduces a technique we will use for the rest of the book: simulating data directly from choice probabilities when utility-plus-error simulation is awkward.
10.1 When Proportional Substitution Fails
Recall the diagnosis from Chapter 7: the MNL’s independent errors give it no way to know that two alternatives are similar in unobserved ways, so improvement in one alternative draws share from all others in proportion to their shares — never disproportionately from its close substitutes.
Put that in laptop terms. Suppose unobserved utility has a large operating-system component: some shoppers, on some occasions, are simply in a macOS mood or locked into Windows workflows we don’t measure. Then a Dell and an Acer (both Windows) share unobserved appeal that neither shares with an Apple. If the Dell cuts its price, common sense says it should poach mostly Windows shoppers — Acer’s customers — and pry away few committed Mac buyers. The MNL is structurally incapable of that asymmetry; the nested logit is built for it.
10.2 The Nested Logit Model
Partition the alternatives in a choice situation into \(M\) nests \(B_1, \ldots, B_M\). Our study: \(B_{\text{win}} = \{\text{Acer}, \text{Dell}\}\), \(B_{\text{mac}} = \{\text{Apple}\}\) (a nest may be a singleton). Utility keeps its form \(U_{ntj} = V_{ntj} + \varepsilon_{ntj}\), but the joint distribution of the errors is generalized so that errors within a nest are correlated while errors across nests remain independent. The degree of within-nest correlation is governed by the dissimilarity parameter \(\lambda_m \in (0, 1]\): the within-nest error correlation is approximately \(1 - \lambda_m^2\) — a standard approximation we take from the literature rather than derive, since the exact correlation has no simple closed form (Train 2009, sec. 4.2) — so
- \(\lambda_m = 1\): no correlation — the nest is cosmetic, and the model collapses to the MNL;
- \(\lambda_m \to 0\): perfect correlation — alternatives in the nest are near-clones on unobservables, and competition within the nest becomes ferocious while the nest as a whole competes with others as a unit.
Values outside \((0,1]\) break consistency with random utility maximization,1 a fact with practical bite: an estimated \(\hat\lambda\) above 1 is a specification alarm (the nesting is wrong), and our estimation code below will constrain the search to the valid interval. For this chapter we give both nests a common \(\lambda\) — the natural choice when one nest is a singleton, whose own \(\lambda\) would be unidentified.
The choice probability (derived in Train 2009, sec. 4.2) has a structure that teaches the model’s behavioral content. Define the inclusive value of nest \(m\):
\[ I_{ntm} = \log \sum_{j' \in B_m} e^{V_{ntj'} / \lambda} \tag{10.1}\]
Then
\[ P_{ntj} = \underbrace{\frac{e^{\lambda I_{ntm(j)}}}{\sum_{m'=1}^{M} e^{\lambda I_{ntm'}}}}_{\Pr(\text{nest } m(j))} \times \underbrace{\frac{e^{V_{ntj}/\lambda}}{\sum_{j' \in B_{m(j)}} e^{V_{ntj'}/\lambda}}}_{\Pr(j \mid \text{nest } m(j))} \tag{10.2}\]
Read it as a two-stage decision — pick a nest, then pick within the nest — with each stage an MNL. The inclusive value is the bridge between stages: \(I_{ntm}\) is the log-sum of nest \(m\), the same “value of a choice set” object that the Gumbel max-stability property produced in Chapter 7, now measuring how attractive the nest is as a whole when the upper stage compares nests.2 The two-stage reading is a mathematical decomposition, not a claim about sequential cognition; the model is a simultaneous utility maximization with correlated errors that happens to factor this way.
Verify the two limiting cases against Equation 10.2 by inspection: at \(\lambda = 1\), the \(I\)’s become ordinary log-sums and the product collapses to Equation 7.4; as \(\lambda \to 0\), the within-nest stage approaches “pick the nest’s best alternative with certainty” while nests compete by their best members. And confirm the IIA repair by forming the odds for \(j \in B_a\) and \(k \in B_b\) with \(a \neq b\). The across-nest denominator is common to both probabilities and cancels, and using the identity \(\sum_{j' \in B_m} e^{V_{ntj'}/\lambda} = e^{I_{ntm}}\) for each within-nest denominator, \[\frac{P_{ntj}}{P_{ntk}} = \frac{e^{\lambda I_{nta}} \; e^{V_{ntj}/\lambda - I_{nta}}}{e^{\lambda I_{ntb}} \; e^{V_{ntk}/\lambda - I_{ntb}}} = e^{(V_{ntj} - V_{ntk})/\lambda} \; e^{(\lambda - 1)(I_{nta} - I_{ntb})}.\] At \(\lambda = 1\) the second factor is 1 and we are back to Equation 7.5. For \(\lambda < 1\), the odds depend — through \(I_{nta}\) and \(I_{ntb}\) — on every alternative in the two nests involved: improve the Acer and the Dell/Apple odds shift, because \(I_{\text{win}}\) rises. Note the precision of the repair: alternatives in any third nest still cancel exactly, so IIA fails at the boundary between the two nests being compared, survives within nests, and survives with respect to everything else. Choosing the boundary is the researcher’s model specification burden — and the model offers no protection against choosing it wrongly, a limitation the mixed logit (Chapter 13) and probit (Chapter 18) will lift at the cost of the closed form.
10.3 A Labeled Variant of the Study
Nesting is a statement about which alternatives share unobserved appeal, so it attaches most naturally to designs where the alternatives have stable identities. For this chapter we switch the running example to a brand-labeled design: every task shows exactly one Acer, one Dell, and one Apple — alternative 1 is the Acer, and so on — with memory, screen, and price still randomized independently. This is a common real-world CBC variant (Chapter 4 introduced the labeled/unlabeled distinction; here it earns its keep), and it changes the bookkeeping in one satisfying way: the brand dummies of Chapter 7 are now, literally, alternative-specific constants — the ASCs that Section 4.5 said labeled designs require. Same columns, same true values, new interpretation. The utility specification and true parameters carry over from Equation 7.6, joined by the new structural parameter:
\[ \lambda^\ast = 0.5 \]
— substantial within-nest correlation (about \(1 - 0.5^2 = 0.75\)), enough to make the substitution patterns visibly non-MNL.
10.4 Simulating Nested Logit Data
Here we hit an instructive wall. Our simulators so far follow the utility recipe: draw errors, add, argmax (Chapter 3, Section 7.6). For the nested logit, the required error draws — marginally Gumbel, dependent within nests in the precise GEV pattern — have no simple constructive recipe like rlogis() or the Cholesky transform.3 But we don’t need one: any model whose choice probabilities we can compute can be simulated directly from those probabilities. Compute \(P_{ntj}\) for the task via Equation 10.2, then draw the choice from the corresponding categorical distribution — in R, sample() with prob =. The utility framework earns its keep in deriving the probabilities; once derived, they are the model.
(Yes, the simulator loops over tasks — and by Chapter 9’s own doctrine that is fine: it runs once, not inside an optimizer. Measure before optimizing includes knowing when not to bother.)
10.4.1 Seeing the repaired substitution
Before estimating, let’s confirm the model does what it was hired to do. Take a representative task — one laptop per brand, identical specs, $1,500 each — cut the Dell’s price by $300, and ask where its gained share comes from, under the nested model (\(\lambda = 0.5\)) versus the MNL (\(\lambda = 1\)):
nested_probs <- function(V, nest, lambda) {
Iv <- tapply(V / lambda, nest, function(z) log(sum(exp(z))))
Pnest <- exp(lambda * Iv) / sum(exp(lambda * Iv))
Pwithin <- exp(V / lambda) / ave(exp(V / lambda), nest, FUN = sum)
as.vector(Pnest[nest] * Pwithin)
}
X_task <- rbind(acer = c(0, 0, 1, 0, 1, 1.5),
dell = c(1, 0, 1, 0, 1, 1.5),
apple = c(0, 1, 1, 0, 1, 1.5))
nest <- c("win", "win", "mac")
V0 <- as.vector(X_task %*% beta_true)
V1 <- V0 + c(0, beta_true["price"] * -0.3, 0) # Dell $300 cheaper: price falls by 0.3 ($000s)
shares <- rbind(
nested_before = nested_probs(V0, nest, 0.5),
nested_after = nested_probs(V1, nest, 0.5),
mnl_before = nested_probs(V0, nest, 1),
mnl_after = nested_probs(V1, nest, 1)
)
colnames(shares) <- rownames(X_task)
round(shares, 3) acer dell apple
nested_before 0.112 0.303 0.585
nested_after 0.074 0.412 0.514
mnl_before 0.186 0.307 0.506
mnl_after 0.164 0.389 0.447Do the accounting. Under the MNL, the Dell’s gain comes from Acer and Apple in proportion to their prior shares — proportional substitution, as Chapter 7 proved it must. Under the nested model, Acer surrenders share at a far higher rate relative to its starting point than Apple does: the price cut plays out mostly as within-nest warfare among the Windows machines, while the Mac buyer barely glances up. That asymmetry — impossible in the MNL — is the entire point of the \(\lambda\) parameter.
10.5 The Nested Logit Log-Likelihood
Estimation needs \(\log P_{ntj}\) for the chosen alternatives. One identity does all the work: by the definition Equation 10.1, the within-nest denominator is \(\sum_{j' \in B_{m(j)}} e^{V_{ntj'}/\lambda} = e^{I_{ntm(j)}}\). So the log of Equation 10.2’s within-nest factor is \(V_{ntj}/\lambda - I_{ntm(j)}\), the log of its nest-choice factor is \(\lambda I_{ntm(j)} - \log \sum_{m'} e^{\lambda I_{ntm'}}\), and (writing \(m(j)\) for the chosen row’s nest) their sum is:
\[ \log P_{ntj} = \underbrace{\left[ \frac{V_{ntj}}{\lambda} - I_{ntm(j)} \right]}_{\text{within-nest}} \; + \; \underbrace{\left[ \lambda I_{ntm(j)} - \log \textstyle\sum_{m'} e^{\lambda I_{ntm'}} \right]}_{\text{across-nest}} \tag{10.3}\]
Both bracketed terms are “utility minus log-sum” shapes — the MNL pattern of Equation 8.2 at two levels. So the vectorized implementation is the Chapter 8 machinery applied twice: once with rows grouped by (task, nest) to build inclusive values, once with nests grouped by task. The parameter vector stacks the two: \(\boldsymbol{\theta}= (\boldsymbol{\beta}, \lambda)\) — theta in the code, unpacked by beta <- theta[1:K] and lambda <- theta[K + 1]; each level gets its own max-shift for numerical stability (the log-sum-exp guard of Section 8.3, applied twice):
loglik_nested <- function(theta, y, X, task_id, nest_id) {
K <- ncol(X)
beta <- theta[1:K]
lambda <- theta[K + 1]
Vs <- as.vector(X %*% beta) / lambda
## level 1: inclusive value per (task, nest) group
g <- as.integer(interaction(task_id, nest_id, drop = TRUE))
m1_row <- ave(Vs, g, FUN = max)
m1 <- as.vector(tapply(m1_row, g, `[`, 1))
I_g <- log(as.vector(rowsum(exp(Vs - m1_row), g))) + m1
task_of_g <- as.vector(tapply(task_id, g, `[`, 1))
## level 2: log-sum of lambda * I across nests within task
lI <- lambda * I_g
m2_g <- ave(lI, task_of_g, FUN = max)
m2 <- as.vector(tapply(m2_g, task_of_g, `[`, 1))
lse2 <- log(as.vector(rowsum(exp(lI - m2_g), task_of_g))) + m2
## eq-nested-loglik-term, chosen rows only
logP <- (Vs - I_g[g]) + (lambda * I_g[g] - lse2[task_id])
sum(logP[y == 1])
}The new idiom is interaction(task_id, nest_id, drop = TRUE): it mints one group label per (task, nest) pair, letting rowsum() and ave() aggregate at the intermediate level, after which I_g[g] and lse2[task_id] broadcast group-level quantities back down — the same label-indexing trick as always, used at two altitudes. Check it the honest way, against a brute-force task-by-task computation of Equation 10.2 at the true parameters:
X <- model.matrix(~ alt + ram + screen + price, data = nested)[, -1]
y <- nested$choice
task_id <- cumsum(!duplicated(nested[, c("id", "task")]))
ll_fast <- loglik_nested(c(beta_true, lambda_true), y, X, task_id, nested$nest)
V <- as.vector(X %*% beta_true); ll_bf <- 0
for (i in unique(task_id)) {
rows <- which(task_id == i)
P <- nested_probs(V[rows], nested$nest[rows], lambda_true)
ll_bf <- ll_bf + log(P[y[rows] == 1])
}
c(vectorized = ll_fast, brute_force = ll_bf) vectorized brute_force
-3216.614 -3216.614 Identical. Note in passing what we just did: the simulator and the likelihood were implemented independently (probability formula task-by-task versus grouped log-sum algebra), so their agreement is a genuine cross-check, not a tautology — the same reason Chapter 8 validated against mlogit rather than against itself.
10.6 Estimation
One wrinkle separates this fit from Chapter 8: \(\lambda\) must stay in \((0, 1]\), and BFGS knows nothing of bounds. Three standard options: reparameterize (estimate \(\alpha\) with \(\lambda = 1/(1+e^{-\alpha})\), unconstrained), penalize, or use a bounded optimizer. We take the third — optim(method = "L-BFGS-B") accepts box constraints directly, and meeting it here matters because bounded parameters recur (class shares in Chapter 11, variance parameters in Chapter 13, where we will choose reparameterization instead and compare notes):
fit <- optim(par = c(rep(0, 6), 0.8),
fn = loglik_nested,
y = y, X = X, task_id = task_id, nest_id = nested$nest,
method = "L-BFGS-B",
lower = c(rep(-Inf, 6), 0.05),
upper = c(rep( Inf, 6), 1),
control = list(fnscale = -1),
hessian = TRUE)
se <- sqrt(diag(solve(-fit$hessian)))
results <- data.frame(truth = c(beta_true, lambda = lambda_true),
estimate = fit$par, se = se, z = fit$par / se)
rownames(results) <- c(colnames(X), "lambda")
round(results, 3) truth estimate se z
altDell 0.5 0.543 0.040 13.674
altApple 1.0 1.042 0.059 17.653
ram16 0.6 0.642 0.049 13.220
ram32 0.9 0.927 0.053 17.502
screen15 0.3 0.288 0.036 7.982
price -1.2 -1.225 0.059 -20.829
lambda 0.5 0.510 0.033 15.238Recovery across the board, including the structural parameter: \(\hat\lambda \approx 0.51\) with a tight standard error, decisively inside \((0,1)\) and decisively away from the MNL’s \(\lambda = 1\). Make that last comparison formal with the likelihood-ratio machinery of Chapter 8 — the MNL is this model with one parameter pinned, so:
loglik_mnl <- function(beta, y, X, task_id) { # from @sec-mnl-mle
V <- as.vector(X %*% beta)
Vmax <- ave(V, task_id, FUN = max)
eV <- exp(V - Vmax)
sum((V - Vmax - log(rowsum(eV, task_id)[task_id]))[y == 1])
}
fit_mnl <- optim(rep(0, 6), loglik_mnl, y = y, X = X, task_id = task_id,
method = "BFGS", control = list(fnscale = -1))
lr <- 2 * (fit$value - fit_mnl$value)
c(LR = round(lr, 1), p = signif(pchisq(lr, df = 1, lower.tail = FALSE), 3)) LR p
1.014e+02 7.660e-24 The data, simulated with real nest correlation, reject the MNL emphatically.4 Estimating a wrong-by-construction MNL on these data is itself instructive — compare fit_mnl$par to the truth and you’ll find the attribute coefficients survive respectably while the substitution predictions (the model’s whole purpose in a competitive analysis) would be badly off: misspecification often hides in the quantities you don’t tabulate.
10.7 Validation Against mlogit
mlogit estimates nested logits via its nests argument; un.nest.el = TRUE imposes our common-\(\lambda\) assumption. With labeled alternatives, mlogit supplies the ASCs itself (its (Intercept):Dell and (Intercept):Apple are our dell and apple dummies — same reference brand, Acer):
library(dfidx)
library(mlogit)
md <- nested
md$chid <- task_id
md$choice <- as.logical(md$choice)
md <- dfidx(md, idx = list(c("chid", "id"), "alt"))
fit_ml <- mlogit(choice ~ ram + screen + price, data = md,
nests = list(win = c("Acer", "Dell"), mac = "Apple"),
un.nest.el = TRUE)
round(rbind(ours = fit$par, mlogit = coef(fit_ml)[c(1:2, 3:6, 7)]), 4) (Intercept):Dell (Intercept):Apple ram16 ram32 screen15 price iv
ours 0.5434 1.0421 0.642 0.9267 0.2878 -1.2252 0.5098
mlogit 0.5434 1.0421 0.642 0.9267 0.2878 -1.2252 0.5098c(ours = fit$value, mlogit = as.numeric(logLik(fit_ml))) ours mlogit
-3215.509 -3215.509 Coefficient-for-coefficient agreement (their iv is our \(\lambda\)), identical log-likelihoods. Two implementations, one summit, again.
10.8 The GEV Family
The nested logit is one member of a family McFadden constructed with a single generating device (McFadden 2001; Train 2009, Ch. 4). Take any function \(G(e^{V_1}, \ldots, e^{V_J})\) satisfying a short list of properties (non-negativity, homogeneity, sign-alternating derivatives); then
\[ P_j = \frac{e^{V_j} \, G_j(e^{V_1}, \ldots, e^{V_J})}{G(e^{V_1}, \ldots, e^{V_J})} \tag{10.4}\]
(with \(G_j = \partial G / \partial e^{V_j}\)) is the choice probability of a well-defined random utility model with generalized extreme value errors — closed-form, RUM-consistent, guaranteed. The family is a probability factory: \(G = \sum_j e^{V_j}\) yields the MNL; \(G = \sum_m \big( \sum_{j \in B_m} e^{V_j/\lambda_m} \big)^{\lambda_m}\) yields the nested logit (differentiate it and Equation 10.2 falls out); overlapping-nest choices of \(G\) yield the cross-nested and paired combinatorial logits, where an alternative belongs partially to several nests — a Dell desktop-replacement laptop that competes in both the “Windows” nest and a “large-screen” nest.
We won’t estimate these siblings — the workflow is exactly the one this chapter just executed: derive \(P\) from \(G\), take logs, write the two-(or three-)level grouped likelihood, constrain the \(\lambda\)’s, climb — and one worked example plus a recipe beats three worked examples. The strategic lesson to carry forward instead: within the GEV family, every substitution pattern must be specified in advance through the structure of \(G\). The researcher declares who competes with whom; the data only measure how strongly. When you do not know the competitive structure — the usual situation — you want models where flexible substitution emerges from estimated quantities instead: random coefficients (Chapter 13) or free error covariance (Chapter 18). Both cost the closed form. The next four chapters pay that cost and build the simulation machinery that covers it.
10.9 Key Learnings
- The nested logit groups alternatives into nests with correlated unobservables; the dissimilarity parameter \(\lambda \in (0,1]\) indexes the correlation (\(\lambda = 1\) recovers the MNL), and estimates outside the interval indict the nesting structure itself.
- The probability (Equation 10.2) factors into two MNL stages linked by the inclusive value — the log-sum of Chapter 7 in its second starring role — and IIA now fails across nests while surviving within them, exactly as designed: our $300 Dell price cut savaged Acer and barely touched Apple.
- New simulation technique: when errors have no convenient constructive form, simulate from the choice probabilities directly (categorical draws). Computable probabilities are all a simulator needs.
- The log-likelihood is the MNL’s “utility minus log-sum” pattern applied at two levels;
interaction()group labels plusrowsum()/broadcast implement it vectorized, with max-shifts at both levels. Simulator and likelihood, implemented independently, agreed at the truth — andmlogitagreed with both. - Bounded parameters met bounded optimization:
L-BFGS-Bwith box constraints. The LR test rejected \(\lambda = 1\) decisively (mind the boundary caveat), and recovery included the structural parameter. - The GEV generating function (Equation 10.4) mass-produces closed-form RUM models — but all of them require the researcher to pre-specify who competes with whom. Flexible, data-driven substitution requires integrals without closed forms: enter simulation-assisted estimation.
Train (2009) Section 4.2 discusses the bounds and the (limited) sense in which \(\lambda\) slightly above 1 can be tolerated. The formal foundation — that the nested logit is the GEV member generated by a particular function \(G\) — appears at the chapter’s end.↩︎
The log-sum’s third appearance is scheduled for Chapter 21, where (divided by the price coefficient) it becomes consumer surplus. One object, three jobs: max-stability byproduct, nest attractiveness, welfare measure.↩︎
Constructions exist (Train 2009, sec. 4.2 notes the possibilities), but they are genuinely awkward, and the probability route is both exact and simpler. The probability route also has a conceptual payoff: it separates “the model” (a distribution over choices given attributes) from “the story that motivated it” (utilities and errors).↩︎
A technical footnote for the careful: the null \(\lambda = 1\) sits on the boundary of the parameter space, which makes the usual \(\chi^2_1\) reference conservative in general. At an LR statistic this size, no correction changes the verdict.↩︎