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.
69,931 characters · 23 sections · 41 citation commands
Fast variational Bayes methods for multinomial probit models
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 \fi
\if10 {
} \fi
{\it Keywords:} Multinomial probit model, Variational inference, Large choice data sets
\spacingset{1.8}
The multinomial probit (MNP) model is a popular tool for analyzing choice behavior, with recent applications including brand choices miyazaki2021dynamic, employment choices mishkin2021gender, and car parking choices paleti2018generalized. The main advantage of the MNP model is the relaxation of the independence of irrelevant alternatives assumption made by multinomial logit models. The MNP achieves this by specifying the conditional covariance matrix of the latent utilities of the choice alternatives. However, estimation of the MNP model is computationally costly. The evaluation of the likelihood function involves high-dimensional integrals, which can be solved by simulation methods. In this paper we propose a variational Bayes (VB) method for estimation in the MNP model, which is accurate and fast even when applied to choice sets that have a large number of choice alternatives and a large number of observations.
Bayesian analysis of the MNP model has reduced the computational complexity of parameter estimation, but its applicability to modern choice data sets is still limited. Bayesian estimation of the MNP model augments the likelihood function with a set of latent utilities, which are then generated inside a Markov chain Monte Carlo (MCMC) scheme albert1993bayesian. This approach avoids the computationally costly step of directly calculating the choice probabilities in the likelihood function via numerical integration. Instead, each MCMC iteration draws a vector of latent utilities from a truncated normal distribution for each observation mcculloch1994exact. Since these draws are highly auto-correlated, a large amount of iterations are required to achieve convergence. Therefore, this approach still has a substantial computational burden, especially when the number of observations or the number of choice alternatives is large.
VB is a computationally scalable alternative to MCMC. Instead of sampling from the posterior, VB calibrates a parametric approximating density by minimizing a divergence function to the posterior. Applying VB to the MNP model poses two main challenges. First, parameter identification in the MNP model requires a restriction on the covariance matrix bunch1991estimability. burgette2012trace show that a trace restriction has the best performance, and introduce a “working parameter” that rescales the trace within the MCMC algorithm. This parameter is not part of the model specification and not identified given the data. Hence an analytical expression for the unnormalized posterior density, which is required for the implementation of VB, is not available. Second, to reduce its computational complexity, the posterior density of the MNP model has to be augmented with a large number of latent variables. Existing VB methods often make strong assumptions on the approximating densities for the latent variables westling2019beyond.
Because of these challenges, existing VB approaches are designed for restrictive specifications of the MNP model. For instance, girolami2006variational model the latent utilities as independent Gaussian processes, that only allow for correlations that are a function of regressors. fasano2022class conduct VB for the coefficients conditional on a fixed covariance matrix. Moreover, both papers impose strong independence assumptions on the family of variational approximations to the posterior density.
This paper proposes a VB method for the MNP model that overcomes the two challenges. First, we construct an analytical expression of the unnormalized augmented posterior density by using the model specification as proposed by loaiza2021scalable. They transform the covariance parameters into a spherical coordinate system. The spherical transformation naturally imposes the trace restriction on the covariance matrix, and therefore the parameters are identified within the model specification. Second, we use the accurate variational approximation for models with multiple latent variables proposed in loaiza2021fast. The approximation for the latent utilities is the exact conditional posterior distribution for the latent utilities, and the approximation for the coefficients and the parameters in the spherical transformation is Gaussian. The combination of these two innovations results in an accurate VB approach that is substantially faster than MCMC.
Additionally, we demonstrate that our VB approach is scalable to data sets with a large number of observations. The VB optimization problem is solved with stochastic gradient ascent (SGA), where each iteration takes a draw from the conditional posterior distribution for the latent utilities. Sampling from the conditional posterior distribution of the latent utilities is computationally costly and hence takes the majority of the computation time in VB. Since SGA allows for subsampling, which means that in each iteration only a subsample of the latent utilities have to be generated, the computational complexity of the proposed methods can be further reduced.
We formulate our VB approach for the general multivariate multinomial probit (MVMNP) model. This means that we provide one common method for fast and accurate inference for a variety of different models, such as the MNP and multivariate probit (MVP) models, that are currently estimated with different identification strategies and different MCMC methods. For instance, zhang2006sampling and talhouk2012efficient fix the covariance matrix to be a correlation matrix in an MVP model. chib1998mcmc and mcculloch2000bayesian fix one element of the covariance matrix in MNP models, and zhang2008bayesian extend this to the MVMNP model. burgette2012trace introduce the trace restriction in the MNP model, and richard2012sparse fix the scale of a factor structure in the covariance matrix in the MVP model. Moreover, our method extends the spherical transformation on the covariance parameters to any MVMNP model, which allows for a factor structure in the covariance matrix that can substantially reduce the number of parameters to be estimated. Hence, our method can be applied to data sets with a large number of choices and a large number of choice alternatives.
Numerical experiments show that our VB method provides accurate parameter estimates and choice probabilities, while it only takes a fraction of the computational cost of MCMC. A numerical experiment with a small data set of 10,000 observations, in which MCMC is feasible, shows a minimal loss in predictive accuracy of VB relative to MCMC. An experiment with one million observations, in which MCMC is infeasible, shows that VB applied to a large data set can improve predictive accuracy relative to MCMC applied to only a subset of the observations.
We illustrate the practical relevance of our method with two empirical applications. First, VB produces similar results as MCMC around 10% of the computation time, in a small real data set of laundry detergent purchases. The second application considers a large-scale choice data set. Data on large choice sets with a large number of observations is widely available nowadays. We estimate the MNP model with the proposed VB method on more than a million pasta purchases in less than 1.5 hours, while MCMC takes more than 90 hours.
The outline of the remainder of this paper is as follows. Section (ref) discusses the model specification and Section (ref) develops our VB method. Section (ref) conducts numerical experiments to evaluate its accuracy and computational costs, and Section (ref) applies the proposed methods to real choice data sets. Section (ref) concludes.
We observe $K$ multinomial choices, where each choice $k=1,\dots,K$ has $J_k+1$ choice alternatives, for individual $i=1,\dots,N$. Let $\boldsymbol{Y}_i=(Y_{i1},\dots,Y_{iK})^\top$ denote the $K$-dimensional random variable describing the joint set of choices for individual $i$, where $Y_{ik}=j$ if individual $i$ chooses $j=0,1,\dots,J_k$ for the $k$-th choice. The number of potential outcomes of $\boldsymbol{Y}_i$ is $\prod_{k=1}^K(J_k+1)$.
Assume that for the $k$-th choice there is a $J_k$-dimensional vector $\boldsymbol{Z}_{ik}=(Z_{ik1},\dots,Z_{ikJ_k})^\top$ of continuous random variables representing the latent utilities for the choice alternatives, and which excludes the base category latent utility $Z_{ik0}$. The multinomial outcome $Y_{ik}$ is determined by the maximum value of $\boldsymbol{Z}_{ik}$ as follows:
where $\max(\boldsymbol{Z}_{ik})$ is the largest element of $\boldsymbol{Z}_{ik}$. The latent utilities corresponding to the choice alternatives in choice $k$ are modeled as
where $X_{ik}$ is a $J_k \times r_k$ regressor matrix, $\boldsymbol{\beta}_{k}$ is an $r_k$-dimensional vector of coefficients and $\boldsymbol{\varepsilon}_{ik}=(\varepsilon_{ik1},\dots,\varepsilon_{ikJ_k})^\top$ is a $J_k$-dimensional disturbance vector with mean zero.
The regressor matrix $X_{ik}$ typically includes $J_k$ choice alternative-specific intercepts, an $n_d$-dimensional vector $\boldsymbol{x}_{i}^d$ of individual-specific characteristics, and a $(J_k+1) \times n_a$ matrix $X_{ik}^a$ of $n_a$ alternative-specific covariates, such that $r_k = J_k+J_kn_d+n_a$ and
with transformation matrix $T_k=[-\boldsymbol{\iota}_{J_k} \quad I_{J_k}]$, where $\boldsymbol{\iota}_{J_k}$ denotes the $J_k$-dimensional vector of ones and $I_{J_k}$ the $J_k\times J_k$ identity matrix.
The model specified in (ref) and (ref) only considers the $J_k$ utilities in $\boldsymbol{Z}_{ik}$. Therefore, the covariates in the regressor matrix $X_{ik}$ in (ref) are transformed to match the dimensions of $\boldsymbol{Z}_{ik}$. The transformation matrix $T_k$ subtracts the covariates corresponding to choice category $j=0$ in choice $k$ from the covariates corresponding to the remaining choice alternatives in choice $k$. This model specification addresses the first parameter identification problem in the MVMNP model: additive redundancy arises if a unique utility is specified for each choice alternative, as discussed in bunch1991estimability. The second identification problem is caused by multiplicative redundancy: multiplying both sides of (ref) by a positive scalar does not change $Y_{ik}$ in (ref). We address this identification problem in Section (ref).
\sloppy By specifying a joint distribution for $\boldsymbol{\varepsilon}_{i} = (\boldsymbol{\varepsilon}_{i1}^\top,\dots,\boldsymbol{\varepsilon}_{iK}^\top)^\top$, the MVMNP model can allow for within-choice and between-choice correlation in the latent utilities. The $K$ models implied by (ref) can be stacked to obtain the multivariate utility model
where $\boldsymbol{Z}_i=(\boldsymbol{Z}_{i1}^\top,\dots,\boldsymbol{Z}_{iK}^\top)^\top$ is a $J-$dimensional vector with $J=\sum_{k=1}^K J_k$, $X_i=\text{blockdiag}(X_{i1},\dots,X_{iK})$ is a $J\times r$ block diagonal regressor matrix with $r=\sum_{k=1}^K r_k$, and $\boldsymbol{\beta}=(\boldsymbol{\beta}_1^\top,\dots,\boldsymbol{\beta}_K^\top)^\top$ is an $r$-dimensional vector of coefficients. The $J \times J$ covariance matrix of the disturbance vector $\boldsymbol{\varepsilon}_{i}$ can be represented as
where $\Sigma_{kl}=\Sigma_{lk}$. The $J_k \times J_l$ covariance matrix $\Sigma_{kl} = \text{cov}(\boldsymbol{\varepsilon}_{ik},\boldsymbol{\varepsilon}_{il})$ captures the correlation across the utilities within-choice $k$ if $k=l$, and the correlation across the utilities between choices $k$ and $l$ if $k \neq l$, with $k=1,\dots,K$ and $l=1,\dots,K$.
The MVMNP model easily simplifies to other commonly used choice models. First, for $K=1$ the MVMNP boils down to a multinomial probit (MNP) model, with $J+1$ choice alternatives corresponding to potentially correlated latent utilities. When $K>1$ and the elements of $\Sigma_{lk}$ equal zero for all $k\neq l$ in (ref), we have $K$ independent MNP models. Second, with $J_k=1$ for $k = 1,\dots,K$, we have a multivariate probit (MVP) model, with $K$ potentially correlated binary choices. When $K=1$ and $J=1$, the model boils down to a simple binary probit model.
For the remainder of the paper we stack the regressor matrices for all individuals in the $(NJ)\times r$ matrix $X = \left[X_1^\top |\dots| X_N^\top\right]^\top$ , the random choice vectors in $\boldsymbol{Y} = \left(\boldsymbol{Y}_1^\top,\dots,\boldsymbol{Y}_N^\top\right)^\top$, and the random latent utility vectors in $\boldsymbol{Z} = \left(\boldsymbol{Z}_1^\top,\dots,\boldsymbol{Z}_N^\top\right)^\top$.
The total number of unique parameters in the covariance matrix $\Sigma$ equals $J(J+1)/2$, which grows quadratically in the number of choices and the number of choice alternatives. Since these parameters have to be estimated from a single multinomial variable $Y_i$, accurate parameter estimation is challenging if either $J$ or $K$ is large, or both, even when a relatively large number of observations $N$ is available.
To reduce the dimension of the parameter space, we specify a factor structure for $\Sigma$. Define the $J\times p$ matrix $B$ with $p \leq J$ and the $J\times J$ diagonal matrix $D$. We model $\Sigma$ as
The total number of parameters in $B$ and $D$ is $n=J(p+1)$. This implies that for a given value of $p$, the number of parameters grows linearly with $J$, instead of quadratically.
To understand the implications of this factor structure on the within- and between-choice correlations, the matrices $B$ and $D$ are partitioned as $B = \left[B_1^\top |\dots| B_K^\top\right]^\top$ and $D = \text{diag}\left(\boldsymbol{d}_1^\top,\dots,\boldsymbol{d}_K^\top\right)$, where $B_k$ is a $J_k\times p$ matrix and $\boldsymbol{d}_k$ a $J_k$-dimensional vector corresponding to choice $k$. The within-choice covariance matrix $\Sigma_{kk}$ is expressed as
which shows that the matrices $\{\Sigma_{kk}\}_{k=1}^K$, and hence the within-choice correlations for each choice, are characterised by disjoint sets of model parameters. The between-choice covariance matrix $\Sigma_{kl}$ with $k \neq l$ is given by
which is only a function of the matrices $B_k$ and $B_l$ corresponding to the within-choice correlations in choices $k$ and $l$. In sum, the matrices $\{B_{k}\}_{k=1}^K$ determine both the within- and between-choice covariances in the latent utilities.
A factor structure with a small number of factors $p$ is correctly specified if the variability in the $J$ latent utilities in the data generating process can be captured by the $p$ latent factors. This is the case if, for instance, choice behavior is driven by a small number of unobserved features of the choice alternatives. Another example is the MVMNP model with uncorrelated choices, which can be estimated as separate univariate MNP models in which the factors only need to capture within-choice correlations. Issues of model misspecification may arise when the total number of unknown underlying factors is greater than $p$. This may be the case if choice behavior is driven by a large number of underlying choice features, or if the correlation pattern cannot be captured by a few factors. For instance, if choice alternatives are spatially related, the covariance matrix may have a banded pattern.
Even when the covariance matrix is misspecified, a low dimensional factor structure could potentially be favoured to trade-off flexibility for parsimony in the model, provided that key outputs from the model, such as predictive performance, remain accurate. The factor analysis literature has proposed several approaches to selecting $p$ optimally, including the use of information criteria, marginal likelihoods and cross-validation fruhwirth2018sparse, and these approaches may also be applied to the MVMNP model. We set $p=K$ in this paper, which allows for different covariance structures within each choice while reducing the number of covariance parameters to be estimated, and hence reducing parameter uncertainty and computation time.
As discussed in Section (ref), the scale of the latent utilities is unidentified. To identify the parameters, we extend the approach of Loaiza-Maya and Nibbering (2021) from an MNP model to the MVMNP model. For each choice $k$, we fix the scale using $\text{trace}(\Sigma_{kk})=J_k$, by transforming the elements of $B_k$ and $\boldsymbol{d}_k$ into a spherical coordinate system.
This parameter identification strategy has three advantages. First, in contrast to alternative identification restrictions, the trace restriction identifies the model parameters without fixing specific elements in the covariance matrix $\Sigma$. burgette2012trace show that Bayesian estimation in the MNP model is sensitive to which elements in the covariance matrix are fixed. Second, the spherical transformation on the factor covariance structure simplifies parameter estimation, as it naturally satisfies the trace restriction. Instead of performing inference on a parameter space with a joint restriction on all elements in each $\Sigma_{kk}$, we perform inference on the angle parameter space for which no joint parameter restriction is required. Third, since the spherical transformation naturally imposes the trace restriction, our approach does not require rescaling of the covariance matrix. Therefore, an analytical expression for the gradient of the likelihood function is available, which allows us to apply VB. Even with $p=J$, in which case there is no dimension reduction, writing $\Sigma$ as in (ref) has the benefit that it can be transformed by a spherical transformation and VB can be applied.
The spherical transformation is applied to the vector $\boldsymbol{\psi}_k$, which is constructed from the elements of $B_k$ and $\boldsymbol{d}_k$ as
where $n_k=J_k(p+1)$, $\text{vec}(\cdot)$ denotes the vectorization operator, and $\text{trace}(\Sigma_{kk})=\sum_{l=1}^{n_k}\psi_{kl}^2$. We transform $\boldsymbol{\psi}_{k}$ into a spherical coordinate system that is defined by a radius, which we set to $\sqrt{J_k}$, and an $(n_k-1)$-dimensional vector of angles $\boldsymbol{\kappa}_k = \left(\kappa_{k1},\dots,\kappa_{k,n_k-1}\right)^\top$.
The spherical transformation reparameterises $\boldsymbol{\psi}_k$ in terms of $\boldsymbol{\kappa}_k$ as
where $\kappa_{kl} \in [0,\pi)$ for $l<n_k-J_k+1$. The remaining angle bounds, $\kappa_{kl} \in [0,\frac{\pi}{2})$ for $n_k-J_k+1\le l \le n_k-1$, ensure that the elements of $\boldsymbol{d}_k$ are strictly greater than zero. This ensures that the map from $\boldsymbol{d}_k$ to $\Sigma_{kk}$ in (ref) is bijective.
The transformation in (ref) satisfies $\sum_{l=1}^{n_k} \psi_{kl}({\boldsymbol{\kappa}_k})^2 = J_k$ for any value of $\boldsymbol{\kappa}_k$. This reparametrization is applied to all $\Sigma_{kk}$, which results in a covariance matrix $\Sigma$ that is characterised by the $n$-dimensional vector $\boldsymbol{\kappa} = \left(\boldsymbol{\kappa}_1^\top,\dots,\boldsymbol{\kappa}_K^\top\right)^\top$, where $n=\sum_{k=1}^Kn_k$. The vector $\boldsymbol{\kappa}$ imposes $K$ trace restrictions simultaneously: one for each choice.
The inverse function of the spherical transformation in (ref) is
where $\kappa_{kl}=0$ if $\psi_{kl}>0$ and $\psi_{k,l+1}=\dots=\psi_{kn_k}=0$, and $\kappa_{kl}=\pi$ if $\psi_{kl}<0$ and $\psi_{k,l+1}=\dots=\psi_{kn_k}=0$.
This section develops a Bayesian method for estimating the parameters in the $m$-dimensional vector $\boldsymbol{\theta} = \left(\boldsymbol{\beta}^\top,\boldsymbol{\kappa}^\top\right)^\top$, with $m = r+n$. We conduct inference of the augmented posterior density
where we use lower case letters to denote realised values of the corresponding random vectors. For instance, $\boldsymbol{y}$ is the realised vector of $\boldsymbol{Y}$.
The augmented likelihood function is given by
where $\phi_{J}\left(\boldsymbol{z}_i;X_i\boldsymbol{\beta},\Sigma(\boldsymbol{\kappa})\right)$ denotes the $J$-variate normal density with mean $X_i\boldsymbol{\beta}$ and covariance matrix $\Sigma(\boldsymbol{\kappa})$, with $\Sigma(\boldsymbol{\kappa})$ the covariance matrix constructed from the vector of angles $\boldsymbol{\kappa}$, $p(\boldsymbol{y}_i|\boldsymbol{z}_i) = \prod_{k=1}^{K}p(y_{ik}|\boldsymbol{z}_{ik})$ and
where $I[A]$ is an indicator function that equals one if $A$ is true and zero otherwise.
We set the prior density as $p(\boldsymbol{\theta}) = p(\boldsymbol\beta)\prod_{k=1}^{K}p(\boldsymbol\kappa_k)$, with $\boldsymbol{\beta}\sim N(\boldsymbol{0}_r, \frac{1}{10}I_r)$ and $p(\boldsymbol\kappa_k)$ specified in online appendix A, with an implied prior mean for $\Sigma$ that equals the equicorrelated covariance matrix $\frac{1}{2}(I_J + \iota_J \iota_J^\top)$.
The posterior density in (ref) can be computed using MCMC sampling. For each individual, the latent utility of each choice alternative is sampled conditional on all the other choice alternatives from a truncated normal. This process induces a sequence of latent utility draws that is highly auto-correlated. Therefore MCMC requires a large number of iterations such that convergence is achieved. Since each iteration involves $N \times J$ draws from a truncated normal, MCMC is computationally costly. Online appendix B describes the MCMC sampling scheme for the MVMNP model.
The computational costs of an MCMC sampling scheme increase in the number of choice alternatives $J_k$ in each choice $k$, the number of choices $K$, and the number of observations $N$. As a result, when the total number of choice alternatives $J$ is large, MCMC is considered to be computationally practical as long as the number of observations is small. There are two empirical settings where this is the case. First, univariate choice sets ($K=1$) that have a large number of choice alternatives. Second, applications that consider multiple choices and for which the overall number of choice alternatives $J$ is large. However, it is precisely in these type of settings where having a large number of observations is key for accurate estimation of the high-dimensional covariance matrix of the latent utilities.
To circumvent the computational challenges of MCMC, we utilize variational Bayes. VB approximates the posterior density in (ref) by a parametric density $q_{\widehat{\lambda}}\left(\boldsymbol{\theta},\boldsymbol{z}\right)\in\mathcal{Q}$ from the class of density functions $\mathcal{Q} = \{q_\lambda\left(\boldsymbol{\theta},\boldsymbol{z}\right): \boldsymbol{\lambda}\in\Lambda\}$, where $q_\lambda\left(\boldsymbol{\theta},\boldsymbol{z}\right)$ is indexed by the variational parameter vector $\boldsymbol{\lambda}\in\Lambda$. The optimal variational parameter $\widehat{\boldsymbol{\lambda}}$ is obtained by maximizing the evidence lower bound (ELBO) function $\mathcal{L}\left(\boldsymbol{\lambda}\right)=E_{q_{\lambda}}\left[\log g(\boldsymbol{\theta},\boldsymbol{z})-\log q_\lambda(\boldsymbol{\theta},\boldsymbol{z})\right]$:
where $g(\boldsymbol{\theta},\boldsymbol{z}) = p(\boldsymbol{y}|\boldsymbol{z})p(\boldsymbol{z}|X,\boldsymbol{\theta})p(\boldsymbol{\theta})$ is the unnormalized posterior density. ormerod2010explaining show that the optimization problem in (ref) is equivalent, but computationally more efficient, to minimizing the Kullback-Leibler (KL) divergence between $q_{{\lambda}}\left(\boldsymbol{\theta},\boldsymbol{z}\right)$ and the exact posterior density $p(\boldsymbol{\theta},\boldsymbol{z}|\boldsymbol{y},X)$.
While MCMC generates draws from the exact posterior distribution, VB can only construct an approximation to it. However, VB has three main advantages over MCMC for the MVMNP model. First, MCMC may show high auto-correlation in its chain for this model, leading to substantial computational costs. VB relies on optimization rather than sampling, and therefore reduces the computation time. Second, VB requires much less storage memory as the output from VB is the calibrated parameter vector of the approximation, rather than a large number of parameter draws. Third, VB can readily incorporate subsampling of the latent utilities in the optimization routine, which can further reduce the computational burden.
Key to the implementation of VB is the choice of the variational family $\mathcal{Q}$. We set $q_\lambda\left(\boldsymbol{\theta},\boldsymbol{z}\right) = p(\boldsymbol{z}|\boldsymbol{\theta},\boldsymbol{y},X)q_\lambda(\boldsymbol{\theta})$, where $p(\boldsymbol{z}|\boldsymbol{\theta},\boldsymbol{y},X)$ is the conditional posterior of the latent utilities defined in online appendix B, and we define $q_\lambda(\boldsymbol{\theta})$ below. loaiza2021fast show that the optimization problem in (ref) with this variational family is equivalent to the optimization problem that considers the KL divergence between the intractable posterior $p(\boldsymbol{\theta}|\boldsymbol{y},X)$ and $q_\lambda(\boldsymbol{\theta})$. Alternative specifications for the variational family may result in approximating errors for the latent utilities, which can provide inconsistent estimates, as shown in westling2019beyond.
For the choice of $q_\lambda(\boldsymbol{\theta})$, we follow ong2018gaussian and employ a Gaussian density with mean $\boldsymbol{\mu}$ and covariance matrix $\Omega = CC^\top+E^2$, where $C$ is a matrix of dimension $m\times s$ for $s<m$, $E=\text{diag}(\boldsymbol{e})$ and $\boldsymbol{e}$ an $m$-dimensional vector. The variational parameter vector for this approximating class is $\boldsymbol{\lambda} = \left(\boldsymbol{\mu}^\top,\text{vech}(C)^\top,\boldsymbol{e}^\top\right)^\top$, where the operator vech denotes the half vectorization of a rectangular matrix such that $\text{vech}(C)= \left(C_{1:m,1}^\top,\dots,C_{s:m,s}^\top\right)^\top$ with $C_{j:m,j} = \left(C_{jj},\dots,C_{mj}\right)^\top$ for $j = 1,\dots,s$.
We solve the optimization problem in (ref) using SGA methods. SGA calibrates the variational parameter by iterating over
until convergence is achieved. The vector $\boldsymbol{\rho}^{[j]}$ contains the so called “learning parameters”, which we set according to the ADADELTA approach in zeiler2012adadelta. The vector $\widehat{\nabla_\lambda \mathcal{L}\left(\boldsymbol{\lambda}^{[j]}\right)}$ is an unbiased estimate of the gradient of the ELBO evaluated at $\boldsymbol{\lambda}^{[j]}$.
We construct $\widehat{\nabla_\lambda \mathcal{L}\left(\boldsymbol{\lambda}\right)}$ using the following expression of the gradient
where $\boldsymbol{\theta}(\boldsymbol{\zeta},\boldsymbol{\lambda}) = \boldsymbol{\mu}+C\boldsymbol{w}+\boldsymbol{e}\circ\boldsymbol{\epsilon}$, ${\boldsymbol{w}}\sim N(\boldsymbol{0}_{s},I_{s})$, $\boldsymbol{\epsilon}\sim N(\boldsymbol{0}_{m},I_{m})$, and $\boldsymbol{\zeta} = \left({\boldsymbol{w}}^\top,{\boldsymbol{\epsilon}}^\top\right)^\top$. This expression is derived in loaiza2021fast using the “re-parametrization trick” in kingma2013auto. We derive $\nabla_{\theta}\log g\left[\boldsymbol{\theta}(\boldsymbol{\zeta},\boldsymbol{\lambda}),\boldsymbol{z}\right]$ for the MVMNP model in online appendix C, and $\frac{\partial\boldsymbol{\theta}(\boldsymbol{\zeta},\boldsymbol{\lambda}) }{\partial\boldsymbol{\lambda}}$ and $\nabla_\theta \log q_\lambda\left[\boldsymbol{\theta}(\boldsymbol{\zeta},\boldsymbol{\lambda}) \right]$ are provided in ong2018gaussian.
An unbiased estimate of $\nabla_\lambda \mathcal{L}\left(\boldsymbol{\lambda}\right)$ is constructed by using a sample estimate of the expectation. At each SGA iteration $[j]$, we calculate a sample estimate of (ref) based on only one draw for both $\boldsymbol\zeta$ and $\boldsymbol{z}$: $\boldsymbol\zeta^{[j]}\sim N(\boldsymbol{0}_{s+m},I_{s+m})$ and $\boldsymbol{z}^{[j]}\sim p(\boldsymbol{z}|\boldsymbol{\theta}(\boldsymbol{\zeta}^{[j]},\boldsymbol{\lambda}^{[j]}),\boldsymbol{y},X)$. Since sampling directly from $p(\boldsymbol{z}|\boldsymbol{\theta}(\boldsymbol{\zeta},\boldsymbol{\lambda}),\boldsymbol{y},X)$ is infeasible, we generate the latent utility vector draw $\boldsymbol{z}^{[j]}$ via $G$ Gibbs steps of the truncated normal algorithm proposed in mcculloch1994exact. The Gibbs sampling algorithm is started at the last iterate value $\boldsymbol{z}^{[j-1]}$. Since these Gibbs draws are highly correlated, a larger value for $G$ increases the accuracy of the estimate for the gradient. However, each Gibbs step includes $N\times J$ draws from a truncated normal distribution, which is computationally costly. We find that $G=10$ balances well accuracy and computational speed .
The objective of VB is to compute the optimal variational parameter vector, so it suffices to run enough SGA iterations until convergence is reached for all the elements of $\boldsymbol{\lambda}^{[j]}$. SGA generally requires a small number of iterations, making it much faster than MCMC. However, the majority of the computation time is still spent on the generation of the latent utilities. Since SGA allows for subsampling of the observations, the computational burden of the latent utilities can be substantially reduced in VB.
Instead of sampling the latent utilities $\boldsymbol{z}_i$ at each iteration for all individuals, SGA can estimate the gradient unbiasedly using only a subsample of the latent utilities. The ELBO gradient can be rewritten in terms of the variable $A\subset\{1,\dots,N\}\sim f(A)$, where a draw from $f(A)$ is a random subsample of indexes without replacement. Define the subsample of latent utilities as $\boldsymbol{z}_A = \{\boldsymbol{z}_i\}_{i\in A}$. Since $E_A\left[{\nabla_{\theta}\log g}(\boldsymbol{\theta},\boldsymbol{z}_A)\right]={\nabla_{\theta}\log g}(\boldsymbol{\theta},\boldsymbol{z})$, it holds that
Online appendix C provides the expression for ${\nabla_{\theta}\log g}(\boldsymbol{\theta},\boldsymbol{z}_A)$ required to compute the subsampling gradient estimate for the MNP model. An unbiased estimate of (ref) is constructed using a sample estimate of the expectation, using the draws $A^{[j]}\sim f(A)$, $\boldsymbol{\zeta}^{[j]}\sim N(\boldsymbol{0}_{s+m},I_{s+m})$, and $\boldsymbol{z}_{A^{[j]}}^{[j]}\sim p(\boldsymbol{z}_{A^{[j]}}|\boldsymbol{\theta}(\boldsymbol{\zeta}^{[j]},\boldsymbol{\lambda}^{[j]}),\boldsymbol{y},X)$.
The VB predictive probability mass function for $\boldsymbol{Y}_{i}$ is given by
where $X_i$ denotes the attributes of the observation $i$ to be predicted. We construct an estimate $\hat{p}_{\hat{\lambda}}(\boldsymbol{Y}_i|X_i)$ for (ref) as the empirical probability mass implied by the draws $\{\boldsymbol{y}_i^{[m]}\}_{m=1}^M$, obtained by drawing from $\boldsymbol{\theta}^{[m]}\sim q_{\hat{\lambda}}(\boldsymbol{\theta})$, $\boldsymbol{z}_i^{[m]}\sim p(\boldsymbol{z}_i|\boldsymbol{\theta}^{[m]},X_i)$ and $\boldsymbol{y}_i^{[m]}\sim p(\boldsymbol{Y}_i|\boldsymbol{z}_i^{[m]})$.
To evaluate predictive performance, we employ the logarithmic score (log-score), which is a probabilistic measure of predictive accuracy. The choice-specific log-score is given as
and the total model fit can be assessed by the average log-scores across choices.
The point forecast $\hat{\boldsymbol{Y}}_i$ for $\boldsymbol{Y}_i$ is constructed as the mode of $\hat{p}_{\hat{\lambda}}(\boldsymbol{Y}_i|X_i)$. The point prediction accuracy can be measured in terms of the hit-rate given as
For both the hit-rate and the log-score large values are preferred.
This section presents two numerical experiments to assess the accuracy and the computational costs of the proposed VB approach. The first experiment compares VB to MCMC in a moderately sized data set in which MCMC is computationally feasible. The second experiment is on a large dataset for which VB estimation is feasible, but MCMC is not.
We generate a data set from the data generating process specified in (ref) and (ref) with $K=2$ and $J_1=J_2=10$. The elements of the matrices $X_{i1}^a$ and $X_{i2}^a$ are independently generated from normal distributions with corresponding mean $\mu=0$ and variance $\sigma^2=1$. These elements can be interpreted as the logarithm of the prices of the choice categories. We do not include individual-specific characteristics $\boldsymbol x_{i}^d$.
The true parameter vector $\boldsymbol\beta_0$ consists of $J$ intercepts drawn independently from uniform distributions $U(-0.5,0)$, and the coefficients for $X_{i1}^a$ and $X_{i2}^a$ are fixed at -0.3 and -0.6, respectively. The true covariance matrix $\Sigma_0$ is set as a draw from the inverse Wishart distribution with equicorrelated scale matrix $\frac{1}{2}(I_J + \iota_J \iota_J^\top)$ and degrees of freedom $J+3$.
We apply our VB method to two generated data sets, one with $N=10,000$ and the second one with $N=1,000,000$. For both settings we generate an additional 10,000 observations for out-of-sample evaluation. VB with subsampling is denoted by VB($\frac{M}{N}100\%$), where $\frac{M}{N}100\%$ denotes the percentage of the total estimation sample $N$ used in each VB iteration step. For the purpose of comparison, we also estimate an MVMNP model with the covariance matrix fixed at the identity matrix (VB-I), as described in online appendix D.
The VB methods estimate the model with 5000 iterations of SGA, with 10 Gibbs sampling steps in each SGA iteration. We take 10,000 draws from the variational posterior and predictive distribution to construct the results. The results from MCMC sampling are based on 200,000 iterations, from which the first 100,000 are discarded and we use a thinning value of 10. This results in 10,000 draws from the posterior and predictive distribution. All methods use the prior specification as discussed in Section (ref). The methods are implemented in a HP Z240 SFF Workstation with an Intel i7-7700 CPU \@ 3.6GHz.
First, we assess convergence of SGA in our VB methods. The ELBO in (ref), which is typically used as convergence measure in VB, is not available in closed-form. The hit-rate defined in (ref), with $\hat{\boldsymbol{Y}}_i$ as the mode of $\hat{p}_{\lambda^{[j]}}(\boldsymbol{Y}_i|X_i)$ in iteration $[j]$, can be used instead. To reduce the computational costs, we evaluate the conditional hit-rate for a fixed random subsample of 500 observations, using 200 draws from $\hat{p}_{\lambda^{[j]}}(\boldsymbol{Y}_i|X_i)$, in each tenth iteration.
Figure (ref) shows the conditional hit-rate for choice 1 and 2, in VB and VB(1%) by a yellow and black line, respectively. The figure indicates that 5000 iterations are sufficient for convergence. The conditional hit-rates remain wiggly because they are constructed using an estimate for $\hat{p}_{\lambda^{[j]}}(\boldsymbol{Y}_i|X_i)$ that is based on a $\boldsymbol\lambda^{[j]}$ that is updated using an estimate for the gradient. Therefore, the final estimate for the variational parameter $\boldsymbol\lambda$ is constructed as the average over the $\boldsymbol{\lambda}^{[j]}$ in the final 100 iterations. We find that VB converges in less iterations than VB(1%). However, VB(1%) has a faster convergence, as an average iteration takes 0.018 seconds, compared to 0.562 seconds per iteration in VB.
The computation time of VB(1%), VB(10%), and VB is 0.02, 0.11, and 0.77 hours respectively. Since MCMC takes 6.1 hours, this means that VB uses less than 14% of the time required for MCMC, and the time can be further decreased with subsampling.
Second, we assess the accuracy of the posterior distribution for the parameters. Panels (a) to (c) in Figure (ref) show the VB against the MCMC posterior means. The closer the circles lie to the 45 degree line, the closer the VB posterior means are to those of MCMC. The VB estimates are scattered around the 45 degree line, which indicates that they are close to the exact posterior means. Panels (d) to (f) in Figure (ref) show that VB tends to underestimate the posterior standard deviation, which is a well-documented property of variational approximations blei2017variational,yu2021assessment.
Online appendix E compares the posterior means and standard deviations of VB(10%) and VB(1%) to those of MCMC. The posterior means of VB with subsampling are still scattered around the 45 degree line, but the deviations from this line slightly increase with smaller subsamples. This suggests that the reduction in computational costs induced by subsampling comes at the cost of a small loss in accuracy. VB with subsampling does not seem to underestimate the posterior standard deviation.
Although we have demonstrated the accuracy of VB at estimating the posterior of the parameters of the model, these parameter estimates themselves are hard to interpret and as such are not the key output from the model. Instead, the posterior choice probabilities are the quantity of interest in most empirical applications. Figure (ref) shows the choice probability of one of the categories for choice 1 (Panel (a)) and for choice 2 (Panel (b)) as a function of their price, with the prices of the other categories fixed at their mean. The solid yellow lines correspond to the posterior probabilities of MCMC and the dashed black lines to VB. These lines are almost identical, and we find the same result for the other categories, and when comparing VB(1%) to MCMC. Hence we conclude that VB and VB with subsampling accurately estimate the posterior choice probabilities.
Third, we examine the predictive accuracy of our VB approach. The log-score can be used to assess the impact of subsampling in VB on the estimated model fit. Figure (ref) shows the estimation time and in- and out-of-sample log-score across choices corresponding to VB with different subsampling sizes and MCMC. The log-score increases in the size of the subsample, with a big increase corresponding to subsampling with 1% to 10%. After 10%, an increase in the subsample results in small gains in the log-scores.
Figure (ref) shows that the impact of subsampling in VB on the estimated model fit is small compared to the gains in computational efficiency. For instance, VB(10%) has an in- and out-of-sample log-score close to VB and MCMC. The differences between these three methods are small relative to the difference between MCMC and the log-score of the oracle: the log-score computed using the true parameter values. This is a striking result, as VB(10%) is estimated in 7 minutes, VB takes 47 minutes, and MCMC more than 6 hours.
Table (ref) provides a more detailed overview on the estimated model fit. The upper panel shows the in- and out-of-sample log-score and hit-rate for the first choice, and the lower panel for the second choice. These measures are higher in the second choice, indicating a stronger signal. MCMC performs better than VB on the log-scores, but does not always outperform the hit-rates of the different VB approaches. The average log-score of VB across choices declines with smaller subsamples.
VB in the MVMNP model can also be compared to VB in a choice model with an identity covariance matrix, VB-I. We find that VB-I is outperformed on all measures by VB. In general, all models perform better than the naive method on all metrics: the naive method sets the forecast equal the most frequently observed category in the data. The oracle, that uses the true parameter values to construct a forecast, corresponds to the highest log-scores, but does not always attain the highest hit-rates.
To study the robustness of the predictive performance of the model to the choice of the number of factors $p$, we repeat the predictive exercise for an increasing number of factors $p = 0,\dots,10$. Because the true DGP is generated from a full covariance matrix, smaller values of $p$ indicate a higher level of model misspecification. We find that the increase in predictive performance is modest beyond $p=2$ factors, while the estimation time grows linearly with $p$. Details are deferred to online appendix E.
The second experiment illustrates that VB with subsampling makes the estimation of the MVMNP model computationally feasible on big data sets. Figure (ref) shows that MCMC takes more than six hours with 10,000 observations. In practice, choice data sets may have much larger samples as we illustrate in Section (ref). With one million observations, MCMC takes around eighteen days. These computational costs make MCMC impractical in many choice applications. VB takes almost three days, which is still a substantial computational cost. On the other hand, VB(1%) takes less than 50 minutes to be implemented.
An equally time efficient approach that could be implemented instead of subsampling VB, would be to consider VB on a random subsample of the data. This approach does not make use of the complete data set, and as such can be suboptimal in terms of predictive accuracy. To show this, we apply VB to a subsample of 10,000 observations, which requires approximately the same computation time as VB(1%) on a million observations. Table (ref) shows the log-scores and hit-rates of the two approaches on the same out-of-sample observations, and the average log-score of VB(1%) is indeed higher than that of VB.
To illustrate our VB method with real data, we fit a multinomial probit model to two consumer choice data sets with different dimensions. First, Section (ref) employs a commonly used traditional data set on laundry detergent brand purchases with a few thousand observations. Second, Section (ref) uses a modern data set on pasta brand purchases with more than one million observations. We discuss the posterior choice probabilities and the predictive performance of VB and MCMC, and defer the results on the posterior parameter distributions to online appendix F. The implementation and prior settings of the proposed methods are discussed in Section (ref). We randomly allocate 80% of the observations for estimation of the model, and the remaining 20% are employed for out-of-sample evaluation.
This section uses a small choice data set to illustrate on real data that VB is several times faster than MCMC, yet produces similar choice probabilities and predictive accuracy. The data contains 2657 purchases of six brands of laundry detergents and the log price per ounce of each brand. The data set is described in detail by chintagunta1998empirical and available in imai2005mnp. We follow imai2005bayesian, burgette2019symmetric, and loaiza2021scalable by fitting multinomial probit models with an intercept and the log price for each brand.
We find that the posterior purchase probabilities of VB and MCMC are similar. Panel (a) in Figure (ref) shows the probability of buying the most popular brand as a function of its price, with the prices of the other brands fixed at their mean. The solid yellow line corresponds to the posterior probabilities of MCMC, the dashed black line to VB, and the dotted red line to VB(1%). VB produces posterior probabilities that are almost identical to MCMC, and VB(1%) only shows small differences compared to MCMC. We also find negligible differences between VB and MCMC, and small differences between VB(1%) and MCMC, for the purchase probabilities for the other five brands.
VB also attains similar predictive accuracy to MCMC. Table (ref) shows that the in- and out-of-sample log-scores of VB and MCMC are almost identical. The log-scores slowly decrease in the order of subsampling in VB. The hit rates do not seem to be very sensitive to the approximations by VB, or in the VB with subsampling. Both the log-scores and hit-rates of VB(10%) are relatively close to MCMC, especially when we consider the difference in these metrics between MCMC and the naive forecasting method, in which the forecast equals the most frequently observed category in the data. The fact that VB with an identity covariance matrix results in lower log-scores than VB with a full covariance matrix, indicates that the correlations across the latent utilities matter in this application.
Although the differences in parameter estimates and predictive accuracy are minimal, VB is more than eight times faster than MCMC. The final row of Table (ref) shows that MCMC takes 720 seconds while VB only takes 81 seconds. This computation time can be further reduced by subsampling. VB(10%) takes 38 seconds, as it required 10,000 rather than 5,000 SGA iterations to converge, and VB(1%) only 16 seconds.
Due to its low computational costs, VB is well suited to estimate multiple prior specifications to study the robustness of the posterior results to the choice of base category. For instance, burgette2021symmetric show that the posterior choice probabilities of Bayesian MNP models can depend on the base category specification. Panel (b) in Figure (ref) shows the probability of buying the least popular detergent brand `All' estimated with VB. The dotted red line uses brand `All' as the base category, and the solid yellow lines correspond to the five other base category specifications. The lines are different, and `All' as base category results in substantially higher purchase probabilities than with other base categories. The differences are less pronounced for choice probabilities of more popular brands. Setting the number of factors equal to the number of brands shows the same base category sensitivities, which is in line with the findings in burgette2021symmetric who specify a full covariance matrix. As a robust alternative for estimating the choice probabilities, the posterior probabilities corresponding to different base category specifications can be pooled. The dashed black line in Panel (b) in Figure (ref) shows the average purchase probability across all specifications.
This section shows that our approach can be scaled to real data with many observations. We use more than one million purchases from ten pasta brands in a consumer choice data set made available by Dunnhumby\footnote{https://www.dunnhumby.com/source-files/} as “Carbo-Loading: A Relational Database". From this data set, we select the purchases of pasta brands, excluding the private labels, that do not involve coupons. Since the brands with a small purchase volume are of less interest to a marketing manager, we focus on the 96.782% of purchases that corresponds to the ten top-selling pasta brands. The final sample contains 1,070,436 observations and the purchase frequencies vary from 6,280 to 316,018.
We consider the same MNP models as with the small data set on laundry detergent purchases, also including an intercept and the log price for each brand. We define the log price per ounce of each brand in the same way as, for instance, allenby1991quality and loaiza2021scalable. The Dunnhumby data set only contains the amount of dollar spent on a product at purchase dates. We impute the prices for brands that are not sold on a certain purchase date by taking the mean of the observed prices of a specific brand on the nearest date in the same week. In cases where there is no purchase record in the same week, we take the most recent observed price.
Figure (ref) shows the purchase probabilities as a function of its price for the most popular pasta brand in Panel (a) and the least popular pasta brand in Panel (b). The solid yellow line corresponds to the posterior probabilities of MCMC, the dashed black line to VB, and the dotted red line to VB(1%). These lines are almost identical for the most purchased pasta brand. For the least purchased pasta brand, the probabilities of MCMC are different from the probabilities of VB(1%) and VB, although the difference with the latter is small. In general, VB produces posterior probabilities that are accurate for the relatively large pasta brands, and lose some accuracy for brands with a small number of purchases.
Table (ref) reports the predictive performance measures. Based on the in-sample and out-of-sample log-scores and hit-rates, VB and VB(10%) show almost no loss in accuracy relative to MCMC. The table also indicates that the difference in predictive performance between MCMC and VB(1%) is small compared to the difference between MCMC and the naive method. We also find that, irrespective of the subsampling size, VB outperforms VB-I, which emphasises the importance of taking correlations into account.
Table (ref) also shows that VB makes the estimation of multinomial probit models feasible on real choice data sets with many observations. The final row shows that implementation of VB takes more than fifteen hours, which is a fraction of the 92 hours of MCMC. VB can even further reduce the computation time with subsampling, with VB(10%) only taking around 1.5 hours, and VB(1%) around eleven minutes.
Multinomial probit models are widely used for analyzing choice behavior. The main benefit of the model is the specification of the covariance matrix of the latent utilities. To accurately estimate the covariance parameters from a single categorical dependent variable, a large number of observations is required. Choice data sets with many observations are nowadays widely available. For instance, scanner data as used in the empirical application in this paper have records of millions of transactions. However, MCMC methods that are currently used for parameter estimation are computationally costly.
This paper proposes a variational Bayes method that employs the conditional posterior of the latent utilities as a part of the variational family. This allows for accurate approximations to the exact posterior. The method is faster than MCMC with moderately sized data sets, and is scalable to large-scale data in which MCMC estimation is infeasible.
Numerical experiments and an empirical application to a laundry detergent choice set demonstrate that our approach produces accurate approximating densities to the MCMC exact posterior densities, while only requiring a small fraction of the MCMC computation time. The computational cost for our approach can be further reduced by considering subsampling methods inside the stochastic gradient ascent algorithm, with small impact in its predictive accuracy relative to MCMC.
The new method improves the applicability of the multinomial probit model to modern choice data sets. We illustrate the potential of the new approach in large samples by applying it to a pasta choice data set that consists of more than one million observations.
\setcounter{page}{1}
\setcounter{figure}{0} \setcounter{table}{0} \setcounter{section}{0}
This online appendix has six parts: