On the Bayesian Two-Sample Problem for Ranking Data

Abstract

We consider the two-sample problem involving a new class of angle-based models for ranking data. These models are functions of the cosine of the angle between a ranking and a consensus vector. A Bayesian approach is employed to determine the corresponding predictive densities. Two competing hypotheses are considered, and we compute the Bayes factor to quantify the evidence provided by the observed data under each hypothesis. We apply the results to a real data set.

1. Introduction

Ranking data have been frequently collected and analyzed in a variety of different contexts, including problems of trend, market research, politics and judging situations. A comprehensive discussion of various probability ranking models may be found in Alvo and Yu [1], whereas Fagin et al. [2] considered models for the distribution of rankings for the top-k items. In the present context, scores are assigned to items to reflect their relative preference. The probability of observing a ranking is modeled to be a function of the cosine of the angle between the vector of scores and a consensus score vector denoted by $\mathit{\theta }$. A Bayesian approach is then taken wherein a joint von Mises–Fisher prior distribution is placed on the parameters. The estimation is subsequently conducted by making use of the variational inference method. The model is successfully applied to the Sushi data set.

In the present article, our attention is focused on the two-sample problem for ranking data, and the interest is in quantifying the evidence for two competing hypotheses. This is achieved by computing Bayes factors. In Section 2, we first describe the one-sample angle-based model, as well as the Bayesian MCMC approach. In Section 3, we describe the Bayesian two-sample problem and propose the use of Bayes factors to quantify the evidence provided by data in favor of one hypothesis compared to another. In Section 4, we describe an application of the two-sample model to the well-known Sushi data set. We conclude with a discussion in Section 5.

2. One-Sample Angle-Based Model

In this section, we recall the one-sample problem. A ranking for t items labeled $1,\dots ,t$, $\mathit{R}={(R\left(1\right),\dots ,R\left(t\right))}^T$ may be viewed as a permutation of the integers $1,\dots ,t$, which may be conveniently expressed in standardized form as

$\mathit{y}={\displaystyle \frac{\mathit{R}-\frac{t+1}{2}}{\sqrt{\frac{t(t^2-1)}{12}}}},$

where $\mathit{y}$ is the $t\times 1$ vector with $∥\mathit{y}∥=1$. The probability distribution of the ranking y may be modeled as follows:

(1)$p\left(\mathit{y}\right|\kappa ,\mathit{\theta })=C(\kappa ,\mathit{\theta })exp\left\{\kappa {\mathit{\theta }}^T\mathit{y}\right\}$

where the parameter $\mathit{\theta }$ may be interpreted as a $t\times 1$ “consensus” vector with $∥\mathit{\theta }∥=1$, parameter $\kappa \ge 0$, and $C(\kappa ,\mathit{\theta })$ is the normalizing constant. Since both $\mathit{\theta }$ and y are standardized geometrically, the cosine of the angle between the consensus score vector $\mathit{\theta }$ and the observation $\mathit{y}$ is given by ${\mathit{\theta }}^T\mathit{y}$. The parameter $\kappa$ can be viewed as a concentration parameter. A small value of $\kappa$ implies that the distribution of rankings is uniform, whereas a large value of $\kappa$ points to a distribution that is more concentrated around the consensus score vector. The normalizing constant is computed by summing $exp\left\{\kappa {\mathit{\theta }}^T\mathit{y}\right\}$ over all $t!$ possible permutations of the integers $1,\dots ,t$. For large values of t, this becomes prohibitive.

The model in (1) closely resembles the continuous von Mises–Fisher distribution, abbreviated as $vMF\left(\mathit{x}\right|\mathit{m},\kappa )$, which is defined on a $\left(p-1\right)$-dimensional unit sphere with mean direction $\mathit{m}$ and concentration parameter $\kappa$:

$p\left(\mathit{x}\right|\kappa ,\mathit{m})=V_p\left(\kappa \right)exp\left(\kappa {\mathit{m}}^T\mathit{x}\right),$

where

$V_p\left(\kappa \right)={\displaystyle \frac{{\kappa }^{\frac{p}{2}-1}}{{\left(2\pi \right)}^{\frac{p}{2}}{I}_{\frac{p}{2}-1}\left(\kappa \right)}},$

and $I_{{\displaystyle \frac{p}{2}}-1}\left(\kappa \right)$ is the modified Bessel function of the first kind with order ${\displaystyle \frac{p}{2}}-1.$ Consequently, as shown in Xu et al. [3], we may approximate the normalizing constant in the discrete model by

$C(\kappa ,\mathit{\theta })\simeq C_t\left(\kappa \right)={\displaystyle \frac{{\kappa }^{\frac{t-3}{2}}}{2^{\frac{t-3}{2}}t!I_{\frac{t-3}{2}}\left(\kappa \right)\Gamma \left(\frac{t-1}{2}\right)}},$

where $\Gamma (·)$ is the gamma function.

2.1. Maximum Likelihood Estimation (MLE)

The estimation of the parameters $\left(\kappa ,\mathit{\theta }\right)$ in (1) proceeds by making use of the method of maximum likelihood. Let $\mathit{Y}=\left\{{\mathit{y}}_1,\dots ,{\mathit{y}}_N\right\}$ be a random sample of N standardized rankings drawn from $p\left(\mathit{y}\right|\kappa ,\mathit{\theta })$. Then, the log-likelihood of $\left(\kappa ,\mathit{\theta }\right)$ is given by

(2)$L\left(\mathit{y}\right|\kappa ,\mathit{\theta })=Nln{C}_t\left(\kappa \right)+\sum _{i=1}^N\kappa {\mathit{\theta }}^T{\mathit{y}}_i.$

The maximum likelihood estimator of $\mathit{\theta }$ subject to $∥\mathit{\theta }∥=1$ and $\kappa \ge 0$ is given by

${\widehat{\mathit{\theta }}}_{MLE}={\displaystyle \frac{{\sum }_{i=1}^N{\mathit{y}}_i}{∥{\sum }_{i=1}^N{\mathit{y}}_i∥}},$

whereas $\widehat{\kappa }$ is calculated from the relation

$A_t\left(\kappa \right)\equiv {\displaystyle \frac{-C_t^{\prime }\left(\kappa \right)}{C_t\left(\kappa \right)}}={\displaystyle \frac{I_{\frac{t-1}{2}}\left(\kappa \right)}{I_{\frac{t-3}{2}}\left(\kappa \right)}}={\displaystyle \frac{∥{\sum }_{i=1}^N{\mathit{y}}_i∥}{N}}\equiv r$

Whereas Banerjee et al. [4] proposed

${\widehat{\kappa }}_{MLE}={\displaystyle \frac{r(t-1-r^2)}{1-r^2}}$

a more precise approximation can be obtained recursively as

${\kappa }_{i+1}={\kappa }_i-{\displaystyle \frac{A_t\left({\kappa }_i\right)-r}{1-A_t{\left({\kappa }_i\right)}^2-\frac{t-2}{{\kappa }_i}{A}_t\left({\kappa }_i\right)}},i=0,1,2,\dots .$

2.2. Bayesian Method with Conjugate Prior

Define a conjugate prior for $(\kappa ,\mathit{\theta })$ as follows:

$\begin{array}{cc}\hfill p(\kappa ,\mathit{\theta })& =p\left(\mathit{\theta }\right|\kappa \left)p\right(\kappa )\hfill \\ \hfill & =vMF\left(\mathit{\theta }\right|{\mathit{m}}_0,{\beta }_0\kappa )G\left(\kappa \right|a_0,b_0),\hfill \end{array}$

where $G\left(\kappa \right|a_0,b_0)$ is the gamma density function with shape parameter $a_0$ and rate parameter $b_0$ and a mean equal to ${\displaystyle \frac{a_0}{b_0}}$. Given $\mathit{Y}\mathbf{=}\mathit{y}$, the posterior density of $(\kappa ,\mathit{\theta })$ can be expressed by

$p(\kappa ,\mathit{\theta }|\mathit{y})\propto exp\left\{\beta \kappa {\mathit{m}}^T\mathit{\theta }\right\}{V}_t\left(\beta \kappa \right){\displaystyle \frac{{\left[C_t\left(\kappa \right)\right]}^{N+{\nu }_0}}{V_t\left(\beta \kappa \right)}},$

where $\mathit{m}=\left({\beta }_0{\mathit{m}}_{\mathbf{0}}+{\sum }_{i=1}^N{\mathit{y}}_i\right){\beta }^{-1},$ $\beta =∥{\beta }_0{\mathit{m}}_0+{\sum }_{i=1}^N{\mathit{y}}_i∥$. The posterior density can be factorized as

(3)$p(\kappa ,\mathit{\theta }|\mathit{y})=p(\mathit{\theta }|\kappa ,\mathit{y})p\left(\kappa \right|\mathit{y}),$

where $p\left(\mathit{\theta }\right|\kappa ,\mathit{y})\sim vMF(\mathit{\theta }|\mathit{m},\beta \kappa )$ and

$\begin{array}{ccc}\hfill p\left(\kappa \right|\mathit{y})& \propto & {\displaystyle \frac{{\left[C_t\left(\kappa \right)\right]}^{N+{\nu }_0}}{V_t\left(\beta \kappa \right)}}\hfill \\ & =& {\displaystyle \frac{{\kappa }^{\frac{t-3}{2}({\upsilon }_0+N)}{I}_{\frac{t-2}{2}}\left(\beta \kappa \right)}{{\left[I_{\frac{t-3}{2}}\left(\kappa \right)\right]}^{{\nu }_0+N}{\left(\beta \kappa \right)}^{\frac{t-2}{2}}}}\hfill \end{array}$

It is not possible to compute the normalizing constant for $p\left(\kappa \right|Y\left)\right)$ in closed form. Nunez-Antonio and Gutiérrez-Pena [5] exploited a sampling importance resampling (SIR) procedure with a proposal density chosen to be the gamma density with variance equal to some pre-specified number such as 50 or 100. However, this choice for the variance is crucially related to the performance of SIR. An improper choice of variance may lead to slow or unsuccessful convergence. Moreover, the MCMC method led to intensive computational complexity. Furthermore, $\beta \kappa$ can be very large if N is large, which complicates the computation of the term $V_t\left(\beta \kappa \right)$. Thus, it would not be possible to calculate the weights in the SIR method when N is large. As a result, the posterior distribution was approximated using the method of variational inference (abbreviated as VI henceforth).

As shown in (Xu et al. [3]), $p\left(\kappa \right|Y\left)\right)$ can be approximated by $Gamma\left(\kappa \right|a,b)$ with shape a and rate b, where the posterior mode $\overline{\kappa }$ is

$\overline{\kappa }=\left\{\begin{array}{cc}{\displaystyle \frac{a-1}{b}}\hfill & \mathrm{if}a>1,\hfill \\ {\displaystyle \frac{a}{b}}\hfill & \mathrm{otherwise}\hfill \end{array}\right.$

and

$a=a_0+N\left({\displaystyle \frac{t-3}{2}}\right)+\beta \overline{\kappa }\left[{\displaystyle \frac{\partial }{\partial \beta \kappa }}ln{I}_{{\displaystyle \frac{t-2}{2}}}\left(\beta \overline{\kappa }\right)\right],$

$b=b_0+N{\displaystyle \frac{\partial }{\partial \kappa }}{I}_{{\displaystyle \frac{t-3}{2}}}\left(\overline{\kappa }\right)+{\beta }_0\left[{\displaystyle \frac{\partial }{\partial {\beta }_0\kappa }}ln{I}_{{\displaystyle \frac{t-2}{2}}}\left({\beta }_0\overline{\kappa }\right)\right].$

Here, we approximated the posterior distribution by making use of variational inference, which has been used in many applications. In the variational inference approach, we first propose a candidate family of densities and then select the member of that family that comes closest to the target posterior density in the Kullback–Leibler sense. Specifically, let $q\left(\mathbf{Z}\right)$ represent the candidate family and suppose that $p\left(\mathbf{Z}|\mathit{Y}\right)$ represents the target posterior density. Then, the Kullback–Leibler divergence is given by

$KL\left(q|p\right)=E_q\left[ln{\displaystyle \frac{q\left(\mathbf{Z}\right)}{p\left(Z|\mathit{Y}\right)}}\right].$

For further background on Bayesian methods, see Gelman et al. [6].

3. Two-Sample Problem

In the two-sample problem, there are two independent random samples ${\mathit{Y}}_i=\left\{{\mathit{y}}_{i1},\dots ,{\mathit{y}}_{i{N}_i}\right\}$ for $i=1,2$ of standardized rankings each drawn, respectively, from $p\left({\mathit{y}}_{\mathit{i}}\right|\kappa ,{\mathit{\theta }}_{\mathit{i}}).$ Conditional on $\kappa$, we assume that there are independent von Mises conjugate priors, respectively, for $({\mathit{\theta }}_1,{\mathit{\theta }}_2)$ as

$p\left({\mathit{\theta }}_{\mathit{i}}\right|\kappa )\propto {\left[C_t\left(\kappa \right)\right]}^{{\nu }_0}exp\left\{\kappa {\mathit{m}}_0^T{\mathit{\theta }}_{\mathit{i}}\right\},$

where $∥{\mathit{m}}_0∥=1$, ${\nu }_0,\ge 0$. To motivate interest in the two-sample problem, we shall investigate the difference in food preference patterns with respect to Sushi between Eastern and Western Japan. Individual models are first fitted to each sample. Then, we are interested in comparing the parameters ${\mathit{\theta }}_{\mathbf{1}},{\mathit{\theta }}_{\mathbf{2}}$.

The Bayes factor is a statistical measure used in Bayesian analysis to quantify the evidence in favor of one hypothesis compared to another. It compares the likelihood of the observed data under two competing hypotheses, typically the null hypothesis $H_0$ vs. the alternative hypothesis $H_1$. The Bayes factor is defined as the ratio of the likelihoods under the alternative and the null hypothesis, respectively:

$BF={\displaystyle \frac{P(Data\mid H_1)}{P(Data\mid H_0)}}$

Bayes factors are widely used for model selection, hypothesis testing, and decision-making within the Bayesian framework; see Kass and Raftery [7]. They provide a ratio of how much more (or less) likely the data are under one hypothesis compared to the other. Unlike p-values, they provide a continuous measure of evidence and allow for a direct comparison between hypotheses that incorporate prior beliefs. They penalize overfitting by integrating over parameters, thereby discouraging overly complex models.

Bayes factors have diverse applications across various domains in statistics, science, and decision-making—for example, in testing the effectiveness of a new drug compared to a placebo; in model selection; to evaluate multiple competing models in psychology and cognitive science; in forensic science, genetics, and bioinformatics; in econometrics; and, finally, in artificial intelligence and machine learning. By quantifying evidence in a probabilistic manner, Bayes factors offer an intuitive and rigorous alternative to approaches that make use of p-values. A somewhat practical interpretation scale is as follows.

Criterion	Interpretation
1 ≤ BF < 3	Weak evidence for $H_1$
3 ≤ BF < 10	Moderate evidence for $H_1$
BF ≥ 10	Strong evidence for $H_1$
BF < 1/3	Moderate evidence for $H_0$
BF < 1/10	Strong evidence for $H_0$

The Bayes factor was computed in order to compare the following two models. Under model 1, denoted $M_1,$${\mathit{\theta }}_{\mathbf{1}}={\mathit{\theta }}_{\mathbf{2}}$, whereas, under model 2, denoted $M_2$, equality is not assumed. We are concerned with the conditional Bayes factor given $\kappa$ between the two models. The Bayes factor enables us to compute the posterior odds of $M_2$ to $M_1$. It is defined as the ratio of the predictive densities of $M_2$ to $M_1$, namely

$\begin{array}{ccc}\hfill B_{21}\left(\kappa \right)& =& {\displaystyle \frac{\int p\left({\mathit{y}}_{\mathbf{1}}\right|\kappa ,{\mathit{\theta }}_{\mathbf{1}})p\left({\mathit{y}}_{\mathbf{2}}\right|\kappa ,{\mathit{\theta }}_{\mathbf{2}})p({\mathit{\theta }}_{\mathbf{1}},{\mathit{\theta }}_{\mathbf{2}}\left|\kappa \right|y_1{y}_2)d{\mathit{\theta }}_{\mathbf{1}}d{\mathit{\theta }}_{\mathbf{2}}}{\int p\left({\mathit{y}}_{\mathbf{1}}\right|\kappa ,\mathit{\theta })p\left({\mathit{y}}_{\mathbf{2}}\right|\kappa ,\mathit{\theta })p\left(\mathit{\theta }\right|\kappa ,y_1{y}_2)d\mathit{\theta }}}\hfill \end{array}$

Consider

$\begin{array}{ccc}\hfill \int p\left({\mathit{y}}_{\mathbf{1}}\right|\kappa ,{\mathit{\theta }}_{\mathbf{1}})p({\mathit{\theta }}_{\mathbf{1}}\left|\kappa \right|y_1)d{\mathit{\theta }}_{\mathbf{1}}& =& \int C_t^{N_1}\left(\kappa \right)V_t\left(\kappa \right)exp\left\{\kappa {\mathit{\theta }}_1^T\left(\sum _i{y}_{1i}+{\beta }_0{m}_0\right)\right\}d{\mathit{\theta }}_{\mathbf{1}}\hfill \\ \\ & =& C_t^{N_1}\left(\kappa \right)V_t\left(\kappa \right)V_t^{-1}\left(\kappa {\beta }_1\right)\hfill \end{array}$

where

${\beta }_1=∥\sum _i{y}_{1i}+{\beta }_0{m}_0∥$

Similarly,

$\begin{array}{ccc}\hfill \int p\left({\mathit{y}}_{\mathbf{1}}\right|\kappa ,{\mathit{\theta }}_{\mathbf{2}})p({\mathit{\theta }}_{\mathbf{2}}\left|\kappa \right|Y_2)d{\mathit{\theta }}_{\mathbf{2}}& =& C_t^{N_2}\left(\kappa \right)V_t\left(\kappa \right)V_t^{-1}\left(\kappa {\beta }_2\right)\hfill \end{array}$

where

${\beta }_2=∥\sum _i{y}_{2i}+{\beta }_0{m}_0∥$

Consequently, the numerator under model $M_2$ is given by

$\left[C_t^{N_1}\left(\kappa \right)V_t\left(\kappa \right)V_t^{-1}\left(\kappa {\beta }_1\right)\right]\left[C_t^{N_2}\left(\kappa \right)V_t\left(\kappa \right)V_t^{-1}\left(\kappa {\beta }_2\right)\right]$

As for the denominator in $B_{21}$,

$\begin{array}{ccc}\hfill \int p\left({\mathit{y}}_{\mathbf{1}}\right|\kappa ,\mathit{\theta })p\left({\mathit{y}}_{\mathbf{2}}\right|\kappa ,\mathit{\theta })p\left(\mathit{\theta }\right|\kappa \left|\right)d\mathit{\theta }& =& \int C_t^{N_1+N_2}\left(\kappa \right)V_t\left(\kappa \right)exp\left\{\kappa {\mathit{\theta }}^T\left(\sum _i{y}_{1i}+\sum _i{y}_{2i}+{\beta }_0{m}_0\right)\right\}d\mathit{\theta }\hfill \\ \\ & =& C_t^{N_1+N_2}\left(\kappa \right)V_t\left(\kappa \right)V_t^{-1}\left(\kappa \beta \right)\hfill \end{array}$

where

$\beta =∥\sum _i{y}_{1i}+\sum _i{y}_{2i}+{\beta }_0{m}_0∥$

It follows that the conditional Bayes factor given $\kappa$ is the ratio

$\begin{array}{ccc}\hfill B_{21}\left(\kappa \right)& =& {\displaystyle \frac{C_t^{N_1}\left(\kappa \right)V_t\left(\kappa \right)V_t^{-1}\left(\kappa {\beta }_1\right)C_t^{N_2}\left(\kappa \right)V_t\left(\kappa \right)V_t^{-1}\left(\kappa {\beta }_2\right)}{C_t^{N_1+N_2}\left(\kappa \right)V_t\left(\kappa \right)V_t^{-1}\left(\kappa \beta \right)}}\hfill \\ \\ & =& {\displaystyle \frac{V_t\left(\kappa \beta \right)V_t\left(\kappa \right)}{V_t\left(\kappa {\beta }_1\right)V_t\left(\kappa {\beta }_2\right)}}\hfill \\ \\ & =& {\displaystyle \frac{\frac{{\left(\beta \kappa \right)}^{\frac{p}{2}-1}}{{\left(2\pi \right)}^{\frac{p}{2}}{I}_{\frac{p}{2}-1}\left(\beta \kappa \right)}\frac{{\kappa }^{\frac{p}{2}-1}}{{\left(2\pi \right)}^{\frac{p}{2}}{I}_{\frac{p}{2}-1}\left(\kappa \right)}}{\frac{{\left({\beta }_1\kappa \right)}^{\frac{p}{2}-1}}{{\left(2\pi \right)}^{\frac{p}{2}}{I}_{\frac{p}{2}-1}\left({\beta }_1\kappa \right)}\frac{{\left({\beta }_2\kappa \right)}^{\frac{p}{2}-1}}{{\left(2\pi \right)}^{\frac{p}{2}}{I}_{\frac{p}{2}-1}\left({\beta }_2\kappa \right)}}}\hfill \\ \\ & =& {\left({\displaystyle \frac{\beta }{{\beta }_1{\beta }_2}}\right)}^{{\displaystyle \frac{p}{2}}-1}{\displaystyle \frac{I_{\frac{p}{2}-1}\left({\beta }_1\kappa \right)I_{\frac{p}{2}-1}\left({\beta }_2\kappa \right)}{I_{\frac{p}{2}-1}\left(\beta \kappa \right)I_{\frac{p}{2}-1}\left(\kappa \right)}}\hfill \end{array}$

where $p=t-1.$

There is an approximation to the Bessel functions for a large x:

$I_{{\displaystyle \frac{p}{2}}-1}\left(x\right)\simeq {\displaystyle \frac{e^x}{\sqrt{2\pi x}}}$

In this case, for a large $\kappa$, the conditional Bayes factor is approximately

$B_{21}\left(\kappa \right)={\left({\displaystyle \frac{\beta }{{\beta }_1{\beta }_2}}\right)}^{{\displaystyle \frac{p}{2}}-1}{\displaystyle \frac{I_{\frac{p}{2}-1}\left({\beta }_1\kappa \right)I_{\frac{p}{2}-1}\left({\beta }_2\kappa \right)}{I_{\frac{p}{2}-1}\left(\beta \kappa \right)I_{\frac{p}{2}-1}\left(\kappa \right)}}\simeq {\left({\displaystyle \frac{\beta }{{\beta }_1{\beta }_2}}\right)}^{{\displaystyle \frac{p-1}{2}}}{e}^{\left({\beta }_1+{\beta }_2-\beta -1\right)\kappa }$

4. Application

Using the mixture ranking model, we were interested in uncovering differences in the food preference patterns between Eastern and Western Japan. We made use of the two Sushi data sets of Kamishima and Akaho [8]. Historically, Western Japan has been mainly affected by the culture of the Mikado emperor and nobles, while Eastern Japan was the home of the Shogun and Samurai warriors. Therefore, the preference patterns in food are different between these two regions [8].

We made use of the complete rankings of $t=10$ different types of sushi given by 5000 respondents according to their preferences as two samples distinguished by region. There were $N_1=3285$ respondents from Eastern Japan and $N_2=1715$ from Western Japan. We applied the approach described earlier to both the Eastern and Western Japan data. The sample sizes are quite large compared to t, and, consequently, the estimated models for all three methods are quite similar. In Figure 1, we present the posterior means of $\mathit{\theta }$ between Eastern Japan (blue bar) and Western Japan (red bar). A more negative value of ${\mathit{\theta }}_i$ implies that sushi i is more preferred. Eastern and Western Japan differ mostly in their preferences for salmon roe, squid, sea eel, shrimp, and tuna. Respondents from Eastern Japan prefer salmon roe and tuna more than those from Western Japan, who exhibit a greater preference for squid, shrimp, and sea eel. From Table 1, we see the larger posterior mean of $\kappa$, which shows that the Eastern Japanese are slightly more cohesive in their preferences than the Western Japanese in their preferences.

The boxplots of the top smallest posterior means are displayed in Figure 2. Eastern and Western Japan appear to differ mostly in their preferences for sea eel, salmon roe, tuna, sea urchin, and sea bream. The Eastern Japanese exhibit a greater preference for salmon roe, tuna, and sea urchin sushi, whereas the Western Japanese like sea eel and sea bream. Generally, tuna and sea urchin are oilier foods, whereas salmon roe and tuna are more seasonal. So, from the analysis of both data sets, we can conclude that the Eastern Japanese prefer more oily and seasonal food compared to the Western Japanese [8]. The boxplots in Figure 2 also provide a measure of the uncertainty in the values of the ${\mathit{\theta }}_{i}s.$

We computed the Bayes factor (3) using ${\beta }_0=0,\beta =2209.22,\kappa =4.7.$ We have that

$B_{21}\left(\kappa \right)\cong 3.6\times {10}^{-11}$

This leads to favoring the $M_1$ model, namely that ${\mathit{\theta }}_1={\mathit{\theta }}_2$. The calculation of the Bayes factor comparing the two models supports the conclusion reached in the original paper.

5. Conclusions and Discussion

We considered the two-sample ranking data problem using angle-based ranking models. Each model is itself a function of a consensus score vector $\mathit{\theta }$, which exhibits detailed information on item preferences, unlike a distance-based model, which only provides an equally spaced modal ranking. We proceeded with a Bayesian approach and placed a conjugate prior on the parameters. We then postulated two different hypotheses and proceeded with the calculation of Bayes factors to compare them. Unlike p-values, Bayes factors provide a continuous measure of evidence between hypotheses that incorporate prior beliefs. We made use of Bayesian variational inference to approximate both the posterior density and the predictive density. We applied our results to the popular sushi data and compared the food preferences of Western and Eastern Japanese consumers.

In further research, it would be of interest to incorporate additional arguments or covariates into the model. These could include judge-specific covariates such as age, gender, and income, whereas item-specific covariates could be prices, weights, and brands. The combination of judge–item covariates could include personal experience on using each phone or brand. The addition of these covariates in the model may greatly enhance the predictive power of the model. Finally, research could explore score functions based on the Kendall distance to enhance ranking comparisons.

Data Availability Statement

Data is contained within the article.

Acknowledgments

The author would like to thank the Natural Sciences and Engineering Research Council of Canada OGP0009068.

Conflicts of Interest

The author declares no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1 Posterior means of $\mathit{\theta }$ for the sushi complete ranking data ($t=10$) in Eastern Japan (blue bar) and Western Japan (red bar) obtained by Bayesian-VI.

Figure 2 Boxplots of the top 10 smallest posterior means of $\mathit{\theta }$ for the sushi incomplete rankings ($t=100$). Eastern Japan is shown in blue, where circles indicate outliers. Western Japan is shown in red, where red pluses indicate outliers.

Table 1

Posterior parameters for the sushi complete ranking data ($t=10$) in Eastern Japan and Western Japan obtained by Bayesian-VI.

Posterior Parameter	Eastern Japan	Western Japan
$\beta$	1458.85	741.61
a	18,509.84	9462.70
b	3801.57	2087.37
Posterior Mean of $\kappa$	4.87	4.53

References

1. Alvo, M.; Yu, P. Statistical Methods for Ranking Data; Springer: New Yrok, NY, USA, 2014.

2. Fagin, R.; Kumar, R.; Sivakumar, D. Comparing top k lists. SIAM J. Discret.; 2003; 17, pp. 134-160. [DOI: https://dx.doi.org/10.1137/S0895480102412856]

3. Xu, H.; Alvo, M.; Yu, P. Angle-based models for ranking data. Comput. Stat. Data Anal.; 2018; 121, pp. 113-136.

4. Banerjee, A.; Dhillon, I.; Ghosh, J.; Sra, S. Clustering on the unit hypersphere using von Mises-Fisher distributions. J. Mach. Learn. Res.; 2005; 6, pp. 1345-1382.

5. Nunez-Antonio, G.; Gutiérrez-Pena, E. A bayesian analysis of directional data using the von misesfisher distribution. Commun.-Stat.-Simul. Comput.; 2005; 34, pp. 989-999. [DOI: https://dx.doi.org/10.1080/03610910500308495]

6. Gelman, A.; Carlin, J.B.; Stern, H.S.; Rubin, D.B. Bayesian Data Analysis; CRC Press: Boca Raton, FL, USA, 2013.

7. Kass, R.E.; Raftery, A.E. Bayes factors. J. Am. Stat. Assoc.; 1995; 90, pp. 773-795.

8. Kamishima, T.; Akaho, S. Efficient clustering for orders. Proceedings of the 2nd International Workshop on Mining Complex Data; Hong Kong, China, 18–22 December 2006; pp. 274-278.