Discrete Choice Models: Zero to Hero Tutorial

This comprehensive tutorial takes you from the foundational concepts of Discrete Choice modelling all the way through advanced extensions, assumption testing, heterogeneity analysis, and practical usage within the DataStatPro application. Whether you are a complete beginner or an experienced analyst, this guide is structured to build your understanding step by step.

Table of Contents

1. Prerequisites and Background Concepts
2. What are Discrete Choice Models?
3. The Mathematical Framework
4. Key Assumptions
5. Identification and Causal Inference
6. Binary Choice Models: Logit and Probit
7. Hypothesis Testing and Inference
8. Effect Size Measures
9. Model Fit and Evaluation
10. Diagnostics and Assumption Testing
11. Extensions: Multinomial and Conditional Logit
12. Extensions: Ordered Choice Models
13. Extensions: Nested Logit and Mixed Logit
14. Extensions: Panel Data Discrete Choice
15. Using the Discrete Choice Component
16. Computational and Formula Details
17. Worked Examples
18. Common Mistakes and How to Avoid Them
19. Troubleshooting
20. Quick Reference Cheat Sheet

1. Prerequisites and Background Concepts

Before diving into Discrete Choice Models, it is helpful to be familiar with the following foundational concepts. Do not worry if you are not — each concept is briefly explained here.

1.1 Random Variables and Probability Distributions

A random variable $Y$ is a variable whose value is determined by a random process. In discrete choice modelling, the outcome variable $Y$ takes on a finite set of values representing alternative choices (e.g., $Y \in \{0, 1\}$ for binary outcomes, or $Y \in \{1, 2, \dots, J\}$ for multiple alternatives).

Key distributions used in discrete choice models:

• Bernoulli distribution: $Y \sim \text{Bernoulli}(p)$, where $P(Y = 1) = p$ and $P(Y = 0) = 1 - p$. Used for binary outcomes.
• Logistic (Gumbel) distribution: The basis of the logit model. The standard logistic CDF is $F(x) = \Lambda(x) = \frac{e^x}{1 + e^x} = \frac{1}{1 + e^{-x}}$.
• Standard Normal distribution: The basis of the probit model. The CDF is $F(x) = \Phi(x) = \int_{-\infty}^x \frac{1}{\sqrt{2\pi}}e^{-t^2/2}dt$.
• Type I Extreme Value (Gumbel) distribution: The basis of the multinomial logit model. Its CDF is $F(x) = e^{-e^{-x}}$.

1.2 Likelihood and Maximum Likelihood Estimation

Maximum Likelihood Estimation (MLE) is the primary estimation method for discrete choice models. For a sample of $n$ independent observations, the likelihood function is:

$\mathcal{L}(\boldsymbol{\theta}; \mathbf{y}, \mathbf{X}) = \prod_{i=1}^n f(y_i \mid \mathbf{x}_i; \boldsymbol{\theta})$

The log-likelihood (more convenient for optimisation) is:

$\ell(\boldsymbol{\theta}) = \ln \mathcal{L}(\boldsymbol{\theta}) = \sum_{i=1}^n \ln f(y_i \mid \mathbf{x}_i; \boldsymbol{\theta})$

The MLE $\hat{\boldsymbol{\theta}}$ maximises $\ell(\boldsymbol{\theta})$:

$\hat{\boldsymbol{\theta}}_{MLE} = \arg\max_{\boldsymbol{\theta}} \ell(\boldsymbol{\theta})$

MLE has attractive properties under regularity conditions: it is consistent, asymptotically normal, and asymptotically efficient.

1.3 Latent Variable Models

Many discrete choice models are derived from an underlying latent variable (unobserved continuous variable). Define a latent utility or propensity:

$Y_i^* = \mathbf{x}_i^T\boldsymbol{\beta} + \epsilon_i$

The observed discrete outcome is a deterministic function of the latent variable:

$Y_i = \begin{cases} 1 & \text{if } Y_i^* > 0 \\ 0 & \text{if } Y_i^* \leq 0 \end{cases}$

The distribution assumed for $\epsilon_i$ determines the model family:

• $\epsilon_i \sim \text{Logistic}(0,1)$ → Logit model
• $\epsilon_i \sim \mathcal{N}(0, 1)$ → Probit model

1.4 Random Utility Maximisation (RUM)

The Random Utility Model (McFadden, 1974) provides the economic foundation for discrete choice. Each decision-maker $i$ assigns a utility $U_{ij}$ to each alternative $j$:

$U_{ij} = V_{ij} + \epsilon_{ij}$

Where:

• $V_{ij} = \mathbf{x}_{ij}^T\boldsymbol{\beta}$ is the systematic (deterministic) utility component.
• $\epsilon_{ij}$ is the random utility component — unobserved factors affecting utility.

The decision-maker chooses the alternative that maximises utility:

$Y_i = j \iff U_{ij} > U_{ik} \quad \forall k \neq j$

Different distributional assumptions for $\epsilon_{ij}$ yield different discrete choice model families.

1.5 Ordinary Least Squares and Its Limitations for Discrete Outcomes

OLS regression applied to a binary outcome ($Y \in \{0,1\}$) produces the Linear Probability Model (LPM):

$E[Y_i \mid \mathbf{x}_i] = \mathbf{x}_i^T\boldsymbol{\beta}$

The LPM has well-known limitations:

• Predicted probabilities outside $[0,1]$: OLS can predict probabilities less than 0 or greater than 1.
• Heteroscedastic errors: $\text{Var}(\epsilon_i) = p_i(1-p_i)$ varies with $\mathbf{x}_i$, violating the OLS homoscedasticity assumption.
• Non-constant marginal effects: The true relationship between covariates and $P(Y=1)$ is typically non-linear (sigmoid-shaped), not linear.

Discrete choice models address all three limitations by modelling probabilities through a monotone transformation that maps the real line to $(0,1)$.

1.6 Multinomial Outcomes and Ordinal Data

Multinomial outcomes have more than two unordered categories: $Y \in \{1, 2, \dots, J\}$ where no natural ordering exists (e.g., mode of transport: car, bus, train, bicycle).

Ordinal outcomes have more than two categories with a natural ordering: $Y \in \{1, 2, \dots, J\}$ where $1 < 2 < \dots < J$ (e.g., satisfaction: low, medium, high).

Different model families are designed for each type of outcome.

2. What are Discrete Choice Models?

2.1 The Core Idea

Discrete Choice Models (DCMs) are statistical models designed to explain and predict the choices made by individuals (or firms, households, or other decision-making units) when they face a finite set of mutually exclusive alternatives.

The core modelling challenge: the outcome is not a continuous variable but a discrete category. Standard regression is misspecified for such outcomes. DCMs model the probability of choosing each alternative as a function of:

• Decision-maker characteristics: Age, income, gender, education of the chooser.
• Alternative-specific attributes: Price, quality, travel time, distance of each option.
• Contextual factors: Market conditions, constraints, availability.

The general DCM probability structure:

$P(Y_i = j \mid \mathbf{x}_i) = F_j(\mathbf{x}_i, \boldsymbol{\theta})$

Where $F_j(\cdot)$ is a function mapping covariates and parameters to probabilities, with:

$\sum_{j=1}^J P(Y_i = j \mid \mathbf{x}_i) = 1 \quad \text{and} \quad P(Y_i = j \mid \mathbf{x}_i) \in (0,1)$

2.2 A Taxonomy of Discrete Choice Models

2.3 Real-World Applications

Discrete choice models are applied across virtually every field involving individual decision-making:

• Transportation: Modal choice (car vs. bus vs. train vs. bicycle), route choice, departure time choice.
• Health Economics: Insurance plan choice, treatment adoption, healthcare provider selection, smoking/drinking behaviour.
• Marketing: Brand choice, product adoption, willingness to pay for product attributes.
• Labour Economics: Occupational choice, labour force participation, migration decisions, unionisation.
• Environmental Economics: Willingness to pay for environmental goods, habitat choice, recreational site choice.
• Political Science: Voting behaviour, party affiliation, referendum outcomes.
• Housing: Residential location choice, tenure choice (rent vs. own), housing type choice.
• Finance: Portfolio allocation categories, credit/default behaviour, investment vehicle choice.
• Education: School choice, field of study selection, dropout decisions.

2.4 Discrete Choice Models vs. Other Regression Methods

3. The Mathematical Framework

3.1 The Binary Logit Model

The logit model specifies the probability of the outcome $Y_i = 1$ as:

$P(Y_i = 1 \mid \mathbf{x}_i) = \Lambda(\mathbf{x}_i^T\boldsymbol{\beta}) = \frac{e^{\mathbf{x}_i^T\boldsymbol{\beta}}}{1 + e^{\mathbf{x}_i^T\boldsymbol{\beta}}} = \frac{1}{1 + e^{-\mathbf{x}_i^T\boldsymbol{\beta}}}$

And:

$P(Y_i = 0 \mid \mathbf{x}_i) = 1 - \Lambda(\mathbf{x}_i^T\boldsymbol{\beta}) = \frac{1}{1 + e^{\mathbf{x}_i^T\boldsymbol{\beta}}}$

The log-odds (logit) transformation linearises the model:

$\ln\left(\frac{P(Y_i=1)}{P(Y_i=0)}\right) = \ln\left(\frac{P_i}{1-P_i}\right) = \mathbf{x}_i^T\boldsymbol{\beta}$

This is the log-odds ratio (or logit), and $\boldsymbol{\beta}$ are the log-odds coefficients.

3.2 The Binary Probit Model

The probit model specifies:

$P(Y_i = 1 \mid \mathbf{x}_i) = \Phi(\mathbf{x}_i^T\boldsymbol{\beta})$

Where $\Phi(\cdot)$ is the standard normal CDF. From the latent variable representation:

$Y_i^* = \mathbf{x}_i^T\boldsymbol{\beta} + \epsilon_i, \quad \epsilon_i \sim \mathcal{N}(0,1)$

$P(Y_i = 1) = P(Y_i^* > 0) = P(\epsilon_i > -\mathbf{x}_i^T\boldsymbol{\beta}) = 1 - \Phi(-\mathbf{x}_i^T\boldsymbol{\beta}) = \Phi(\mathbf{x}_i^T\boldsymbol{\beta})$

3.3 The Logit vs. Probit Comparison

The logistic and standard normal CDFs are very similar in shape. Key differences:

The rule of thumb for converting coefficients: $\hat{\beta}_{logit} \approx \frac{\pi}{\sqrt{3}} \hat{\beta}_{probit} \approx 1.81 \hat{\beta}_{probit}$.

3.4 The Multinomial Logit Model

For $J > 2$ unordered alternatives, the Multinomial Logit (MNL) specifies, for alternative $j \in \{1, \dots, J\}$ with reference category $j = 1$:

$P(Y_i = j \mid \mathbf{x}_i) = \frac{e^{\mathbf{x}_i^T\boldsymbol{\beta}_j}}{\sum_{k=1}^J e^{\mathbf{x}_i^T\boldsymbol{\beta}_k}}$

With the normalisation $\boldsymbol{\beta}_1 = \mathbf{0}$ (reference category), so:

$P(Y_i = j \mid \mathbf{x}_i) = \frac{e^{\mathbf{x}_i^T\boldsymbol{\beta}_j}}{1 + \sum_{k=2}^J e^{\mathbf{x}_i^T\boldsymbol{\beta}_k}}, \quad j = 2, \dots, J$

$P(Y_i = 1 \mid \mathbf{x}_i) = \frac{1}{1 + \sum_{k=2}^J e^{\mathbf{x}_i^T\boldsymbol{\beta}_k}}$

The log-odds ratio relative to the reference category:

$\ln\left(\frac{P(Y_i = j)}{P(Y_i = 1)}\right) = \mathbf{x}_i^T\boldsymbol{\beta}_j$

3.5 The Conditional Logit Model

The Conditional Logit (CL) model (McFadden, 1974) allows attributes to vary across alternatives. The utility of alternative $j$ for individual $i$:

$U_{ij} = \mathbf{z}_{ij}^T\boldsymbol{\gamma} + \mathbf{x}_i^T\boldsymbol{\beta}_j + \epsilon_{ij}$

Where $\mathbf{z}_{ij}$ are alternative-specific attributes (e.g., price of alternative $j$, travel time of option $j$) with a common coefficient $\boldsymbol{\gamma}$, and $\mathbf{x}_i$ are individual-specific characteristics with alternative-specific coefficients $\boldsymbol{\beta}_j$.

The choice probability:

$P(Y_i = j \mid \mathbf{Z}_i) = \frac{e^{\mathbf{z}_{ij}^T\boldsymbol{\gamma} + \mathbf{x}_i^T\boldsymbol{\beta}_j}}{\sum_{k=1}^J e^{\mathbf{z}_{ik}^T\boldsymbol{\gamma} + \mathbf{x}_i^T\boldsymbol{\beta}_k}}$

3.6 The Ordered Logit Model

For an ordinal outcome $Y_i \in \{1, 2, \dots, J\}$, the Ordered Logit (Proportional Odds) model uses a single latent variable:

$Y_i^* = \mathbf{x}_i^T\boldsymbol{\beta} + \epsilon_i, \quad \epsilon_i \sim \text{Logistic}(0,1)$

With $J-1$ threshold (cut-point) parameters $\tau_1 < \tau_2 < \dots < \tau_{J-1}$:

$Y_i = j \iff \tau_{j-1} < Y_i^* \leq \tau_j$

Where $\tau_0 = -\infty$ and $\tau_J = +\infty$. The choice probabilities are:

$P(Y_i = j \mid \mathbf{x}_i) = \Lambda(\tau_j - \mathbf{x}_i^T\boldsymbol{\beta}) - \Lambda(\tau_{j-1} - \mathbf{x}_i^T\boldsymbol{\beta})$

$P(Y_i \leq j \mid \mathbf{x}_i) = \Lambda(\tau_j - \mathbf{x}_i^T\boldsymbol{\beta})$

3.7 The Nested Logit Model

The Nested Logit partitions the $J$ alternatives into $M$ mutually exclusive nests $B_1, B_2, \dots, B_M$. The choice probability for alternative $j$ in nest $m$:

$P(Y_i = j) = P(\text{nest } m) \times P(j \mid \text{nest } m)$

$P(j \mid \text{nest } m) = \frac{e^{V_{ij}/\lambda_m}}{\sum_{k \in B_m} e^{V_{ik}/\lambda_m}}$

$P(\text{nest } m) = \frac{e^{W_{im} + \lambda_m I_{im}}}{\sum_{l=1}^M e^{W_{il} + \lambda_l I_{il}}}$

Where $I_{im} = \ln\sum_{k \in B_m} e^{V_{ik}/\lambda_m}$ is the inclusive value (log-sum), $\lambda_m \in (0,1]$ is the dissimilarity parameter for nest $m$, and $W_{im}$ contains nest-level attributes.

3.8 The Mixed Logit Model

The Mixed Logit (also called the Random Parameters Logit) allows coefficients to vary across individuals:

$P(Y_i = j \mid \mathbf{x}_i, \boldsymbol{\beta}_i) = \frac{e^{\mathbf{x}_{ij}^T\boldsymbol{\beta}_i}}{\sum_{k=1}^J e^{\mathbf{x}_{ik}^T\boldsymbol{\beta}_i}}$

Where $\boldsymbol{\beta}_i \sim f(\boldsymbol{\beta} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma})$ — individual-specific random coefficients drawn from a mixing distribution (typically normal or log-normal).

The unconditional choice probability integrates over the random coefficient distribution:

$P(Y_i = j \mid \mathbf{x}_i) = \int \frac{e^{\mathbf{x}_{ij}^T\boldsymbol{\beta}}}{\sum_{k=1}^J e^{\mathbf{x}_{ik}^T\boldsymbol{\beta}}} f(\boldsymbol{\beta} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}) \, d\boldsymbol{\beta}$

This integral has no closed form and is evaluated by simulation (see Section 16.7).

4. Key Assumptions

4.1 Independence of Irrelevant Alternatives (IIA)

The most important and controversial assumption in multinomial logit models is the Independence of Irrelevant Alternatives (IIA):

The ratio of probabilities for any two alternatives $j$ and $k$ is independent of all other alternatives in the choice set.

Formally, for the MNL model:

$\frac{P(Y_i = j)}{P(Y_i = k)} = \frac{e^{\mathbf{x}_i^T\boldsymbol{\beta}_j}}{e^{\mathbf{x}_i^T\boldsymbol{\beta}_k}} = e^{\mathbf{x}_i^T(\boldsymbol{\beta}_j - \boldsymbol{\beta}_k)}$

This ratio depends only on $j$ and $k$, not on any other alternative $l$.

The Red Bus / Blue Bus Problem: A classic IIA failure. Suppose individuals choose between Car and Red Bus with equal probability (50/50). If a Blue Bus (identical to Red Bus except colour) is added, IIA predicts all three have 1/3 probability each. But intuitively, the split should be 50% car and 50% bus (25% red + 25% blue). IIA allocates "competition" uniformly across all alternatives rather than within similar alternatives.

IIA is implied by: The Type I Extreme Value (Gumbel) distributional assumption on $\epsilon_{ij}$ and the independence across alternatives.

IIA fails when: Alternatives are correlated substitutes — i.e., some alternatives are more similar to each other than to others. In such cases, the error terms $\epsilon_{ij}$ are correlated across alternatives.

4.2 The Proportional Odds Assumption

For the Ordered Logit model, the proportional odds assumption (also called parallel lines assumption) requires that the effect of covariates on the log-odds is constant across all thresholds:

$\ln\left(\frac{P(Y_i > j)}{P(Y_i \leq j)}\right) = \mathbf{x}_i^T\boldsymbol{\beta} - \tau_j \quad \forall j$

The coefficient vector $\boldsymbol{\beta}$ is the same for all outcome categories — only the intercept $\tau_j$ changes. This is a strong assumption that should be explicitly tested (see Section 10.3).

4.3 Random Utility Consistency

For the RUM foundation to be valid:

1. Completeness: Decision-makers can rank all alternatives.
2. Transitivity: Preferences are transitive (if $A \succ B$ and $B \succ C$, then $A \succ C$).
3. Utility maximisation: Decision-makers always choose the alternative with the highest utility.
4. Stable preferences: Preferences do not change during the observation period.

4.4 Independence of Observations

Standard discrete choice models assume independent observations across individuals. In panel data (Section 14), this assumption is relaxed by allowing within-individual correlation across repeated choices.

4.5 Correct Specification of the Choice Set

The model assumes:

• All relevant alternatives are included in the choice set.
• The choice set is the same for all individuals (or, in some extensions, individual-specific choice sets are correctly specified).
• No irrelevant alternatives contaminate the model.

5. Identification and Causal Inference

5.1 What Discrete Choice Models Identify

Identification in discrete choice models means the ability to recover the structural parameters $\boldsymbol{\beta}$ from the data. Key identification conditions:

• Scale normalisation: The scale of the latent utility is unidentified. In the binary probit, $\text{Var}(\epsilon) = 1$ is imposed; in logit, $\text{Var}(\epsilon) = \pi^2/3$ is imposed.
• Location normalisation: Only utility differences are identified, not absolute utility levels. In the MNL, one alternative's coefficient vector is normalised to zero.
• No perfect multicollinearity: Covariates must not be perfectly linearly dependent.
• Exclusion of constants: In models with alternative-specific attributes only, a constant cannot be separately identified from the normalisation.

5.2 Endogeneity in Discrete Choice Models

Endogeneity arises when a regressor $X_{ij}$ is correlated with the unobserved utility component $\epsilon_{ij}$. Common sources:

• Omitted variables: Unobserved factors affecting both the covariate and the choice.
• Simultaneity: The choice affects the covariate (e.g., price determined by anticipated demand).
• Measurement error: Classical measurement error attenuates estimates toward zero.

Consequences: MLE estimates are inconsistent under endogeneity — standard corrections are required.

Remedies:

• Control Function Approach: Add the residuals from an auxiliary regression of the endogenous variable on instruments as an additional regressor in the discrete choice model.
• IV Probit / IV Logit: Two-stage estimation using valid instruments.
• Berry-Levinsohn-Pakes (BLP): For market-level discrete choice with endogenous prices, uses product-level instrumental variables.

5.3 Average Partial Effects vs. Structural Parameters

In discrete choice models, the raw coefficients $\boldsymbol{\beta}$ are not directly interpretable as marginal effects. The Average Partial Effect (APE) of continuous covariate $X_k$ on $P(Y_i = j)$ is:

$APE_k = \frac{1}{n}\sum_{i=1}^n \frac{\partial P(Y_i = j \mid \mathbf{x}_i)}{\partial X_{ik}}$

For binary logit:

$APE_k^{logit} = \frac{1}{n}\sum_{i=1}^n \hat{p}_i(1-\hat{p}_i)\hat{\beta}_k$

For binary probit:

$APE_k^{probit} = \frac{1}{n}\sum_{i=1}^n \phi(\mathbf{x}_i^T\hat{\boldsymbol{\beta}})\hat{\beta}_k$

The APE (also called the Average Marginal Effect, AME) is the primary reported quantity of interest in discrete choice models — analogous to the regression coefficient in OLS.

5.4 Partial Effects at the Mean (PEM) and Partial Effects at Representative Values

Partial Effect at the Mean (PEM): Evaluate the marginal effect at the sample mean $\bar{\mathbf{x}}$:

$PEM_k = \frac{\partial P(Y_i = j \mid \mathbf{x}_i = \bar{\mathbf{x}})}{\partial X_k}$

For binary logit: $PEM_k = \Lambda(\bar{\mathbf{x}}^T\hat{\boldsymbol{\beta}})[1 - \Lambda(\bar{\mathbf{x}}^T\hat{\boldsymbol{\beta}})]\hat{\beta}_k$

⚠️ The PEM evaluates the marginal effect at a potentially non-existent "average individual." The APE is generally preferred because it averages the marginal effect across actual observations, accounting for the non-linearity of the model.

5.5 Willingness to Pay (WTP) in Choice Models

In models with a cost or price attribute (e.g., transport cost, product price), the Willingness to Pay (WTP) for a change in attribute $k$ is:

$WTP_k = -\frac{\hat{\gamma}_k}{\hat{\gamma}_{cost}}$

Where $\hat{\gamma}_k$ is the coefficient on attribute $k$ and $\hat{\gamma}_{cost}$ is the coefficient on cost. This ratio gives the marginal rate of substitution between attribute $k$ and money — a central output of stated preference and transport choice studies.

6. Binary Choice Models: Logit and Probit

6.1 The Log-Likelihood for Binary Models

For a binary outcome $Y_i \in \{0, 1\}$, the log-likelihood is:

$\ell(\boldsymbol{\beta}) = \sum_{i=1}^n \left[Y_i \ln P_i + (1 - Y_i) \ln(1 - P_i)\right]$

Where $P_i = P(Y_i = 1 \mid \mathbf{x}_i)$.

For logit: $P_i = \Lambda(\mathbf{x}_i^T\boldsymbol{\beta})$, giving:

$\ell_{logit}(\boldsymbol{\beta}) = \sum_{i=1}^n \left[Y_i \mathbf{x}_i^T\boldsymbol{\beta} - \ln(1 + e^{\mathbf{x}_i^T\boldsymbol{\beta}})\right]$

For probit: $P_i = \Phi(\mathbf{x}_i^T\boldsymbol{\beta})$, giving:

$\ell_{probit}(\boldsymbol{\beta}) = \sum_{i=1}^n \left[Y_i \ln\Phi(\mathbf{x}_i^T\boldsymbol{\beta}) + (1-Y_i)\ln\Phi(-\mathbf{x}_i^T\boldsymbol{\beta})\right]$

Both log-likelihoods are globally concave, ensuring a unique maximum.

6.2 The Score and Hessian

The score vector (gradient of the log-likelihood):

$\mathbf{s}(\boldsymbol{\beta}) = \frac{\partial \ell}{\partial \boldsymbol{\beta}} = \sum_{i=1}^n \left(Y_i - P_i\right)\mathbf{x}_i$

For logit, this simplifies elegantly because $\partial P_i / \partial \eta_i = P_i(1-P_i)$ where $\eta_i = \mathbf{x}_i^T\boldsymbol{\beta}$.

The Hessian matrix (second derivative):

$\mathbf{H}(\boldsymbol{\beta}) = \frac{\partial^2 \ell}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^T} = -\sum_{i=1}^n w_i \mathbf{x}_i\mathbf{x}_i^T$

Where $w_i = P_i(1-P_i)$ for logit and $w_i = \phi(\mathbf{x}_i^T\boldsymbol{\beta})^2 / [P_i(1-P_i)]$ for probit. The negative Hessian is positive definite, confirming concavity.

6.3 Newton-Raphson and IRLS Estimation

The MLE is obtained iteratively using Newton-Raphson (or equivalently, Iteratively Reweighted Least Squares — IRLS):

$\boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} - \left[\mathbf{H}(\boldsymbol{\beta}^{(t)})\right]^{-1}\mathbf{s}(\boldsymbol{\beta}^{(t)})$

IRLS Interpretation: At each iteration, solve a weighted OLS problem:

$\boldsymbol{\beta}^{(t+1)} = \left(\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X}\right)^{-1}\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{z}^{(t)}$

Where $\mathbf{W}^{(t)} = \text{diag}(w_1^{(t)}, \dots, w_n^{(t)})$ is a diagonal weight matrix and $\mathbf{z}^{(t)} = \boldsymbol{\eta}^{(t)} + (\mathbf{W}^{(t)})^{-1}(\mathbf{y} - \hat{\mathbf{p}}^{(t)})$ is the adjusted dependent variable.

6.4 Asymptotic Properties of MLE

Under regularity conditions, the MLE is:

$\sqrt{n}(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}_0) \xrightarrow{d} \mathcal{N}(\mathbf{0}, \mathcal{I}(\boldsymbol{\beta}_0)^{-1})$

Where the Fisher information matrix is:

$\mathcal{I}(\boldsymbol{\beta}) = -E\left[\frac{\partial^2 \ell_i}{\partial \boldsymbol{\beta}\partial\boldsymbol{\beta}^T}\right] = \sum_{i=1}^n w_i \mathbf{x}_i\mathbf{x}_i^T$

The variance-covariance matrix of $\hat{\boldsymbol{\beta}}$ is estimated by the inverse of the observed information matrix:

$\hat{\mathbf{V}}(\hat{\boldsymbol{\beta}}) = \left[-\mathbf{H}(\hat{\boldsymbol{\beta}})\right]^{-1} = \left(\mathbf{X}^T\hat{\mathbf{W}}\mathbf{X}\right)^{-1}$

6.5 Interpreting Logit Coefficients as Odds Ratios

For the logit model, exponentiating the coefficient gives the odds ratio:

$OR_k = e^{\hat{\beta}_k}$

Interpretation: A one-unit increase in $X_k$ multiplies the odds of $Y=1$ by $e^{\hat{\beta}_k}$:

$\frac{P(Y=1 \mid X_k + 1) / P(Y=0 \mid X_k + 1)}{P(Y=1 \mid X_k) / P(Y=0 \mid X_k)} = e^{\hat{\beta}_k}$

• If $OR_k > 1$: $X_k$ increases the odds of $Y=1$.
• If $OR_k < 1$: $X_k$ decreases the odds of $Y=1$.
• If $OR_k = 1$: No effect on odds.

⚠️ Odds ratios are not the same as probability ratios (relative risks). Do not interpret the odds ratio as "X% more likely." Convert to marginal probabilities via the APE for clearer communication.

6.6 Marginal Effects in Binary Models

For a continuous covariate $X_k$, the marginal effect of $X_k$ on $P(Y=1)$ for individual $i$:

Logit: $\frac{\partial P_i}{\partial X_{ik}} = \hat{p}_i(1-\hat{p}_i)\hat{\beta}_k = \Lambda(\mathbf{x}_i^T\hat{\boldsymbol{\beta}})[1-\Lambda(\mathbf{x}_i^T\hat{\boldsymbol{\beta}})]\hat{\beta}_k$

Probit: $\frac{\partial P_i}{\partial X_{ik}} = \phi(\mathbf{x}_i^T\hat{\boldsymbol{\beta}})\hat{\beta}_k$

For a discrete/binary covariate $X_k \in \{0,1\}$, the marginal effect is the discrete change in predicted probability:

$\Delta P_i = P(Y_i=1 \mid X_{ik}=1, \mathbf{x}_{i,-k}) - P(Y_i=1 \mid X_{ik}=0, \mathbf{x}_{i,-k})$

6.7 Standard Errors for Average Partial Effects (Delta Method)

The APE is a nonlinear function of $\hat{\boldsymbol{\beta}}$. Standard errors are obtained via the delta method:

$\widehat{SE}(APE_k) = \sqrt{\mathbf{g}_k^T \hat{\mathbf{V}}(\hat{\boldsymbol{\beta}}) \mathbf{g}_k}$

Where $\mathbf{g}_k = \partial \widehat{APE}_k / \partial \hat{\boldsymbol{\beta}}$ is the gradient of the APE with respect to the coefficient vector. Alternatively, use the bootstrap for more reliable inference with small samples.

7. Hypothesis Testing and Inference

7.1 The Wald Test

The Wald test for $H_0: \beta_k = 0$ uses the asymptotic normality of the MLE:

$z = \frac{\hat{\beta}_k}{SE(\hat{\beta}_k)} \sim \mathcal{N}(0,1)$

Or equivalently, $z^2 \sim \chi^2_1$. For a vector of $q$ restrictions $H_0: \mathbf{R}\boldsymbol{\beta} = \mathbf{r}$:

$W = (\mathbf{R}\hat{\boldsymbol{\beta}} - \mathbf{r})^T\left[\mathbf{R}\hat{\mathbf{V}}(\hat{\boldsymbol{\beta}})\mathbf{R}^T\right]^{-1}(\mathbf{R}\hat{\boldsymbol{\beta}} - \mathbf{r}) \sim \chi^2_q$

7.2 The Likelihood Ratio Test

The Likelihood Ratio (LR) test compares a restricted model (imposing $H_0$) to an unrestricted model:

$LR = -2\left[\ell(\hat{\boldsymbol{\beta}}_{restricted}) - \ell(\hat{\boldsymbol{\beta}}_{unrestricted})\right] \sim \chi^2_q$

Where $q$ is the number of restrictions. The LR test is generally preferred over the Wald test because it is invariant to reparameterisation and often has better finite-sample properties.

Special case: The LR test comparing a model with covariates to an intercept-only model:

$LR = -2[\ell_0 - \ell_{full}] \sim \chi^2_k$

Where $\ell_0 = n[p_0\ln p_0 + (1-p_0)\ln(1-p_0)]$ and $p_0 = \bar{Y}$ is the sample proportion.

7.3 The Score (Lagrange Multiplier) Test

The Score test (Rao test) only requires estimating the restricted model:

$S = \mathbf{s}(\hat{\boldsymbol{\beta}}_{restricted})^T \hat{\mathcal{I}}(\hat{\boldsymbol{\beta}}_{restricted})^{-1} \mathbf{s}(\hat{\boldsymbol{\beta}}_{restricted}) \sim \chi^2_q$

Useful when the unrestricted model is computationally expensive to estimate.

7.4 Test Equivalences and Recommendations

The three tests are asymptotically equivalent but differ in finite samples. The LR test is generally most reliable.

7.5 Confidence Intervals

A $(1-\alpha)\times100\%$ Wald confidence interval for $\beta_k$:

$\hat{\beta}_k \pm z_{\alpha/2} \times SE(\hat{\beta}_k)$

A profile likelihood confidence interval (more reliable for small samples):

$CI_{PL} = \left\{\beta_k : -2[\ell(\hat{\boldsymbol{\beta}}_{-k}, \beta_k) - \ell(\hat{\boldsymbol{\beta}})] \leq \chi^2_{\alpha,1}\right\}$

7.6 Testing IIA with the Hausman-McFadden Test

The Hausman-McFadden test for IIA in the MNL compares the full-sample MNL estimates to estimates obtained after removing one alternative from the choice set:

$H_{IIA} = (\hat{\boldsymbol{\beta}}_s - \hat{\boldsymbol{\beta}}_f)^T\left[\hat{\mathbf{V}}_s - \hat{\mathbf{V}}_f\right]^{-1}(\hat{\boldsymbol{\beta}}_s - \hat{\boldsymbol{\beta}}_f) \sim \chi^2_k$

Where $\hat{\boldsymbol{\beta}}_s$ are estimates from the restricted choice set and $\hat{\boldsymbol{\beta}}_f$ are estimates from the full choice set. Rejection suggests IIA violation.

⚠️ The Hausman-McFadden test has poor finite-sample properties and can produce negative test statistics. The Small-Hsiao test offers an alternative. Neither test is definitive. Subject-matter knowledge about alternative similarity remains essential.

7.7 Testing the Proportional Odds Assumption

The Brant test for the proportional odds assumption in ordered logit estimates a separate binary logit for each cumulative split and tests whether the coefficients are equal across splits:

$H_0: \boldsymbol{\beta}^{(1)} = \boldsymbol{\beta}^{(2)} = \dots = \boldsymbol{\beta}^{(J-1)}$

A chi-squared test statistic is formed from the sum of squared differences in estimates across cumulative splits, weighted by their precision. Rejection indicates the proportional odds assumption is violated, and a generalised ordered logit or multinomial logit should be considered.

7.8 Robust Standard Errors in Discrete Choice Models

While MLE standard errors are derived from the information matrix, misspecification-robust (sandwich) standard errors are available:

$\hat{\mathbf{V}}_{sandwich}(\hat{\boldsymbol{\beta}}) = \hat{\mathbf{H}}^{-1}\hat{\mathbf{B}}\hat{\mathbf{H}}^{-1}$

Where $\hat{\mathbf{H}} = -\sum_i \partial^2\ell_i/\partial\boldsymbol{\beta}\partial\boldsymbol{\beta}^T$ and $\hat{\mathbf{B}} = \sum_i (\partial\ell_i/\partial\boldsymbol{\beta})(\partial\ell_i/\partial\boldsymbol{\beta})^T$ (the outer product of scores).

• Use heteroscedasticity-robust (sandwich) SEs when the distributional assumption may be misspecified.
• Use cluster-robust SEs when observations are grouped (e.g., individuals within firms or households within regions).

8. Effect Size Measures

8.1 Average Partial Effects (APE / AME)

The primary effect size in discrete choice models is the Average Partial Effect (APE), also called the Average Marginal Effect (AME):

$APE_k = \frac{1}{n}\sum_{i=1}^n \frac{\partial P(Y_i = j \mid \mathbf{x}_i)}{\partial X_{ik}}$

Interpretation: The average change in the probability of outcome $j$ associated with a one-unit increase in $X_k$, averaging over all individuals in the sample.

For a binary covariate $D_i \in \{0,1\}$:

$APE_k = \frac{1}{n}\sum_{i=1}^n \left[P(Y_i=1 \mid D_i=1, \mathbf{x}_{i,-k}) - P(Y_i=1 \mid D_i=0, \mathbf{x}_{i,-k})\right]$

8.2 Odds Ratios and Relative Risk

8.3 Predicted Probability Changes

For practical communication, report predicted probabilities at meaningful covariate values:

$\Delta\hat{P} = \hat{P}(Y_i=1 \mid \mathbf{x}_{high}) - \hat{P}(Y_i=1 \mid \mathbf{x}_{low})$

Where $\mathbf{x}_{high}$ and $\mathbf{x}_{low}$ represent two substantively meaningful covariate profiles (e.g., high-income vs. low-income; treated vs. untreated).

8.4 Standardised Coefficients in Discrete Choice Models

To compare the relative importance of different covariates, standardise the APE by the standard deviation of the outcome:

$\beta^*_k = APE_k \times \frac{s_{X_k}}{s_Y}$

Where $s_{X_k}$ is the standard deviation of $X_k$ and $s_Y = \sqrt{\bar{p}(1-\bar{p})}$ for binary outcomes. This produces an effect size interpretable as the change in probability (in units of the outcome SD) per SD change in $X_k$.

8.5 McFadden's Pseudo-$R^2$ as Effect Size

McFadden's pseudo-$R^2$ measures the proportional improvement in log-likelihood:

$\rho^2 = 1 - \frac{\ell(\hat{\boldsymbol{\beta}})}{\ell_0}$

Where $\ell_0 = \ell(\beta_0 \text{ only})$ is the log-likelihood of the intercept-only model. While not a pure effect size, it provides a scale for comparing model fit improvement:

8.6 Willingness to Pay (WTP) as Effect Size in Choice Experiments

In stated or revealed preference studies, WTP contextualises effect sizes economically:

$WTP_k = -\frac{APE_k}{APE_{cost}} = -\frac{\hat{\gamma}_k}{\hat{\gamma}_{cost}}$

Report WTP with confidence intervals obtained via the delta method or Krinsky-Robb simulation.

9. Model Fit and Evaluation

9.1 Goodness-of-Fit Statistics

9.2 Pseudo-R² Measures

Multiple pseudo-$R^2$ measures exist; they capture different aspects of fit:

McFadden (1974): $\rho^2_{McFadden} = 1 - \frac{\ell(\hat{\boldsymbol{\beta}})}{\ell_0}$

Cox-Snell: $R^2_{CS} = 1 - \left(\frac{\mathcal{L}_0}{\mathcal{L}(\hat{\boldsymbol{\beta}})}\right)^{2/n}$

Nagelkerke (normalised Cox-Snell): $R^2_{N} = \frac{R^2_{CS}}{1 - \mathcal{L}_0^{2/n}}$

⚠️ No single pseudo-$R^2$ is universally "correct." Report multiple, and always prefer out-of-sample predictive performance metrics (AUC, Brier score) for evaluating predictive models.

9.3 Classification Metrics for Binary Models

For binary models, at a threshold $c$ (default $c = 0.5$):

$\hat{Y}_i = \mathbf{1}[\hat{p}_i \geq c]$

The Receiver Operating Characteristic (ROC) curve plots sensitivity vs. $(1-\text{specificity})$ across all classification thresholds $c \in [0,1]$. The Area Under the Curve (AUC) summarises discrimination:

9.4 Calibration

Calibration assesses whether predicted probabilities match observed outcome rates.

Hosmer-Lemeshow test: Partition observations into $G$ (typically 10) quantile groups by predicted probability. Compare observed and expected counts in each group:

$H = \sum_{g=1}^G \frac{(O_g^1 - E_g^1)^2}{E_g^1} + \frac{(O_g^0 - E_g^0)^2}{E_g^0} \sim \chi^2_{G-2}$

Where $O_g^1$ and $E_g^1$ are observed and expected counts of $Y=1$ in group $g$. Rejection suggests poor calibration.

Calibration plot: Plot mean predicted probability vs. observed proportion in each decile group. A well-calibrated model lies along the 45° diagonal.

9.5 Information Criteria for Model Comparison

When comparing non-nested models (e.g., logit vs. probit; different covariate sets):

$AIC = -2\ell(\hat{\boldsymbol{\beta}}) + 2k$ $BIC = -2\ell(\hat{\boldsymbol{\beta}}) + k\ln(n)$

Lower values indicate better fit. BIC imposes a heavier penalty on model complexity, favouring parsimony. AIC and BIC are only directly comparable for models fitted to the same dataset with the same outcome variable.

9.6 Out-of-Sample Validation

For predictive models, always assess performance on held-out data:

• $k$-fold cross-validation: Partition data into $k$ folds; train on $k-1$ and test on 1; rotate and average performance metrics.
• Train-test split: Randomly assign 70-80% to training and 20-30% to test.
• Temporal split: For time-indexed data, train on earlier periods and test on later periods.
• Brier score: Mean squared error for probability predictions: $BS = n^{-1}\sum_i (\hat{p}_i - Y_i)^2$.

10. Diagnostics and Assumption Testing

10.1 Residuals in Discrete Choice Models

Unlike OLS, residuals in discrete choice models require careful definition.

Pearson residuals: $r_i^P = \frac{Y_i - \hat{p}_i}{\sqrt{\hat{p}_i(1-\hat{p}_i)}}$

Deviance residuals: $r_i^D = \text{sign}(Y_i - \hat{p}_i)\sqrt{-2\left[Y_i\ln\hat{p}_i + (1-Y_i)\ln(1-\hat{p}_i)\right]}$

The deviance (sum of squared deviance residuals) is:

$D = \sum_{i=1}^n (r_i^D)^2 = -2\ell(\hat{\boldsymbol{\beta}})$

Standardised residuals for outlier detection: $r_i^{std} = r_i^D / \sqrt{1 - h_{ii}}$, where $h_{ii}$ is the hat-value (leverage).

10.2 Influence and Leverage

Leverage in logit/probit: $h_{ii} = \hat{w}_i \mathbf{x}_i^T(\mathbf{X}^T\hat{\mathbf{W}}\mathbf{X})^{-1}\mathbf{x}_i$

Cook's distance analogue: $CD_i = \frac{(r_i^P)^2 h_{ii}}{k(1-h_{ii})^2}$

DFFITS and DFBETAS analogues are available for identifying influential observations. Flag observations with $|CD_i| > 4/n$ or $|r_i^{std}| > 2$ for inspection.

10.3 Testing the Proportional Odds Assumption (Ordered Logit)

Brant (1990) test: Estimates a binary logit for each of the $J-1$ cumulative dichotomisations and tests whether coefficients are equal. Available both as a global test (all covariates) and variable-specific tests:

$\chi^2_{Brant} = \sum_{j=1}^{J-2} (\hat{\boldsymbol{\beta}}_j - \hat{\boldsymbol{\beta}}_{J-1})^T \left[\hat{\mathbf{V}}_j + \hat{\mathbf{V}}_{J-1}\right]^{-1} (\hat{\boldsymbol{\beta}}_j - \hat{\boldsymbol{\beta}}_{J-1})$

Graphical check: Plot ordered logit coefficients estimated separately for each binary cumulative split. Coefficients that vary substantially suggest a violation.

Remedy if violated:

• Generalised Ordered Logit (partial proportional odds): Allow some (but not all) coefficients to vary across thresholds.
• Multinomial Logit: Drop the ordinal structure entirely; less efficient but unrestricted.
• Stereotype Logit (reduced-rank MNL): Intermediate model allowing partial ordering.

10.4 Testing IIA

Multiple tests for IIA are available, each with limitations:

Remedy if IIA fails:

• Nested Logit: Group correlated alternatives into nests.
• Mixed Logit: Allow correlation across alternatives through random coefficients.
• Multinomial Probit: Directly models correlated errors; computationally intensive.

10.5 Checking for Complete Separation

Complete separation occurs when a covariate or linear combination of covariates perfectly predicts the outcome — the MLE does not exist (the log-likelihood has no finite maximum):

• Perfect separation: $P(\hat{Y}_i = Y_i) = 1$ for some linear predictor.
• Quasi-complete separation: The outcome is perfectly predicted for a subset of observations.

Detection: MLE algorithm fails to converge; extremely large coefficient estimates with very large standard errors; implausible predicted probabilities near 0 or 1.

Remedies:

• Firth penalised MLE: Modifies the score equations by a Jeffreys prior penalty — reduces bias in small samples and resolves separation.
• Ridge-penalised logit: $\ell_{penalised} = \ell(\boldsymbol{\beta}) - \lambda\|\boldsymbol{\beta}\|^2$.
• Bayesian logit/probit: Place weakly informative priors on coefficients.
• Drop/combine categories: Merge sparse categories causing separation.

10.6 Heteroscedasticity in Probit (Heteroscedastic Probit)

In the standard probit, $\text{Var}(\epsilon_i) = 1$ for all $i$. If the true error variance is heteroscedastic:

$\text{Var}(\epsilon_i) = \sigma_i^2 = [e^{\mathbf{h}_i^T\boldsymbol{\delta}}]^2$

The heteroscedastic probit models:

$P(Y_i = 1 \mid \mathbf{x}_i, \mathbf{h}_i) = \Phi\left(\frac{\mathbf{x}_i^T\boldsymbol{\beta}}{e^{\mathbf{h}_i^T\boldsymbol{\delta}}}\right)$

Standard probit estimates are inconsistent under heteroscedasticity (unlike OLS which remains consistent, only losing efficiency). The linktest (adding the squared predicted index as a covariate) checks for systematic misspecification.

10.7 Goodness-of-Link Tests

The linktest (Pregibon, 1980) adds $\hat{\eta}_i^2 = (\mathbf{x}_i^T\hat{\boldsymbol{\beta}})^2$ as an additional regressor to the fitted model:

$P(Y_i=1) = F(\hat{\eta}_i \cdot \delta_1 + \hat{\eta}_i^2 \cdot \delta_2)$

Under correct specification, $\hat{\delta}_2 = 0$ (the squared term should not be significant). A significant $\hat{\delta}_2$ indicates link function misspecification or omitted non-linear terms.

11. Extensions: Multinomial and Conditional Logit

11.1 MNL Log-Likelihood

For outcome $Y_i \in \{1, \dots, J\}$ with reference category $j=1$, the MNL log-likelihood:

$\ell(\{\boldsymbol{\beta}_j\}_{j=2}^J) = \sum_{i=1}^n \sum_{j=1}^J \mathbf{1}[Y_i=j] \ln P(Y_i=j \mid \mathbf{x}_i)$

$= \sum_{i=1}^n \left[\sum_{j=2}^J \mathbf{1}[Y_i=j]\mathbf{x}_i^T\boldsymbol{\beta}_j - \ln\left(1 + \sum_{k=2}^J e^{\mathbf{x}_i^T\boldsymbol{\beta}_k}\right)\right]$

11.2 Marginal Effects in MNL

For the MNL, the marginal effect of $X_k$ on $P(Y_i = j)$:

$\frac{\partial P(Y_i = j)}{\partial X_k} = P(Y_i=j)\left[\beta_{jk} - \sum_{l=1}^J P(Y_i=l)\beta_{lk}\right]$

Note: Cross-effects — the effect of a covariate on a different category's probability — may be positive or negative, depending on model parameters.

Average Partial Effect: $APE_{jk} = \frac{1}{n}\sum_{i=1}^n \hat{P}_{ij}\left[\hat{\beta}_{jk} - \sum_{l=1}^J \hat{P}_{il}\hat{\beta}_{lk}\right]$

11.3 The Conditional Logit and Mixed-Effects Specification

The full Mixed Logit specification that includes both individual-varying and alternative-varying attributes:

$V_{ij} = \underbrace{\mathbf{z}_{ij}^T\boldsymbol{\gamma}}_{\text{alternative attributes}} + \underbrace{\mathbf{x}_i^T\boldsymbol{\beta}_j}_{\text{individual chars. × alt. FE}}$

Where:

• $\mathbf{z}_{ij}$: Attributes that vary across both individuals and alternatives (e.g., travel time from $i$'s origin to alternative $j$).
• $\mathbf{x}_i$: Characteristics of the individual (e.g., income), interacted with alternative-specific dummy variables to allow the effect to vary across alternatives.