Bayesian Estimation of Mixed Multinomial Logit Models: Advances and Simulation-Based Evaluations

\thispagestyle{empty}

\thispagestyle{empty}

Abstract

Variational Bayes (VB) methods have emerged as a fast and computationally-efficient alternative to Markov chain Monte Carlo (MCMC) methods for scalable Bayesian estimation of mixed multinomial logit (MMNL) models. It has been established that VB is substantially faster than MCMC at practically no compromises in predictive accuracy. In this paper, we address two critical gaps concerning the usage and understanding of VB for MMNL. First, extant VB methods are limited to utility specifications involving only individual-specific taste parameters. Second, the finite-sample properties of VB estimators and the relative performance of VB, MCMC and maximum simulated likelihood estimation (MSLE) are not known. To address the former, this study extends several VB methods for MMNL to admit utility specifications including both fixed and random utility parameters. To address the latter, we conduct an extensive simulation-based evaluation to benchmark the extended VB methods against MCMC and MSLE in terms of estimation times, parameter recovery and predictive accuracy. The results suggest that all VB variants with the exception of the ones relying on an alternative variational lower bound constructed with the help of the modified Jensen’s inequality perform as well as MCMC and MSLE at prediction and parameter recovery. In particular, VB with nonconjugate variational message passing and the delta-method (VB-NCVMP-$\Delta$) is up to 16 times faster than MCMC and MSLE. Thus, VB-NCVMP-$\Delta$ can be an attractive alternative to MCMC and MSLE for fast, scalable and accurate estimation of MMNL models. \\ \\ Keywords: Variational Bayes; Bayesian inference; mixed logit; nonconjugate variational message passing.

\pagenumbering{arabic}

Introduction

The mixed multinomial logit (MMNL) model mcfadden2000mixed is the workhorse model in many disciplines---such as economics, health, marketing and transportation---that are concerned with the analysis and prediction of individual choice behavior. While maximum simulated likelihood estimation train2009discrete is the predominant estimation strategy for MMNL models, the Bayesian approach represents an alternative estimation strategy, which entails the key benefit that the whole posterior distribution of all model parameters including the individual-specific parameters can be obtained. Posterior inference in MMNL models is typically performed with the help of Markov chain Monte Carlo (MCMC) methods, which approximate the posterior distribution of the MMNL model parameters through samples from a Markov chain whose stationary distribution is the posterior distribution of interest rossi2012bayesian, train2009discrete. While MCMC methods constitute a powerful framework for posterior inference in complex probabilistic models gelman2013bayesian, these methods are subject to several bottlenecks, which inhibit their scalability to large datasets, namely i) long computation times, ii) high costs for the storage of the posterior draws and iii) difficulties in assessing convergence blei2017variational, braun2010variational, depraetere2017comparison, tan2017stochastic.

Variational Bayes (VB) methods blei2017variational, jordan1999introduction, ormerod2010explaining have emerged as an alternative to MCMC and promise to address the shortcomings of MCMC methods. The basic intuition behind VB is to view approximate Bayesian inference as an optimization problem rather than a sampling problem. VB aims at finding a parametric variational distribution over the unknown model parameters, whereby the parameters of the variational distribution are optimized such that the probability distance (typically measured in terms of the Kullback-Leibler divergence) between the exact posterior distribution and the variational distribution is minimal. A key challenge in the application of VB to posterior inference in MMNL models is that the expectation of the logarithm of the choice probabilities---or, to be precise, the expectation of the log-sum of exponentials (E-LSE) term---cannot be expressed in closed form, because there is no general conjugate prior for the multinomial logit model. As a consequence, updates for variational factors pertaining to utility parameters require special treatment. The literature proposes different methods to facilitate VB for posterior inference in MMNL models braun2010variational, depraetere2017comparison, tan2017stochastic. In essence, these approaches proceed as follows: The E-LSE term is approximated either analytically or by simulation, or an alternative variational lower bound is defined. Then, updates for the nonconjugate variational factors are performed with the help of either quasi-Newton (QN) methods nocedal2006numerical or the nonconjugate variational message passing (NCVMP) approach knowles2011non.

Extant studies of VB methods for posterior inference in MMNL models establish that VB is substantially faster than MCMC at negligible compromises in predictive accuracy braun2010variational, depraetere2017comparison, tan2017stochastic. However, these studies find wanting in several important ways. First, the QN and NCVMP updating strategies have been studied in isolation from each other and their relative performance is not known. Second, none of these studies compare VB to the widely-used MSLE method. Third, the performance of the considered estimation approaches has only been evaluated in terms of predictive accuracy, while the finite sample properties, i.e. the ability to recover true parameters, of the estimators are not known. Fourth, VB methods have only been implemented and tested for posterior inference in MMNL models with only individual-specific utility parameters despite the practical relevance of fixed utility parameters in discrete choice modeling applications.

Consequently, the objective of this paper is twofold: First, we extend several VB methods to allow for posterior inference in MMNL models with a more general utility specification including both fixed and random utility parameters.\footnote{Strictly, all model parameters are random quantities in Bayesian estimation. Here, we adopt the nomenclature used by train2009discrete and refer to utility parameters that are invariant across decision-makers as fixed utility parameters and to utility parameters that are individual-specific and (normally) distributed across decision-makers as random utility parameters.} Then, we carry out a comprehensive simulation-based evaluation, in which we contrast the performance of different VB methods, MCMC and MSLE in terms of estimation times, parameter recovery and predictive accuracy.\footnote{In this paper, we compare VB and MCMC with MSLE, as MSLE continues to represent the most widely used estimation strategy for MMNL models. We acknowledge that Bhat and co-authors have developed the Maximum Approximate Composite Marginal Likelihood (MACML) approach bhat2011simulation for frequentist estimation of mixed multinomial probit (MMNP) models. MACML has been shown to be faster and more accurate than MSLE patil2017simulation. In addition, the approach is flexible, as it has been used for the estimation of integrated choice and latent variable models bhat2014new and MMNP models with non-normal parametric mixing distributions bhat2018new. Despite its limitation to MMNP, MACML thus represents an attractive alternative to MSLE for frequentist estimation of mixed random utility models. However, we concur with bhat2018new that MMNP is no more or less general than MMNL. Comparisons between Bayesian estimation methods for MMNL and MACML for MMNP are admittedly intriguing but are beyond the scope of the current paper.}

We emphasize that the inclusion of fixed utility parameters, in addition to individual-specific utility parameters, is important in practice bansal2018extending: First, alternative-specific fixed effects can be introduced by including alternative-specific constants (ASCs) in the utility specification. Assuming ASCs to be random may result in empirical identification issues, especially if their distribution is similar to that of the error term train2009discrete. Second, utility parameters corresponding to individual-specific characteristics (e.g. age, gender etc.) are typically assumed to be fixed. Treating these alternative-specific parameters as random may not provide substantive behavioral insights and may unnecessarily inflate the number of random parameters so that the “curse of dimensionality” becomes a concern cherchi2012monte. Third, systematic taste variation can be parsimoniously represented through the inclusion of additional fixed parameters that pertain to interactions of the alternative-specific attribute (e.g. cost or travel time) and relevant individual-specific attributes bhat1998accommodating.

In the case of MSLE, one can easily accommodate fixed utility parameters by specifying them as random utility parameters with a constrained variance, because the individual-specific parameters are integrated out so that that the fixed parameters can be jointly updated with the parameters of the mixing distribution. This approach is not feasible for Bayesian estimation methods, because the individual-specific parameters are directly estimated train2009discrete, rossi2012bayesian. If the fixed utility parameters were specified as random with a constrained variance in VB estimation, the respective variational factors would have to be identical across decision-makers. However, it is not straightforward to impose this restriction in the existing VB methods. This is because the variational factors of the individual-specific parameters are updated independently for each individual, while updates for the variational factors of the fixed parameters necessarily depend on all observations.

We organize the remainder of this paper as follows: First, we provide a fully Bayesian formulation of the MMNL model (Section (ref)). To be self-contained, we present the default MCMC method for posterior inference in MMNL models (Section (ref)). Then, we describe different VB methods for posterior inference in MMNL models with a more general utility specification including a combination of both fixed and individual-specific utility parameters (Section (ref)). Next, we present the simulation-based evaluation (Section (ref)) and finally, we conclude (Section (ref)).

Mixed multinomial logit model

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_{n} \}$, 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 the set $C_{nt}$. 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 \in C_{nt}$ on choice occasion $t$ is

where $y_{nt} \in C_{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}$:

In this paper, 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 huang2013Simple,akinc2018Bayesian. 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 generative process implies the following joint distribution of data and model parameters:

where $\omega = \nu + K - 1$, $\boldsymbol{B} = 2\nu \text{diag}(\boldsymbol{a})$, $s = \frac{1}{2}$ and $r_{k} = A_{k}^{-2}$.\footnote{To be clear, the following forms of the Gamma and inverse Wishart distributions are considered:

whereby $\boldsymbol{\Omega}$ and $\boldsymbol{B}$ are $K \times K$ positive-definite matrices.} By Bayes' rule, the posterior distribution of interest is then given by

Exact inference of this posterior distribution is not possible, because the model evidence $\int P (\boldsymbol{y}_{1:N}, \boldsymbol{\theta}) d \boldsymbol{\theta}$ is not tractable. In the following sections, we discuss different strategies to approximate the posterior distribution of the MMNL model parameters and provide our extensions to some of these strategies under the more general linear-in-parameters utility specification including both fixed and and random taste parameters.

Markov chain Monte Carlo

The general idea of Markov chain Monte Carlo (MCMC) methods is to approximate a difficult-to-compute posterior distribution through samples from a Markov chain whose stationary distribution is the posterior distribution of interest robert2004monte.

In the present application, a Markov chain for the posterior distribution of the MMNL model parameters $\boldsymbol{\theta}$ can be constructed by taking samples from the conditional distributions of $\boldsymbol{\theta}$. Direct sampling from the conditional distributions of $\boldsymbol{\zeta}$, $\boldsymbol{\Omega}$ and $\boldsymbol{a}$ is possible, because the conditional distributions belong to known families of distributions. However, updates for $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}_{1:N}$ need to be generated with the help of random-walk (RW) Metropolis algorithms, because the nonconjugacy of the multinomial logit kernel and the normal priors leads to unrecognizable conditional distributions. The resulting MCMC algorithm is a blocked Gibbs sampler with two embedded Metropolis steps. A pseudo-code representation of the sampler is shown in Algorithm (ref). Here, $\rho_{\boldsymbol{\alpha}}$ and $\rho_{\boldsymbol{\beta}}$ denote step sizes, which need to be tuned.\footnote{In the subsequent applications of the sampling scheme, we apply the same tuning mechanism as train2009discrete, i.e. we let $\rho_{\boldsymbol{\alpha}} = 0.01$ and set $\rho_{\boldsymbol{\beta}}$ to an initial value of 0.1. After each iteration, $\rho_{\boldsymbol{\beta}}$ is decreased by 0.001, if the average acceptance rate across all decision-makers is less than 0.3; $\rho_{\boldsymbol{\beta}}$ is increased by 0.001, if the average acceptance rate across all decision-makers is more than 0.3.} The sampling scheme outlined in Algorithm (ref) is identical to the one studied by akinc2018Bayesian with the only difference that updates for the fixed parameters $\boldsymbol{\alpha}$ are incorporated. It is also known as the Allenby-Train procedure rossi2012bayesian, train2009discrete.

A bottleneck of Algorithm (ref) is its reliance on two RW Metropolis steps for the fixed and the individual-specific parameters, respectively. Notwithstanding that these steps are easy to implement and to vectorize, the RW Metropolis algorithm can be inefficient when it is tuned suboptimally rossi2012bayesian. If the step size is too small, the chain moves too quickly and the draws exhibit high serial correlation. If the step size is too large, the posterior is not properly explored and the algorithm can get stuck. The RW Metropolis algorithm can be replaced by an independence Metropolis algorithm rossi2012bayesian, which takes draws around the posterior mode. However, a complication of this approach is that at each iteration, a maximization needs to be performed to find the posterior mode, which is particularly challenging for the individual-specific parameters.

An emerging method to generate samples from a Markov chain is Hamiltonian Monte Carlo neal2011mcmc. HMC uses information contained in the gradient of the log target density to efficiently explore the posterior distribution of interest and to reduce the amount of serial correlation in the chains. A variant of HMC is the No-U-Turn sampler hoffman2014no-u-turn, which automatically adapts the number of leapfrog steps required for the discretization of the Hamiltonian dynamics underlying HMC. NUTS is interfaced by Stan carpenter2017stan, a probabilistic programming language that enables posterior inference on a wide variety of user-defined models. However, the generality of Stan comes at an immense computational cost, which is further aggravated when the model of interest depends on many parameters as is the case for MMNL.\footnote{ We also explored the use of Stan as part of the current research study but found that estimation times were prohibitive for the sample sizes considered in the simulation evaluation presented in Section (ref). Our experiences with Stan are generally consistent with the literature. ben2019foundations contrast NUTS with the Allenby-Train procedure and find that both methods perform equally well at recovering the true parameter values. However, whereas the reported estimation time for the Allenby-Train procedure is 12 minutes, NUTS had to be run “overnight”. vij2017random attempted to use Stan to estimate a MMNL model on a large dataset containing 30,166 observations from 17,700 individuals but were unable to do so due to memory constraints. A possible avenue for future research is to custom-code a NUTS procedure with analytical gradients to enable fast and scalable posterior inference for MMNL.}

Variational Bayes

Background

Variational Bayes blei2017variational, jordan1999introduction, ormerod2010explaining differs from MCMC in that approximate Bayesian inference is viewed as optimization problem rather than a sampling problem. Figure (ref) illustrates the conceptual differences between MCMC and VB. In MCMC, the posterior distribution of interest $P(\theta \vert \boldsymbol{y})$ is approximated through samples from a Markov chain whose stationary distribution is the posterior distribution of interest. In VB, the posterior distribution of interest is approximated through a parametric variational distribution $q(\theta \vert \boldsymbol{\nu})$ whose parameters $\boldsymbol{\nu}$ are fit such that the $P(\theta \vert \boldsymbol{y})$ and the approximating variational distribution are close in probability distance.

Casting approximate Bayesian inference as an optimization problem comes with several benefits which enable scaling Bayesian estimation to large datasets. First, the memory issues of MCMC are overcome, as only the variational parameters rather than the posterior draws need to be stored. Second, convergence can be straightforwardly assessed by evaluating the change in the variational lower bound (an alternative measure for the distance between the posterior distribution of interest and the approximating variational distribution) or the change in the estimates of the variational parameters from one iteration to another. Third, serial correlation is no longer a concern, as no samples are taken.

To build further intuition about the fundamental principles of VB, we consider a generative model $P(\boldsymbol{y}, \boldsymbol{\theta})$ consisting of observed data $\boldsymbol{y}$ and unknown parameters $\boldsymbol{\theta}$. Our goal is to find an approximation of the posterior distribution $P(\boldsymbol{\theta} \vert \boldsymbol{y})$. VB aims at finding a variational distribution $q(\boldsymbol{\theta})$ over the unknown parameters that is close to the actual posterior distribution $P(\boldsymbol{\theta} \vert \boldsymbol{y})$. A computationally-convenient way to measure the distance between two probability distributions is the Kullback-Leibler (KL) divergence kullback1951information. The KL divergence between $q(\boldsymbol{\theta})$ and $P(\boldsymbol{\theta} \vert \boldsymbol{y})$ is given by

The goal of VB is to minimize this divergence, i.e.

However, the expectation $\mathbb{E}_{q} \left \{ \ln P(\boldsymbol{\theta} \vert \boldsymbol{y}) \right \} = \mathbb{E}_{q} \left \{ \ln P(\boldsymbol{y}, \boldsymbol{\theta}) \right \} - \ln P(\boldsymbol{y})$ in expression (ref) is not analytically tractable, because there is not closed-form expression for $\ln P(\boldsymbol{y})$. Therefore, we consider the following alternative objective function:

The term $\mathbb{E}_{q} \left \{ \ln P(\boldsymbol{y}, \boldsymbol{\theta}) \right \}- \mathbb{E}_{q} \left \{ \ln q(\boldsymbol{\theta}) \right \}$ is referred to as the evidence lower bound (ELBO). Maximizing the ELBO is equivalent to minimizing the KL divergence between the approximate variational distribution and the intractable exact posterior distribution. Consequently, the goal of VB can be re-formulated as

The functional form of the variational distribution $q(\boldsymbol{\theta})$ remains to be chosen. In principle, the complexity of the variational distribution determines the quality of the approximation of the posterior and the difficulty of the optimisation problem blei2017variational. Here, we appeal to the mean-field family of distributions jordan1999introduction, under which the variational distribution factorises as

where $j \in \{1, \ldots, J\}$ indexes the model parameters collected in $\boldsymbol{\theta}$. The mean-field assumption breaks the dependence between the model parameters by imposing mutual independence of the variational factors. It can be shown that the optimal density of each variational factor is given by

i.e. the optimal density of each variational factor is proportional to the exponentiated expectation of the logarithm of the joint distribution of $\boldsymbol{y}$ and $\boldsymbol{\theta}$, where the expectation is taken with respect to all parameters other than $\boldsymbol{\theta}_{j}$ ormerod2010explaining,blei2017variational. Provided that the model of interest is conditionally conjugate, the optimal densities of all variational factors belong to recognizable families of distributions blei2017variational. Due to the implicit nature of the expectation operator $\mathbb{E}_{- \boldsymbol{\theta}_{j}}$, the ELBO can then be maximized via a simple iterative coordinate ascent algorithm bishop2006pattern, in which the variational factors are updated one at a time conditional on the current estimates of the other variational factors. With this algorithm, iterative updates with respect to each variational factor are performed by equating each of the variational factors to its respective optimal density, i.e. we set $q(\boldsymbol{\theta}_{j}) = q^{*}(\boldsymbol{\theta}_{j})$ for $j = 1, \ldots, J$. Because the ELBO is convex with respect to each of the variational factors, the ELBO is guaranteed to converge to a local optimum boyd2004convex. Moreover, an important result from the frequentist perspective is the variational Bernstein-von Mises theorem, which states that under benign conditions, the mean-field variational Bayes estimate $\check{\boldsymbol{\theta}} = \int \boldsymbol{\theta} q^{*}(\boldsymbol{\theta}) d \boldsymbol{\theta}$ is consistent wang2018frequentist.

Finally, we observe that VB can be viewed as a tractable approximation of the expectation-maximization (EM) algorithm dempster1977maximum. To make this analogy clear, we partition the model parameters into global parameters $\boldsymbol{\theta}_{G} = \{ \boldsymbol{\alpha}, \boldsymbol{\zeta}, \boldsymbol{\Omega}, \boldsymbol{a}\}$ and local parameters (latent variables) $\boldsymbol{\theta}_{L} = \boldsymbol{\beta}_{1:N}$. Since the EM algorithm is a frequentist estimation procedure, point estimates (instead of the posterior distribution) of the global parameters $\boldsymbol{\theta}_{G}$ are of interest and are obtained by maximizing the log-likelihood via a two-step iterative procedure. In the expectation step (E-step), the distribution of local parameters conditional on the current estimates of the global parameters is calculated. In the maximization step (M-step), the conditional expectation (i.e. the lower bound on the log-likelihood) is maximized over the unknown global parameters. In Bayesian estimation, the global parameters are also treated as random variables and the posterior distribution of both the local and the global parameters is estimated. VB becomes useful when the conditional expectation relative to these parameters is intractable. Whereas the EM algorithm works with the exact conditional distribution on the local parameters, VB approximates the intractable conditional distributions of the parameters of interest with the help of a simpler, parametric variational distribution. In a similar way as the EM algorithm, VB updates the parameters of the variational distribution by iteratively maximizing the ELBO (which is analogous to the lower bound of the log-likelihood in EM); each VB iteration tightens the gap between the variational distribution and the actual posterior distribution. For more details on the connection between VB and the EM algorithm, we refer to beal2003variational.

Variational Bayes for posterior inference in mixed multinomial logit models

General strategy

In the present application, we are interested in approximating the posterior distribution of the MMNL model parameters $\{\boldsymbol{\alpha}, \boldsymbol{\zeta}, \boldsymbol{\Omega}, a_{1:K}, \boldsymbol{\beta}_{1:N} \}$ (see expression (ref)) through a fitted variational distribution. We posit a variational distribution from the mean-field family, i.e. the variational distribution factorizes as follows:

Recall that the optimal densities of the variational factors are given by $q^{*}(\theta_{i}) \propto \exp \mathbb{E}_{- \theta_{i}} \left \{ \ln P(\boldsymbol{y}, \boldsymbol{\theta}) \right \}$. We find that $q^{*}(\boldsymbol{\zeta} \lvert \boldsymbol{\mu}_{\boldsymbol_{\zeta}}, \boldsymbol{\Sigma}_{\boldsymbol{\zeta}})$, $q^{*}(\boldsymbol{\Omega} \lvert w, \boldsymbol{\Theta})$ and $q^{*}(a_{k} \lvert c, d_k)$ are common probability distributions (see Appendix (ref)). However, $q^{*}(\boldsymbol{\alpha})$ and $q^{*}(\boldsymbol{\beta}_{n})$ are not members of recognizable families of distributions, because the MNL kernel does not have a general conjugate prior. For simplicity and computational convenience, we assume that $q(\boldsymbol{\alpha}) = \text{Normal}(\boldsymbol{\mu}_{\boldsymbol{\alpha}}, \boldsymbol{\Sigma}_{\boldsymbol{\alpha}})$ and $q(\boldsymbol{\beta}_{n}) = \text{Normal}(\boldsymbol{\mu}_{\boldsymbol{\beta}_{n}}, \boldsymbol{\Sigma}_{\boldsymbol{\beta}_{n}})$ for all $n \in \{1, \ldots, N \}$. For notational convenience, we can combine the variational factors such that $q(\boldsymbol{\alpha}) q(\boldsymbol{\beta}_{n}) = q(\boldsymbol{\Gamma}_{n}) = \text{Normal}( \boldsymbol{\Gamma}_{n0}, \boldsymbol{V}_{\boldsymbol{\Gamma}_{n0}})$ with $\boldsymbol{\Gamma}_{n0} =

^{\top}$ and $\boldsymbol{V}_{\boldsymbol{\Gamma}_{n0}} =

$ for $n = 1, \dots, N$. The negative entropy of the variational distribution is given by

Moreover, the logarithm of the joint distribution of the data and the unknown model parameters is given by

Taking expectations, we obtain

Hence, the ELBO of MMNL is:

The ELBO is maximized using an iterative coordinate ascent algorithm. Iterative updates of $q(\boldsymbol{\zeta})$, $q(\boldsymbol{\Omega})$, and $q(a_{k})$ are performed by equating each variational factor to its respective optimal distribution $q^{*}(\boldsymbol{\zeta})$, $q^{*}(\boldsymbol{\Omega})$ and $q^{*}(a_{k})$, respectively. However, updates of $q(\boldsymbol{\alpha})$ and $q(\boldsymbol{\beta}_{n})$ require special treatment, because there is no closed-form expression for the expectation of the log-sum of exponentials (LSE) in equation (ref). To be precise, the LSE term is given by

and $\mathbb{E}_{q} \left \{ g_{nt}(\boldsymbol{\Gamma}_{n}) \right \}$ (henceforth, E-LSE) is not tractable.

Approximations, bounds and updating strategies

The literature proposes different methods for enabling VB for posterior inference in MMNL models with only individual-specific utility parameters (i.e. $\boldsymbol{\Gamma}_{n} = \boldsymbol{\beta}_{n}$) braun2010variational, depraetere2017comparison, tan2017stochastic. In essence, these methods proceed as follows: The E-LSE term is approximated either analytically or by simulation, or an alternative variational lower bound is defined. Then, updates for the nonconjugate variational factors are performed with the help of either quasi-Newton (QN) methods nocedal2006numerical or nonconjugate variational message passing knowles2011non.

Table (ref) provides an overview of relevant instances of VB methods for posterior inference in MMNL models and classifies these approaches according to their E-LSE approximation method or lower bound and their updating strategy. Table (ref) also shows which methods are extended in the current paper to allow for posterior inference in MMNL models with both fixed and random utility parameters. In this study, we consider one analytical approximation method, namely the Delta ($\Delta$) method bickel2015mathematical, one simulation-based approximation method, namely quasi-Monte Carlo (QMC) integration dick2010digital, as well as an alternative variational lower bound of E-LSE defined with the help of the modified Jensen's inequality knowles2011non in combination with QN- and NCVMP-based updates.\footnote{QMC methods are widely used in statistics and related areas to approximate intractable integrals by simulation. For a general treatment of QMC methods, we refer to dick2010digital. For in-depth treatments of QMC methods in the context of simulation-assisted estimation of discrete choice models, the reader is directed to bhat2001quasi, sivakumar2005simulation and train2009discrete.}\textsuperscript{,}\footnote{ In this study, we do not consider NCVMP in combination with QMC integration (henceforth, VB-NCVMP-QMC), as the calculations of the gradients of the expectations of the logarithm of the joint distribution involve inversions of large matrices. As a consequence, VB-NCVMP-QMC becomes numerically unstable and positive-definiteness of the updates of the covariance matrices $\boldsymbol{\Sigma}_{\boldsymbol{\alpha}}$ and $\boldsymbol{\Sigma}_{\boldsymbol{\beta}_{n}}$ cannot be guaranteed. The updates of the nonconjugate variational factors in VB-NCVMP-QMC can be made available upon request.}

We select the analytical and simulation-based E-LSE approximation methods and the alternative variational lower bound as well as the updating strategies for the nonconjugate variational factors based on the findings of earlier studies: tan2017stochastic also adopts the stochastic linear regression (SLR) approach salimans2013fixed for posterior inference in MMNL models with only individual-specific utility parameters. SLR is a VB variant, which involves stochastic simulations to update the variational distributions in non-conjugate models. In this paper, we do not extend VB-SLR for posterior inference in MMNL model with a more general utility specification involving a combination of fixed and random utility parameters, because it is computationally expensive to condition the iterative and simulation-based updates of one set of parameters on the approximate posterior distribution of the other set of parameters. tan2017stochastic further uses Laplace's method to approximate E-LSE and then employs QN methods to update $q(\boldsymbol{\beta}_{n})$ (henceforth, VB-QN-L). However, VB-QN-L is found to provide inferior predictive accuracy in comparison with MCMC, VB-NCVMP-$\Delta$ and VB-SLR. Moreover, braun2010variational also consider the original version of Jensen's inequality to define an alternative variational lower bound and then use QN methods to update $q(\boldsymbol{\beta}_{n})$. However, the modified Jensen's inequality proposed by knowles2011non provides a tighter lower bound. depraetere2017comparison study a variety of other quadratic lower bounds but find that these bounds are outperformed by the modified Jensen's inequality. From Table (ref), it can further be seen that the relative performance the QN- and NCVMP-based updating strategies are not known, as these updating strategies have been studied in isolation from each other.

In what follows, we describe the considered methods to approximate E-LSE and the alternative variational lower bound:

Next, we outline the updating strategies for the nonconjugate variational factors:

Algorithm (ref) succinctly summarizes the considered VB methods for posterior inference in MMNL models with a linear-in-parameters utility specification including both fixed and random utility parameters

Simulation evaluation

Data and experimental setup

For the simulation study, we devise a semi-synthetic data generating process (DGP), under which the choice sets and population parameters are based on real data from a stated choice study on consumer preferences for alternative fuel vehicles in Germany achtnicht2012german. The real data comprise 3,588 observations from 598 individuals. In the original study, respondents were presented with six choice sets, each of which consisted of seven alternatives, which in turn were characterized by six attributes, namely fuel type and propulsion technology (gasoline, diesel, hybrid, LPG/CNG, biofuel, hydrogen, electric), purchase price, operating costs, engine power, CO\textsubscript{2} emissions and fuel availability.

We generate the semi-synthetic choice data as follows: Decision-makers are assumed to be utility maximizers and to evaluate alternatives based on the utility specification $U_{ntj} = \boldsymbol{X}_{ntj,F} \boldsymbol{\alpha} + \boldsymbol{X}_{ntj,R} \boldsymbol{\beta}_{n} + \epsilon_{ntj}$. Here, $n \in \{ 1, \ldots, N \}$ indexes decision-makers, $t \in \{1, \ldots, T \}$ indexes choice occasions, and $j \in \{1, \ldots, 7 \}$ indexes alternatives. $\boldsymbol{X}_{ntj,F}$ is a row-vector of attributes for which tastes $\boldsymbol{\alpha}$ are invariant across decision-makers (gasoline, hybrid, LPG/CNG, biofuel, hydrogen, electric, purchase price); $\boldsymbol{X}_{ntj,R}$ is a row-vector of attributes for which tastes $\boldsymbol{\beta}_{n}$ are individual-specific (operating costs, engine power, CO\textsubscript{2} emissions, fuel availability). The choice sets $\boldsymbol{X}_{nt,1:7}$ with $\boldsymbol{X}_{ntj} =

$ are drawn from the real data with equal probability and with replacement. $\epsilon_{ntj}$ is a stochastic disturbance sampled from $Gumbel(0,1)$. The individual-specific taste parameters are drawn from a multivariate normal distribution, i.e. $\boldsymbol{\beta}_{n} \sim N(\boldsymbol{\zeta}, \boldsymbol{\Omega})$ for $n = 1, \ldots, N$ with $\boldsymbol{\Omega} = diag(\boldsymbol{\sigma}) \boldsymbol{\Psi} diag(\boldsymbol{\sigma})$, where $\boldsymbol{\sigma}$ is a standard deviation vector, and $\boldsymbol{\Psi}$ is a correlation matrix. The values of $\boldsymbol{\alpha}$, $\boldsymbol{\zeta}$, and $\boldsymbol{\sigma}$ are based on maximum simulated likelihood point estimates of the parameters of a mixed multinomial logit model fit to the real data. The scale of the population-level parameters is set such that the error rate is approximately 50%, i.e. in 50% of the cases decision-makers deviate from the deterministically-best alternative due to the stochastic utility component.

We consider four experimental scenarios: In scenarios 1 and 2, the fixed taste parameters and their corresponding attributes are omitted from the utility specification in the DGP, and only the individual-specific parameters are estimated. In scenarios 3 and 4, the full utility specification is used in the DGP, and both sets of taste parameters are estimated. Furthermore, the degree of correlation among individual-specific taste parameters is relatively low in scenarios 1 and 3, whereas it is relatively high in scenarios 2 and 4. In Appendix (ref), we enumerate the values of $\boldsymbol{\alpha}$, $\boldsymbol{\zeta}$, $\boldsymbol{\sigma}$, and $\boldsymbol{\Psi}$ for each experimental scenario. In each scenario, $N$ takes a value in $\{500, 2000\}$, and $T$ takes a value in $\{5, 10\}$. For each experimental scenario and combination of $N$ and $T$, we consider 20 replications, whereby the data for each replication are generated based on a different random seed.

Accuracy assessment

We employ two performance metrics to assess the accuracy of the estimation approaches:

Implementation details

We implement all estimation approaches described above by writing our own Python code\footnote{The Python code is publicly available at \url{https://github.com/RicoKrueger/bayes_mxl}.} and make an effort that the implementations of the different estimators are as similar as possible to allow for fair comparisons of estimation times. The computation of the simulated log-likelihood for MSLE and all sampling steps of the MCMC algorithm can be fully vectorized. However, VB estimation necessarily involves loops to update the variational factors pertaining to the individual-specific taste parameters. For MSLE, choice probabilities are simulated using 1,000 simulation draws generated via the Modified Latin Hypercube Sampling method hess2006use. For VB-QN-QMC, we use 64 simulation draws generated via the same method; we also explored larger numbers of simulation draws (128, 256) for VB-QN-QMC but found that increases in the number of simulation draws resulted in prohibitive estimation times. For MSLE and VB-QN, we employ the Broyden-Fletcher-Goldfarb-Shanno algorithm nocedal2006numerical included in Python's SciPy library jones2001open to carry out the numerical optimizations; the default settings of the algorithm are used and analytical or simulated gradients are supplied. To assure positive-definiteness of the covariance matrices, all numerical optimizations are in fact performed with respect to the Cholesky factors of the covariance matrices. For MCMC, the sampler is executed with two parallel Markov chains and 100,000 iterations for each chain, whereby the initial 50,000 iterations of each chain are discarded for burn-in. After burn-in, every fifth draw is retained to reduce the amount of autocorrelation in the chains. For the VB methods, we apply the same stopping criterion as tan2017stochastic: We define $\boldsymbol{\vartheta} =

^{\top}$ and let $\vartheta_{i}^{(\tau)}$ denote the $i$th element of $\boldsymbol{\vartheta}$ at iteration $\tau$. We terminate the iterative coordinate ascent algorithm, when $\delta^{(\tau)} = \operatorname*{arg\,max}_{i} \frac{\vert \vartheta_{i}^{(\tau + 1)} - \vartheta_{i}^{(\tau)}\vert }{\vert \vartheta_{i}^{(\tau)} \vert} < 0.005$. As $\delta^{(\tau)}$ can fluctuate, $\boldsymbol{\vartheta}^{(\tau)}$ is substituted by its average over the last five iterations. The simulation experiments are conducted on the Katana high performance computing cluster at the Faculty of Science, UNSW Australia.

Results

Tables (ref) to (ref) enumerate the results for scenarios 1 to 4, respectively. Each table gives the means and the standard errors of the considered performance metrics for 20 replications under different combinations of sample sizes $N \in \{500,2000\}$ and choice occasions per decision-maker $T \in \{5,10\}$. In principle, a statistical testing procedure such as ANOVA could be used to compare the performance metrics of the different methods. Here, we will simply compare mean estimates, as the standard errors are generally small.

First, we examine the impact of the sample size $N$ and the number of choice occasions per decision-maker $T$ on the performance of the estimation methods. For all methods, the mean RMSE of $\boldsymbol{\alpha}$, $\boldsymbol{\zeta}$, $\boldsymbol{\Omega}_{U}$ and $\boldsymbol{\beta}_{1:N}$ as well as the mean TVD decrease with the sample size $N$ and the number of occasions $T$. These findings numerically validate the consistency of VB methods wang2018frequentist. In our subsequent discussion, we only make explicit mention of numerical results for $\{N=2000,T=10\}$, as the comparative performance of the estimation methods is generally consistent across all combinations of $N$ and $T$.

All methods recover the mean vector $\boldsymbol{\zeta}$ and the individual-specific parameters $\boldsymbol{\beta}_{1:N}$ equally well in the considered scenarios. For example, the mean RMSE values of $\boldsymbol{\zeta}$ fall into tight intervals of $[0.0181, 0.0281]$, $[0.0191, 0.0252]$, $[0.0246, 0.0286]$, and $[0.0246,0.0291]$. Likewise, the corresponding ranges for $\boldsymbol{\beta}_{1:N}$ are $[0.7198, 0.7257]$, $[0.6929, 0.6959]$, $[0.7217, 0.7254]$, and $[0.6955, 0.6982]$. Furthermore, the results of scenarios 3 and 4 show that the fixed parameters $\boldsymbol{\alpha}$ are also recovered equally well by the considered methods. Narrow ranges of the corresponding mean RMSE values across all methods in both scenarios support this observation: $[0.0269, 0.0277]$, $[0.0298, 0.0307]$.

With the exception of the VB methods relying on the MJI-based alternative variational lower bound, all methods perform equally well at recovering the covariance matrix $\boldsymbol{\Omega}$. Excluding VB-QN-MJI and VB-NCVMP-MJI, the mean RMSE values of $\boldsymbol{\Omega}_{U}$ lie in narrow ranges of $[0.0568, 0.0800]$, $[0.0570, 0.0665]$, $[0.0711, 0.0736]$ and $[0.0572, 0.0692]$, whereas the mean RMSE values of $\boldsymbol{\Omega}_{U}$ for VB-QN-MJI and VB-NCVMP-MJI are substantially larger. Upon close inspection of the simulation results, it can be seen that that the magnitudes of the relative differences in the mean RMSE value of $\boldsymbol{\Omega}_{U}$ between the MJI-based VB methods and the other methods increase, as $N$ rises. For all methods, the recovery of $\boldsymbol{\Omega}_{U}$ ameliorates, as the number of choice occasions per decision-maker increases. Furthermore, we observe that the degree of correlation does not affect the quality of the estimation of all methods.

Next, we compare the predictive accuracy of the estimation methods. With the exception of VB-QN-MJI and VB-NCVMP-MJI, the estimation approaches perform equally well at prediction. The lower predictive accuracy of MJI-based methods can be attributed to a less accurate recovery of the covariance matrix $\boldsymbol{\Omega}$. In the majority of the considered experimental conditions, the MJI-based VB methods perform noticeably worse than the competing methods, which implies that the alternative variational lower bound defined with the help of the modified Jensen's inequality affords less accurate inferences than the analytical and simulation-based E-LSE approximations. This finding is consistent with depraetere2017comparison. We also observe that the TVD proxy for MSLE is comparable to the actual TVD calculated for the Bayesian methods.

Finally, we contrast the computational efficiency of the estimation methods. For VB, we observe that NCVMP updates are substantially faster than QN updates at virtually no compromises in parameter recovery and predictive accuracy. In contrast to earlier studies braun2010variational, depraetere2017comparison, we do not find that the QN-based VB methods are faster than MCMC, even though we use similar numbers of draws for the posterior simulations. A possible explanation for this discrepancy is that earlier studies rely on the bayesm rossi2012bayesian package for R to carry out the MCMC estimations, whereas we develop our own efficient Python implementation. Of the considered VB methods, VB-NCVMP-$\Delta$ performs best at balancing fast estimation times, acceptable parameter recovery and good predictive accuracy. Across the considered experimental conditions, VB-NCVMP-$\Delta$ is on average between 1.7 to 16.2 times faster than MCMC and MSLE, while performing nearly as well at prediction and parameter recovery. Whereas earlier studies reported occasional convergence issues for the delta-method-based E-LSE approximation depraetere2017comparison, tan2017stochastic, we encountered no such issues in the current simulation study.

Conclusions

This study extends several variational Bayes (VB) methods to allow for posterior inference in mixed multinomial logit (MMNL) models with a linear-in-parameters utility specification involving both taste parameters that vary normally across decision-makers as well as taste parameters that are invariant across decision-makers. In addition, extensive simulation-based evaluations provide new evidence into the finite-sample properties and the predictive accuracy of VB methods for MMNL in comparison with Markov chain Monte Carlo (MCMC) methods and maximum simulated likelihood estimation (MSLE). Our findings suggest that VB with nonconjugate variational message passing and a delta-method-based approximation of the expectation of log-sum of exponential (E-LSE) term (VB-NCVMP-$\Delta$) is an attractive alternative to MCMC and MSLE for fast and scalable estimation of MMNL models. The substantial gains in computational efficiency come at practically no compromises in parameter recovery and predictive accuracy.

There are several directions in which future research can build on the work presented in the current paper. First, VB methods for posterior inference in MMNL models are currently limited to MMNL models with normal mixing distributions and utility specifications in preference space. Extending VB methods to accommodate more flexible parametric, nonparametric, and semiparametric mixing distributions as well as utility specifications in willingness-to-pay space is an immediate step to support the use of VB methods in empirical applications. Second, VB methods can be devised for extended discrete choice models walker2001extended such as the integrated choice and latent variable model. As excessive estimation times continue to represent a bottleneck in empirical applications of such advanced discrete choice models, VB methods could facilitate the use of these and other behaviourally-rich models in novel contexts and applications. Third, we have shown that VB methods perform reasonably well at recovering individual-level parameters and lend themselves well to applications in which fast predictions are paramount. Thus, our analysis may inform the development of online estimation procedures that could enable near real time learning and prediction of individual preferences. Fourth, to further accelerate VB estimation for large datasets, stochastic variational inference methods can be leveraged hoffman2013stochastic, tan2017stochastic.

Adaptations of VB to other discrete choice models, new contexts and applications may benefit from fundamental advancements in the underlying VB procedure. First, in this paper, we have considered extensions to a standard VB approach, which relies on the KL divergence and the mean-field assumption. While computationally-convenient, the KL divergence is known to be a relatively loose bound, which may in turn lead to an underestimation of posterior variances zhang2018advances. Thus, other probability divergences such as $\alpha$- and $f$-divergences zhang2018advances may be explored in future work. The mean-field assumption is computationally convenient, but it restricts the flexibility of the variational distribution to an extent that the exact posterior can never be assumed by its variational approximation zhang2018advances. The quality of the variational distribution may be improved by injecting structure into the formulation of the variational distribution. This may be achieved by explicitly recognizing that some parameters are hierarchically dependent ranganath2016hierarchical. Second, Markov chain variational inference salimans2015markov, wolf2016variational seeks to combine the conceptual benefits of MCMC and VB, i.e. i.e. accurate inferences and fast estimation, respectively. Developing an MCVI method for MMNL is another potential direction for future research. Third, enhancements in the analytical and simulation-based approximation of E-LSE could lead to further improvements in the computational efficiency and quality of the VB methods. Improvements in computational efficiency may also be realized by leveraging advancements in technical computing soft- and hardware.

Finally, another avenue for future research is to contrast the VB methods considered in the current study with other emerging analytical approximation methods proposed in the frequentist context such as the Maximum Approximate Composite Marginal Likelihood (MACML) approach bhat2014new, bhat2018new, bhat2011simulation, patil2017simulation.

Acknowledgements

We are grateful to Chandra Bhat, Abdul Pinjari and two anonymous reviewers for their critical assessment of our work. We would also like to thank Martin Achtnicht for sharing the stated choice data, Tim Hillel for help with the Python implementation, and Naveen Sunder for many helpful comments on the first draft. PB and RAD are thankful to the National Science Foundation CAREER Award CBET-1253475 for financially supporting this research. PB is also thankful to Prof. Joan Walker and Prof. Kenneth Train for their guidance during his visit to UC Berkeley under the Exchange Scholar Program. RK and THR acknowledge financial support from the Australian Research Council (DE170101346). This research includes computations using the Linux computational cluster Katana supported by the Faculty of Science, UNSW Australia.

Author contribution statement

PB: conception and design, model extension and implementation (VB-NCVMP), data preparation and analysis, manuscript writing and editing. \\ RK: conception and design, model extension and implementation (VB-QN), model implementation (MCMC, MSLE, VB-NCVMP), data preparation and analysis, manuscript writing and editing. \\ MB: conception and design, manuscript editing, supervision. \\ RAD: conception and design, manuscript editing, supervision. \\ THR: conception and design, manuscript editing, supervision. \\