{'fraction_non_alphanumeric': 0.034266528415543375, 'fraction_numerical': 0.015158248160766306, 'mean_word_length': 4.824992809893587, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Robust discrete choice models with t-distributed kernel errors

Outliers in discrete choice response data may result from misclassification and misreporting of the response variable and from choice behaviour that is inconsistent with modelling assumptions (e.g. random utility maximisation). In the presence of outliers, standard discrete choice models produce biased estimates and suffer from compromised predictive accuracy. Robust statistical models are less sensitive to outliers than standard non-robust models. This paper analyses two robust alternatives to the multinomial probit (MNP) model. The two models are robit models whose kernel error distributions are heavy-tailed t-distributions to moderate the influence of outliers. The first model is the multinomial robit (MNR) model, in which a generic degrees of freedom parameter controls the heavy-tailedness of the kernel error distribution. The second model, the generalised multinomial robit (Gen-MNR) model, is more flexible than MNR, as it allows for distinct heavy-tailedness in each dimension of the kernel error distribution. For both models, we derive Gibbs samplers for posterior inference. In a simulation study, we illustrate the excellent finite sample properties of the proposed Bayes estimators and show that MNR and Gen-MNR produce more accurate estimates if the choice data contain outliers through the lens of the non-robust MNP model. In a case study on transport mode choice behaviour, MNR and Gen-MNR outperform MNP by substantial margins in terms of in-sample fit and out-of-sample predictive accuracy. The case study also highlights differences in elasticity estimates across models.


Introduction
Discrete choice models are widely used to explain and predict the demand for discrete alternatives such as retail products, medical treatments and modes of transport. Discrete choice models are usually formulated based on random utility theory (McFadden, 1981) which posits that a rational decisionmaker chooses the alternative with the highest utility from a set of mutually exclusive alternatives. In principle, utility is not identifiable in absolute terms, and only differences in utility matter. Utility differences are typically assumed to follow independent and identical logistic or multivariate Gaussian distributions. Whereas the former assumption leads to a logit kernel, the latter assumption leads to a probit kernel.
Logit is popular in practice due to its closed-form choice probabilities. However, logit suffers from two drawbacks. First, the independence of irrelevant alternatives property implies restrictive substitution patterns. Second, logit assumes that the random error terms are homoskedastic across alternatives, which implies that the same level of decision uncertainty applies to all alternatives in a choice set. Probit overcomes the limitations of logit by allowing for the estimation of a full error covariance, subject to identification restrictions (Train, 2009). Probit choice probabilities lack closed-form expressions. However, advances in computational power combined with progress in simulation techniques (Burgette and Nordheim, 2012;Hajivassiliou et al., 1996;Imai and Van Dyk, 2005;McCulloch and Rossi, 1994;Train, 2009) and analytical approximation approaches (Bhat, 2011) have attenuated the relevance of this drawback.
Both logit and probit lack robustness to outliers in the response data due to their light-tailed kernel error distributions (Benoit et al., 2016;Hausman et al., 1998). Outliers in discrete choice response data may result from various contamination processes, including misclassification and misreporting of the response variable as well as choice behaviour that is inconsistent with the postulated decision rule (e.g. random utility maximisation). Ignoring outliers in logit and probit applications can lead to biased and inconsistent parameter estimates (Hausman et al., 1998). Even though response data contamination is a well acknowledged problem in discrete choice analysis (Hausman et al., 1998;Paleti and Balan, 2019), the development of robust multinomial choice models is understudied.
Extant approaches lack flexibility and suffer from computational drawbacks.
In this paper, we address this research gap by analysing two robust alternatives to the multinomial probit (MNP) model. Both models belong to the family of robit models whose kernel error distributions are heavy-tailed t-distributions to moderate the influence of outliers. In the first model, the multinomial robit (MNR) model, a generic degrees of freedom parameter controls the heavy-tailedness of the kernel error distribution. The second model, the generalised multinomial robit (Gen-MNR), is more flexible than MNR, as it allows for different marginal heavy-tailedness of the kernel error distribution.
For both models, we devise Gibbs samplers for posterior inference. We first use simulated data to investigate the properties of the proposed models and their estimation methods in terms of parameter recovery and elasticity estimates. Then, we compare MNP, MNR and Gen-MNR in a case study on transport mode choice behaviour in London, UK.
The remainder of this paper is structured as follows: In Section 2, we elaborate on the research gap and discuss connections to related work. In Section 3, we exhibit the mathematical formulations of MNP, MNR and Gen-MNR. In Section 4, we describe inference and implementation details. In Section 5, we present the simulated examples; in Section 6, we present the case study. Finally, we conclude in Section 7.

Related work
Robust statistical models are less sensitive to outliers than standard "non-robust" models (Gelman et al., 2013). Heavy-tailed distributions such as the t-distributions are an appealing choice for enhancing the robustness of statistical models (Gelman and Hill, 2006;Lange et al., 1989). Compared to the Gaussian distribution, the t-distribution has one more parameter, called the degrees of freedom (DOF) parameter, which controls the heavy-tailedness of the distribution and moderates outliers.
There are several pertinent examples of outlier robust models that are based on t-distributions.
In the context of generalised linear models, Liu (2004) proposes the binary robit model, which is based on a t-distribution with unknown degrees of freedom, as a robust alternative to logistic and probit regression models. Furthermore, Ding (2014) constructs a robust Heckman selection model using the t-distribution as kernel error distribution. Jiang and Ding (2016) formulate Heckman selection and multivariate robit models based on t-distributions with different marginal DOF. However, robust multinomial choice models with heavy-tailed kernel error distributions have received limited attention.
To the best of our knowledge, only three studies formulate robust multinomial choice models using heavy-tailed distributions. Dubey et al. (2020) present the first multinomial robit (MNR) model, i.e. a multinomial choice model defined through a t-distributed kernel error with an estimable DOF. Dubey et al. (2020) make a strong empirical case to adopt the MNR model over the multinomial probit (MNP) model. First, the estimates of the MNP model are inconsistent, if the kernel errors in the data generating process are heavy-tailed. Second, the robustness of MNR results in superior in-sample fit and out-of-sample predictive ability for class-imbalanced datasets. Furthermore, Peyhardi (2020) formulates a MNR model in the context of generalised linear models and shows that MNR can help in identifying artificial aspects in the design of stated preference experiments. Finally, Benoit et al. (2016) devise a multinomial choice model, in which utility differences follow a symmetric and heavy-tailed multivariate Laplace distribution. Using simulated and real data, the authors show that their proposed model is less sensitive to outliers than MNP.
We identify two research gaps in the formulation and estimation of existing robust multinomial choice models. First, the kernel error distributions of existing models lack flexibility. Dubey et al. (2020) and Peyhardi (2020) constrain the flexibility of the kernel error distribution by assuming that a single, generic DOF parameter controls the heavy-tailedness of the kernel error distribution.
This modelling assumption implies that all utility differences between different pairs of alternatives exhibit the same level of aberrance. Unlike the t-distribution, the Laplace distribution underlying the formulation of the multinomial choice model proposed by Benoit et al. (2016) cannot exhibit varying levels of heavy-tailedness, since the Laplace distribution does not have a third parameter controlling the heavy-tailedness of the distribution. Second, the estimation approaches employed in the studies by Dubey et al. (2020) and Peyhardi (2020) are not scalable. Dubey et al. (2020) are unable to derive analytical gradients of the MNR model and thus rely on computationally-expensive numerical gradient approximations during the maximisation of the simulated log-likelihood of the model. Peyhardi (2020) estimates the DOF parameter by performing a grid search, which requires the model to be estimated at multiple values of the DOF parameter. Clearly, this approach suffers from the curse of dimensionality if the underlying kernel distribution has multiple DOF parameters.
Importantly, in both approaches, hierarchical model extensions such as random parameters, latent variables and shrinkage priors are computational expensive, since these extensions necessitate an additional layer of simulation in the computation of the log-likelihood.
In this paper, we address the first limitation of existing MNR models (i.e. generic heavy-tailedness) by formulating a generalised MNR (Gen-MNR) model with alternative-specific DOF parameters. To that end, we adopt the non-elliptical contoured t-distribution (Jiang and Ding, 2016) as kernel error distribution. To address the second limitation (i.e. computationally-expensive estimation), we devise gradient-free Bayesian estimation approaches for both the MNR and the Gen-MNR models. In the construction of the Bayesian estimation approaches, we exploit the hierarchical normal mixture representation of the t-distribution. To circumvent complex likelihood computations in the estimation of the MNR and the Gen-MNR models, we employ a combination of Bayesian data augmentation techniques used in the estimation of MNP models (Albert and Chib, 1993;McCulloch and Rossi, 1994) as well as of non-multinomial robit models (Ding, 2014;Jiang and Ding, 2016). Bayesian estimation also facilitates potential hierarchical model extensions.

Multinomial probit (MNP)
We consider a standard random utility model in which an agent i = 1, . . . , N chooses from a set of J mutually exclusive alternatives. In principle, utility is not identified at an absolute level.
Therefore, the MNP model is defined through a J − 1-dimensional Gaussian latent variable vec- McCulloch and Rossi, 1994). The elements of w i correspond to the utility differences with respect to the base alternative J. The observed choice y i ∈ {1, . . . , J} is assumed to arise from The latent variable w i is represented as Here, where X obs i j is the observed attribute vector of alternative j for agent i. β is a K vector of taste parameters. Σ is a (J − 1) × (J − 1) covariance matrix. The latent variable representation (2) is not identified, because w i can be multiplied by any positive scalar c without changing the likelihood (1), i.e. y i (w i ) = y i (c w i ). Therefore, we must set the scale of the model. Following Burgette and Nordheim (2012), we impose a trace restriction on Σ with tr(Σ) = J − 1. We implement the trace restriction using a constrained inverse Wishart prior (Imai and Van Dyk, 2005;Burgette and Nordheim, 2012). To construct this prior, we introduce a working parameter α, which is not identified given the observed data y, but is identified given {y, w }. We transformw i = αw i for i = 1, . . . , N so thatw i ∼ N (X iβ ,Σ) withβ = αβ andΣ = α 2 Σ. We then put an inverse Wishart prior on the intermediate, unidentified quantityΣ so thatΣ ∼ IW (ρ,Λ). After the transformation Σ =Σ/α 2 with α 2 = tr(Σ)/(J − 1), the implied prior on the tuple {Σ, α 2 } is where α 0 is some positive constant such thatΛ = α 0 Λ. The conditional distribution of α 2 is To complete the specification of the MNP model, we specify the prior β ∼ N (0, B −1 0 ). It is known that predictions under the Bayesian formulation of the MNP model can be sensitive to the selection of the base alternative J (Burgette and Nordheim, 2012).

Multinomial robit (MNR)
The MNR model assumes a t-distributed kernel error for the latent variable w i , i.e.
where Σ is a (J − 1) × (J − 1) covariance matrix and ν is a scalar DOF parameter. The t-distribution has the following normal mixture representation (Ding, 2014): The latent variables q = {q 1 , . . . , q N } allow for heavy-tailedness in the distribution of the kernel error by increasing the variability of i across different i. Figure 1 illustrates the relationship between the χ 2 -distribution (which controls the distribution of q ) and a t-distribution (which controls the distribution of ) with unit variance for different DOF ν. For small ν < 30, the t-distribution exhibits heavy tails. As ν approaches ∞, the t-distribution converges to the normal distribution. We use the same priors for β and Σ as in MNP and also introduce the working parameter α. In addition, we put a Gamma prior on ν with ν ∼ Gamma(α 0 , β 0 ). Predictions under the Bayesian formulation of the MNR model can be sensitive to the selection of the base alternative in the same way as predictions under the Bayesian formulation of the MNP model.

Generalised multinomial robit (Gen-MNR)
We generalise MNR by allowing for different marginal heavy-tailedness in the distribution of the latent variable w i . The Gen-MNR model assumes that the kernel error of w i is drawn from a non-elliptically contoured t-distribution (NECT; Jiang and Ding, 2016). We have where Σ is a (J − 1) × (J − 1) covariance matrix and ν = {ν 1 , . . . , ν M } is a M vector of DOF with The NECT distribution has the following normal mixture representation (Jiang and Ding, 2016): where where m( j) maps dimension j onto its associated DOF. In the rest of this work, we assume that M = J − 1 without loss of generality. The Gen-MNR model uses the same prior distributions as the MNR model and also includes the working parameter α. We let ν j ∼ Gamma(α 0 , β 0 ) for j = 1, . . . , J − 1. Predictions of the Gen-MNR model can be sensitive to the selection of the base alternative in the same way as predictions of the MNP and MNR models.

Inference and implementation details
For the estimation of the MNP, MNR and Gen-MNR models, we employ Markov chain Monte Carlo methods in the form of Gibbs sampling (Robert and Casella, 2013). The sampling schemes for MNP, MNR and Gen-MNR are presented in Algorithms 1, 2 and 3, respectively. Algorithm 1 is based on the Gibbs samplers proposed by Burgette and Nordheim (2012) and Imai and Van Dyk (2005). Algorithms 2 and 3 are extensions of Algorithm 1. Whereas Algorithm 2 incorporates elements of the Gibbs sampler proposed by Ding (2014) for the robust Heckman selection model, Algorithm 3 incorporates elements of the Gibbs sampler proposed by Jiang and Ding (2016) for the multivariate robit model.
Data augmentation (Tanner and Wong, 1987) facilitates the construction of the samplers. The central idea of Bayesian data augmentation is to treat latent variables as unknown model parameters, which are imputed in additional sampling steps. Each of the samplers uses the data augmentation scheme developed by Albert and Chib (1993) and McCulloch and Rossi (1994) to impute the latent variable w (see Appendix A.1 for details). The samplers for the MNR and Gen-MNR models additionally incorporate the data augmentation schemes devised by Ding (2014) and Jiang and Ding (2016), respectively, to impute the latent variable q . Data augmentation circumvents complex likelihood calculations in the estimation of the MNP, MNR and Gen-MNR models. This is because conditional on w and q (if applicable), the models reduce to standard Bayesian linear models.
Furthermore, following Imai and Van Dyk (2005), we leverage marginal data augmentation (MDA; Van Dyk and Meng, 2001;Van Dyk, 2010) to increase the rate of convergence of the MCMC algorithms.
The central idea of MDA is to marginalise out the working parameter α, which is not identified given the observed data y, but is identified given the augmented data {y, w }, in some of the conditional updates of the model parameters in order to improve mixing of the Markov chains. Imai and Van Dyk (2005) argue that because α 2 is assigned a distribution with positive variance in (4), the marginal model P(y, w |θ , α)p(α|θ )dα is more diffuse and thus mixes faster than the conditional model The full conditional distribution of ν in the MNR model as well as the full conditional distributions of ν j and q i j in the Gen-MNR model are nonstandard. To draw from these intricate distributions, we use Metropolised Independence samplers (Liu, 2008) with approximate Gamma proposals, as devised by Ding (2014) and Jiang and Ding (2016) (see Appendices A.2 and A.3 for details).
Draws from the univariate truncated normal distribution are generated using a combination of inverse transform sampling and rejection sampling with a Rayleigh proposal (Botev and l'Ecuyer, 2016). This hybrid approach is adopted to warrant reliable and efficient posterior simulations. On its own, the inverse transform sampling method breaks down due to numerical issues when the truncation interval is too far in the tail. Rejection sampling remedies this problem by allowing for reliable simulation in the tail. The Rayleigh proposal has been shown to be faster than other proposal distributions. For a detailed analysis of sampling methods for the truncated normal distribution, the reader is directed to Botev and l'Ecuyer (2016).
We implement the MCMC algorithms for all models considered in this paper in Julia (Bezanson et al., 2017). 1 In the subsequent applications, the Gibbs samplers are executed with a single chain consisting of 100,000 draws including a warm-up period of 50,000 draws. A thinning factor of 10 is applied to the post warm-up draws. Convergence is assessed with the help of the potential scale reduction factor (Gelman et al., 1992).
Predictive choice distributions are obtained using simulation since the considered models lack closed-form choice probabilities. Choice probabilities in MNP are calculated using the GHK simulator (Hajivassiliou et al., 1996;Train, 2009); choice probabilities in MNR and Gen-MNR obtained using frequency simulators (Lerman and Manski, 1981;Geweke et al., 1994). Frequency simulators are fast and easy to implement and provide unbiased approximations of choice probabilities (Geweke et al., 1994). However, in some applications, frequency simulators may have a positive probability that the 1 The estimation code is available at https://github.com/RicoKrueger/robit. simulated choice probability of an alternative is zero for any finite number of draws (Train, 2009).
Zero simulated choice probabilities are especially likely when the underlying true choice probability is small (Train, 2009). Zero probabilities create issues in the calculation of the log-likelihood because the logarithm of zero is undefined (Train, 2009). In the applications presented in this paper, we only encountered the problem of zero simulated choice probabilities in the MNP model. Therefore, we adopt the GHK simulator for the calculation of predictive choice distributions in MNP. The recursive GHK simulator does not suffer from the issue of zero simulated choice probabilities and has been shown to provide accurate and reliable approximations of probit choice probabilities (Geweke et al., 1994).
Update ν j as shown for ν in Algorithm 1, Step 2. end for Step 3: Update w , α 2 as described in Algorithm 1, Step 1.

Simulation study
We consider two simulated examples to i) assess the ability of the proposed Gibbs samplers to recover parameters in finite samples and ii) illustrate the effect of ignoring the non-normality of the kernel error distribution on the sensitivities of the choice distributions to changes in an exogenous variable.

Example I: Data generated according to MNR model
In the first example, we simulate data according to the generative process underlying the MNR model.
We let N = 10, 000 and J = 4. We set β = (−0.5, 0.5, −0.5, 1, −1) , where the first three entries are bias terms pertaining to alternative-specific constants and the last two entries are utility parameters pertaining to alternative-specific attributes. Furthermore, we let Σ = Since the first three predictors are alternative-specific constants, we have X obs i jk = 1 for i = 1, . . . , N , j = 1, 2, 3 and k = 1, 2, 3. The remaining predictors are alternative-specific attributes. We draw . . , J and k = 4, 5. The fourth alternative is set as reference alternative both in data generation and model estimation. For the sake of simplicity, we do not perform a search over the specification of the reference alternative. We employ the synthetic data to estimate MNP, MNR and Gen-MNR using the methods described in Section 4.   Figure 3 shows the posterior means of the predicted choice distributions as a function of attribute X obs jk with j = j and k = 4. To summarise the results shown in Figure 3 and to quantify the posterior uncertainty in the predictions, we also compute the posterior distributions of the Euclidean distances between the true and the predictive choice distribution sensitivities. We summarise these posterior distributions with mean, standard deviation and 95%-credible intervals in Table 1.
From Figure 3, it can be seen that MNR and Gen-MNR are able to recover the choice distributions of the true model in all considered scenarios. By contrast, the MNP model, which is unable to capture the heavy-tailedness of the kernel error distribution, produces biased choice distributions.
These insights are corroborated by Table 1. In all considered scenarios, the posterior means of the Euclidean distances between the true and predicted choice distribution sensitivities are substantially smaller for MNR and Gen-MNR than for MNP. In nearly all considered scenarios, the credible intervals of the posterior distribution of the Euclidean distance between the true and predicted choice distribution sensitivity for MNP does not overlap with the credible intervals for MNR and Gen-MNR.
In addition, the credible intervals of MNR and Gen-MNR overlap in all scenarios, which indicates that the differences between MNR and Gen-MNR are not statistically significant at a 5%-significance level.
For example, for j = 4 and alternative 2, we observe substantial differences between MNP and the true model Figure 3, while MNR and Gen-MNR appear to coincide with the true model. According to Table 1, the posterior mean of the Euclidean distance between the true model and MNP is 0.310, while the posterior mean of the Euclidean distance between the true model and MNR (Gen-MNR) is 0.048 (0.080). Table 1 further suggests that for the same scenario, the credible interval of the posterior distribution of the Euclidean distance between the true and predicted choice distribution sensitivity for MNP does not overlap with the credible intervals for MNR and Gen-MNR, while the credible intervals for MNR and Gen-MNR coincide with each other.
In sum, these observations suggest that MNR and Gen-MNR perform significantly better than MNP at recovering the true choice distribution sensitivities.

Example II: Data generated according to Gen-MNR model
In the second example, we consider data generated according to the Gen-MNR model. The data generating process is essentially same as in Example I, with the only difference we allow for different marginal heavy-tailedness by setting ν = (5, 3, 1) .    Table 2).
Overall, Figure 5 indicates that in Example II with distinct marginal heavy-tailedness, only Gen-MNR is able to recover the choice distributions of the true model in all scenarios. This is because MNR allows for distinct marginal heavy-tailedness, whereas MNR only allows for generic heavy-tailedness and MNP does not capture any heavy-tailedness.
Upon closer inspection, Figure 5 and Table 2 highlight interesting differences in the abilities of the three models to recover choice probabilities. The differences between the three models and the true model are relatively minor for alternative 1 (see the first column of Figure 5). Table 2 largely confirms this observation. For alternative 1, the differences between the three models are comparatively minor, as the credible intervals of the posterior distributions of the Euclidean distances between the true and predicted choice distribution sensitivities of the three models overlap in all scenarios except when j = 1.
For alternative 2, MNP is biased, while MNR and Gen-MNR perform well at recovering the choice distribution of the true model (see the second column of Figure 5). For alternative 2, Table 2 shows that the posterior means of the Euclidean distances between the true and predictive choice distribution sensitivities are considerably larger for MNP than for MNR and Gen-MNR; also, the credible intervals of the posterior distribution of the Euclidean distances for MNP do not overlap with the credible intervals for MNR and Gen-MNR.
For alternative 3, only Gen-MNR is able to recover the choice distribution of the true model, whereas MNP and MNR produce biased choice distributions (see the third column of Figure 5). Table 2  For alternative 4, Figure 5 and Table 2 suggest that the differences between the three models are comparatively minor when j ∈ {1, 3}. However, when j ∈ {2, 4}, only MNR and Gen-MNR perform well at recovering the true choice distributions, while MNP is biased.
In sum, these patterns of differences are reflective of the distinct marginal heavy-tailedness of the kernel error distribution induced by the different elements of the degrees of freedom parameter ν = (ν 1 , ν 2 , ν 3 ) = (5, 3, 1) . Only Gen-MNR is able to account for the extreme heavy-tailedness caused by ν 3 = 1 in the true model.

. Data and utility specification
Revealed preference data for the case study are sourced from the London Passenger Mode Choice (LPMC) dataset (Hillel et al., 2018). The LPMC dataset consists of trip records from the London Travel Demand Survey, which was conducted from 2012 to 2015. For each trip record, Hillel et al. (2018) imputed tailored choice sets including the attributes of the chosen and the non-chosen alternatives using an online directions application programming interface. As the LPMC dataset is a revealed preference mode choice dataset derived from a household travel survey, we conjecture that the observed response data may be contaminated due to misreporting and misclassification. In addition, it is possible that random utility maximisation does not accurately represent the decision rules governing some of the observed choices.
In this case study, we restrict our analysis to home-based trips reported by individuals who are at least 12 years old. The original dataset contains 58,584 of such observations from 26,904 individuals.
For our analysis, we consider 10,820 observations from 5,000 randomly selected individuals for model training. In addition, we randomly select one observation each from an additional set of 1,250 randomly selected individuals for out-of-sample validation. There are four mode choice alternatives, namely walking, cycling, transit and driving. In the training data, the observed market shares are 16.4%, 3.0%, 37.8% and 42.8%, respectively. In the test data, the observed market shares are 14.5%, 3.9%, 41.8% and 39.9%, respectively.
Each of the three models uses the same specification of the systematic utility, as shown in Table 3.
The variable "traffic variability" is a measure of the driving travel time uncertainty for the given origin-destination pair. It is defined as the difference between the travel times in a pessimistic traffic scenario and in an optimistic one divided by the travel time in a typical, best-guess traffic scenario (see Hillel et al., 2018).
The drive alternative is set as reference alternative in the estimation of all models. We performed a search over the specification of the reference alternative but found no substantive differences in parameter estimates, in-sample fit and out-of-sample predictive accuracy for different specifications of the reference alternative.

Variable Walk Cycle Transit Drive
Alternative-specific constants β asc, cycle β asc, transit β asc, drive    with just one more parameter. We also observe that MNR and Gen-MNR provide a slightly better out-of-sample accuracy than MNP. In practice, Gen-MNR should be selected over MNP and MNR

In-sample fit and out-of-sample predictive accuracy
given that Gen-MNR provides the best in-sample fit and out-of-sample predictive accuracy out of the three considered models. Regardless of this recommendation, we contrast the three models in terms of parameter and elasticity estimates in the following subsections in order to highlight practical differences between the three models.  Table 4.: In-sample fit and out-of-sample predictive ability

Parameter estimates
In Table 5, we summarise the posterior distributions of the parameters of MNP, MNR and Gen-MNR.
For each parameter, we report the posterior mean, posterior standard deviation and bounds of the 95% credible interval.
The posterior distributions of the DOF parameters provide information about the heavy-tailedness of the kernel error distributions of MNR and Gen-MNR. The posterior mean of the generic DOF parameter ν in MNR is 2.132, which indicates that the heavy-tailedness of the respective kernel error distribution is substantial. The posterior distributions of the DOF parameters ν walk−drive , ν cycle−drive , ν transit−drive in Gen-MNR suggest that the heavy-tailedness of the respective kernel error distribution is also substantial, while the extent of heavy-tailedness varies significantly across the dimensions of the kernel error distribution. Heavy-tailedness is most pronounced for the utility differences All three models suggest that the considered alternative-specific attributes (see Table 3) have a statistically significant and negative influence on mode choice propensities. Since the scales of the estimated utility parameters β are not guaranteed to be identical in the three models, the estimated effects cannot be compared in absolute terms. Therefore, we contrast the sensitivities to alternativespecific attributes in terms of their implied willingness to pay (WTP), which is scale free and allows for a money-metric representation sensitivities. In general, WTP for an attribute k is defined as the ratio of the marginal utility with respect to attribute k and the marginal utility with respect to price.
If the utility is specified as linear-in-parameters with only main effects as in the current application, WTP for an attribute k is the given by the ratio β k /β p of the non-price coefficient β k pertaining to the attribute of interest k and the price coefficient β p .  The models also provide insights into the influence of individual-and context-specific attributes on mode choice propensities. The differences between the three models are minor. For example, all models suggest that female travellers are relatively less likely to cycle and relatively more likely to use transit. While MNP and Gen-MNR suggest that women are more likely to select the driving mode, MNR does not suggest any gender differences in the propensity to use the driving mode. In all three models, old age reduces the propensities to use transit and the driving mode. Furthermore, travel during winter months is found to reduce the propensity to cycle in all models. As expected, higher levels of car ownership increase the propensity to select the driving mode. measures the responsiveness of the quantity demanded of an alternative j by a group of observations relative to a change in attribute k of alternative j .
An aggregate arc elasticity as the ratio of the percentage change in the quantity demanded and the percentage change of the attribute. Thus, we have where W ( j) and W * ( j) denote the predicted market shares of alternative j in the group under consideration before and after the change, respectively. Similarly, X j k and X * j k denote the average of attribute k of alternative j in the group before and after the change, respectively. In the current application, we compute aggregate arc elasticities for the whole sample. However, in principle, it is also possible to calculate aggregate arc elasticities for different exogenously defined subgroups of the sample. Table 7 shows that the aggregate arc elasticity estimates can differ substantially across the three models. Overall, the results illustrate that Gen-MNR is more flexible than MNP and MNR. For example, the elasticities of the demand for cycling with respect to changes in cycling travel times differ quite considerably across the three models. For a 10% decrease in the out-of-vehicle time of cycling, the own-elasticities of cycling are −0.94, −0.76 and −1.03 in MNP, MNR and Gen-MNR, respectively.
Reducing cycling travel times is an important policy-relevant scenario. For instance, the construction of cycling superhighways (e.g. Rayaprolu et al., 2020) and a stimulation of e-bike uptake (e.g. Dill and Rose, 2012) could result in widespread decreases of cycling travel times. Given the differences in the estimated aggregate arc elasticities, Gen-MNR makes a stronger case for the effectiveness of these interventions than MNP and MNR.
The three models also produce different elasticities for changes to walking out-of-vehicle travel times. The demand for walking is estimated to be more elastic to changes in walking time in MNR and Gen-MNR than in MNP. More specifically, the aggregate arc elasticity of walking demand for a 10% reduction in walking time is −1.78 in both MNR and Gen-MNR but is −1.57 in MNP. Besides, the demand for cycling is estimated to be more elastic to changes in walking time in Gen-MNR than in MNP and MNR. The aggregate arc elasticity of cycling demand for a 10% reduction in walking time is 0.27 in Gen-MNR, while it is 0.02 and 0.04 in MNP and MNR. One possible way to reduce walking travel times are technological innovations such as fast-moving walkways (e.g. Scarinci et al., 2017). Policies to support such walkways are more compelling under MNR and Gen-MNR estimates.
Furthermore, we observe that the cross-elasticities of the demand for cycling with respect to changes in driving traffic variability, transit in-vehicle travel time and transit out-of-vehicle travel time differ markedly across the three models. For example, the aggregate arc elasticity of the demand for cycling with respect to a change in transit out-of-vehicle travel time is 0.31 in MNP and Gen-MNR but is 0.12 in MNR. Interestingly, in each of these three scenarios, the elasticity estimates of MNR are quite different from the elasticity estimates of MNP and Gen-MNR. These differences are reflective of the rigidity imposed by the generic DOF parameter in MNR.

Model Walk Cycle Transit Drive
Cycling

Conclusion
In this paper, we analysed two robust alternatives to the multinomial probit (MNP) model. Both alternatives belong to the family of robit models whose kernel error distributions are heavy-tailed t-distributions that moderate the influence of outliers. The first model is the multinomial robit (MNR) model, in which a single, generic degrees of freedom parameter controls the heavy-tailedness of the kernel error distribution. The second model, the generalised multinomial robit (Gen-MNR) model, is based on the non-elliptically contoured multivariate t-distribution, which allows for distinct heavy-tailedness in each dimension of the kernel error distribution. For both models, we devised Gibbs samplers for posterior inference.
Using simulated data, we demonstrated the ability of the proposed Gibbs samplers for MNR and Gen-MNR to recover parameters in finite samples. We also showed that MNR and Gen-MNR outperform MNP at recovering choice distributions when the response data contain outliers. In a case study on transport mode choice behaviour, we demonstrated that MNR and Gen-MNR outperform MNP by considerable margins in terms of in-sample fit and out-of-sample predictive ability. We also illustrated that Gen-MNR provides the most flexible elasticity estimates.
On the whole, our analysis suggests that Gen-MNR is a useful addition to the choice modelling toolbox due to its robustness properties. In general, Gen-MNR should be preferred over the previously studied MNR model in daily data analysis practice because of its more flexible kernel error distribution.
Gen-MNR is also a parsimonious model. It requires estimating only J − 2 additional parameters compared to MNR (whereby J denotes the number of alternatives). In practice, the non-elliptical contoured t-distribution used in the formulation of Gen-MNR can also be specified in a way such that one DOF parameter controls the heavy-tailedness of more than one marginal of the kernel error distribution. Analysts can exploit this feature of Gen-MNR to achieve more parsimonious model specifications, especially in applications with many choice alternatives.
Our analysis suggests several directions for future research. First, hierarchical model extensions such as random parameters, latent variables and shrinkage priors could be considered. Furthermore, aspects of the systematic parts of the latent utility differences could be represented using Bayesian additive regression trees (BART), which automatically partition large predictor spaces to capture interaction effects and nonlinearities (Chipman et al., 2010;Kindo et al., 2016). As these extensions are rooted in the Bayesian inferential paradigm, they can be incorporated into MNR and Gen-MNR with relative ease. A second direction for future research is to use Bayesian modelling to automate aspects of the specification of MNR and Gen-MNR. For example, Burgette et al. (2020) propose a symmetric prior which obviates the specification of a base alternative for MNP kernel error covariances. Such symmetric priors could also be developed for MNR and Gen-MNR. Similarly, more diffuse prior formulations for the covariance of kernel error distributions of MNP, MNR and Gen-MNR can be investigated along the lines of Huang et al. (2013). In addition, a prior could be designed to automatically infer a parsimonious mapping of degrees of freedom parameters onto utility differences in the Gen-MNR. Finally, a third direction for future research is to formulate new robust discrete choice models using the skew-t-distribution (Kim et al., 2008;Lee and Mclachlan, 2014). The skew-t-distribution is an asymmetric probability distribution which can assume various smooth density shapes with either light or heavy tails and a longer tail either to the left or right.
The skew-t-distribution is tractable in practical applications and parsimonious in the number of parameters controlling skewness (Kim et al., 2008;Lee and Mclachlan, 2014). Hence, discrete choice models based on the skew-t-distribution would afford further improvements in robustness relative to Gen-MNR which relies on a symmetric kernel error distribution. Developing MCMC algorithms for posterior inferences in discrete choice models based on the skew-t-distribution is also a promising direction for future research since efficient data augmentation schemes could be exploited in similar ways as in the current paper.

Author contribution statement
RK: conception and design, method development and implementation, data processing and analysis, manuscript writing and editing, supervision.
MB: conception and design, manuscript editing, supervision.
TG: conception and design, method development and implementation, manuscript writing and editing.
PB: conception and design, manuscript writing and editing.

Conflicts of interest
The authors declare that they have no conflicts of interests.