Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
16,561 characters · 6 sections · 13 citation commands
The standard Gibbs sampler of Mixed Multinomial Logit (MMNL) models involves sampling from conditional densities of utility parameters using Metropolis-Hastings (MH) algorithm due to unavailability of conjugate prior for logit kernel. To address this non-conjugacy concern, we propose the application of P{\'o}lygamma data augmentation (PG-DA) technique for the MMNL estimation. The posterior estimates of the augmented and the default Gibbs sampler are similar for two-alternative scenario (binary choice), but we encounter empirical identification issues in the case of more alternatives ($J \geq 3$).
The mixed multinomial logit (MMNL) model mcfadden2000mixed is established as follows: We consider a standard discrete choice setup, in which on choice occasion $t \in \{1, \ldots T \}$, a decision-maker $n \in \{1, \ldots N \}$ derives utility $U_{ntj} = V(\boldsymbol{X}_{ntj}, \boldsymbol{\Gamma}_{n}) + \epsilon_{ntj}$ from alternative $j \in \{1, \ldots J \}$. Here, $V()$ denotes the representative utility, $\boldsymbol{X}_{ntj}$ is a row-vector of covariates, $\boldsymbol{\Gamma}_{n}$ is a collection of taste parameters, and $\epsilon_{ntj}$ is a stochastic disturbance. The assumption $\epsilon_{ntj} \sim \text{Gumbel}(0,1)$ leads to a multinomial logit (MNL) kernel such that the probability that decision-maker $n$ chooses alternative $j$ on choice occasion $t$ is
where $y_{nt}$ captures the observed choice. The choice probability can be iterated over choice scenarios to obtain the probability of observing a decision-maker's sequence of choices $\boldsymbol{y}_{n}$:
We consider a general utility specification under which tastes $\boldsymbol{\Gamma}_{n}$ are partitioned into fixed taste parameters $\boldsymbol{\alpha}$, which are invariant across decision-makers, and random taste parameters $\boldsymbol{\beta}_{n}$, which are individual-specific, such that $\boldsymbol{\Gamma}_{n} =
^{\top}$, whereby $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}_{n}$ are vectors of lengths $L$ and $K$, respectively. Analogously, the row-vector of covariates $\boldsymbol{X}_{ntj}$ is partitioned into attributes $\boldsymbol{X}_{ntj,F}$, which pertain to the fixed parameters $\boldsymbol{\alpha}$, as well as into attributes $\boldsymbol{X}_{ntj,R}$, which pertain to the individual-specific parameters $\boldsymbol{\beta}_{n}$, such that $\boldsymbol{X}_{ntj} =
$. For simplicity, we assume that the representative utility is linear-in-parameters, i.e.
The distribution of tastes $\boldsymbol{\beta}_{1:N}$ is assumed to be multivariate normal, i.e. $\boldsymbol{\beta}_{n} \sim \text{N}(\boldsymbol{\zeta}, \boldsymbol{\Omega})$ for $n = 1, \dots, N$, where $\boldsymbol{\zeta}$ is a mean vector and $\boldsymbol{\Omega}$ is a covariance matrix. In a fully Bayesian setup, the invariant (across individuals) parameters $\boldsymbol{\alpha}$, $\boldsymbol{\zeta}$, $\boldsymbol{\Omega}$ are also considered to be random parameters and are thus given priors. We use normal priors for the fixed parameters $\boldsymbol{\alpha}$ and for the mean vector $\boldsymbol{\zeta}$. Following tan2017stochastic and akinc2018bayesian, we employ Huang's half-t prior huang2013Simple for covariance matrix $\boldsymbol{\Omega}$, as this prior specification exhibits superior noninformativity properties compared to other prior specifications for covariance matrices. In particular, akinc2018bayesian show that Huang's half-t prior outperforms the inverse Wishart prior, which is often employed in fully Bayesian specifications of MMNL models train2009discrete, in terms of parameter recovery.
Stated succinctly, the generative process of the fully Bayesian MMNL model is:
where ((ref)) and ((ref)) induce Huang's half-t prior huang2013Simple. $\{ \boldsymbol{\lambda}_{0}, \boldsymbol{\Xi}_{0}, \boldsymbol{\mu}_{0}, \boldsymbol{\Sigma}_{0}, \nu, A_{1:K} \}$ are known hyper-parameters, and $\boldsymbol{\theta} = \{ \boldsymbol{\alpha}, \boldsymbol{\zeta}, \boldsymbol{\Omega}, \boldsymbol{a}, \boldsymbol{\beta}_{1:N}\}$ is a collection of model parameters whose posterior distribution we wish to estimate.
The default Gibbs sampler for posterior inference in MMNL models involves Metropolis steps to take draws from conditional densities of the utility parameters ($\bm{\beta}_{n}$ and $\bm{\alpha}$) due to the unavailability of a conjugate prior for the MNL kernel. MCMC estimation of binary and multinomial logistic regression models encounters a similar issue of non-conjugacy. P{\'o}lya-Gamma data augmentation (PG-DA) is the state-of-the-art technique to handle non-conjugacy in MCMC estimation of binary logistic regression models polson2013bayesian. PG-DA augments the Gibbs sampler by introducing an additional P{\'o}lya-Gamma distributed latent variable, which circumvents the need of the Metropolis algorithm by ensuring conjugate updates. polson2013bayesian also extend PG-DA to the multinomial logistic regression model. Yet, this extension requires all utility (or link function) parameters to be alternative-specific. We use the same idea in deriving a PG-DA-based Gibbs sampler for MMNL, but we have to consider the same restriction on utility specification, i.e. replace $\bm{\Gamma}_{n}$ by $\bm{\Gamma}_{nj}$.
The representative utility is: $V_{ntj} = \bm{X}_{ntj} \bm{\Gamma}_{nj} =\bm{X}_{ntj,F} \bm{\alpha}_{j} + \bm{X}_{ntj,R} \bm{\beta}_{nj}$, where $\bm{\beta}_{nj} \sim \text{N}(\bm{\zeta}_j,\bm{\Omega})$. The hyper-parameters remain the same, but the model parameters are $\bm{\theta} = \{ \bm{\alpha}_{1:J}, \bm{\zeta}_{1:J}, \bm{\Omega}, a_{1:K}, \bm{\beta}_{\{1:N,1:J\}} \}$. Adhering to the original notation, we can write the joint distribution of the data and the model parameters:
Algorithm (ref) presents the augmented Gibbs sampler for the MMNL model. The conditional densities of $a_{1:K}$, $\bm{\Omega}$, and $\bm{\zeta}_{1:J}$ are similar to those of the Allenby-Train procedure akinc2018bayesian. The next subsection details the derivation of conditional densities of $\bm{\beta}_{\{1:N,1:J\}}$ and $\bm{\alpha}_{1:J}$.
Using holmes2006bayesian, we can convert the multinomial logit likelihood expression to the binary logit likelihood:
where $\bm{\theta}_{-\bm{\beta}_{nj}}$ is a resulting parameter vector after removing $\bm{\beta}_{nj}$ and
We now introduce a P{\'o}lya--Gamma distributed auxiliary variable $\phi_{ntk} \sim \text{PG}(1,0) \; \forall n,t,k$ and $\kappa_{ntk} = y_{ntk}- \frac{1}{2}$. Now consider the identity derived by polson2013bayesian:
The conditional density of $\bm{\beta}_{nj}$ is:
The conditional distribution of $\bm{\beta}_{nj}$ is:
The conditional density of $\bm{\alpha}_{j}$ can be derived similarly:
We test the performance of the PG-DA-based Gibbs sampler against the Metropolis-based Gibbs sampler in a Monte Carlo study. We first consider the MNL model ($\bm{\Gamma}_{nj} = \bm{\alpha}_{j}$), where both samplers perform equally well. For MMNL model, the posterior estimates of the proposed PG-DA approach and the default Gibbs sampler are similar for $J=2$, but we encounter an explosion of the conditional distribution parameters in the case of more alternatives ($J \geq 3$). Results for $J=2$ and MATLAB code is available upon request.
This appears to be an issue of empirical identification because of too many model parameters. Before the PG-DA sampler diverges, representative utilities are either very small or very large in magnitude for all alternatives across all observations. Therefore, instead of the actual magnitude of utilities, their comparative scales determine the choice probabilities. Thus, the algorithm might have a tendency to increase the relative scale of the latent utilities by increasing the scale of the parameters.
In fact, prior to divergence the probability estimates of the chosen and non-chosen alternatives are close to one and zero, respectively. We speculate that such behavior might be a consequence of too many model parameters, which might allow the algorithm to find a parameter configuration that can fit the data very well (in terms of choice probabilities). Once the algorithm finds that configuration, it starts to increase the relative scale between the utilities (thus allowing the chosen alternatives to have probability close to one), causing the parameter explosion.
As future research, stick-breaking constructions can be explored to adopt PG-DA in MCMC estimation of MMNL while keeping a parsimonious model specification, i.e. with generic utility parameters linderman2015dependent,zhang2017permuted. However, before adopting these constructions, consistency with microeconomic theory needs to be established first.