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.
68,941 characters · 15 sections · 53 citation commands
\thispagestyle{empty}
\thispagestyle{empty}
Variational Bayes (VB), a method originating from machine learning, enables fast and scalable estimation of complex probabilistic models. Thus far, applications of VB in discrete choice analysis have been limited to mixed logit models with unobserved inter-individual taste heterogeneity. However, such a model formulation may be too restrictive in panel data settings, since tastes may vary both between individuals as well as across choice tasks encountered by the same individual. In this paper, we derive a VB method for posterior inference in mixed logit models with unobserved inter- and intra-individual heterogeneity. In a simulation study, we benchmark the performance of the proposed VB method against maximum simulated likelihood (MSL) and Markov chain Monte Carlo (MCMC) methods in terms of parameter recovery, predictive accuracy and computational efficiency. The simulation study shows that VB can be a fast, scalable and accurate alternative to MSL and MCMC estimation, especially in applications in which fast predictions are paramount. VB is observed to be between 2.8 and 17.7 times faster than the two competing methods, while affording comparable or superior accuracy. Besides, the simulation study demonstrates that a parallelised implementation of the MSL estimator with analytical gradients is a viable alternative to MCMC in terms of both estimation accuracy and computational efficiency, as the MSL estimator is observed to be between 0.9 and 2.1 times faster than MCMC. \\ \\ Keywords: Variational Bayes; Bayesian inference; mixed logit; inter- and intra-individual heterogeneity.
\pagenumbering{arabic}
The representation of taste heterogeneity is a principal concern of discrete choice analysis, as information on the distribution of tastes is critical for demand forecasting, welfare analysis and market segmentation allenby1998marketing, ben2019foundations. From the analyst's perspective, taste variation is often random, as differences in sensitivities cannot be related to observed or observable characteristics of the decision-maker or features of the choice context bhat1998accommodating, bhat2000incorporating.
Mixed random utility models such as mixed logit mcfadden2000mixed provide a powerful framework to account for unobserved taste heterogeneity in discrete choice models. When longitudinal choice data are analysed with the help of mixed random utility models, it is standard practice to assume that tastes vary randomly across decision-makers but not across choice occasions encountered by the same individual revelt1998mixed. The implicit assumption underlying this treatment of unobserved heterogeneity is that an individual's tastes are unique and stable stigler1977gustibus. Contrasting views of preference formation postulate that preferences are constructed in an ad-hoc manner at the moment of choice bettman1998constructive or learnt and discovered through experience kivetz2008synthesis.
From the perspective of discrete choice analysis, these alternative views of preference formation justify accounting for both inter- and intra-individual random heterogeneity hess2015intra. A straightforward way to accommodate unobserved inter- and intra-individual heterogeneity in mixed random utility models is to augment a mixed logit model with a multivariate normal mixing distribution in a hierarchical fashion such that the case-specific taste parameters are generated as normal perturbations around the individual-specific taste parameters becker2018bayesian, bhat2002unified, bhat2006impact, danaf2019online, hess2009allowing, hess2011recovery, hess2015intra, yanez2011treatment. Besides, bhat2011simulation accommodate random parameters with unobserved inter- and inter-individual within the multinomial probit framework and employ the maximum approximate composite marginal likelihood bhat2011maximum approach for model estimation.
Mixed logit models with unobserved inter- and intra-individual heterogeneity can be estimated with the help of maximum simulated likelihood (MSL) estimation methods hess2009allowing, hess2011recovery. However, this estimation strategy is computationally expensive, as it involves the simulation of iterated integrals. becker2018bayesian propose a Markov chain Monte Carlo (MCMC) method, which builds on the Allenby-Train procedure train2009discrete for mixed logit models with only inter-individual heterogeneity. Notwithstanding that MCMC methods constitute a powerful framework for posterior inference in complex probabilistic models gelman2013bayesian, they are subject to several bottlenecks, which inhibit their scalability to large datasets, namely long computation times, serial correlation, high storage costs for the posterior draws and difficulties in assessing convergence bansal2020bayesian, depraetere2017comparison.
Variational Bayes methods blei2017variational, jordan1999introduction, ormerod2010explaining have emerged as a fast and computationally-efficient alternative to MCMC methods for posterior inference in discrete choice models. VB addresses the shortcomings of MCMC by recasting Bayesian inference as an optimisation problem in lieu of a sampling problem. Whilst in MCMC, the posterior distribution of interest is approximated through samples from a Markov chain, VB approximates the posterior distribution of interest through a variational distribution whose parameters are optimised such that the probability distance between the posterior distribution of interest and the approximating variational distribution is minimal. Several studies bansal2020bayesian, braun2010variational, depraetere2017comparison, tan2017stochastic derive and assess VB methods for mixed logit models with only inter-individual heterogeneity. All of these studies establish that VB is substantially faster than MCMC at practically no compromises in predictive accuracy.
Motivated by these recent advances in Bayesian estimation of discrete choice models, this paper has two objectives: First, we derive a VB method for posterior inference in mixed logit models with unobserved inter- and intra-individual heterogeneity. Second, we benchmark the proposed VB method against MSL and MCMC in a simulation study in terms of parameter recovery, predictive accuracy and computational efficiency.
The remainder of this paper is organised as follows. First, we present the mathematical formulation of mixed logit with unobserved inter- and intra-individual heterogeneity (Section (ref)). Then, we develop a VB method for this model (Section (ref)) and contrast the performance of this method against MSL and MCMC in a simulation study (Section (ref)). Last, we conclude with a summary and an outline of directions for future research (Section (ref)).
The mixed logit model with unobserved inter- and intra-individual heterogeneity hess2009allowing, hess2011recovery is established as follows: On choice occasion $t \in \{1, \ldots T \}$, a decision-maker $n \in \{1, \ldots N \}$ derives utility
from alternative $j$ in the set $C_{nt}$. Here, $V()$ denotes the representative utility, $\boldsymbol{X}_{ntj}$ is a row-vector of covariates, $\boldsymbol{\beta}_{nt}$ is a collection of taste parameters, and $\epsilon_{ntj}$ is a stochastic disturbance. The assumption $\epsilon_{ntj} \sim \text{Gumbel}(0,1)$ leads to the multinomial logit (MNL) model such that the probability that decision-maker $n$ chooses alternative $j \in C_{nt}$ on choice occasion $t$ can be expressed as
where $y_{ntj}$ equals $1$ if alternative $j \in C_{nt}$ is chosen and zero otherwise.
In equation (ref), the taste parameters $\boldsymbol{\beta}_{nt}$ are specified as observation-specific. To allow for dependence between repeated observations for the same individual and to accommodate inter-individual taste heterogeneity, it has become standard practice to adopt Revelt's and Train's (revelt1998mixed) panel estimator for the mixed logit model. Under this specification, taste homogeneity across replications is assumed such that $\boldsymbol{\beta}_{nt} = \boldsymbol{\beta}_{n}$ $\forall t \in \{ 1, \ldots, T \}$. To accommodate intra-individual taste heterogeneity in addition to inter-individual taste heterogeneity, the taste vector $\boldsymbol{\beta}_{nt}$ can be defined as a normal perturbation around an individual-specific parameter $\boldsymbol{\mu}_{n}$, i.e. $\boldsymbol{\beta}_{nt} \sim \mbox{N}(\boldsymbol{\mu}_{n}, \boldsymbol{\Sigma}_{W})$ for $t = 1, \ldots, T$, where $\boldsymbol{\Sigma}_{W}$ is a covariance matrix. The distribution of individual-specific parameters $\boldsymbol{\mu}_{1:N}$ is then also assumed to be multivariate normal, i.e. $\boldsymbol{\mu}_{n} \sim \text{N}(\boldsymbol{\zeta}, \boldsymbol{\Sigma}_{B})$ for $n = 1, \dots, N$, where $\boldsymbol{\zeta}$ is a mean vector and $\boldsymbol{\Sigma}_{B}$ is a covariance matrix.
Under a fully Bayesian approach, the parameters $\boldsymbol{\zeta}$, $\boldsymbol{\Sigma}_{B}$, $\boldsymbol{\Sigma}_{W}$ are considered to be random, unknown quantities and are thus given priors. Here, we use a normal prior for mean vector $\boldsymbol{\zeta}$ and a marginally-noninformative Huang's half-t prior akinc2018Bayesian, huang2013Simple for the covariance matrices $\boldsymbol{\Sigma}_{B}$ and $\boldsymbol{\Sigma}_{W}$.
Stated succinctly, the generative process of mixed logit with unobserved inter- and intra-individual heterogeneity is as follows:
where $\{ \boldsymbol{\xi}_{0},\boldsymbol{\Xi}_{0}, \nu_{B}, \nu_{W}, A_{B,1:K}, A_{W,1:K} \}$ are known hyper-parameters, and $\boldsymbol{\theta} = \{\boldsymbol{a}_{B}, \boldsymbol{a}_{W},\boldsymbol{\Sigma}_{B}, \boldsymbol{\Sigma}_{W}, \boldsymbol{\zeta}, \allowbreak \boldsymbol{\mu}_{1:N}, \allowbreak \boldsymbol{\beta}_{1:N,1:T_n}\}$ is a collection of model parameters whose posterior distribution we wish to estimate.
The generative process (expressions (ref)--(ref)) implies the following joint distribution of the data and the model parameters:
where $\omega_{B} = \nu_{B} + K - 1$, $\boldsymbol{B}_{B} = 2\nu_{B} \text{diag}(\boldsymbol{a}_{B})$, $\omega_{W} = \nu_{W} + K - 1$, $\boldsymbol{B}_{W} = 2\nu_{W} \text{diag}(\boldsymbol{a}_{W})$, $s = \frac{1}{2}$, $r_{B,k} = A_{B,k}^{-2}$ and $r_{W,k} = A_{W,k}^{-2}$. By Bayes' rule, the posterior distribution of interest is given by
Exact inference of this posterior distribution is not possible, because the model evidence $\int P (\boldsymbol{y}_{1:N}, \boldsymbol{\theta}) \allowbreak d \boldsymbol{\theta}$ is not tractable. Hence, we resort to approximate inference methods. becker2018bayesian propose a Markov chain Monte Carlo (MCMC method for posterior inference in the described model. This method is presented in Appendix (ref). In the subsequent section, we derive a variational Bayes (VB) method for scalable inference. Mixed logit with unobserved inter- and intra-individual heterogeneity can also be estimated in a frequentist way using maximum simulated likelihood (MSL) estimation hess2009allowing, hess2011recovery. The MSL method is described in Appendix (ref).
This section is divided into two subsections: In Section (ref), we outline foundational principles of posterior inference with variational Bayes (VB), building on related work by bansal2020bayesian.\footnote{We also direct the reader to the statistics and machine learning literature blei2017variational, ormerod2010explaining, zhang2018advances for more technical reviews of VB estimation.} In Section (ref), we derive a VB method for posterior inference in mixed logit with unobserved inter- and intra-individual heterogeneity.
The general idea underlying VB estimation is to recast posterior inference as an optimisation problem. VB approximates an intractable posterior distribution $P(\boldsymbol{\theta} \vert \boldsymbol{y}) = \frac{P(\boldsymbol{\theta} , \boldsymbol{y})}{\int P(\boldsymbol{\theta} , \boldsymbol{y}) d \boldsymbol{\theta}}$ through a parametric variational distribution $q(\boldsymbol{\theta} \vert \boldsymbol{\nu})$. The optimisation problem consists of manipulating the parameters $\boldsymbol{\nu}$ of the variational distribution such that the probability distance between the posterior distribution of interest and the approximating variational distribution is minimal. In contradistinction, MCMC treats posterior inference as a sampling problem, i.e. the posterior distribution of interest is approximated through samples from a Markov chain whose stationary distribution is the posterior distribution of interest.
Translating posterior inference into an optimisation problem effectively addresses the shortcomings of MCMC bansal2020bayesian. First, only the current estimates of the variational parameters rather than thousands of posterior draws from the Markov chains need to be stored. Second, convergence can be easily assessed by considering the change of the variational lower bound or the estimates of the variational parameters between successive iterations. Last, serial correlation becomes irrelevant, because no samples are taken.
The Kullback-Leibler (KL) divergence kullback1951information provides a measure of the probability distance between the variational distribution $q(\boldsymbol{\theta})$ and the posterior distribution $P(\boldsymbol{\theta} \vert \boldsymbol{y})$. It is defined as
VB seeks to minimise this divergence. Yet, $D_{\text{KL}} \left (q(\boldsymbol{\theta}) \vert \vert P(\boldsymbol{\theta} \vert \boldsymbol{y}) \right )$ is not analytically tractable, because the term $\ln P(\boldsymbol{y})$ lacks a closed-form expression. For this reason, we define an alternative variational lower bound, which is referred to as the evidence lower bound (ELBO):
The ELBO is $\ln P(\boldsymbol{y})$ (which does not depend on $\boldsymbol{\theta}$) minus the KL divergence. Minimising the KL divergence between the variational distribution and the posterior distribution of interest (expression (ref)) is equivalent to maximising the ELBO (expression (ref)). Thus, the goal of VB can be re-formulated as
The functional form of the variational distribution $q(\boldsymbol{\theta})$ needs to be configured by the analyst. The complexity of the variational distribution influences both the expressiveness of the variational distribution (and thus the quality of the approximation to the posterior) as well as the difficulty of the optimisation problem. A computationally-convenient family of variational distributions is the mean-field family of distributions jordan1999introduction. The mean-field assumption imposes independence between blocks of model parameters such that the variational distribution exhibits the form $q(\boldsymbol{\theta}_{1:B}) = \prod_{b=1}^{B} q(\boldsymbol{\theta}_{b})$, where $b \in \{1, \ldots, B\}$ indexes the independent blocks of model parameters. Under the mean-field assumption, the optimal density of each variational factor is $q^{*}(\boldsymbol{\theta}_{b}) \propto \exp \mathbb{E}_{- \boldsymbol{\theta}_{b}} \left \{ \ln P(\boldsymbol{y}, \boldsymbol{\theta}) \right \} $, 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}$, whereby expectation is calculated with respect to all parameters other than $\boldsymbol{\theta}_{b}$ blei2017variational, ormerod2010explaining. If the model of interest has a conditionally-conjugate structure, the optimal densities of all variational factors are known probability distributions. In this case, the variational objective can be maximised with the help of a simple iterative coordinate ascent algorithm bishop2006pattern, which involves updating the variational factors one at a time conditionally on the current estimates of the other variational factors.
To apply VB to posterior inference in mixed logit with unobserved inter- and intra-individual heterogeneity, we re-formulate the generative process of the model (see expressions (ref)--(ref)) such that the hierarchical dependence between the individual- and the observation-specific taste parameters $\boldsymbol{\mu}_{n}$ and $\boldsymbol{\beta}_{nt}$ is removed. We let
and define
This modification does not change the model but alters the updates in the VB procedure. If the hierarchical dependence between the individual- and the observation-specific taste parameters is preserved, we find that the non-informative structure of $q^{*}(\boldsymbol{\mu}_{n})$ results in a severe underestimation of the inter-individual covariance $\boldsymbol{\Sigma}_{B}$.\footnote{If $\boldsymbol{\mu}_{n}$ and $\boldsymbol{\beta}_{nt}$ are treated as hierarchically dependent, the optimal density of the conjugate variational factors pertaining to the individual-specific parameters $\boldsymbol{\mu}_{n}$ is $q^{*}(\boldsymbol{\mu}_{n}) \propto \text{Normal}(\boldsymbol{\mu}_{\boldsymbol{\mu}_{n}}, \boldsymbol{\Sigma}_{\boldsymbol{\mu}_{n}})$ with $\boldsymbol{\Sigma}_{\boldsymbol{\mu}_{n}} = \left ( \mathbb{E}_{- \boldsymbol{\mu}_{n}} \left \{ \boldsymbol{\Sigma}_{B}^{-1} \right \} + T \mathbb{E}_{- \boldsymbol{\mu}_{n}} \left \{ \boldsymbol{\Sigma}_{W}^{-1} \right \} \right )^{-1}$ and $\boldsymbol{\mu}_{\boldsymbol{\mu}_{n}} = \boldsymbol{\Sigma}_{\boldsymbol{\mu}_{n}} \Big ( \mathbb{E}_{- \boldsymbol{\mu}_{n}} \big \{ \boldsymbol{\Sigma}_{B}^{-1} \big \} \mathbb{E}_{- \boldsymbol{\mu}_{n}} \{\boldsymbol{\zeta} \} + \allowbreak \mathbb{E}_{- \boldsymbol{\mu}_{n}} \big \{ \boldsymbol{\Sigma}_{W}^{-1} \big \} \allowbreak \sum_{t = 1}^{T} \mathbb{E}_{- \boldsymbol{\mu}_{n}} \boldsymbol{\beta}_{nt} \Big )$, where $\boldsymbol{\Sigma}_{\boldsymbol{\mu}_{n}}$ only varies across individuals, if the number of choice occasions differs across individuals.}
The modified posterior distribution is
with $\boldsymbol{\theta} = \{\boldsymbol{a}_{B}, \boldsymbol{a}_{W},\boldsymbol{\Sigma}_{B}, \boldsymbol{\Sigma}_{W}, \boldsymbol{\zeta}, \boldsymbol{\mu}_{1:N}, \allowbreak \boldsymbol{\gamma}_{1:N,1:T_n}\}$ and
We wish to approximate this posterior distribution through a fitted variational distribution. To this end, we posit a variational distribution from the mean-field family such that
Under the mean-field assumption, the optimal densities of the variational factors are given by $ q^{*}(\boldsymbol{\theta}_{j}) \propto \exp \mathbb{E}_{- \boldsymbol{\theta}_{i}} \left \{ \ln P(\boldsymbol{y}, \boldsymbol{\theta}) \right \} $. We find that $q^{*}(a_{B,k} \vert c_{B}, d_{B,k})$, $q^{*}(a_{W,k} \vert c_{W}, d_{W,k})$, $q^{*}(\boldsymbol{\Sigma}_{B} \lvert w_{B}, \boldsymbol{\Theta}_{B})$, $q^{*}(\boldsymbol{\Sigma}_{W} \lvert w_{W}, \boldsymbol{\Theta}_{W})$ and $q^{*}(\boldsymbol{\zeta} \lvert \boldsymbol{\mu}_{\boldsymbol_{\zeta}}, \allowbreak \boldsymbol{\Sigma}_{\boldsymbol{\zeta}})$ are common probability distributions (see Appendix (ref)). However, $q^{*}(\boldsymbol{\mu}_{n})$ and $q^{*}(\boldsymbol{\gamma}_{nt})$ are not members of recognisable families of distributions, because the multinomial logit model lacks a general conjugate prior. For computational convenience, we assume $q(\boldsymbol{\mu}_{n}) = \text{Normal}(\boldsymbol{\mu}_{\boldsymbol{\mu}_{n}}, \boldsymbol{\Sigma}_{\boldsymbol{\mu}_{n}})$ and $q(\boldsymbol{\gamma}_{nt}) = \text{Normal}(\boldsymbol{\mu}_{\boldsymbol{\gamma}_{nt}}, \boldsymbol{\Sigma}_{\boldsymbol{\gamma}_{nt}})$.
The evidence lower bound (ELBO) of mixed logit with unobserved inter- and intra-individual heterogeneity is presented in Appendix (ref). Because of the mean-field assumption, the ELBO can maximised via an iterative coordinate ascent algorithm. Iterative updates of $q(a_{B,k})$, $q(a_{W,k})$, $(\boldsymbol{\Sigma}_{B})$, $(\boldsymbol{\Sigma}_{W})$ and $(\boldsymbol{\zeta})$ are performed by equating each variational factor to its respective optimal distribution $q^{*}(a_{B,k})$, $q^{*}(a_{W,k})$, $q^{*}(\boldsymbol{\Sigma}_{B})$, $q^{*}(\boldsymbol{\Sigma}_{W})$, $q^{*}(\boldsymbol{\zeta})$.
Yet, updates of $q(\boldsymbol{\mu}_{n})$ and $q(\boldsymbol{\gamma}_{nt})$ demand special treatment, because there is no closed-form expression for the expectation of the log-sum of exponentials (E-LSE) term in equation (ref). The intractable E-LSE term is given by
with
bansal2020bayesian analyse several methods to approximate the E-LSE term and to update the variational factors corresponding to utility parameters in the context of mixed logit with unobserved inter-individual heterogeneity. In this paper, we use quasi-Monte Carlo (QMC) integration bhat2001quasi, dick2010digital, sivakumar2005simulation, train2009discrete to approximate the E-LSE term and then use quasi-Newton (QN) methods nocedal2006numerical to update $q(\boldsymbol{\mu}_{n})$ and $q(\boldsymbol{\gamma}_{nt})$, as we find that the alternative methods discussed in bansal2020bayesian perform poorly for mixed logit with unobserved inter- and intra-individual heterogeneity.
With QMC integration, the intractable E-LSE terms are approximated by simulation. We have
where $\boldsymbol{\mu}_{n,d} = \boldsymbol{\mu}_{\boldsymbol{\mu}_{n}} + \mbox{chol}(\boldsymbol{\Sigma}_{\boldsymbol{\mu}_{n}}) \boldsymbol{\xi}^{(\boldsymbol{\mu})}_{n,d}$ and $\boldsymbol{\gamma}_{nt,d} = \boldsymbol{\mu}_{\boldsymbol{\gamma}_{nt}} + \mbox{chol}(\boldsymbol{\Sigma}_{\boldsymbol{\gamma}_{nt}}) \boldsymbol{\xi}^{(\boldsymbol{\gamma})}_{nt,d}$. Moreover, $\boldsymbol{\xi}^{(\boldsymbol{\mu})}_{n,d}$ and $\boldsymbol{\xi}^{(\boldsymbol{\gamma})}_{nt,d}$ denote standard normal simulation draws for the inter- and the intra-individual taste parameters, respectively. Note that the simulation approximation given in expression (ref) is much simpler than the simulation approximation required for MSL estimation (see Appendix (ref)), because the intra-individual draws need not be conditioned on the inter-individual draws, which is an immediate consequence of the mean-field assumption.
Updates of $q(\boldsymbol{\mu}_{n})$ and $q(\boldsymbol{\gamma}_{nt})$ can then be performed with the help of quasi-Newton methods. We have
and
whereby the intractable E-LSE terms $\mathbb{E}_{q} \left \{ g_{nt}(\boldsymbol{\mu}_{n}, \boldsymbol{\gamma}_{nt}) \right \}$ are replaced by the approximation given in expression (ref).
The algorithm presented in Figure (ref) succinctly summarises the VB method for posterior inference in mixed logit models with unobserved inter- and intra-individual heterogeneity.
For the simulation study, we rely on synthetic choice data, which we generate as follows: The choice sets comprise five unlabelled alternatives, which are characterised by four attributes. Decision-makers are assumed to be utility maximisers and to evaluate the alternatives based on the utility specification
Here, $n \in \{ 1, \ldots, N \}$ indexes decision-makers, $t \in \{1, \ldots, T \}$ indexes choice occasions, and $j \in \{1, \ldots, 5 \}$ indexes alternatives. Furthermore, $\boldsymbol{X}_{ntj}$ is a row-vector of uniformly distributed alternative-specific attributes, and $\epsilon_{ntj}$ is a stochastic disturbance sampled from $\text{Gumbel}(0,1)$. We consider two experimental scenarios for the generation of the case-specific taste parameters $\boldsymbol{\beta}_{nt}$. In both scenarios, $\boldsymbol{\beta}_{nt}$ are drawn via the following process:
where $\boldsymbol{\Sigma}_{B} = \text{diag}(\boldsymbol{\sigma}_{B}) \boldsymbol{\Omega}_{B} \text{diag}(\boldsymbol{\sigma}_{B})$ and $\boldsymbol{\Sigma}_{W} = \text{diag}(\boldsymbol{\sigma}_{W}) \boldsymbol{\Omega}_{W} \text{diag}(\boldsymbol{\sigma}_{W})$. Here, $\{ \boldsymbol{\sigma}_{B}, \boldsymbol{\sigma}_{W} \}$ represent standard deviation vectors and $\{ \boldsymbol{\Omega}_{B}, \boldsymbol{\Omega}_{W} \}$ are correlation matrices. In each scenario, we vary the degree of correlation across the inter- and intra-individual taste parameters. In scenario 1, the degree of correlation is relatively low, whereas it is relatively high in scenario 2. In both scenarios, we let $\boldsymbol{\sigma}_{B}^{2} = 2 \cdot \frac{2}{3} \cdot \vert \boldsymbol{\zeta} \vert $ and $\boldsymbol{\sigma}_{W}^{2} = 2 \cdot \frac{1}{3} \cdot \vert \boldsymbol{\zeta} \vert$, i.e. the total variance of each random parameter is twice the absolute value of its mean, whereby two thirds of the total variance are due to inter-individual taste variation, and one third of the total variance is due to intra-individual taste variation. The assumed values of $\boldsymbol{\zeta}$, $\boldsymbol{\Omega}_{B}$ and $\boldsymbol{\Omega}_{W}$ for each scenario are enumerated in Appendix (ref). In both scenarios, the alternative-specific attributes $\boldsymbol{X}_{ntj}$ are drawn from $\text{Uniform}(0, 2)$, which implies an error rate of approximately 50%, i.e. in 50% of the cases decision-makers deviate from the deterministically-best alternative due to the stochastic utility component. In each scenario, $N$ takes a value in $\{250, 1000 \}$, and $T$ takes a value in $\{8, 16\}$. For each experimental scenario and for each combination of $N$ and $T$, we consider 30 replications, whereby the data for each replication are generated using a different random seed.
We evaluate the accuracy of the estimation approaches in terms of their ability to recover parameters in finite sample and in terms of their predictive accuracy.
To assess how the estimation approaches perform at recovering parameters, we calculate the root mean square error (RMSE) for selected parameters, namely the mean vector $\boldsymbol{\zeta}$ and the unique elements $\{ \boldsymbol{\Sigma}_{B,U}, \boldsymbol{\Sigma}_{W,U} \}$ of the covariance matrices $\{ \boldsymbol{\Sigma}_{B}, \boldsymbol{\Sigma}_{W} \}$. Given a collection of parameters $\boldsymbol{\theta}$ and its estimate $\hat{\boldsymbol{\theta}}$, RMSE is defined as
where $J$ denotes the total number of scalar parameters collected in $\boldsymbol{\theta}$. For MSL, point estimates of $\boldsymbol{\zeta}$, $\boldsymbol{\Sigma}_{B}$ and $\boldsymbol{\Sigma}_{W}$ are directly obtained. For MCMC, estimates of the parameters of interest are given by the means of the respective posterior draws. For VB, we have $\hat{\boldsymbol{\zeta}} = \boldsymbol{\mu}_{\boldsymbol{\zeta}}$, $\hat{\boldsymbol{\Sigma}}_{B} = \frac{\boldsymbol{\Theta}_{B}}{w_{B} - K - 1}$, and $\hat{\boldsymbol{\Sigma}}_{W} = \frac{\boldsymbol{\Theta}_{W}}{w_{W} - K - 1}$. As our aim is to evaluate how well the estimation methods perform at recovering the distributions of the realised individual- and observation-specific parameters $\{ \boldsymbol{\mu}_{1:N}, \boldsymbol{\beta}_{1:N,1:T} \}$, we use the sample mean $\boldsymbol{\zeta}_{0} = \frac{1}{N} \sum_{n = 1}^{N} \boldsymbol{\mu}_{n}$ and the sample covariances $\boldsymbol{\Sigma}_{B,0} = \frac{1}{N} \sum_{n = 1}^{N} (\boldsymbol{\mu}_{n} - \boldsymbol{\zeta}_{0}) (\boldsymbol{\mu}_{n} - \boldsymbol{\zeta}_{0})^{\top}$ and $\boldsymbol{\Sigma}_{W,0} = \frac{1}{NT} \sum_{n = 1}^{N} \sum_{t = 1}^{T} (\boldsymbol{\beta}_{nt} - \boldsymbol{\mu}_{n}) (\boldsymbol{\beta}_{nt} - \boldsymbol{\mu}_{n})^{\top}$ as true parameter values for $\boldsymbol{\zeta}$, $\boldsymbol{\Sigma}_{B}$ and $\boldsymbol{\Sigma}_{W}$, respectively.
To assess the predictive accuracy of the Bayesian methods, we consider two out-of-sample prediction scenarios. In the first scenario, we predict choice probabilities for a new set of individuals, i.e. we predict between individuals. In the second scenario, we predict choice probabilities for new choice sets for individuals who are already in the sample, i.e. we predict within individuals. For each of these scenarios, we calculate the total variation distance braun2010variational between the true and the estimated predictive choice distributions.
We highlight three important advantages of using TVD as a measure of predictive accuracy. First, TVD is a strictly proper scoring rule gneiting2007strictly in that TVD for a given choice set is exclusively minimised by the true predictive choice distribution as opposed to alternative measures of predictive accuracy such as the hit rate. Second, TVD allows for a robust assessment of predictive accuracy, as the estimated predictive choice distribution is compared against the true predictive choice distribution rather than against a single realised choice. Third, TVD fully accounts for the posterior uncertainty captured by the estimation methods. Disregarding posterior uncertainty leads to tighter predictions but also gives a false sense of precision, in particular when the posterior variance of relevant model parameters is large.
We proceed as follows to calculate TVD in the two prediction scenarios:
To quantify the benefits of accommodating intra-individual heterogeneity in addition to inter-individual heterogeneity, we also estimate mixed logit models with only inter-individual heterogeneity via MCMC bansal2020bayesian and compute $\text{TVD}_{B}$ and $\text{TVD}_{W}$ for these models.
We implement the MSL, MCMC and VB estimators by writing our own Python code.\footnote{The code is publicly available at \url{https://github.com/RicoKrueger/mixed logit_inter_intra}.}
For MSL and VB, the numerical optimisations are carried out with the help of the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm nocedal2006numerical contained in Python's SciPy library jones2001open and analytical gradients are supplied. Details concerning the implementation of MSL including the required gradient expressions are provided in Appendix (ref). For MSL, we use 200 inter-individual simulation draws per decision-maker and 200 intra-individual simulation draws per observation. For VB, we use 100 inter-individual and 100 intra-individual simulation draws. For both methods, the simulation draws are generated via the Modified Latin Hypercube sampling (MLHS) approach hess2006use. To assure that the covariance matrices maintain positive-definiteness, the optimisations are performed with respect to the Cholesky factors of the covariance matrices. For VB, 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 replaced by its average over the last five iterations.
For MSL and VB, we also take advantage of Python's parallel processing capacities to improve the computational efficiency of the methods. For MSL, we process the likelihood computations in eight parallel batches, each of which corresponds to 25 inter-individual simulation draws. For VB, the updates corresponding to the individual- and the observation-specific parameters are processed in two or eight parallel batches, each of which comprises the observations of 125 individuals.
For MCMC, the sampler is executed with two parallel Markov chains and 400,000 iterations for each chain, whereby the initial 200,000 iterations of each chain are discarded for burn-in. After burn-in, only every tenth draw is retained to moderate storage requirements and to facilitate post-simulation computations. For mixed logit with only inter-individual heterogeneity, the MCMC 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 kept.
The simulation experiments are conducted on the Katana high performance computing cluster at the Faculty of Science, UNSW Australia.
In Tables (ref) and (ref), we enumerate the simulation results for scenarios 1 and 2, respectively. In each table, we report the means and standard errors of the considered performance metrics across 30 replications for different combinations of sample sizes $N \in \{250,1000\}$ and choice occasions per decision-maker $T \in \{8, 16\}$.
First, we examine the finite-sample properties of the estimators for mixed logit with unobserved inter- and intra-individual heterogeneity. By and large, MSL, MCMC and VB perform equally well at recovering the mean vector $\boldsymbol{\zeta}$ and the covariance matrix $\boldsymbol{\Sigma}_{B}$. For large samples ($N = 1000$), MSL and MCMC perform slightly better than VB at recovering $\boldsymbol{\Sigma}_{B}$, when the number of choice occasions per decision-maker is small ($T = 8$). However, the differences between the methods diminish, as $T$ increases. Furthermore, VB outperforms MCMC and MSL at recovering $\boldsymbol{\Sigma}_{W}$ in the majority of the considered experimental settings. The differences between VB and the two competing methods in the recovery of $\boldsymbol{\Sigma}_{W}$ are particularly pronounced, when the sample size is small ($N = 250$). Overall, the MSL and MCMC estimators perform equally well at recovering the parameters of interest in the considered experimental scenarios. This observation is largely consistent with becker2018bayesian.
Next, we contrast the between-individuals predictive accuracy of the Bayesian methods. Note that in Tables (ref) and (ref), $\text{TVD}_{B}$ is reported in percent. Overall, the differences in between-individuals predictive accuracy between MCMC and VB for mixed logit with unobserved inter-and intra-individual heterogeneity are negligibly small, as the absolute difference in mean $\text{TVD}_{B}$ does not exceed one permille in any of the considered experimental settings. In small samples ($N = 250$), VB performs slightly better than MCMC, but the differences between the methods diminish, as $T$ increases. By contrast, in large samples ($N = 1000$), MCMC performs marginally better than VB. Furthermore, we observe that in all of the considered experimental settings MCMC for mixed logit with only inter-individual heterogeneity is outperformed by the competing methods in term of between-individuals predictive accuracy. This suggests that for the considered data generating process, under which two thirds of the total variance are due to between-individual heterogeneity and the remaining third of the total variance is due to intra-individual heterogeneity, more accurate between-individuals predictions can be produced, if both inter- and intra-individual heterogeneity are accounted for.
Interestingly, accounting for intra-individual heterogeneity in addition to inter-individual heterogeneity does not result in more accurate within-individuals predictions in any of the considered experimental scenarios. This is because in a parametric hierarchical model, the estimates of the unit-specific parameters are biased towards the population mean and thus do not accurately reflect the true value of the underlying parameter danaf2019online. From the results presented in Tables (ref) and (ref), it can be seen that this bias persists, even if the number of choice tasks per decision-maker is relatively large ($T = 16$).
Finally, we compare the computational efficiency of the estimation methods for mixed logit with unobserved inter- and intra-individual heterogeneity. Across all of the considered experimental settings, VB is on average between 4.6 and 15.1 times faster than MCMC and between 2.8 and 17.7 times faster than MSL. Interestingly, the MSL method (which uses analytical gradients in combination with parallel processing) is either substantially faster than MCMC or only marginally slower. In small samples ($N = 250$), MSL estimation outperforms MCMC and is between 1.4 and 2.1 times faster than MCMC. In large samples ($N = 1000$), MSL is on par with MCMC and is between 0.9 and 1.1 times faster than MCMC. By contrast, becker2018bayesian observe that MSL estimation is substantially slower than MCMC. However, their implementation of MSL estimation relies on numerical gradients and does not use parallel processing. We also highlight that during MCMC estimation, large files for the storage of the posterior draws of $\boldsymbol{\zeta}$, $\boldsymbol{\Sigma}_{B}$, $\boldsymbol{\Sigma}_{W}$ and $\boldsymbol{\mu}_{n}$ for $n = 1, \ldots N$ are produced. For the large samples ($N = 1000$), the posterior draws of both Markov chains require approximately 1.3 Gigabyte of disk space. This is in spite of the fact that only every tenth posterior draw is retained after burn in and also in spite of the fact that the efficient Hierarchical Data Format hdf5 is used for the storage of the posterior draws. Both VB and MSL estimation do not consume any disk space during estimation, as no posterior draws need to be stored.
Motivated by recent advances in scalable Bayesian inference for discrete choice models, this paper derives a variational Bayes (VB) method for posterior inference in mixed logit models with unobserved inter- and intra-individual heterogeneity. The proposed VB method uses quasi-Monte Carlo (QMC) integration to approximate the intractable expectations of the log-sum of exponentials term and relies on quasi-Newton methods to update the variational factors of the utility parameters. In a simulation study, we benchmark the performance of the proposed VB method against maximum simulated likelihood (MSL) estimation and Markov chain Monte Carlo (MCMC) in terms of parameter recovery, predictive accuracy and computational efficiency.
In summary, the simulation study demonstrates that VB can be an attractive alternative to MSL and MCMC for the computationally-efficient estimation of mixed logit models with unobserved inter- and intra-individual heterogeneity. In the majority of the considered experimental settings, VB performs at least as well as MSL and MCMC at parameter recovery and out-of-sample prediction, while being between 2.8 and 17.7 times faster than the competing methods. The simulation study also shows that MSL with analytical gradients and parallelised log-likelihood and gradient computations represents a viable alternative to MCMC in terms of both parameter recovery and computational efficiency. This finding stands in contrast to the study by becker2018bayesian. In their study, the authors observe that MSL estimation is substantially slower than MCMC, Yet, their implementation of MSL estimation is not parallelised and does not use analytical gradients.
There are four main directions in which future research can build on the current paper. One direction for future research is to further improve the accuracy of VB. As highlighted by bansal2020bayesian, a promising research direction in this regard is to combine the specific benefits of MCMC and VB by marrying the two approaches in an integrated framework salimans2015markov, wolf2016variational. Another direction for future work is to contrast the performance of the VB, MCMC and MSL estimators for mixed logit with unobserved inter- and intra-individual heterogeneity with the maximum approximate composite marginal likelihood bhat2011maximum estimator for an equivalent mixed probit model bhat2011simulation. A third avenue for future research is to explore how discrete choice models can be formulated so that they can provide accurate individualised recommendations based on personalised predictions. Our analysis suggests that mixed logit models that account for both unobserved inter- and intra-individual heterogeneity do not provide benefits over simpler mixed logit models that only account for inter-individual heterogeneity in regard to within-individual predictions. In the machine learning literature, filtering approaches are used to generate personalised recommendations in recommender systems koren2009matrix. Thus, an intriguing avenue for future research is to incorporate such filtering approaches into discrete choice models donnelly2019counterfactual, zhu2019personalized. Finally, a fourth avenue for future research is to investigate the behavioural sources of what can be detected as intra-individual taste variation and to incorporate explicit representations of the underlying behavioural processes into discrete choice models balbontin2019better.
RK: conception and design, method derivation and implementation, data preparation and analysis, manuscript writing and editing. \\ PB: conception and design, method derivation, data analysis, manuscript editing. \\ MB: conception and design, manuscript editing, supervision. \\ RAD: conception and design, manuscript editing, supervision. \\ THR: conception and design, manuscript editing, supervision, funding acquisition, project administration, resources. \\