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,778 characters · 25 sections · 46 citation commands
Scalable Bayesian estimation in the multinomial probit model
{\bf Keywords:} Multinomial probit model, Factor analysis, Parameter identification, Spherical coordinates \\ {\bf JEL Classification:} C11, C25, C35, C38
\thispagestyle{empty} \setcounter{page}{1}
The multinomial probit (MNP) is an important model for analysing choice behavior, because it allows the latent utilities of the choice alternatives to be correlated. These correlations capture general substitution patterns among choice alternatives, in contrast to the case with models that impose the independence of irrelevant alternatives property hausman1984specification.
However, current specifications of the multinomial probit model are not scalable to discrete choice problems with a large number of choice alternatives, as the number of parameters in the covariance matrix of the latent utitilities grows quadratically in the number of choice alternatives burgette2013multiple. This curse of dimensionality is exacerbated by the fact that, contrary to standard covariance matrix estimation settings, where multiple continuous variables are observed, all parameters in the covariance matrix have to be estimated from a single categorical variable.
Standard dimension reduction techniques, such as factor analysis, cannot straightforwardly be applied to the covariance matrix of the latent utitilities. Since the scale of the latent utilities is not identified bunch1991estimability, the identification of the model parameters requires a restriction on the covariance matrix. The main challenge is to reduce the number of parameters that characterize the covariance matrix, while imposing an identifying restriction.
This paper proposes a multinomial probit model specification that is scalable to modern choice data with many choice alternatives. Specifically, we employ a factor structure on the covariance matrix, where the number of parameters scales linearly, rather than quadratically, with the number of choice alternatives. The parameters are identified by imposing a trace restriction on the covariance matrix. To impose the trace restriction on the factor representation, we transform the covariance parameters to a spherical coordinate system of angles and a spherical radius. The radius is, by construction, equal to the square root of the trace of the covariance matrix. Therefore, the trace restriction is readily imposed by setting the squared radius equal to the number of choice alternatives.
To conduct Bayesian estimation, prior densities on the angles in the reparameterization must be selected. We elicit the priors on the angles from well understood prior assumptions popularly used in the Bayesian factor analysis literature. The process of elicitation can be performed using a fast algorithm provided in this paper. The computation of the posterior distribution involves a Markov Chain Monte Carlo (MCMC) sampler with Gibbs sampling steps for the coefficients and latent utilities, and a Metropolis-Hastings step for the angle parameters. A numerical experiment confirms that the MCMC sampler succeeds in accurately estimating the model parameters, in similar computation time as existing MNP specifications.
An application to real consumer choice data illustrates the empirical relevance of the scalable multinomial probit model. We construct consumer choice data with 50 alternatives as a modern counterpart of commonly used laundry detergent and margarine purchase data sets with only six alternatives. In these large choice sets, our approach produces better predictive performance than existing methods. The model is able to identify correlations across a large set of products. The results show that limiting the analysis to only a few products may, for instance, severely bias price elasticity estimates. The proposed model has similar performance to existing multinomial probit specifications when applied to the traditional laundry detergent and margarine choice data with six alternatives.
This paper makes three important contributions to the multinomial probit literature. First, the proposed approach addresses the scalability of the multinomial probit model directly. piatek2017multinomial specify a factor structure for the covariance matrix under the assumption that the factor loadings are known. This assumption might be unrealistic in many choice problems, especially when the number of alternatives is large. cripps2009parsimonious propose a covariance selection prior that permits elements of the inverse of the covariance matrix to be zero. This approach allows for a sparse representation of the model, but does not reduce the number of parameters to be estimated.
Second, this paper contributes to the literature on parameter identification in multinomial probit models. burgette2012trace show that fixing the trace of the covariance matrix should be preferred over fixing a diagonal element, as in mcculloch2000bayesian and imai2005bayesian. Our model reparametrization satisfies the trace restriction and also naturally imposes parsimony, which is not embedded in the marginal data augmentation approach used by burgette2012trace.
Third, this is the first study that applies the multinomial probit model to real choice data with a large number of choice alternatives. Empirical applications of multinomial probit models have been limited to only a few choice alternatives. For instance, imai2005bayesian consider six clothing detergent brands, mcculloch1994exact and burgette2012trace six margarine brands, piatek2017multinomial two education levels and three occupation categories, and cripps2009parsimonious five tests for cervical cancer. This limitation is of particular concern today, with the widespread availability of data on large choice sets.
The outline of the remainder of this paper is as follows. Section (ref) discusses the model specification and Section (ref) introduces a scalable covariance matrix specification. Section (ref) discusses prior specifications and the MCMC sampler. Section (ref) conducts a numerical experiment to evaluate estimation accuracy, and Section (ref) applies the proposed methods to real consumer choice data sets. Section (ref) concludes.
Let $Y_i$ be an observable unordered random categorical variable with support on the set $A_J = \{0,1,2,\dots,J\}$, with $J+1$ the number of choice alternatives, and $i=1,\dots,N$, with $N$ the number of individuals. Let $\tilde{Z}_i = (\tilde{z}_{i0},\dots,\tilde{z}_{iJ})^\top$ be a $(J+1)\times1$ vector of continuous random variables that can be interpreted as latent utilities, with
The latent utilities are modeled as
where $\tilde{X}_i$ is a $(J+1)\times \tilde{K}$ matrix of observed regressors, $\tilde{\beta}$ is a vector of coefficients, and $\tilde{\varepsilon}_i$ is an independent normally distributed disturbance vector with covariance matrix $\tilde{\Sigma}$. The regressor matrix typically includes an intercept, a $k_d$-dimensional vector $x_{i,d}$ of individual-specific characteristics, and a $(J+1)\times k_a$ matrix $x_{i,a}$ of $k_a$ alternative-specific covariates, such that
where $I_{J+1}$, denotes the identity matrix of dimensions $(J+1)\times(J+1)$.
The parameters $\tilde{\beta}$ and $\tilde{\Sigma}$ in the multinomial probit model specified in (ref) are not identified bunch1991estimability. There are two parameter identification problems. First, the location of the latent utilities is unidentified, since $Y_i(\tilde{Z}_i+c) = Y_i(\tilde{Z}_i)$ for $c\in\mathbb{R}$. Second, the scale of the latent utilities is also unidentified, as $Y_i(c\tilde{Z}_i) = Y_i(\tilde{Z}_i)$ for $c\in\mathbb{R}^+$.
A standard solution to the identification problem in the location is to difference the utilities with respect to a baseline category. Define choice category $j=0$ as the base category and define the differences in utilities as $Z_{ij}=\tilde{z}_{ij}-\tilde{z}_{i0}$. The dependent variable $Y_i$ equals
where $\max(Z_i)$ is the largest element of $Z_i$.
The utility model in (ref) is transformed to differences in utilities by
with transformation matrix $T=\left[-\iota_J \quad I_J\right]$, the $J \times K$ transformed regression matrix $X_i=[I_{J} \quad x_{i,d}^\top \otimes I_{J} \quad Tx_{i,a}]$, and the $J \times J$ transformed covariance matrix $\Sigma=T\tilde{\Sigma}T^\top$. For the remainder of this paper we employ this location identification approach, and refer to $Z_i$ as utilities. Moreover, we define $Y=(Y_1,\dots,Y_N)^\top$, $Z=(Z_1^\top,\dots,Z_N^\top)^\top$, and $X=(X_1^\top,\dots,X_N^\top)^\top$.
There are multiple solutions to the unidentified scale in the latent utilities, but they all impose a constraint on the covariance matrix $\Sigma$. mcculloch2000bayesian develop a prior that fixes the (1,1) element of $\Sigma$ to be equal to one. burgette2012trace argue that the assignment of a choice alternative to the unit variance can have a large effect on the posterior choice probabilities. They propose to fix the trace of the covariance matrix instead of fixing one of its elements. We follow this approach and restrict the trace of the covariance matrix $\Sigma$ in (ref) to be equal to $J$.
The total number of unique parameters in the covariance matrix, $J(J+1)/2$, grows quadratically with $J$. To reduce the dimension of the parameter space, Section (ref) specifies a factor structure for $\Sigma$. Section (ref) introduces a reparametrization of the factor structure that imposes a trace restriction on the covariance matrix.
Denote as ${\gamma}$ a $J\times q$ matrix with $q<J$, and as $D$ a diagonal matrix with positive diagonal elements $d=(d_1,\dots,d_J)$. We model $\Sigma$ as
For the purpose of identification, the upper triangular elements of $\gamma$ are fixed at zero geweke1996measuring. In this factor covariance structure, the total number of parameters that characterise the covariance matrix is $n = J(q+1)-q(q-1)/2$, which implies that for a given value of $q$, the number of parameters grows linearly with $J$.
A major challenge in the factor covariance structure is the implementation of the identifying restriction $\text{trace}(\Sigma) = J$, which implies
where the scalar $\gamma_{jk}$ denotes the element in row $j$ and column $k$ in $\gamma$, while the scalar $d_j$ is the $j^{\text{th}}$ element in $d$.
The trace restriction implies that by construction, the elements of $\gamma$ and $d$ need to be constrained to the surface of an $n-$sphere of radius $\sqrt{J}$. To show this, define the $n$-dimensional vector $\psi$ as
where vech$(\gamma) =\left(\gamma_{1:J,1}^\top,\dots,\gamma_{q:J,q}^\top\right)^\top$ and $\gamma_{k:J,k}=\left(\gamma_{kk},\dots,\gamma_{Jk}\right)^\top$. The trace restriction in (ref) is equivalent to the spherical restriction $\sum_{l=1}^n\psi_l^2 =J$ on $\psi$, from which follows that the elements of $\psi$ are restricted to the $n-$dimensional spherical space $\mathbb{S}^n = \{\psi: \sum_{l=1}^n\psi_l^2 =J\}$.
As one must guarantee the {spherical} restriction, which involves all the elements in $\psi$, estimation of $\psi$ is challenging. To avoid direct implementation of this restriction, we exploit the fact that $\psi$ is restricted to $\mathbb{S}^n$, which naturally allows for a spherical transformation from $\psi$ into an $(n-1)$-dimensional vector of angles $\kappa = \left(\kappa_1,\dots,\kappa_{n-1}\right)^\top$ and the radius $\sqrt{J}$, with $\kappa\in\mathbb{A}$ and $\mathbb{A} = [0,\pi)^{n-2}\times[0,2\pi)$.
The spherical transformation reparametrises $\psi_l$ as
The main advantage of this transformation is that inference on the parameter space $\mathbb{A}$ is a more accessible problem, as no {spherical} restriction is required for $\kappa$. Moreover, the inverse function of the transformation is available in closed-form,
Notice from (ref) that $\kappa_l(c{\psi}) = \kappa_l({\psi})$, for any positive scalar $c$. Setting $c=\sum_{l=1}^n\psi_l^2=\text{trace}(\Sigma)$, shows that $\kappa_l$ is a function of the scale of $\psi_l$ relative to the other elements in $\psi$, rather than the trace of $\Sigma$.
This section develops a Bayesian method for estimating the identified parameters $(\beta,\kappa)$, subject to the trace restriction in (ref). The {density} of interest is the augmented posterior
where we assume the priors of $\beta$ and $\kappa$ to be independent, with corresponding hyperparametes $B$ and $\theta$. The prior on the coefficients $\beta$ is standard and specified as
We develop the prior choice for $p(\kappa|\theta)$ in the next section.
Developing prior beliefs directly on the angular coordinates $\kappa$ is challenging, because these coordinates lack interpretation in relation to the covariance matrix $\Sigma$. On the other hand, because $\psi$ directly determines the covariance matrix $\Sigma$ via (ref) and (ref), the implications that a prior on $\psi$ has on the prior beliefs on $\Sigma$ are well understood. At the same time, choosing a prior on $\psi$ imposes the prior on $\kappa$ required in (ref).
Therefore, we first select a prior $p\left(\psi|\theta\right)$. This prior implies $p(\kappa|\theta)$, which is required for parameter estimation, and $p(\Sigma|\theta)$. Because the elements of $\kappa$ and $\psi$ do not have a monotonically increasing relationship, analytical derivation of $p(\kappa|\theta)$ is challenging. Instead, we construct a parametric approximation to $p\left(\kappa|\theta\right)$ from a flexible parametric density class $\mathcal{P}$ with elements $\tilde{p}\left({\kappa}|{\lambda}\right)$ indexed by $\lambda\in\Lambda$.
The approximating prior density $\tilde{p}({\kappa}|\hat{\lambda})$ is calibrated by minimizing an estimate of the Kullback-Leibler divergence $\text{KL}\left[p\left(\kappa|\theta\right)||\tilde{p}\left({\kappa}|\lambda\right)\right]$, with respect to $\lambda$. Specifically, we minimize
where the first term can be ignored in the minimization problem, and the second term is computed using $M$ draws $\{\kappa^{[m]}\}_{m=1}^{M}$. These draws are produced by generating from the prior distribution $p({\psi}|\theta)$, and then transforming into draws from the prior $p\left({\kappa}|\theta\right)$. Algorithm (ref) outlines the steps of the optimization process.
Once calibrated, we can use the prior $\tilde{p}(\kappa|\hat{\lambda})$ for inference, by plugging it into (ref) instead of $p\left({\kappa}|\theta\right)$. Key to accurate approximation of $p\left({\kappa}|\theta\right)$ is the selection of a flexible parametric density class $\mathcal{P}$. Here we use $\tilde{p}(\kappa|{\lambda})=\prod_{l=1}^{n-1}\tilde{p}(\kappa_l|{\lambda}_l)$ with
where $G\left(\kappa_l\right) = \Phi_1^{-1}\left(\frac{\kappa_l}{\pi}\right)$ for $l<n-1$, $G\left(\kappa_l\right) = \Phi_1^{-1}\left(\frac{\kappa_l}{2\pi}\right)$ for $l=n-1$, $G'()$ is the derivative of $G()$, $t_{\eta}()$ is the yeo2000new transformation, $t_{\eta}'()$ its first derivative, while $\phi_1()$ and $\Phi^{-1}_1()$ denote the density and inverse distribution function of a standard normal variable, respectively.
The density function in (ref) is capable of accurately approximating the margins of the prior $p(\kappa|\theta)$, which in turn results in accurate approximation to the implied prior $p(\Sigma|\theta)$, as we will show in Section (ref). Moreover, the evaluation of the density is computationally efficient, which increases the speed of the sampling algorithm discussed in Section (ref). Therefore, we select this density class over generally more computationally involved non-parametric alternatives. For details on the properties and construction of this distribution we refer to Appendix (ref).
To set the prior $p(\kappa|\theta)$, the practitioner needs to specify a prior on the vector ${\psi}$. Although simulation from $p(\psi|\theta)$ is required for Algorithm (ref), this prior density does not have to be available in closed-form.
We propose the following characterisation of the prior distribution $p(\psi|\theta)$,
where $\nu$ and $s$ denote the shape and rate parameters of the Inverse-Gamma distribution.
This choice of prior links the trace-restricted parameters in $\psi$ to an unrestricted parameter vector $\ddot{\psi}=(\ddot{d}^\top,\text{vech}(\ddot{\gamma})^\top)^\top$. As a result, we induce a prior on the restricted parameters $\psi$ by selecting a prior for the unrestricted parameters $\ddot{\psi}$. For $\ddot{\psi}$ any prior that fits the type of factor structure in (ref) may be employed. We use a particularly well-established prior in the factor literature, that assumes a normal distribution on $\ddot{\gamma}_{jk}$ and an Inverse-Gamma distribution for $\ddot{d}^2_j$ (see for instance lopes2014modern and references therein).
Since $\psi$ is a function of the scale of the elements of $\ddot{\psi}$ relative to the norm of $\ddot{\psi}$, the location of $p(\psi|\theta)$ is unidentified. To solve this, we anchor the mean of the Inverse-Gamma prior in (ref) at one by setting $s = \nu-1$. Thus, the hyperparameters for the prior are $\theta = \left(\mu_{\gamma},\sigma_\gamma,\nu\right)^\top$.
This section assesses the accuracy of the approximating density for the choice of prior in the previous section. We focus on the parameter space of most interest to the practitioner, the covariance matrix $\Sigma$, which is implied by $\tilde{p}(\kappa|\hat{\lambda})$. Specifically, we focus on assessing the accuracy of $\tilde{p}(\Sigma_{2,2}|\hat{\lambda})$ at replicating ${p}(\Sigma_{2,2}|\theta)$, and the accuracy of $\tilde{p}(\rho_{2,3}|\hat{\lambda})$ at replicating ${p}(\rho_{2,3}|\theta)$ with $\rho_{2,3} =\frac{\Sigma_{2,3}}{\sqrt{\Sigma_{2,2}\Sigma_{3,3}}}$, for $J=6$, $q=1$, and $\theta=(0,1,5)^\top$. All the results in this section also hold for the other elements in $\Sigma$ and for any number of choice alternatives $J$.
The prior ${p}(\Sigma_{2,2}|\theta)$ is constructed via simulation. First, we generate draws from the prior $p(\kappa|\theta)$, then transform them into draws for $\Sigma$, and finally construct a kernel density estimator of $p(\Sigma_{2,2}|\theta)$. Similarly, we use Algorithm (ref) to calibrate the approximating prior $\tilde{p}(\kappa|\hat{\lambda})$, from which we also obtain a kernel density estimator of $\tilde{p}(\Sigma_{2,2}|\hat{\lambda})$ via simulation. We construct $p(\rho_{2,3}|\theta)$ and $\tilde{p}(\rho_{2,3}|\hat{\lambda})$ in the same way.
In Panel (a) of Figure (ref), the yellow solid line and the black dashed line represent the implied prior densities $p(\Sigma_{2,2}|\theta)$ and $\tilde{p}(\Sigma_{2,2}|\hat{\lambda})$, respectively. The approximating prior is an accurate representation of the prior. The remaining panels in Figure (ref) demonstrate that the approximating prior remains accurate for alternative hyperparameter values. A similar result is obtained when we compare the implied priors $p(\rho_{2,3}|\theta)$ and $\tilde{p}(\rho_{2,3}|\hat{\lambda})$, as shown in Figure (ref).
Additionally, Figure (ref) shows that also for a larger number of factors $q=4$ the approximating prior is accurate.
The user of the prior proposed in (ref)-(ref) only has to set the hyperparameters $\theta$. Here we discuss the impact that these hyperparameters have on the implied prior for $\Sigma$. The parameters $\sigma_{\gamma}$ and $\nu$ jointly control the dispersion of the variances around their prior mean of one, and $\sigma_{\gamma}$ also governs the variance of the correlations around their prior mean.
To illustrate this, Panels (b) and (d) in Figure (ref) show that, for the small value of $\sigma_{\gamma}=0.1$, a larger value of $\nu$ makes the prior on the diagonal elements of $\Sigma$ tighter around one. In Panels (a) and (c), we observe that for large values of $\sigma_{\gamma}$ the hyperparameter $\nu$ has little effect on the prior. On the other hand, comparing Panel (a) to (b), and Panel (c) to (d), indicates that for smaller values of $\sigma_{\gamma}$ the impact of $\nu$ on the prior is more pronounced.
Panels (a) and (b) in Figure (ref) show that $\sigma_\gamma$ governs the prior variance on the underlying correlations of $\Sigma$, with larger values for $\sigma_\gamma$ associated with a prior with larger variance. The parameter $\nu$ does not affect the implied prior on the correlations. So a small value for $\sigma_{\gamma}$ shrinks the correlations towards their prior means. A comparison of Panels (a) and (b) to Panels (c) and (d) indicates that the prior mean of $\Sigma$ equals $I_J$ for $\mu_{\gamma}=0$, because the correlations have a prior mean of zero, and it equals an equicorrelated matrix when $\mu_{\gamma}\ne0$.
The implied prior for $\Sigma$ is also sensitive to the number of factors considered. Figure (ref) shows how these implied prior densities change when the number of factors is set to $q=4$. Increasing the number of factors has a small effect in the priors for both the variance and correlation parameters.
For the remainder of this paper we consider two choices for $\mu_{\gamma}$. The first choice, $\mu_{\gamma}=0$, sets the prior mean of the covariance matrix $\Sigma$ equal to the identity matrix, which is a standard MNP prior choice as we will discuss in the next section. However, $\Sigma$ is the covariance matrix of the differences in utilities $Z$, for which an identity covariance matrix does not necessarily imply a symmetric correlation structure for the untransformed utilities $\tilde{Z}$. Therefore, the second choice, $\mu_{\gamma}=\mu_\gamma^*$, sets the prior mean of $\Sigma$ equal to the equicorrelated covariance matrix $\frac{1}{2}(I_J + \iota_J \iota_J^\top)$. This follows geweke1994alternative, who shrink the covariance matrix of the untransformed utilities $\tilde{\Sigma}$ to an identity matrix, which is equivalent to shrinking $\Sigma$ to this equicorrelated matrix. The value $\mu_\gamma^*$ that produces the equicorrelated matrix above can be computed as in Appendix (ref), for any given values of $\sigma_{\gamma}$, $\nu$ and $q$.
For the remaining hyperparameters, we use the values $\sigma_{\gamma}=1$, $\nu=5$, and $q=1$. These hyperparameter values produce prior densities for the variances that have close to zero probability mass at zero. They also imply prior densities for the correlations that have low probability mass at one and minus one. These properties make the method computationally stable, and guarantee that extreme values for the variances and correlations are the result of a strong signal in the data instead of highly uninformative priors.
This section compares our prior assumptions on $\Sigma$ to two well-known alternatives in the MNP literature. The first, from here on in referred to as MNP-MPR, is proposed by mcculloch2000bayesian, who set
with degrees of freedom $\delta=J+3$ and scale matrix $C=(\delta-J)(1-\tau)I_{J-1}$, where $\tau=\frac{1}{8}$. The second alternative, from here on in MNP-BN, was proposed in burgette2012trace, who specify
with degrees of freedom $s=J+3$ and scale matrix $S=I_J$.
Panel (a) in Figure (ref) compares the implied prior on $\Sigma_{2,2}$ for our proposed multinomial probit model with factor structure, from here on in referred to as MNP-FS, to MNP-BN and MNP-MPR. For MNP-FS we consider $q=1$ and $\theta=(0,1,5)^\top$. The three densities are positively asymmetric and most of their probability mass lies between zero and four. Aside from the fact that the MNP-FS density has low mass at zero, all three prior densities have a similar shape.
Panel (b) in Figure (ref) compares the implied priors on the correlation element $\rho_{2,3}$. The MNP-BN and MNP-MPR priors assign more probability mass to the edges of the support, while the MNP-FS prior is {slightly} tighter and has a spike at zero. Although not visible in the figure, the MNP-FS prior still assigns enough probability mass at extreme values of the support, to allow for accurate estimation of correlations whose true parameter values are extreme.
To construct the posterior in (ref), we employ the following MCMC sampling scheme:\\ \ \\ $\underline{\text{Sampling Scheme}}$\\ \ \ Step 1: Generate from $\beta|Z,\kappa,X,B$. \\ \ \ Step 2: Generate from $Z|\beta,\kappa,Y,X$.\\ \ \ Step 3: Generate from $\kappa|Z,\beta,X,\theta$.\\ \ \\ Steps 1 and 2 are standard Gibbs sampling steps, see for instance mcculloch1994exact. These steps require one to transform from $\kappa$ to $\Sigma$, which can be easily achieved following the instructions outlined in Table (ref).
For step 3 we employ a random walk Metropolis-Hastings sampler. At the beginning of each iteration, the elements of $\kappa$ are randomly assigned to groups of five elements. The groups are then sampled, one group conditional on the other, with the $5$-dimensional proposal density equal to the product of five independent truncated univariate normals. The variances of the proposal densities are set adaptively to target acceptance rates between $15\%$ and $30\%$.
The random assignment into groups, plus the parameter-specific adaptive steps, allow one to target parameter-specific acceptance rates without having to sample each element of $\kappa$ one at a time. roberts2009examples provide more details on adaptive MCMC, and smith2015copula provides an illustration on random allocation within MCMC. Appendix (ref) discusses the sampling steps in more detail.
{Although the reparametrization of $\Sigma$ into $\kappa$ results in a non-conjugate sampling step, the added computational cost is small. This is because Step 2 - also required for the competing MNP estimation approaches - is the most time consuming of the sampling scheme, especially for large $J$. Section (ref) illustrates this in the simulation experiment.}
{Potential alternatives to our reparametrization that produce a conjugate Gibbs sampling scheme may complicate the inference of the factor structure. For instance, the factor structure in (ref) could potentially be sampled by augmenting the parameter space with a set of latent factors as in lopes2014modern, in combination with the marginal data augmentation approach in burgette2012trace to impose the trace restriction. However, as discussed in piatek2017multinomial, measurement of latent factors in an MNP model without extra data or a priori knowledge of the factor loadings $\gamma$ is challenging.}
This section presents a numerical experiment to assess the accuracy of the parameter estimates in the proposed multinomial probit model. First we describe the simulation design, second we discuss the results, and finally we examine how sensitive the fitted probabilities are to different base category specifications.
We generate a data set from the data generating process specified in (ref) and (ref). The data closely matches the empirical application in Section (ref), with the number of discrete choices $J+1 = 50$ and number of observations $N = 5000$. The elements of the vector $x_{i,a}$ are independently generated from normal distributions with corresponding mean $\mu=0$ and variance $\sigma^2=1$, and can be interpreted as the logarithm of the prices of the choice categories. We do not include individual-specific characteristics $x_{i,d}$.
The true parameter vector $\beta_0$ consists of $J$ intercepts drawn independently from normal distributions with $\mu=0$ and $\sigma^2=\sqrt{0.5}$, and the coefficient for $x_{i,a}$ {which is} fixed at -0.7. The true covariance matrix $\Sigma_0$ is set by drawing $\tilde{\Sigma}_0$ from the $\text{Inverse-Wishart}(S,J+3)$ distribution, where the scale matrix $S$ has ones on the diagonal and 0.5 as the common off-diagonal element. We set $\Sigma_0=J\tilde{\Sigma}_0/\text{trace}(\tilde{\Sigma}_0)$.
We apply our method, MNP-FS, to the generated dataset. For the purpose of comparison we also implement the MNP-BN and MNP-MPR approaches. For our model, we set $\theta=(0,1,5)^\top$ with $q=1$ in the prior for $\kappa$. Setting the number of factors to one substantially reduces the number of covariance parameters to be estimated, especially with 50 choice alternatives. The implied prior densities at these parameter values are discussed in detail in Section (ref).
We use the prior for the coefficients specified in (ref) with $B^{-1}=0.1I_K$ for all models. This is a rather uninformative prior, given that the covariates are scaled to a variance of one in the sampler. Since settings with large choice sets are vulnerable to numerical instabilities, we avoid the use of improper prior specifications.
The posterior results are based on 200,000 iterations of the MCMC samplers, from which the first 100,000 are discarded.
Figure (ref) compares the true parameter values of $\beta_0$ and $\Sigma_0$ ($x-$axis) against the corresponding posterior mean estimates ($y-$axis). The yellow and black circles correspond to MNP-FS and MNP-BN, respectively. The closer the circles lie to the 45 degree diagonal line, the closer the posterior means are to the true parameter values.
Panel (a) in Figure (ref) {shows} that the MNP-FS approach provides more accurate estimates of $\beta_0$ than MNP-BN. This is also the case for the diagonal elements of $\Sigma_0$ in Panel (b). Finally, the results for the posterior mean estimates of the underlying correlations of $\Sigma_0$ are presented in panel (c) of Figure (ref). {The correlation estimates from MNP-FS are clearly the most tightly scattered around the diagonal line, and as such the most accurate.} {This result is particularly striking given the fact that the MNP-FS is the one imposing the most restrictive covariance structure.} Appendix (ref) shows that the comparison between MNP-FS and MNP-MPR, result in similar conclusions.
{The error measures reported in} Table (ref) confirm the {conclusions above}. The smallest root mean squared error (RMSE) of the posterior mean estimates for the variances and correlations correspond to the {MNP-FS} specification. The same conclusion is reached when comparison is conducted in terms of the mean absolute error (MAE), confirming that our approach allows for more accurate estimation of the elements in the covariance matrix. For the coefficients, the RMSE and MAE values from MNP-MPR are the smallest, closely followed by MNP-FS.
{The computational costs of the methods are similar. Since MNP-FS requires several MH steps for the generation of $\kappa$, it takes $0.06$ seconds more per sample iteration than MNP-BN and MNP-MPR, which need on average $0.36$ and $0.35$ seconds for each sample iteration in this simulation exercise.}
burgette2019symmetric show that the estimated probabilities from Bayesian MNP models can depend on the base category specification. This section analyses the sensitivity of the MNP-FS results to the base category specification under the identity prior $\theta=(0,${$1$}$,5)^\top$ and the equicorrelated prior $\theta=(${$1.525$}$,${$1$}$,5)^\top$ specifications and $q=1$.
First, we consider the identity prior. Panel (a) in Figure (ref) presents the estimated probabilities as a function of price for the least popular category, which has 16 observations. The solid black line corresponds to the true probabilities, the dashed black line is estimated using the same base category as in Section (ref) which has 85 observations, while the yellow line corresponds to the specification that has the largest category with 440 observations as the base category. In all panels, the price of the other categories is fixed at the mean across all observations.
Panel (a) in Figure (ref) shows that the yellow line and the dashed black line are different from each other. Panel (b) presents the equivalent plot for the most popular category. In this panel, there is a slightly bigger difference between the yellow and dashed black line, indicating that the estimated probabilities of the least popular choice are less sensitive to the base category specification than the most popular category.
As discussed in Section (ref), the identity prior does not necessarily imply an identity covariance matrix for the untransformed utilities. Since the utilities are in differences with the base category, this prior specification may result in estimates that are sensitive to the base category specification. In contrast, the equicorrelated prior follows from an identity covariance matrix for the untransformed utilities.
Panels (c) and (d) show the estimated probabilities for the least and most popular categories, respectively, when using an equicorrelated prior. Panel (c) indicates that the results for the least popular category are not as sensitive as those when using an identity prior. The estimated probabilities from the different base category specifications are now almost identical. We find the same results for the most popular category in Panel (d).
In sum, we conclude that the estimated probabilities are sensitive to the base category specification when a prior specification that shrinks $\Sigma$ to an identity matrix is employed. However, the impact of the base category can be substantially decreased by specifying an equicorrelated prior. An alternative way to deal with the sensitivity to the base category, is to pool the estimated probabilities across all models with different base category specifications to obtain probabilities that do not depend on one base category. This approach is computationally costly, especially when the number of choice alternatives is large.
This section fits the multinomial probit model to consumer choice data sets of different dimensions. First, Section (ref) considers a traditional consumer choice data set on laundry detergent purchases with only six alternatives. Section (ref) constructs a laundry detergent purchases data set with 50 choice alternatives from a big transaction data set. Section (ref) considers margarine brands in both a widely used small choice set and a newly constructed large choice set.
The proposed multinomial probit model is compared to benchmark specifications. These specifications and prior settings are discussed in Section (ref). Moreover, we estimate a multinomial probit model with the covariance matrix fixed to the identity matrix, referred to as MNP-I.
The in-sample and out-of-sample predictive accuracy of the models are evaluated in terms of the predictive hit-rate and the logarithmic score (log-score). The predictive probability mass function for $Y_i$ is given by
where $X_i$ denotes the attributes of the observation $i$ to be predicted. For ease of notation we refer to $p({Y}_i|X_i,Y,X,B,\theta)$ as $p({Y}_i|X_i)$. An estimate $\hat{p}({Y}_i|X_i)$ for the predictive in (ref) is constructed as the empirical probability mass implied by the draws ${Y}_i^{[m]}$ obtained from $p(Y_i|X_i,\beta^{[m]},\kappa^{[m]})$, where $\{\beta^{[m]}\}_{m=1}^M$ and $\{\kappa^{[m]}\}_{m=1}^M$ denote the MCMC draws.
The point forecast $\hat{Y}_i$ for $Y_i$ is constructed as the mode of $\hat{p}({Y}_i|X_i)$. The hit-rate is defined as
where $I[A]$ is an indicator function. The log-score is defined as
For both the hit-rate and the log-score large values are preferred. As a rough measure of statistical significance, we test the difference of the hit-rates between MNP-FS and the benchmark models by a normal test and the difference of the log-scores by the giacomini2006tests test.
We randomly allocate 80% of the observations for estimation of the model, and the remaining 20% are employed for out-of-sample evaluation.
imai2005bayesian and burgette2019symmetric fit multinomial probit models to purchase data of laundry detergents. This data contains purchases of 2657 households out of six brands of laundry detergents, and the log price of each brand. The data set is described in detail by chintagunta1998empirical and available in imai2005mnp. We follow imai2005bayesian and burgette2019symmetric and fit the multinomial probit models with an intercept and the log price for each brand.
Table (ref) shows that for a small choice set the hit-rates and log-scores of our proposed model are similar to the ones of the benchmark methods. The symbol $(^+)$ indicates that MNP-FS significantly outperforms the benchmark method, while the symbol $(^-)$ indicates that the benchmark significantly outperforms MNP-FS. Only the in-sample log-score of MNP-BN significantly improves upon MNP-FS. All models show substantial improvements over the naive forecast, which is constructed as the observed in-sample frequency of the choice alternatives. However, MNP-FS, MNP-BN and MNP-MPR do not have a better in-sample hit-rate or out-of-sample log-score than the MNP-I.
The posterior parameter estimates from MNP-FS, MNP-BN and MNP-MPR are also similar. For instance, the posterior mean of the correlation matrix of the latent utilities presented in Figure (ref), has similar patterns across all three methods. The small differences in the posterior estimates are also observed for the price coefficient. Figure (ref) shows the posterior densities for the different models. The densities are concentrated around similar values.
Nowadays, almost all real-life consumer choice sets contain many more choice alternatives than six. To illustrate the importance of a scalable multinomial probit model in these settings, we analyse a laundry detergent purchase data set with 50 choice alternatives.
We use the Complete Journey dataset published by Dunnhumby\footnote{https://www.dunnhumby.com/sourcefiles}. This dataset contains all purchases of 92,339 products over two years from a group of 2,500 households at a retailer. We filter the purchases of products with the description “Laundry Detergents", which results in 300 unique products with different brands, sizes and variants such as liquid or powder detergents. Since the unique products with a small purchase volume are of less interest to a marketing manager, we focus on the 50 top-selling products.
We define the log price of each brand in the same way as, for instance, allenby1991quality and wan2017modeling. The Dunnhumby data set only contains records of shelf prices at purchase dates. We impute the prices for products that are not sold on a certain purchase date by taking the mean of the observed prices of a specific product 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. We remove the observations for which we cannot impute a price for each product.
The final sample contains 4839 observations on 50 categories, which contains 64% of the laundry detergent purchases and the purchase frequency varies from 30 to 274 per category. We fit the same models as in the exercise with six choice alternatives, also considering an intercept and log price coefficient.
For this large choice set, Table (ref) shows that the in-sample and out-of-sample log-score of our proposed model are larger than those of the benchmark models, with the out-of-sample improvement also being statistically significant. The hit-rates are not statistically different from those of MNP-BN and MNP-MPR. However, the MNP-FS reports both hit-rates and log-scores significantly larger than that of MNP-I. Comparing this result to the relative performance of MNP-FS and MNP-I with six choice alternatives, suggests that accounting for correlations across utilities is especially important when the choice set is large.
For the large laundry detergent choice set, the posterior parameter estimates show differences across the different models. Figure (ref) presents the posterior means of the elements of the correlation matrix of the latent utilities. While some general patterns are common to all three methods, the posterior correlations of MNP-BN and MNP-MPR show more variation. The factor structure in MNP-FS has fewer parameters, which restrict the patterns in the correlation matrix. Parsimony may become more important in this large choice set, as the number of estimated covariance parameters by MNP-BN and MNP-MPR equals $49\times50/2=1225$ relative to 4839 observations. This might explain the fact that the predictive accuracy of the MNP-FS, which only estimates $49\times2=98$ covariance parameters, significantly improves upon benchmark models in the large choice set, while the differences are not significant in the small choice set.
Figure (ref) reports the posterior densities for the price coefficient. We find that the differences in the specifications of the covariance matrix are also reflected in the posterior of the price coefficient. The posterior mean of the price coefficient in the MNP-FS model equals -0.262, compared to -0.372 and -0.958 in the MNP-BN and MNP-MPR respectively. The corresponding posterior standard deviation is respectively 0.018, 0.025, and 0.031. Hence we conclude that the benchmark specifications result in a larger price effect estimate with more posterior uncertainty than the MNP-FS model.
One might argue that in most practical settings only a small set of high volume products are of interest. To examine the importance of considering a large choice set, we compare the effect of price on the purchase probability of the six most popular products between an MNP model estimated on only those six products and an MNP model estimated on the total choice set. Figure (ref) suggest that, even when only the six top selling products are of interest, including the other products in the model is important for an effective pricing strategy.
Panel (a) of Figure (ref) shows the probability of buying one of the top six selling products as a function of the price of the most popular product. The dashed black line corresponds to MNP-FS with only the top six products included as choice alternatives, and the solid yellow line to MNP-FS with 50 products included. The model that includes all 50 products indicates that the probability of buying in the top six decreases when the price of the top product increases. In other words, consumers substitute away from the top product to products outside the top six. This effect cannot be captured by the model that only includes the top six products, which sets the probability of buying in the top six equal to one by construction.
Panel (b) of Figure (ref) shows the probability of buying the most popular product as a function of its price. The dashed black line shows that demand for laundry detergent is inelastic between a price of zero and four, and highly elastic for prices higher than four. However, the solid yellow line shows that including the purchases of all 50 products, results in purchase probabilities that change smoothly with price. This result suggests that excluding purchases of low volume products from the analysis can bias the estimated price effect of high volume products.
We compare the hit-rates and log-scores of the MNP-FS to the benchmark methods on another commonly used choice data set. mcculloch1994exact, burgette2012trace, and burgette2019symmetric fit multinomial probit models to panel data on purchases of margarine by 516 households. They estimate an intercept and a log price coefficient for six brands. The data is described in detail by allenby1991quality and available in rossi2012bayesian. We follow burgette2012trace and burgette2019symmetric and fit the multinomial probit models to the first purchase of each household.
As an alternative for the small margarine choice set, we construct a large margarine choice data set in the same way as for laundry detergents in Section (ref). We filter the purchases of products with the description “Margarines" from the Complete Journey dataset which results in 178 unique products. The final sample contains 11754 observations on 50 categories, which contains 96% of the margarine purchases and the purchase frequency varies from 32 to 1206 per category.
Table (ref) shows that for the margarine data our proposed model has higher in-sample and out-of-sample hit-rates and log-scores than the benchmark models when the number of choice alternatives is large. The log-scores of the MNP-FS are significantly larger than the log-scores of all the benchmark models. The MNP-FS model does not show significantly different performance from the MNP-I in the small choice set, but shows significant improvements in the large choice set on all metrics. This result is in line with the laundry detergent application and supports the claim that it is important to take correlations into account in large choice sets.
Table (ref) in Appendix (ref) reports the hit-rates and log-scores for MNP-FS when the equicorrelated prior is considered. The results are similar and indicate that employing an equicorrelated prior does not necessarily lead to an increase in predictive performance.
This paper proposes a factor structure on the covariance matrix in the multinomial probit model that makes the model scalable to the dimensions of modern choice sets. The model parameters are identified by a reparamatrization of the factor structure that imposes a trace-restriction on the covariance matrix.
A numerical experiment shows that the model parameters can be accurately estimated on a choice set with 50 alternatives in the proposed multinomial probit specification. On a real data set with 50 choice alternatives, the hit-rates and log-scores demonstrate {significant predictive} improvements relative to benchmark approaches.
The large size of modern assortments and the increasing amount of product differentiation makes the scalable choice model that accounts for correlations across choice alternatives of managerial relevance. The empirical application to retail data of 50 laundry detergents suggests that managers may overestimate the price effect when only analysing top-selling products, relative to including all products in the multinomial probit model.