2021

These authors contributed equally to this work.

These authors contributed equally to this work.

1] DG of Economics, Statistics and Research, Bank of Italy, Rome, Italy

2]Department of Statistical Sciences, Sapienza University of Rome, Italy

# Efficient and accurate inference for mixtures of Mallows models with
Spearman distance

Marta Crispino<EMAIL_ADDRESS>Cristina Mollica
<EMAIL_ADDRESS>Valerio Astuti<EMAIL_ADDRESS>Luca Tardella<EMAIL_ADDRESS>[  [

###### Abstract

The Mallows model occupies a central role in parametric modelling of ranking
data to learn preferences of a population of judges. Despite the wide range of
metrics for rankings that can be considered in the model specification, the
choice is typically limited to the Kendall, Cayley or Hamming distances, due
to the closed-form expression of the related model normalizing constant. This
work instead focuses on the Mallows model with Spearman distance.

An efficient and accurate EM algorithm for estimating finite mixtures of
Mallows models with Spearman distance is developed, by relying on a twofold
data augmentation strategy aimed at i) enlarging the applicability of Mallows
models to samples drawn from heterogeneous populations; ii) dealing with
partial rankings affected by diverse forms of censoring. Additionally, a novel
approximation of the model normalizing constant is introduced to support the
challenging model-based clustering of rankings with a large number of items.
The inferential ability of the EM scheme and the effectiveness of the
approximation are assessed by extensive simulation studies. Finally, we show
that the application to three real-world datasets endorses our proposals also
in the comparison with competing mixtures of ranking models.

###### keywords:

Ranking data, distance-based models, model-based clustering, EM algorithm,
censoring

##  1 Introduction

Ranking data arise when a sample of $N$ people is presented a finite set of
$n$ alternatives, called items, and is asked to rank them according to a
certain criterion, such as personal preferences or attitudes. Thus, a ranking
is the result of a comparative judgment on a set of competing alternatives
expressed in the form of order relation.

The increasing interest in ranked data analysis is motivated by several
research contexts, among which marketing and political surveys, where items
could be consumer goods, political candidates or goals, but also by
psychological and behavioral studies consisting, for instance, in ordering of
words/topics according to the perceived association with a reference subject
(Gormley and Murphy, 2008a; Caron et al., 2014). Another typical context for
ranking data is sports and, more generally, competitions. Some examples are
horse or car races, soccer and basket championships, where the rankings are
determined by the relative ability of the competing players or teams (Henery,
1981).

However, ranking data are not necessarily the result of a human comparative
choice where, typically, only a limited number of options can be submitted to
the judgment of the rankers. The need of handling large ordered lists of
alternatives can be motivated, for example, with the analysis of the output
obtained from search algorithms. This can include the results of web-search
engines, as well as of programs for matching biological sequences or specific
traits with a reference database, that return a list of items ordered
according to the degree of agreement with the input information (DeConde et
al., 2006). A further field of application concerns the widespread digital
platforms of on demand services, such as movies, songs and any kind of product
available from the online marketplaces. In this context, various forms of
feedback from the customers (e.g. clicks, ratings, reviews and purchases) are
elaborated from the e-commerce companies to capture customers’ preferences and
supply them with a personalized ordered list of recommended items (Liu et al.,
2019). All these applications are characterized by the challenging task of
dealing with a large value $n$ of items which are often not completely ranked,
giving rise to partial rankings.

In the wide literature on probability distributions for ranking data, see
Marden (1995) and Alvo and Yu (2014) for extensive reviews, the class of
Mallows models (MM) became very popular for its long history, both as a
successful approach for preference learning and as a starting point for
effective methodological advances. Formally, the MM assumes the existence of a
unique modal ranking representing the common comparative judgment on the
items, usually referred to as consensus ranking. Moreover, the variability of
the population of rankers, that is, the agreement with the consensus ranking,
is calibrated by a concentration parameter. A third key ingredient of the MM
is the distance between rankings, since the probability of a generic sequence
depends on its distance from the consensus; as a consequence the MM is also
referred to as distance-based model (DBM).

Even if these models are intuitive and attractive, especially for the fact
that they represent exponential families for random permutations (Diaconis,
1988), postulating a MM as a generative process of the observed rankings does
not prevent noticeable practical issues. First, one of its peculiarities is
the mixed-type parameter space, combining the continuous concentration
parameter and the discrete consensus ranking. The latter, in fact, takes
values in the space $\mathcal{P}_{n}$ of the $n!$ permutations of the first
$n$ integers, whose adequate exploration, needed for making inference on the
model, is computationally demanding and challenging as the number of items
increases. Moreover, inference is difficult to perform in those cases when the
model normalizing constant is intractable. Second, the estimation task is
further complicated in the presence of partial observed rankings and of
heterogeneity in the sample. Finally, some criticism was raised about the low
flexibility of the MM (Mollica and Tardella, 2014), which can limit its
adequacy in describing many real-world cases, especially when the sample size
is high.

For the above reasons, in the last decades various directions to generalize
the MM and improve its fitting performance have been proposed, among which the
mixture approach plays a leading role. Nevertheless, a review of the existing
MM extensions reveals that the version of the MM with Spearman distance (MMS
henceforth) received very little attention, even if this metric induces a
probability distribution on the permutation space that parallels the Gaussian
model and, similarly, enjoys convenient analytical simplifications. The main
reason for this lack of attention lies in the fact that the MMS is
characterized by an intractable normalizing constant, while the MM with other
distances (such as the Kendall and the Cayley) enjoys favourable closed-form
expressions (Fligner and Verducci, 1986; Irurozki, 2014).

This work focuses on the MMS, originally introduced by Mallows (1957) under
the name $\theta$-model, and its extension into the finite mixture framework
for the analysis of full and partial rankings. Specifically, we exploit the
crucial result that, unlike other DBMs, the MMS admits an analytical solution
for the critical estimation of the consensus ranking to generalize the MMS
within the finite mixture approach and perform model-based clustering of
arbitrary forms of partial rankings from the maximum likelihood estimation
(MLE) perspective. Moreover, we propose an original approximation of the MMS
normalizing constant, aiming at effectively addressing the computational
issues arising when inference is conducted on rankings of a large number $n$
of items. This approximation, which may be interesting in its own right (for
instance for sampling procedures), is then compared with existing
approximating approaches. Inference is addressed with a new Expectation-
Maximization (EM) algorithm which is much more accurate and efficient than
existing versions, and allows to successfully investigate the group structure
of large samples of rankings with a high number of alternatives. The merits of
our inferential approach are illustrated with applications to simulated and
real-world ranking data.

The outline of the paper is the following. In Section 2, after reviewing the
MM and its main generalizations, we focus on the inference on the MMS. We then
move to detailing the finite mixture framework, also sketching the EM
algorithm, and finally we explain the data augmentation strategy to deal with
heterogeneous partial rankings. In Section 3, we present the new approximation
of the model normalizing constant and compare it with available alternatives.
An extensive simulation study is performed in Section 4, by considering a wide
range of scenarios in terms of sample size, number of items, heterogeneity of
the sampling distribution, as well as different patterns of data missingness.
Section 5 deals with real-world data illustrations, then the paper ends with
concluding remarks in Section 6.

##  2 Model and inference

###  2.1 The Mallows model: a brief review

Let $\boldsymbol{r}=(r_{1},\dots,r_{n})$ be a ranking of $n$ items with
generic entry $r_{i}$ indicating the rank assigned to item $i$. We adopt the
usual convention that $r_{i}<r_{i^{\prime}}$ means that item $i$ is preferred
to item $i^{\prime}$ (the lower the rank is, the more preferred the item is).
Both items and ranks are identified with the set $\\{1,\dots,n\\}$, implying
that a generic observation $\boldsymbol{r}$ is a permutation of the first $n$
integers and, thus, a point in the finite discrete space $\mathcal{P}_{n}$.

The MM is a well-established parametric class of ranking distributions whose
paternity is attributed to Mallows (1957). The probability of observing the
ranking $\boldsymbol{r}$ under the MM is

$\mathbb{P}(\boldsymbol{r}\,|\boldsymbol{\rho},\theta)=\frac{e^{-\theta\,d(\boldsymbol{r},\boldsymbol{\rho})}}{Z(\theta,\boldsymbol{\rho})}\qquad\boldsymbol{r}\in\mathcal{P}_{n},$
(1)

where $\boldsymbol{\rho}\in\mathcal{P}_{n}$ is the consensus ranking,
$\theta\in\mathbb{R}^{+}_{0}$ is the concentration parameter, $d(\cdot,\cdot)$
is a distance over $\mathcal{P}_{n}$ and
$Z(\theta,\boldsymbol{\rho})=\sum_{\boldsymbol{r}\in\mathcal{P}_{n}}e^{-\theta\,d(\boldsymbol{r},\boldsymbol{\rho})}$
is the normalizing constant or partition function. From (1) it follows that,
under the MM assumption, the probability of observing a ranking decreases
exponentially as its distance from the consensus $\boldsymbol{\rho}$
increases.

In the attempt to enlarge the flexibility of the MM and improve the goodness-
of-fit, various generalizations appeared in the literature. Fligner and
Verducci (1986) proposed the Generalized Mallows model (GMM), which considers
a $(n-1)$-dimensional vector of concentration parameters, each affecting a
particular position in the ranking. The GMM, however, can be defined only when
the distance is decomposable into the sum of independent components associated
to the each single stage of the ranking process, such as for the Kendall,
Cayley and Hamming metrics. Meilǎ and Bao (2010) extended the GMM to infinite
rankings with the Infinite Generalized Mallows model, to handle those cases
where the number of items is very high, or potentially not completely known,
such as in the output of search engines. Another proposal to enhance the
flexibility of the MM is the weighted distance-based model (WDBM) by Lee and
Yu (2010), which amounts to use weighted versions of the commonly used
distances between rankings, including the Spearman one.

In order to explore the unobserved heterogeneity in the sample, the
aforementioned extensions were considered also in the finite mixture modelling
approach. For example, Murphy and Martin (2003) described the general EM
algorithm to derive the MLE of MM mixtures for an arbitrary metric; Lee and Yu
(2012) detailed the same inferential method for the mixture of WDBMs.
Differently from the one illustrated by Murphy and Martin (2003), Busse et al.
(2007) introduce an EM algorithm for the MM mixtures, working with Kendall
distance only but handling both complete and partial rankings of a larger
number of items.

However, all the above contributions do not specifically focus on the MMS and,
if they handle the Spearman distance, they only work in practice with few
items (less than $n=10$, say). This work, instead, is devoted to the MMS and
its mixture extension, with the aim of developing a learning strategy able to
handle also rankings of a large number of items.

###  2.2 Likelihood and inference for the MMS

The Spearman distance between two permutations
$\boldsymbol{r},\boldsymbol{\rho}\in\mathcal{P}_{n}$ is defined as follows

$d(\boldsymbol{r},\boldsymbol{\rho})=\sum_{i=1}^{n}(r_{i}-\rho_{i})^{2}.$ (2)

It is an unnormalized version of the Spearman rank correlation, used to
measure the statistical correlation between the ranks of two variables but,
when rankings are considered as vectors in $\mathbb{R}^{n}$, it is simply the
squared Euclidean distance, or $L_{2}$-norm. By developing the square in (2)
and setting $c_{n}=n(n+1)(2n+1)/6$ and the scalar product
$\boldsymbol{\rho}^{T}\boldsymbol{r}=\sum_{i=1}^{n}\rho_{i}r_{i}$, where the
symbol T denotes the transposition (row vector), one has
$d(\boldsymbol{r},\boldsymbol{\rho})=2\left(c_{n}-\boldsymbol{\rho}^{T}\boldsymbol{r}\right)$.
Hence, the MMS can be written as

$\mathbb{P}(\boldsymbol{r}\,|\boldsymbol{\rho},\theta)=\frac{e^{-2\theta\,\left(c_{n}-\boldsymbol{\rho}^{T}\boldsymbol{r}\right)}}{Z(\theta)},$
(3)

where the normalizing constant
$Z(\theta)=\sum_{\boldsymbol{r}\in\mathcal{P}_{n}}e^{-2\theta\,\left(c_{n}-\boldsymbol{e}^{T}\boldsymbol{r}\right)}$
with $\boldsymbol{e}=(1,2,...,n)$ does not depend on $\boldsymbol{\rho}$
because of the right-invariance property of the Spearman distance (Diaconis,
1988). Note that the number of terms in the sum grows with $n!$, which makes
direct calculation of $Z(\theta)$ unfeasible for all but very small values of
$n$, typically for at most $n=11$.

Let $\underline{\boldsymbol{r}}=(\boldsymbol{r}_{1},\dots,\boldsymbol{r}_{N})$
be a random sample of size $N$ drawn from the MMS and $N_{l}$ be the frequency
of the $l$-th distinct observed ranking $\boldsymbol{r}_{l}$ for
$l=1,\dots,L$, such that $\sum_{l=1}^{L}N_{l}=N$. The observed-data log-
likelihood of the MMS is

$\begin{split}\ell(\boldsymbol{\rho},\theta|\underline{\boldsymbol{r}})=-N\left[\log{Z(\theta)}+2\theta\left(c_{n}-\boldsymbol{\rho}^{T}{\boldsymbol{\bar{r}}}\right)\right],\end{split}$
(4)

where ${\boldsymbol{\bar{r}}}=(\bar{r}_{1},\ldots,\bar{r}_{n})$ is the sample
mean rank vector whose $i$-th entry is
$\bar{r}_{i}=\frac{1}{N}\sum_{l=1}^{L}N_{l}r_{li}$. When considering the
metric (2) in the MM, it is not surprising that the likelihood of the MMS
resembles that of an $n$-variate normal, except for its finite support. In
other words, the MMS is the restriction of the $n$-dimensional Gaussian on
$\mathcal{P}_{n}$ (Crispino and Antoniano-Villalobos, 2022). It follows that,
if $\boldsymbol{\bar{r}}\in\mathcal{P}_{n}$, then the MLE
$\hat{\boldsymbol{\rho}}$ of the consensus ranking simply coincides with the
sample mean rank vector. However, in general
$\boldsymbol{\bar{r}}\notin\mathcal{P}_{n}$, so a further consideration is
required in order to solve the optimization problem. Specifically, the log-
likelihood expression (4) implies that the MLE of the consensus is

$\hat{\boldsymbol{\rho}}=\operatornamewithlimits{argmax}_{\boldsymbol{\rho}\in\mathcal{P}_{n}}\boldsymbol{\rho}^{T}\boldsymbol{\bar{r}}.$
(5)

We observe that the optimization problem (5) is equal to the one discussed in
Marden (1995) for estimating the central ranking under uniformity of the
random permutations, that the author proved to admit the closed-form solution

$\hat{\boldsymbol{\rho}}=(\hat{\rho}_{1},\ldots,\hat{\rho}_{i},\ldots,\hat{\rho}_{n})\,\,\,\text{with}\,\,\,\hat{\rho}_{i}=\text{rank}(\bar{r}_{i})$
(6)

(Theorem 2.2, page 29). This implies that the simple and fast Borda rank
aggregation method, which amounts to order the items according to their sample
average rank, coincides with the MLE of the consensus ranking of the MMS. Even
if this theoretical result, pointed out for the MMS in Feigin and Cohen
(1978), efficiently solves the critical step of estimating the discrete
parameter $\boldsymbol{\rho}$ regardless of the number $n$ of compared items
and the sample size $N$, to the best of our knowledge it has been completely
disregarded in those few works dealing with the MMS and extensions thereof. In
fact, the optimization (5) has been traditionally addressed as if a closed-
form solution did not exist, that is, with computationally intensive
approaches, such as global or local search methods whose effectiveness and
applicability strongly depend on the values of $n$ and $N$ (see for example
Murphy and Martin, 2003; Busse et al., 2007). These search methods require the
exploration of a discrete space with a fast-growing dimension and, hence, are
inefficient both in terms of execution time and convergence achievement as $n$
and $N$ increase, becoming soon infeasible (typically, NP-hard as proven in
Bartholdi et al., 1989). The omission of this important result is also
testified by the lack of specific and effective software for performing
inference on the MMS, that could take advantage of the related computational
simplifications. In fact, the use of the Borda method has been so far mainly
recommended for the MM with Kendall distance, either as approximation of the
MLE of $\boldsymbol{\rho}$ (Ali and Meilă, 2012) or for initializing the local
search (Irurozki et al., 2016).

Concerning the estimation of the concentration parameter $\theta$, one needs
to maximize the profile log-likelihood
$\ell(\hat{\boldsymbol{\rho}},\theta|\underline{\boldsymbol{r}})$ to get
$\hat{\theta}$ as the solution of the following equation

$-Z^{\prime}(\theta)/Z(\theta)=2(c_{n}-\hat{\boldsymbol{\rho}}^{T}{\boldsymbol{\bar{r}}}).$
(7)

As noticed by Fligner and Verducci (1986), $\hat{\theta}$ is the value
equating the expected Spearman distance $\text{E}_{\theta}[D]$ under the MMS,
that is, the left hand side of (7), to the sample average Spearman distance
$\bar{d}=\frac{1}{N}\sum_{l=1}^{L}N_{l}d(\boldsymbol{r}_{l},\hat{\boldsymbol{\rho}})$,
corresponding to the right-hand side of (7). Actually, $\text{E}_{\theta}[D]$
is a strictly decreasing function of $\theta$, but the root of the above
equation is not available in closed-form. So, the MLE can be found
numerically, for instance, via a Newton-Raphson algorithm (Marden, 1995),
although the calculations are demanding for all but small values of $n$,
because of the intractability of the model normalizing constant.

Inference in the presence of intractable normalizing constants has been
tackled in the literature either with algorithms aimed at finding consistent
estimators for the unknown parameters (see for example Murray et al., 2012;
Andrieu and Roberts, 2009), either with off line approximations of the
normalizing constant (Vitelli et al., 2018; Xu et al., 2018). The first stream
of methods is unfeasible in case of the MMS, because it typically requires the
direct simulation from the model likelihood (Vitelli et al., 2018). We here
propose to solve the inferential issue by introducing a new analytical
approximation of $Z(\theta)$ (see Section 3), whose performance is compared
with alternative ones.

###  2.3 Likelihood and inference for the MMS mixture model

Let $\boldsymbol{z}_{l}=(z_{l1},\dots,z_{lG})$ be the latent variable
collecting the binary group membership indicators, that is, $z_{lg}=1$ when
the $l$-th observation belongs to the $g$-th group and $z_{lg}=0$ otherwise.
The complete-data log-likelihood for the $G$-components mixture of MMS is

$\begin{split}\ell_{c}(\boldsymbol{\rho},\boldsymbol{\theta},\boldsymbol{\omega},\underline{\boldsymbol{z}}|\underline{\boldsymbol{r}})=&\sum_{l=1}^{L}\sum_{g=1}^{G}N_{l}z_{lg}\big{[}\log\omega_{g}+\\\
&-2\theta_{g}\left(c_{n}-{\boldsymbol{\rho}^{T}_{g}}\boldsymbol{r}_{l}\right)-\log
Z(\theta_{g})\big{]},\end{split}$

where $\boldsymbol{\omega}=(\omega_{1},\dots,\omega_{G})$ and
$\boldsymbol{\theta}=(\theta_{1},\dots,\theta_{G})$ are, respectively, the
mixture weights and the component-specific concentration parameters, whereas
$\boldsymbol{\rho}$ is a $G\times n$ matrix, whose rows indicate the consensus
rankings of the mixture components. In order to perform MLE in the presence of
latent variables, we employ the EM algorithm (Dempster et al., 1977). In
particular, at iteration $(t+1)$, the proposed EM scheme consists in the
following steps:

E-step: for $l=1,\dots,L$ and $g=1,\dots,G$, - compute the posterior
membership probabilities
$\hat{z}_{lg}^{(t+1)}=\dfrac{\hat{\omega}_{g}^{(t)}\mathbb{P}\left(\boldsymbol{r}_{l}\big{|}\hat{\boldsymbol{\rho}}_{g}^{(t)},\hat{\theta}_{g}^{(t)}\right)}{\sum_{g^{\prime}=1}^{G}\hat{\omega}_{g^{\prime}}^{(t)}\mathbb{P}\left(\boldsymbol{r}_{l}\big{|}\hat{\boldsymbol{\rho}}_{g^{\prime}}^{(t)},\hat{\theta}_{g^{\prime}}^{(t)}\right)}.$
M-step: for $g=1,\dots,G$, - compute the mixing weights
$\hat{\omega}_{g}^{(t+1)}=\frac{\hat{N}_{g}^{(t+1)}}{N}$ where
$\hat{N}_{g}^{(t+1)}=\sum_{l=1}^{L}N_{l}\hat{z}_{lg}^{(t+1)}$ is the number of
sample units allocated in the $g$-th component.
- compute the consensus rankings $\hat{\boldsymbol{\rho}}_{g}^{(t+1)}=\left(\hat{\rho}_{g1}^{(t+1)},\ldots,\hat{\rho}_{gi}^{(t+1)},\ldots,\hat{\rho}_{gn}^{(t+1)}\right),$ with $\hat{\rho}_{gi}^{(t+1)}=\text{rank}\left(\bar{r}_{gi}^{(t+1)}\right)$ and $\bar{r}_{gi}^{(t+1)}=\sum_{l=1}^{L}N_{l}\hat{z}_{lg}^{(t+1)}r_{li}/\hat{N}_{g}^{(t+1)}$.
- determine the concentration parameter estimates $\hat{\theta}_{g}^{(t+1)}$ as the solution of $\text{E}_{\theta_{g}}[D]=2\left(c_{n}-\boldsymbol{\rho}_{g}^{T(t+1)}{\boldsymbol{\bar{r}}}_{g}^{(t+1)}\right).$

Note that Murphy and Martin (2003) first provided an EM algorithm for the MLE
of MM mixtures with a generic metric $d(\cdot,\cdot)$, including the Spearman.
However, their estimation procedure is not specialized for each considered
metric implying that, in the case of Spearman distance, it completely
disregards the closed-form solution of the M-step for $\boldsymbol{\rho}$
whereas, for the decomposable metrics, the convenient expression of the
normalizing constant is totally ignored. Actually, the applications in their
work were limited to full rankings of only $n=5$ items, for which the
permutation space can be explored exhaustively and, even with a time-consuming
inferential process, the MLE is more likely to be achieved. Instead, our
implementation of the EM is completely specialized for the Spearman distance,
thus exploiting the shortcuts explained in the previous sections.

Concerning publicly available softwares, the R package rankdist (Qian and Yu,
2019) performs MLE for the MM in the presence of full and top-$q$ ranking
data, including the MMS and its finite mixture extension. However, the
analytical simplifications associated to the Spearman distance are not
considered, and the MLE solution for $\boldsymbol{\rho}$ is searched locally,
slowing down the computational time. In addition, the program generates all
the $n!$ permutations to compute the normalizing constant $Z(\theta)$ needed
to estimate the concentration parameter, and hence crashes for $n>11$ by
limiting the range of MMS applications to few items. In case of top-$q$
rankings the program works with at most $n=7$ items.

The pmr R package (Lee and Yu, 2013) performs frequentist inference on the
MMS, but again, it generates all the $n!$ possible permutations, to search the
MLE of the consensus ranking globally. As a result, it quickly runs out of
memory restricting the MMS applicability to datasets with a maximum of $n=11$
items.

The R package BayesMallows (Sørensen et al., 2020) implements Bayesian
inference for the MMS with any right-invariant distance, including the
Spearman, and it also supports partial rankings and mixture models. Even if
the procedure scales well with $n$, the implementation of the MMS does not
consider the analytical simplifications of the Spearman distance, which could
speed up the algorithm considerably. In BayesMallows the normalizing constant
of the MM is approximated by means of an importance sampling scheme. In
Section 3 we compare the performance of our novel approximation of $Z(\theta)$
with different types of approximations, including the one available in
Sørensen et al. (2020).

###  2.4 Mixture model extension for partial rankings

In real applications it could happen that the ranking data are not completely
observed. This occurs, for instance, when the assessors rank only a subset of
the items, typically their top-$q$ preferred ones, or there may be missing
data either at random or by design, giving rise to various forms of partial
observations.

In this section we extend the mixture analysis previously described for
accommodating inference on partial data. In the spirit of Beckett (1993),
whose goal was to estimate the parameters of the basic MM from partial
observations, we generalize its MLE approach based on the EM algorithm by
assuming the mixture of MMS as the sampling distribution. We adopt a data
completion argument that suitably combines the two main sets of latent
variables: the missing ranks of the partial sequences and the unobserved group
memberships.

Suppose that $N_{l}$ respondents provide a partial sequence
$\boldsymbol{r}_{l}$ where only a subset
$\mathcal{I}_{l}\subseteq\\{1,2,\dots,n\\}=\mathcal{I}$ of
$n_{l}=|\mathcal{I}_{l}|$ items are actually ranked. As a consequence, the
position $r_{li}$ is known only if $i\in\mathcal{I}_{l}$. We stress that, in
our setting, arbitrary types of censoring patterns are allowed, implying that
the set of the assigned positions may be any not-empty subset of the first $n$
integers, leading to $n_{l}=1,\dots,n-1$. Moreover, a full ranking
$\tilde{\boldsymbol{r}}_{l}$ is defined as the complete sequence obtained by
combining the $n_{l}$ observed entries $\boldsymbol{r}_{l}$ with the latent
unassigned $n-n_{l}$ ranks for the items
$i^{\prime}\in\mathcal{I}\setminus\mathcal{I}_{l}$. The latent components of
the partially ranked vector are inferred through the EM algorithm, that
exploits the information provided by the entire sample to fill the missing
entries of the sequence coherently with the observed ones.

Let $\mathcal{C}(\boldsymbol{r}_{l})\subset\mathcal{P}_{n}$, be the set of
full rankings which are compatible with the partial sequence
$\boldsymbol{r}_{l}$. This means that $\mathcal{C}(\boldsymbol{r}_{l})$
collects all the complete ranked vectors on which the application of a certain
censoring pattern yields the partial observation $\boldsymbol{r}_{l}$. Under
the mixture model, the conditional probability of observing the full ranking
$\tilde{\boldsymbol{r}}_{m}\in\mathcal{C}(\boldsymbol{r}_{l})$ is

$\begin{split}p_{r_{l}}(\tilde{\boldsymbol{r}}_{m})&\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\mathbb{P}(\tilde{\boldsymbol{r}}_{m}\,|\boldsymbol{r}_{l},\boldsymbol{\rho},\boldsymbol{\theta},\boldsymbol{\omega})\\\
&=\frac{\sum\limits_{g=1}^{G}\omega_{g}e^{-2\theta_{g}\left(c_{n}-{\boldsymbol{\rho}^{T}_{g}}\tilde{\boldsymbol{r}}_{m}\right)-\log
Z(\theta_{g})}}{\sum\limits_{\tilde{\boldsymbol{s}}\in\mathcal{C}(\boldsymbol{r}_{l})}\sum\limits_{g=1}^{G}\omega_{g}e^{-2\theta_{g}\left(c_{n}-{\boldsymbol{\rho}^{T}_{g}}\tilde{\boldsymbol{s}}\right)-\log
Z(\theta_{g})}}.\end{split}$ (8)

Equality (8) extends expression (6) in Beckett (1993) to the finite mixture
case. The observed frequency $N_{l}$ of respondents reporting each partial
ranking $\boldsymbol{r}_{l}$ is then distributed over the full rankings that
may give rise to $\boldsymbol{r}_{l}$, in order to obtain the E-step estimate
of the frequency of $\tilde{\boldsymbol{r}}_{m}$ as

$\hat{N}_{m}=\sum_{l:\,\tilde{\boldsymbol{r}}_{m}\in\mathcal{C}(\boldsymbol{r}_{l})}N_{l}p_{r_{l}}(\tilde{\boldsymbol{r}}_{m}).$
(9)

The above sum accounts for the fact that, under varying censoring patterns,
the same full ranking $\tilde{\boldsymbol{r}}_{m}$ can be compatible with
different partial rankings; for instance, $\tilde{\boldsymbol{r}}=(1,2,4,3)$
is compatible with both $\boldsymbol{r}_{1}=(1,2,NA,NA)$ and
$\boldsymbol{r}_{2}=(1,NA,NA,NA)$. Obviously, if the observed data share the
same censoring pattern (e.g. they are all top-3 implying $n_{l}=3$ for all
$l=1,\dots,L$), the sum in (9) can be avoided and the computation simplifies
greatly.

In practice, the EM scheme detailed in Section 2.3 is adapted for handling
partial rankings by adding a further E-step, needed to account for the latent
ranked entries, and by considering the augmented data
$\tilde{\boldsymbol{r}}_{m}\in\mathcal{C}(\boldsymbol{r}_{1})\cup\cdots\cup\mathcal{C}(\boldsymbol{r}_{L})$
as the sample to draw inference on.

##  3 A new approximation of the partition function

If, on the one hand, the existence of the closed-form for the MLE of
$\boldsymbol{\rho}$ avoids the computationally demanding local search for
inferring the consensus ranking, on the other, considering large values of $n$
affects the calculation of $Z(\theta)$ which is needed in the estimation steps
of (i) the concentration parameter (see eq. (7)), (ii) the posterior
membership probabilities and (iii) the conditional probabilities of observing
the full rankings compatible with the partial information (see eq. (8)). In
fact, the rapidly-increasing size $n!$ makes the exhaustive enumeration of the
sequences in $\mathcal{P}_{n}$, and hence the direct calculation of
$Z(\theta)=\sum_{\boldsymbol{r}\in\mathcal{P}_{n}}e^{-2\theta\,\left(c_{n}-\boldsymbol{e}^{T}\boldsymbol{r}\right)}$,
feasible only for small $n$.

In this section, we propose a novel approximation of the distribution of the
Spearman distance which enables inference on rankings of a large number $n$ of
items. In fact, this approximation, which may be interesting in its own right
(e.g. for sampling procedures), is then used to approximate both $Z(\theta)$
and $\text{E}_{\theta}[D]$.

###  3.1 The new approximation in practice

Regardless of the metric adopted in the MM, the partition function can be more
conveniently expressed as

$Z(\theta)=\sum_{d_{n}\in\mathcal{D}_{n}}N_{d_{n}}\,e^{-\theta d_{n}},$ (10)

where
$N_{d_{n}}=|\\{\boldsymbol{r}\in\mathcal{P}_{n}\,:\,d(\boldsymbol{r},\boldsymbol{e})=d_{n}\\}|$
is the number of rankings at distance $d_{n}$ from $\boldsymbol{e}$ and
$\mathcal{D}_{n}$ is the set of possible distance values among the
$n$-dimensional rankings. In the case of Spearman distance,
$\mathcal{D}_{n}=\\{d_{n}=2h\,:\,h\in\mathbb{N}_{0}\text{ and }0\leq d_{n}\leq
d_{\text{max},n}=2\binom{n+1}{3}\\}$, implying that the number of terms in
(10) is of order $n^{3}$ rather than $n!$ as in the naïve representation of
$Z(\theta)$ related to (3). The values $N_{d_{n}}$ indicate the frequency
distribution of the Spearman distance under the uniform (null) model and
correspond to the sequence A175929 available only up to $n\leq 14$ in the OEIS
(Sloane, 2017; Vitelli et al., 2018). We here propose a strategy to
approximate $N_{d_{n}}$ when $n>14$.

As proven in Kendall (1970), under the null model the Spearman distance is
asymptotically normally distributed (see Figure 1). However, the convergence
is slow and not uniform over the domain of the distribution. In particular,
even with $n$ large enough, for the distribution to be reliably approximated
in the bulk, the tail behaviour will be different from the one predicted by a
normal distribution.

Figure 1: Comparison between the probability mass distribution of the Spearman
distance under the null model (gray) and the limiting normal density (red) for
increasing $n$ values.

To usefully approximate the distribution of extreme values, we could in
principle use large deviation theory (Cramér, 1938; Varadhan, 1966) but, as
explained in Mukherjee (2016), the computation of explicit results is often
not feasible. We can nevertheless exploit general properties of the large
deviation rate function, and what we know on the values of $N_{d_{n}}$ for
$n\leq 14$, to obtain an approximation. As a matter of fact, in Mukherjee
(2016) it is proved the existence and differentiability of the (scaled)
cumulant generating function,111In the original notation, the author proved
the existence and differentiability of the function
$Z\left(f,\theta\right)=\lim_{n\to\infty}\frac{1}{n}\log\text{E}\left[\exp(n\,\theta\,\nu_{\pi}\left[f\right])\right]$.
By setting $f(x,y)=(x-y)^{2}$ and $\theta=k$, we obtain the desired result.

$\lambda\left(k\right)\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\lim_{n\to\infty}\frac{1}{n}\log\text{E}\left[e^{knx_{n}}\right]\quad
k\in\mathbb{R},$

with $x_{n}=\frac{d_{n}}{d_{\max,n}}$. This result allows us to apply the
Gärtner-Ellis theorem (Gärtner, 1977; Ellis, 1984), to obtain the approximate
rate function

$\xi\left(x_{n}\right)\approx\frac{1}{n}\log\left(\frac{N_{d_{n}}}{n!}\right),$
(11)

with $\xi\left(x_{n}\right)$ only weakly dependent on $n$ for $n\gg 1$. Making
use of the known values of $N_{d_{n}}$ for $n\leq 14$, we can derive for
$\xi\left(x_{n}\right)$ the form

$\xi\left(x_{n}\right)=a_{0}+a_{1}\log\left[x_{n}\left(1-x_{n}\right)\right]+a_{2}x_{n}\left(1-x_{n}\right).$
(12)

The formal derivation of (12) and the value of the coefficients
$(a_{0},a_{1},a_{2})$ are detailed in Supplementary Material A. With this
result, our approximation of $N_{d_{n}}$ for $8<d_{n}<d_{\max,n}-8$ is

$\hat{N}_{d_{n}}\propto n!\exp\left[n\,\xi\left(x_{n}\right)\right].$ (13)

The above approximation is used for $n>14$. The values of ${N}_{d_{n}}$ for
$d_{n}=0,2,4,6$ (as well as those for their symmetrical
$d_{n}=d_{\max,n},d_{\max,n}-2,d_{\max,n}-4$ and $d_{\max,n}-6$) follow from
combinatorial facts and are equal respectively to $1,(n-1),\binom{n-2}{2}$ and
$(n-2)^{3}/6-(n-2)^{2}+23(n-2)/6-1$.

Figure 2: Comparison between the limiting normal density (red line) of the
probability distribution of the Spearman distance under the null model (gray),
and the density obtained with our new approximation (blue line). Left:
original scale. Right: logarithmic scale. $n=14$.

Figure 2 displays the probability mass distribution of the Spearman distance
for $n=14$ (in gray), its limiting normal density (red line) and the density
obtained with our new approximation (blue line) in the original scale (left
panel) and in the logarithmic scale (right panel) for better visualization.

Equipped with the approximation (13), we can approximate the partition
function (10) and
$\text{E}_{\theta}[D]=\sum_{d_{n}\in\mathcal{D}_{n}}dN_{d_{n}}\,e^{-\theta
d_{n}}/\sum_{d_{n}\in\mathcal{D}_{n}}N_{d_{n}}\,e^{-\theta d_{n}}$ by simply
replacing the exact $N_{d_{n}}$ with $\hat{N}_{d_{n}}$.

Figure 3: Ratio between alternative approximations of the partition function
and the exact one for $n=14$ and $\theta\in(0,1]$ (left and right) and
$\theta\in(0,0.02]$ (middle).

###  3.2 A comparison with other solutions

In this section we compare our new approximation
$\hat{Z}_{\text{New}}(\theta)=\sum_{d_{n}\in\mathcal{D}_{n}}\hat{N}_{d_{n}}\,e^{-\theta
d_{n}}$ with two others proposed in the literature, one analytical and the
other simulation-based. The former, based on the von Mises-Fisher (vMF)
distribution and originally proposed in McCullagh (1993), has been recently
exploited by Xu et al. (2018) for inference in their angle-based model for
rankings and is given by

$Z(\theta)\approx\hat{Z}_{\text{vMF}}(\theta)=\frac{2^{\frac{n-3}{2}}n!I_{\frac{n-3}{2}}(\kappa)\Gamma(\frac{n-1}{2})}{\kappa^{\frac{n-3}{2}}e^{\kappa}},$

where $\kappa={\theta n(n^{2}-1)}/6$, $I_{x}(\cdot)$ is the modified Bessel
function of the first kind with order $x$, and $\Gamma(\cdot)$ is the Gamma
function. The reader is referred to Supplementary Material B for the formal
derivation of $\hat{Z}_{\text{vMF}}(\theta)$ in the context of (non-
standardized) rankings. The latter $\hat{Z}_{\text{IS}}(\theta)$, based on the
Importance Sampling (IS) technique, is due to Vitelli et al. (2018) to which
we refer for details.

In Figure 3 we plot the ratio $\hat{Z}(\theta)/Z(\theta)$ as a function of
$\theta$ for $n=14$. In the left and middle panels, where the approximations
$\hat{Z}_{\text{New}}$ (in green) and $\hat{Z}_{\text{vMF}}$ (in orange) are
compared over the interval $\theta\in(0,1]$ and $\theta\in(0,0.02]$
respectively, we notice that for small values of $\theta$, the vMF
approximation and our proposal perform both very good (middle panel). However,
the vMF approximation deteriorates fast already for $\theta>0.015$ and its
negative bias rapidly becomes prominent (left panel). This happens because the
vMF approximation is an integral over the entire ($n-1$)-dimensional sphere
enclosing all the points of $\mathcal{P}_{n}$ and, for growing $\theta$, the
permutations are more and more concentrated in a small region of this sphere
under the MMS, implying that $\hat{Z}_{\text{vMF}}$ is strongly biased. On the
contrary, our new approximation performs very well for increasing values of
$\theta$, because it exploits the knowledge of the true values of $N_{d_{n}}$
available for small $d_{n}$, which are the ones that mostly matter in (10).

Instead, the IS approximation has an oscillating behavior around the exact,
and its error is generally small even if, when compared to our new
approximation, it is closer to the exact for few values of $\theta$. Note
that, even if the performance of the IS is generally good, being based on
random sampling, it is computationally very demanding, while our new
approximation is immediately available for any value of $n$.

From the above results we evince that our proposed new approximation performs
very good for a wide range of values of $\theta$, not only for small ones,
like the analytical approximation resorting on the vMF distribution. This is
particularly relevant when dealing with mixture models, where the cluster-
specific concentrations typically grow with the number of clusters detected.

To complete the study of the proposed approximation, in Supplementary Material
C we document that its effect on the inference on $\theta$ is generally
negligible.

##  4 Simulation studies

###  4.1 Experiments with different sample sizes and number of items

The first simulation study considers data drawn from the homogeneous MMS
($G^{*}=1$) under different parameters settings, and then evaluates the
performance of our fitting procedure in terms of model parameters estimates
with multiple criteria. Specifically, for each combination of $(N,n)$ with
sample size $N\in\\{50,200,500,1000\\}$ and number of items
$n\in\\{5,10,14,15,25,50,100\\}$, we independently simulate 100 datasets from
the homogeneous MMS, by randomly choosing the true consensus ranking as
$\boldsymbol{\rho}\sim\text{Unif}\\{\mathcal{P}_{n}\\}$ and sampling uniformly
$\theta\sim\text{Unif}(\theta_{\text{min}},\theta_{\text{max}})$ over
different intervals depending on the value of $n$ (see Table 1, which collects
the bounds of the uniform sampling intervals). In this way, we have a
comparable concentration over the permutation set in the samples with
different values of $n$.

$n$ | 5 | 10 | 14,15 |  | 25 | 50 | 100
---|---|---|---|---|---|---|---
$\theta_{\text{min}}$ | .15 | .05 | .025 |  | .010 | .002 | .0002
$\theta_{\text{max}}$ | .30 | .10 | .050 |  | .020 | .010 | .0020
Table 1: Bounds of the intervals for the uniform sampling of $\theta$ for
different number of items.

The inferential ability of the proposed EM algorithm to recover the actual
parameters values has been evaluated by averaging over the 100 datasets the
following measures:

  * •

$m_{\theta}=|\hat{\theta}-\theta|/\theta$, the relative error between the true
and estimated concentration;

  * •

$m_{\boldsymbol{\rho}}=d(\hat{\boldsymbol{\rho}},\boldsymbol{\rho})/d_{\max,n}$,
the relative distance between the true and estimated consensus;

  * •

$\phi_{\boldsymbol{\rho}}=\mathbb{1}_{(\boldsymbol{\rho}=\hat{\boldsymbol{\rho}})}$,
the matching between the true and the estimated consensus;

  * •

$\phi_{G}=\mathbb{1}_{(\hat{G}={G^{*}})}$, the matching between the true and
estimated number of components. The optimal number $\hat{G}$ of groups is
automatically selected by means of the Bayesian Information Criterion, BIC
(Schwarz, 1978). In particular, we select $\hat{G}$ by applying the elbow rule
to the screeplot of the BIC values (Zhao and Hautamäki, 2008; Shi et al.,
2021), in place of the more popular choice of selecting the minimum BIC,
because of the tendency of the MMS to over-fit the data (Murphy and Martin,
2003; Qian and Yu, 2019).222Note that this criterion is useful only for
automatizing the simulations: when choosing the number of clusters in
applications the researcher typically uses the BIC, as well as other relevant
measures, as a guidance, and decides the number of groups based also on other
heuristic criteria (such as the size of the clusters, their separation, or
prior information depending on the application at hand).

$n$ | $N$ | $m_{\theta}$ | $m_{\boldsymbol{\rho}}$ | $\phi_{\boldsymbol{\rho}}$ | $\phi_{G}$
---|---|---|---|---|---
5 | 50 | .078 | .0000 | 1.00 | 1.00
5 | 200 | .045 | .0000 | 1.00 | 1.00
5 | 500 | .034 | .0000 | 1.00 | 1.00
5 | 1000 | .020 | .0000 | 1.00 | 1.00
10 | 50 | .071 | .0039 | .54 | .87
10 | 200 | .034 | .0002 | .97 | .91
10 | 500 | .024 | .0000 | 1.00 | 1.00
10 | 1000 | .016 | .0000 | 1.00 | 1.00
14 | 50 | .061 | .0063 | .08 | .99
14 | 200 | .031 | .0010 | .63 | 1.00
14 | 500 | .019 | .0001 | .98 | .98
14 | 1000 | .015 | .0000 | 1.00 | 1.00
(a) $n\leq 14$ $N$ | $m_{\theta}$ | $m_{\boldsymbol{\rho}}$ | $\phi_{\boldsymbol{\rho}}$ | $\phi_{G}$ |
---|---|---|---|---|---
15 | 50 | .062 | .0050 | .06 | 1.00
15 | 200 | .030 | .0006 | .72 | 1.00
15 | 500 | .020 | .0000 | .98 | .99
15 | 1000 | .015 | .0000 | 1.00 | 1.00
25 | 50 | .040 | .0052 | .00 | 1.00
25 | 200 | .029 | .0012 | .10 | 1.00
25 | 500 | .022 | .0002 | .58 | 1.00
25 | 1000 | .021 | .0000 | .93 | 1.00
50 | 50 | .053 | .0051 | .00 | 1.00
50 | 200 | .042 | .0013 | .00 | .99
50 | 500 | .041 | .0005 | .00 | 1.00
50 | 1000 | .047 | .0002 | .05 | 1.00
50 | 10000 | .042 | .0000 | .98 | 1.00
100 | 50 | .104 | .0126 | .00 | 1.00
100 | 200 | .083 | .0026 | .00 | 1.00
100 | 500 | .091 | .0018 | .00 | .93
100 | 1000 | .093 | .0006 | .00 | 1.00
(b) $n>14$

Table 2: Inferential performance of the EM algorithm on simulated full
rankings.

In Table 2 we report the results of the fit. For all $n$, we see that all
performance criteria improve with $N$ and are generally very well, at least
for $n\leq 25$. As a matter of fact, the performance of the fitting procedure
is linked to the sparsity of the problem: when the data are too sparse (that
is, with large $n$ and small $N$), and the concentration parameter is small,
recovering the exact true consensus, measured by $\phi_{\boldsymbol{\rho}}$,
becomes very difficult (see the results for $n=50,100$ in Table 2(b)).
However, $m_{\boldsymbol{\rho}}$ is very small even in these extreme cases,
testifying a general good performance of the procedure. In addition, it is
enough to sufficiently enlarge the sample size to recover exactly the true
$\boldsymbol{\rho}$ in almost all samples considered (see Table 2(b) for
$n=50$ and $N=10000$). Finally, from Table 2, we note also that the
inferential ability does not seem to be worsened by the introduction of the
approximation of $Z(\theta)$. In fact, $m_{\theta}$, is comparable for the
cases $n=14,15$, the first estimated with the exact $Z(\theta)$, and the
second with $\hat{Z}_{\text{New}}(\theta)$.

###  4.2 Experiments with heterogeneous data

In this section we study the performance of our procedure in terms of its
ability to correctly recover the cluster structure.

We simulate 100 datasets from different mixtures of MMS with varying number of
clusters $G^{*}=2,3,4$, while keeping fixed the number of items $n=25$, and
the sample size $N=2000$. For each dataset, we sample (i) the mixture weights
from a symmetric Dirichlet distribution with large shape parameter to favour
populated and balanced clusters, that is,
$(\omega_{1},\dots,\omega_{G^{*}})\sim\text{Dirichlet}(\alpha=(2G^{*},\dots,2G^{*}))$;
(ii) the cluster consensus rankings
$\boldsymbol{\rho}_{1},\dots,\boldsymbol{\rho}_{G^{*}}$ in $\mathcal{P}_{n}$
in order to be at least at Spearman distance $d_{\max,n}/(n-1)=(n^{2}-1)/3$
from each other; and (iii) the concentration parameters uniformly in different
intervals, corresponding to three different scenarios concerning the
separation of the clusters from each other (see Table 3).

| Cluster separation
---|---
| High | Medium | Low
$\theta_{\text{min}}$ | .004 | .003 | .002
$\theta_{\text{max}}$ | .006 | .005 | .004
Table 3: Bounds of the intervals for the uniform sampling of $\theta$ in the
three group separation scenarios.

For each sample, we then fit the data with varying number of components
$G=1,\dots,6$ and select the optimal $\hat{G}$ by means of the BIC like in the
previous section.

| Cluster separation
---|---
| High |  | Medium |  | Low
$G^{*}$ | $\phi_{G}$ | $\phi_{\boldsymbol{z}}$ |  | $\phi_{G}$ | $\phi_{\boldsymbol{z}}$ |  | $\phi_{G}$ | $\phi_{\boldsymbol{z}}$
2 | 1 | .00021 |  | 1.00 | .015 |  | .99 | .131
3 | 1 | .00016 |  | 1.00 | .016 |  | 1.00 | .135
4 | 1 | .00010 |  | .98 | .012 |  | .96 | .104
Table 4: Inferential performance of the EM algorithm on rankings simulated
from alternative MMS-mixtures.

In Table 4 we report $\phi_{G}=\mathbb{1}_{(\hat{G}=G^{*})}$ as before, plus
the mis-classification rate, defined as the fraction of the sampled units
classified in the wrong mixing group,
$\phi_{\boldsymbol{z}}=\sum_{j=1}^{N}\mathbb{1}(z_{j}\not=\dot{z}_{j})/N$,
with the group allocation $\dot{z}_{j}$ determined by the maximum a posteriori
(MAP) of the posterior membership probabilities. As expected, the ability to
correctly reconstruct the clusters decreases as the difficulty of the problem
increases. In general, however, we are happy with the results especially
because, in the most difficult scenario, the mis-classification rate is always
below 15%.

| $n=10$ |  | $n=15$
---|---|---|---
Censoring | $m_{\theta}$ |  | $m_{\boldsymbol{\rho}}$ | $\phi_{\boldsymbol{\rho}}$ |  | $\phi_{G}$ |  | $m_{\theta}$ |  | $m_{\boldsymbol{\rho}}$ | $\phi_{\boldsymbol{\rho}}$ |  | $\phi_{G}$
full | .0175 |  | .0000 | 1.00 |  | .99 |  | .0175 |  | .0000 | 1.00 |  | 1.00
A | .0192 |  | .0000 | 1.00 |  | .90 |  | .0176 |  | .0000 | 1.00 |  | .94
B | .0199 |  | .0002 | .96 |  | .83 |  | .0182 |  | .0002 | .90 |  | .82
Table 5: Inferential performance of the EM algorithm on simulated partial
rankings.

###  4.3 Simulations for incomplete rankings

Through an additional simulation study, we explore the impact of censoring on
the inference, by considering different patterns of missingness of partial
rankings drawn from the homogeneous MMS. To this aim, we first simulate 100
datasets of $N=500$ full rankings for the cases $n=10$ and $n=15$.
Specifically, the observations are drawn from the MMS with
$\boldsymbol{\rho}\sim\text{Unif}\\{\mathcal{P}_{n}\\}$ and
$\theta\sim\text{Unif}(0.1,0.15)$ for $n=10$ and
$\theta\sim\text{Unif}(0.5,0.6)$ for $n=15$. Then, partial top-rankings are
obtained by randomly censoring all the complete sequences according to two
patterns of truncation, labelled respectively as A and B and characterized by
an increasing rate of missingness of the items ranked in the bottom positions.
In both censoring scenarios, the random number of censored bottom positions
varied in $\\{2,3,4,5,6\\}$, but it was sampled by using probability masses
$\\{0.1,0.1,0.1,0.1,0.6\\}$ under the censoring scenario A and
$\\{0.1,0.2,0.4,0.2,0.1\\}$ under censoring scenario B.

From Table 5 we note that the ability to recover the exact parameter value, as
well as the proportion of times that the homogeneous MMS is correctly
selected, decreases as the rate of missingness increases. However, for both
censoring scenarios, all the criteria exhibit very good values when compared
to those obtained in the full ranking analysis, testifying the effectiveness
of our EM procedure to draw inference from partial observations.

##  5 Applications

###  5.1 Analysis of the sports data

We revisit the well-known sports benchmark data (Marden, 1995), available in
the R package Rankcluster (Jacques et al., 2014), which collects $N=130$ full
rankings. These represent the preferences of students at the University of
Illinois towards $n=7$ sports, namely Baseball, Football, Basketball, Tennis,
Cycling, Swimming and Jogging. We explored the presence of a cluster structure
in the sample through the estimation and comparison of alternative mixture
models with $G=1,\dots,5$, and components given by: 1) the MM with Spearman
distance (MMS); 2) the MM with Kendall distance (MMK); 3) the MM with Cayley
distance (MMC); 4) the MM with Hamming distance (MMH) and 5) the Plackett-Luce
model (PL).

| MMS | MMK | MMC | MMH | PL
---|---|---|---|---|---
$\hat{G}$ | $2$ | $2$ | $2$ | $1$ | $2$
BIC | 2144.5 | 2155.1 | 2190.2 | 2206.4 | 2195.5
Table 6: BIC and number of groups for the alternative mixtures fitted to the
sports data.

Table 6 shows the best fitting mixture under each parametric scenario, and the
corresponding estimated number $\hat{G}$ of clusters. With the only exception
of the MMH that does not identify the presence of groups ($\hat{G}=1$), all
the other models agree in recognizing the existence of two clusters. However,
one can note that our MMS is the model associated with the best fitting
performance (BIC=2144.5). The corresponding cluster-specific parameter
estimates (mixing weights, concentration and modal orderings) are displayed
below.

$\hat{\omega}_{1}=0.61\qquad\hat{\theta}_{1}=0.039$

$\hat{\boldsymbol{\rho}}_{1}=$

$\hat{\omega}_{2}=0.39\qquad\hat{\theta}_{2}=0.118$

$\hat{\boldsymbol{\rho}}_{2}=$

Both groups share the middle position assigned to Cycling, but they are
clearly distinguished by top and bottom ranks. In fact, the largest group
($\hat{\omega}_{1}=0.61$) is characterized by the preference towards
individual sports, whereas team sports are the most-liked by students of the
smallest group ($\hat{\omega}_{2}=0.39$).

Figure 4: Screeplot of BIC values for alternative number of groups in the MMS-mixture fitted to the reading genres data. Figure 5: Respondents’ characteristics statistically associated to the model-based clustering obtained from the application of the MMS-mixture to the reading genres data. Group | $\hat{\omega}$ | $\hat{\theta}$ | Consensus
---|---|---|---
1 | 0.33 | 0.048 | (Nov., Fan., Cla., Thr., Tee., Poe., Bio., Hor., Com., Hum., Ess.,)
2 | 0.23 | 0.036 | (Fan., Thr., Hor., Com., Nov., Tee., Hum., Cla., Bio., Ess., Poe.,)
3 | 0.44 | 0.054 | (Nov., Cla., Thr., Ess., Bio., Poe., Fan., Com., Hum., Hor., Tee.,)
Table 7: Parameter estimates of the MMS-mixture fitted to the reading genres
data.

###  5.2 Analysis of the reading genres

In this application we illustrate the effectiveness of our method for handling
partial rankings. To this purpose, we consider a brand new dataset that has
not been analyzed earlier in the ranking literature. This was obtained from an
online survey conducted in Italy in 2019 to investigate reading preferences.
It was promoted by the municipal administration of Latina (Latium, Italy), in
collaboration with Sapienza University of Rome and the School of Government of
the University of Tor Vergata. The dataset is composed of $N=507$ partial
top-5 rankings of $n=11$ reading genres ordered by the respondents according
to their personal preference. The genres are Classic, Novel, Thrillers,
Fantasy, Biography, Teenage, Horror, Comics, Poetry, Essay and Humor.

Due to the large sample size and the heterogeneity of the respondents with
respect to factors typically impacting on the reading preferences, such as
gender, age and education level, we expect the existence of different groups.
We explored this hypothesis through the estimation of $G$-component MMS-
mixtures with $G=1,\dots,6$. Figure 4 shows that BIC values remarkably
decrease up to 3 components. Moreover, for $G\geq 5$ a cluster with a
negligible size is always estimated. According to the elbow rule, we
considered $\hat{G}=3$ as the optimal number of groups and reported the
corresponding cluster-specific parameter estimates in Table 7. The estimates
of the modal rankings reveal that more serious and ponderous readings,
involving Classic and Essay, are ranked in high positions by respondents of
the third cluster, whereas all the genres which are typically the most-liked
ones among younger people all occupy bottom positions. This preference pattern
is completely reversed in the modal ranking of the second group, with all soft
readings (Fantasy, Thrillers, Horror and Comics) ranked in higher positions.
The third component represents more hybrid preferences with a mix of the two
aforementioned reading types in top and bottom positions. Interestingly, the
estimated group structure turned out to be significantly associated to several
individual aspects influencing reading attitudes (see Figure 5). In fact, the
third group is mainly characterized by respondents aged 20 or older, who read
many books in the last year and holding a master or doctoral degree. On the
other hand, respondents of the second group are mostly teenagers and reported
a lower propensity to the reading habit during the last year. Finally, we note
that the percentage of males in the second group is remarkably higher than in
the other two.

###  5.3 An illustration based on gene expression data

| Tumour type
---|---
Group | NB | BL | RMS | EWS
1 | 12 | 0 | 0 | 0
2 | 0 | 8 | 0 | 0
3 | 0 | 0 | 21 | 1
4 | 0 | 0 | 0 | 22
Table 8: Cross tabulation between the classification resulting from the MMS-
mixture application and the true tumour types.

Finally, we apply our methodology on a micro-array gene expression dataset, to
illustrate the ability of our method to analyze heterogeneous rankings of a
large number of items.

We use data adapted from Khan et al. (2001), which contains micro-array gene
expression profiles of four types of small round blue cell tumours (SRBCT) of
childhood for $N=64$ patients. These four different cancers, namely
neuroblastoma (NB), Burkitt lymphom (BL), rhabdomyosarcoma (RMS), and the
Ewing family of tumours (EWS), present a similar histology of SRBCT and thus
often leads to misdiagnosis.

In this illustration we aim at exploiting the MMS-mixture to classify patients
into the 4 different cancer types, based on their ranked gene expression
profiles.

We select the $n=60$ most differentially expressed genes throughout the sample
and, for each patient, we convert the individual quantitative expression
profile into a ranking. The rank transformation is sometimes used in genomics
applications, because inference based on rankings is both more robust than
working on the actual continuous measurements (Eliseussen et al., 2022), and
avoids the arbitrary choice of a numerical pre-processing procedure to handle
the raw outcomes (Mollica and Tardella, 2014).

We then estimate a 4-component MMS-mixture and check the agreement between the
model-based classification with the one based on the true tumour type.

Results are reported in Table 8, where we see that we classify all tumours
correctly except one EWS which is erroneously assigned to the RMS group. From
Table 9, where we report the classification implied by the simple hierarchical
clustering, we note that this method is not able to correctly classify a
considerable subset of RMS patients, which are instead erroneously mixed with
the BL ones.

| Tumour type
---|---
Group | NB | BL | RMS | EWS
1 | 12 | 0 | 0 | 0
2 | 0 | 8 | 10 | 0
3 | 0 | 0 | 11 | 1
4 | 0 | 0 | 0 | 22
Table 9: Cross tabulation of the hierarchical clustering groups and true tumour type. | Rank
---|---
Group | 1 | 2 | 3
1 (NB) | AF1Q | MAP1B | CRMP1
2 (BL) | HLA-DPB1 | MYC | CYFIP2
3 (RMS) | INS-IGF2 | IGF2 | COL3A1
4 (EWS) | FCGRT | MYC | OLFM1
Table 10: Top-3 genes in the consensus rankings of the 4 estimated groups.

In Figure 6 we plot the tumours based on their gene expression levels in the
2D space of the top-2 principal components of the gene expression levels of
the sample. In the figure, the tumours are colored based on their
classification obtained with the MMS. It can be seen that two groups
classified by the MMS (groups 2 and 3) are highly separable in this space,
while groups 1 and 4 appear less separable.

Figure 6: Visualization of MMS classification of cancers. The axes are the
first two principal components of the original quantitative gene expression
levels.

The mixture of MMS also provides the consensus rankings of the four
components, which correspond to the estimated modal orderings of genes for the
four cancer types. In Table 10, we report the top-3 expressed genes of the
different tumours, which are in line with the previous literature, that
identified them as those capable of diagnosing and discerning the considered
cancers types (see, for instance, Table 2 in the work by Chen et al., 2007).

##  6 Conclusions and future work

Within the class of MMs for the analysis of rankings, the version based on the
use of the Spearman distance as metric on the permutation space has been
barely explored in the literature, both from an analytical and a computational
point of view. One of the main reasons lies in the fact that, unlike the
Kendall, the Hamming or the Cayley distances, the Spearman cannot be
decomposed into the sum of independent terms, a property which leads to a
closed-form solution for the partition function and to convenient
simplifications in the estimation procedure. Moreover, with decomposable
metrics, more straightforward extensions of the MM to handle partial rankings
can be defined. Nevertheless, the Spearman distance induces a parametric
ranking distribution which parallels the multivariate normal model and
actually, similarly to the latter, admits an analytical solution for the MLE
of the modal consensus ranking given by the ordering of the average rank
vector. Motivated by this theoretical argument and the resulting computational
simplifications, which instead do not hold for the other metrics, we extend
the MMS to the finite mixture context and develop an efficient and accurate EM
algorithm to perform inference on samples of both full and partial rankings
drawn from heterogeneous populations and affected by censoring of arbitrary
positions.

Additionally, in order to address the inferential issues arising with a large
number of ranked items, this work introduces an approximation of the model
normalizing constant and assesses its performance with respect to alternative
methods. Finally, the utility of our approach has been proven through
applications to both synthetic and real-world data.

Our work provides original research that can be further extended in several
directions. From a computational perspective, an R package to promote a wider
use of the MMS and mixtures thereof is currently under construction. From a
methodological perspective, model-based clustering could be enhanced by
integrating the finite mixture framework with the regression analysis, for
identifying significant factors affecting the preferences and better
characterizing the groups of rankers. This idea is currently under study
(Crispino et al., 2022) by developing a mixture of experts model for the MMS,
in the spirit of Gormley and Murphy (2008b).

## References

  * Ali and Meilă (2012) Ali, A. and M. Meilă (2012).  Experiments with Kemeny ranking: What works when?  Mathematical Social Sciences 64(1), 28–40.
  * Alvo and Yu (2014) Alvo, M. and P. L. H. Yu (2014).  Statistical Methods for Ranking Data.  Frontiers in Probability and the Statistical Sciences. New York, USA: Springer.
  * Andrieu and Roberts (2009) Andrieu, C. and G. O. Roberts (2009).  The pseudo-marginal approach for efficient Monte Carlo computations.  The Annals of Statistics 37(2), 697–725.
  * Bartholdi et al. (1989) Bartholdi, J. J., C. A. Tovey, and M. A. Trick (1989).  The computational difficulty of manipulating an election.  Social Choice and Welfare 6(3), 227–241.
  * Beckett (1993) Beckett, L. A. (1993).  Maximum Likelihood Estimation in Mallows’s Model Using Partially Ranked Data.  In M. A. Fligner and J. S. Verducci (Eds.), Probability Models and Statistical Analyses for Ranking Data, New York, pp. 92–107. Springer New York.
  * Busse et al. (2007) Busse, L. M., P. Orbanz, and J. M. Buhmann (2007).  Cluster analysis of heterogeneous rank data.  In Proceedings of the 24th International Conference on Machine Learning, ICML ’07, New York, USA, pp. 113–120. ACM.
  * Caron et al. (2014) Caron, F., Y. W. Teh, and T. B. Murphy (2014).  Bayesian nonparametric Plackett-Luce models for the analysis of preferences for college degree programmes.  The Annals of Applied Statistics 8(2), 1145–1181.
  * Chen et al. (2007) Chen, Q. R., G. Vansant, K. Oades, M. Pickering, J. S. Wei, Y. K. Song, J. Monforte, and J. Khan (2007).  Diagnosis of the small round blue cell tumors using multiplex polymerase chain reaction.  The Journal of molecular diagnostics 9(1), 80–88.
  * Cramér (1938) Cramér, H. (1938).  Sur un nouveau théoreme-limite de la théorie des probabilités.  Actualités scientifiques et industrielles 736, 2–23.
  * Crispino and Antoniano-Villalobos (2022) Crispino, M. and I. Antoniano-Villalobos (2022).  Informative priors for the consensus ranking in the bayesian mallows model.  To appear, Bayesian Analysis.
  * Crispino et al. (2022) Crispino, M., L. Modugno, and C. Mollica (2022+).  The Mallows model with respondents’ covariates.  Preliminary draft available upon request.
  * DeConde et al. (2006) DeConde, R. P., S. Hawley, S. Falcon, N. Clegg, B. Knudsen, and R. Etzioni (2006).  Combining results of microarray experiments: A rank aggregation approach.  Statistical Applications in Genetics and Molecular Biology 5(1).
  * Dempster et al. (1977) Dempster, A. P., N. M. Laird, and D. B. Rubin (1977).  Maximum likelihood from incomplete data via the EM algorithm.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 39(1), 1–38.  with discussion.
  * Diaconis (1988) Diaconis, P. (1988).  Group representations in probability and statistics, Volume 11 of Lecture Notes–Monograph Series.  Hayward, USA: Institute of Mathematical Statistics.
  * Eliseussen et al. (2022) Eliseussen, E., T. Fleischer, and V. Vitelli (2022).  Rank-based Bayesian variable selection for genome-wide transcriptomic analyses.  Statistics in Medicine, 1–22.
  * Ellis (1984) Ellis, R. S. (1984).  Large deviations for a general class of random vectors.  The Annals of Probability 12(1), 1–12.
  * Feigin and Cohen (1978) Feigin, P. and A. Cohen (1978).  On a Model for Concordance Between Judges.  Journal of the Royal Statistical Society. Series B (Methodological) 40(2), 203–213.
  * Fligner and Verducci (1986) Fligner, M. A. and J. S. Verducci (1986).  Distance based ranking models.  Journal of the Royal Statistical Society B 48(3), 359–369.
  * Gärtner (1977) Gärtner, J. (1977).  On large deviations from the invariant measure.  Theory of Probability & Its Applications 22(1), 24–39.
  * Gormley and Murphy (2008a) Gormley, I. C. and T. B. Murphy (2008a).  Exploring voting blocs within the Irish electorate: A mixture modeling approach.  Journal of the American Statistical Association 103(483), 1014–1027.
  * Gormley and Murphy (2008b) Gormley, I. C. and T. B. Murphy (2008b).  A mixture of experts model for rank data with applications in election studies.  The Annals of Applied Statistics 2(4), 1452–1477.
  * Henery (1981) Henery, R. J. (1981).  Permutation probabilities as models for horse races.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 43(1), 86–91.
  * Irurozki (2014) Irurozki, E. (2014).  Sampling and learning distance-based probability models for permutation spaces.  Ph. D. thesis, Department of Computer Science and Artificial Intelligence of the University of the Basque Country.
  * Irurozki et al. (2016) Irurozki, E., B. Calvo, and A. Lozano (2016).  PerMallows: An R package for Mallows and generalized Mallows models.  Journal of Statistical Software 71(12), 1–30.
  * Jacques et al. (2014) Jacques, J., Q. Grimonprez, and C. Biernacki (2014).  Rankcluster: An R package for clustering multivariate partial rankings.  The R Journal 6(1), 101–110.
  * Kendall (1970) Kendall, M. G. (1970).  Rank correlation methods.  4th ed. Griffin London.
  * Khan et al. (2001) Khan, J., J. S. Wei, M. Ringnér, L. H. Saal, M. Ladanyi, F. Westermann, F. Berthold, M. Schwab, C. R. Antonescu, C. Peterson, and P. S. Meltzer (2001).  Classification and diagnostic prediction of cancers using gene expression profiling and artificial neural networks.  Nature Medicine 7(6), 673–679.
  * Lee and Yu (2010) Lee, P. H. and P. L. H. Yu (2010).  Distance-based tree models for ranking data.  Computational Statistics & Data Analysis 54(6), 1672–1682.
  * Lee and Yu (2012) Lee, P. H. and P. L. H. Yu (2012).  Mixtures of weighted distance-based models for ranking data with applications in political studies.  Computational Statistics & Data Analysis 56(8), 2486–2500.
  * Lee and Yu (2013) Lee, P. H. and P. L. H. Yu (2013).  An R package for analyzing and modeling ranking data.  BMC Medical Research Methodology 3(1), 1–11.
  * Liu et al. (2019) Liu, Q., A. Reiner, A. Frigessi, and I. Scheel (2019, 08).  Diverse personalized recommendations with uncertainty from implicit preference data with the Bayesian Mallows model.  Knowledge-Based Systems 186, 104960.
  * Mallows (1957) Mallows, C. L. (1957).  Non-null ranking models. I.  Biometrika 44(1/2), 114–130.
  * Marden (1995) Marden, J. I. (1995).  Analyzing and Modeling Rank Data, Volume 64 of Monographs on Statistics and Applied Probability.  Cambridge, USA: Chapman & Hall.
  * McCullagh (1993) McCullagh, P. (1993).  Models on spheres and models for permutations.  In M. A. Fligner and J. S. Verducci (Eds.), Probability Models and Statistical Analyses for Ranking Data, New York, pp. 278–283. Springer New York.
  * Meilǎ and Bao (2010) Meilǎ, M. and L. Bao (2010).  An exponential model for infinite rankings.  Journal of Machine Learning Research 11, 3481–3518.
  * Mollica and Tardella (2014) Mollica, C. and L. Tardella (2014).  Epitope profiling via mixture modeling for ranked data.  Statistics in Medicine 33(21), 3738–3758.
  * Mukherjee (2016) Mukherjee, S. (2016).  Estimation in exponential families on permutations.  The Annals of Statistics 44(2), 853–875.
  * Murphy and Martin (2003) Murphy, T. B. and D. Martin (2003).  Mixtures of distance-based models for ranking data.  Computational Statistics & Data Analysis 41(3–4), 645–655.
  * Murray et al. (2012) Murray, I., Z. Ghahramani, and D. MacKay (2012).  MCMC for doubly-intractable distributions.  arXiv preprint arXiv:1206.6848.
  * Qian and Yu (2019) Qian, Z. and P. L. H. Yu (2019).  Weighted distance-based models for ranking data using the R package rankdist.  Journal of Statistical Software 90(5), 1–31.
  * Schwarz (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  Ann. Statist. 6(2), 461–464.
  * Shi et al. (2021) Shi, C., B. Wei, S. Wei, W. Wang, H. Liu, and J. Liu (2021).  A quantitative discriminant method of elbow point for the optimal number of clusters in clustering algorithm.  EURASIP Journal on Wireless Communications and Networking 2021(1), 31.
  * Sloane (2017) Sloane, N. J. A. (2017).  The Encyclopedia of Integer Sequences.
  * Sørensen et al. (2020) Sørensen, Ø., M. Crispino, Q. Liu, and V. Vitelli (2020).  BayesMallows: An R Package for the Bayesian Mallows Model.  The R Journal 12(1), 324–342.
  * Varadhan (1966) Varadhan, S. S. (1966).  Asymptotic probabilities and differential equations.  Communications on Pure and Applied Mathematics 19(3), 261–286.
  * Vitelli et al. (2018) Vitelli, V., Ø. Sørensen, M. Crispino, A. Frigessi, and E. Arjas (2018).  Probabilistic preference learning with the Mallows rank model.  Journal of Machine Learning Research 18(158), 1–49.
  * Xu et al. (2018) Xu, H., M. Alvo, and P. L. H. Yu (2018).  Angle-based models for ranking data.  Computational Statistics & Data Analysis 121, 113–136.
  * Zhao and Hautamäki (2008) Zhao, Q. and V. Hautamäki (2008).  Knee Point Detection in BIC for Detecting the Number of Clusters.  In Blanc-Talon, J., Bourennane, S., Philips, W., Popescu, D., Scheunders, P. (eds) Advanced Concepts for Intelligent Vision Systems. ACIVS 2008\. Lecture Notes in Computer Science, vol 5259. Springer, Berlin, Heidelberg, pp. 664–673.

SUPPLEMENTARY MATERIAL

## A. Formal derivation of equation (12)

Here we give the details for the derivation of equation (12). We first define
the rescaled Spearman distance as

$x_{n}=\frac{d_{n}}{d_{\max,n}}=\frac{c_{n}-\boldsymbol{e}^{T}\boldsymbol{\sigma}}{\binom{n+1}{3}}\qquad\boldsymbol{\sigma}\in\mathcal{P}_{n}.$

For $n\to\infty$, we have that $x_{n}$ will be densely distributed in the unit
interval,
$x\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\lim_{n\to\infty}x_{n}\in\left[0,1\right]$.
The result of Gärtner and Ellis (Gärtner, 1977; Ellis, 1984) allows us to
write the rate function as follows

$\xi\left(x\right)=r\left(x\right)+o(1),$ (14)

with

$r\left(x\right)\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.}}}=\lim_{n\to\infty}\frac{1}{n}\log\mathbb{P}\left(x_{n}\right)$
(15)

independent of $n$ and $o(1)$ going to zero when $n\to\infty$. From the
quadratic nature of the Spearman distance, we can derive the property

$d(\boldsymbol{\sigma},\boldsymbol{e})+d(\boldsymbol{\sigma},\bar{\boldsymbol{e}})=d_{\text{max},n}$

with $\bar{\boldsymbol{e}}=\left(n,n-1,\dots,1\right)=(n+1)-\boldsymbol{e}$.
In addition, we can exploit the invariance property of the Spearman distance
to obtain

$d(\boldsymbol{\sigma},\bar{\boldsymbol{e}})=d\left(\left(\boldsymbol{\sigma}\circ\bar{\boldsymbol{e}}\right),\boldsymbol{e}\right),$

where the symbol $\circ$ indicates the composition of permutations. Given the
last two identities and the invertibility of the mapping
$\boldsymbol{\sigma}\circ\bar{\boldsymbol{e}}=(\sigma_{n},\sigma_{n-1},\dots,\sigma_{1})$
, we conclude that to every permutation $\boldsymbol{\sigma}$ at distance
${d}_{n}$ from $\boldsymbol{e}$, it corresponds one permutation at distance
$d_{\text{max},n}-{d}_{n}$ from $\boldsymbol{e}$. This implies
$N_{{d}_{n}}=N_{\left(d_{\text{max},n}-{d}_{n}\right)}$, and similarly for
$r(x)$ we have

$r(x)=r(1-x).$

The distribution $\mathbb{P}\left(x_{n}\right)$ is asymptotically unimodal and
symmetric, hence its maximum is attained for $x=\frac{1}{2}$. By definition of
the rate function (14) this implies $r\left(\frac{1}{2}\right)=0$. For the
tail values $x\to 0$ and $x\to 1$, we make the assumption of a logarithmic
divergence. By considering also the symmetry of $r(x)$, we obtain the tail
behaviour

$r(x)\sim\alpha_{1}\left[\log(x)+\log(1-x)\right]\quad x\to 0^{+},1^{-}.$

For the bulk of the distribution, we use a normal approximation

$r(x)\approx\alpha_{1}\left[\log(x)+\log(1-x)\right]+\alpha_{2}\,x\,(1-x).$
(16)

We can derive the value of $\alpha_{1}$ from the knowledge of the values of
$N_{{d}_{n}}$ for ${d}_{n}=2$, which give
$x=\binom{n+1}{3}^{-1}\approx\frac{1}{n^{3}}$. On the other hand, there is a
fraction of $\frac{n-1}{n!}$ permutations at distance $d_{n}=2$, which implies

$r(n^{-3})=\log(n)+o(1)$

or, equivalently, $\alpha_{1}=\frac{1}{3}$. By exploiting the fact that
$r\left(\frac{1}{2}\right)=0$ and equation (16), we derive the further
constraint $8\log(2)\alpha_{1}-\alpha_{2}=0$, so that only the sub-leading
corrections of order $o(1)$ in equation (14) remain to be estimated. To
perform this estimation, we exploit the known values of $N_{{d}_{n}}$ for
small $n$. First, we rewrite $\xi(x)$ as

$\displaystyle\xi(x)=$ $\displaystyle r(x)+\beta_{0}+$
$\displaystyle+\beta_{1}\left[\log(x)+\log(1-x)\right]+\beta_{2}x\,(1-x)$

with $\beta_{0},\beta_{1},\beta_{2}=o(1)$. Then, we proceed to estimate
$\beta_{0},\beta_{1},\beta_{2}$ with a generalized least squares fit. In
particular, we select $\beta_{0},\beta_{1},\beta_{2}$ in order to minimize the
cost function

$\sum_{n=n_{\text{min}}}^{n_{\text{max}}}n^{2}\sum_{d_{n}}\left[\xi\left(\frac{d_{n}}{d_{\text{max},n}}\right)-\frac{1}{n}\log\left(\frac{N_{d_{n}}}{n!}\right)\right]^{2},$

where $n_{\text{min}}$ and $n_{\text{max}}$ are parameters chosen on the basis
of the known values of $N_{d_{n}}$ in the OEIS sequence A175929 (Sloane,
2017). These values are available only for $n\leq 14$, and the first $n$ for
which there are no null values of $N_{d_{n}}$ in the bulk of the distribution
is $n=4$. Thus we select the values $n_{\text{min}}=4$ and
$n_{\text{max}}=14$. In addition, we modify the least squares minimization
with the $n^{2}$ prefactor to assign more weight to errors in the estimation
for high values of $n$, both because for these values we are closer to the
normal limit distribution, and because ultimately we are interested in
extending the estimate to values of $n$ greater than the ones for which the
$N_{d_{n}}$ series are available.

Finally, we obtain

$\xi(x)=a_{0}+a_{1}\left[\log(x)+\log(1-x)\right]+a_{2}\,x\,(1-x),$

with

$\displaystyle a_{0}=$ $\displaystyle\alpha_{0}\approx-\frac{0.24}{\sqrt{n}},$
$\displaystyle a_{1}=$
$\displaystyle\alpha_{1}+\beta_{1}\approx\frac{1}{3}-\frac{0.1784}{\sqrt{n}},$
$\displaystyle a_{2}=$
$\displaystyle\alpha_{2}+\beta_{2}\approx\frac{8}{3}\log(2)-\frac{5.5241}{\sqrt{n}}.$

The highest value of $n$ for which the series $N_{d_{n}}$ is available is
$n=14$. For this value we obtain a root mean square error equal to

$\sqrt{\frac{\sum_{d_{n}}\left[\xi\left(\frac{d_{n}}{d_{\text{max},n}}\right)-\frac{1}{n}\log\left(\frac{N_{d_{n}}}{n!}\right)\right]^{2}}{\binom{n+1}{3}}}=0.0027.$

## B. Derivation of the approximation $\hat{Z}_{\text{vMF}}(\theta)$

Let $n\in\mathbb{N}^{+}$, $a_{n}=(n+1)/2$, $b_{n}=n(n^{2}-1)/12$ and
$c_{n}=n(n+1)(2n+1)/6$. Let us denote with $\boldsymbol{r}$ a generic rankings
of $n$ items, with $\boldsymbol{e}=(1,\dots,n)$ the identity permutation and
with $\boldsymbol{y}$ and $\boldsymbol{\pi}$ their standardized versions, i.e.
$\boldsymbol{y}=(\boldsymbol{r}-a_{n})/\sqrt{b_{n}}$ and
$\boldsymbol{\pi}=(\boldsymbol{e}-a_{n})/\sqrt{b_{n}}$. Let $\mathcal{P}_{n}$
be the set of all $n!$ rankings and $\mathcal{S}_{n}$ the set of all $n!$
standardized rankings. We want to show that, by relying on in the result in Xu
et al. (2018), the normalizing constant of the MMS over $\mathcal{P}_{n}$ can
be written as

$Z(\theta)=e^{-\kappa}C(\kappa,\theta)^{-1}$

where $\kappa=2\theta b_{n}$ and

$C(\kappa,\theta)^{-1}=\sum_{\boldsymbol{y}\in\mathcal{S}_{n}}\exp\big{\\{}\kappa\boldsymbol{y}^{T}\boldsymbol{\pi}\big{\\}}$

is the normalizing constant of the MMS over $\mathcal{S}_{n}$. Note that

$\boldsymbol{r}^{T}\boldsymbol{e}=b_{n}\boldsymbol{y}^{T}\boldsymbol{\pi}+na_{n}^{2}\qquad\text{and}\qquad
b_{n}=c_{n}-na_{n}^{2}.$

Then, it follows

$\begin{split}Z(\theta)=&\sum_{\boldsymbol{r}\in\mathcal{P}_{n}}\exp\\{-\theta[2c_{n}-2\boldsymbol{r}^{T}\boldsymbol{e}]\\}\\\
=&\sum_{\boldsymbol{y}\in\mathcal{S}_{n}}\exp\\{-\theta[2c_{n}-2b_{n}\boldsymbol{y}^{T}\boldsymbol{\pi}-2na_{n}^{2}\\}\\\
=&\sum_{\boldsymbol{y}\in\mathcal{S}_{n}}\exp\\{-\theta[2b_{n}-2b_{n}\boldsymbol{y}^{T}\boldsymbol{\pi}]\\}\\\
=&\sum_{\boldsymbol{y}\in\mathcal{S}_{n}}\exp\\{-\kappa+\kappa\boldsymbol{y}^{T}\boldsymbol{\pi}\\}=e^{-\kappa}C(\kappa,\theta)^{-1}.\end{split}$

Since Xu et al. (2018, Appendix A) proposed the approximation

$C(\kappa,\theta)^{-1}\approx\frac{2^{\frac{n-3}{2}}n!I_{\frac{n-3}{2}}(\kappa)\Gamma(\frac{n-1}{2})}{\kappa^{\frac{n-3}{2}}},$

we can substitute it in the last expression of $Z(\theta)$, yielding

$Z(\theta)\approx\hat{Z}_{\text{vMF}}(\theta)=\frac{2^{\frac{n-3}{2}}n!I_{\frac{n-3}{2}}(\kappa)\Gamma(\frac{n-1}{2})}{\kappa^{\frac{n-3}{2}}e^{\kappa}}.$
(17)

Note that in our manuscript $\kappa=\theta n(n^{2}-1)/6=2\theta b_{n}$, while
in Xu et al. (2018), page 114 $\kappa=\lambda n(n^{2}-1)/12$. Hence, one can
recover the exact formula in Xu et al. (2018) by simply setting
$\theta=\lambda/2$.

$\theta$ | $\bar{\eta}_{\hat{\theta}_{\text{OEIS}}}$ | $\bar{\eta}_{\hat{\theta}_{\text{new}}}$ | $\bar{\eta}_{\hat{\theta}_{\text{vMF}}}$
---|---|---|---
1e-04 | .00387 | .00384 | .00409
5e-04 | .00349 | .00345 | .00369
.001 | .00309 | .00305 | .00329
.005 | .00159 | .00156 | .00178
.01 | .00173 | .00170 | .00183
.05 | .00209 | .00211 | .00215
.1 | .00217 | .00225 | .00369
.5 | .00885 | .00958 | .16272
1 | .01765 | .01776 | .84306
(a) $n=6$ $\theta$ | $\bar{\eta}_{\hat{\theta}_{\text{OEIS}}}$ | $\bar{\eta}_{\hat{\theta}_{\text{new}}}$ | $\bar{\eta}_{\hat{\theta}_{\text{vMF}}}$
---|---|---|---
1e-04 | .00151 | .00148 | .00158
5e-04 | .00115 | .00112 | .00123
.001 | .00083 | .00080 | .00091
.005 | .00042 | .00041 | .00045
.01 | .00045 | .00047 | .00047
.05 | .00088 | .00140 | .00386
.1 | .00126 | .00366 | .01397
.5 | .00662 | .02223 | .23848
1 | .01347 | .01590 | 1.06776
(b) $n=10$

$\theta$ | $\bar{\eta}_{\hat{\theta}_{\text{OEIS}}}$ | $\bar{\eta}_{\hat{\theta}_{\text{new}}}$ | $\bar{\eta}_{\hat{\theta}_{\text{vMF}}}$
---|---|---|---
1e-04 | .00115 | .00115 | .00119
5e-04 | .00080 | .00080 | .00085
.001 | .00057 | .00057 | .00062
.005 | .00032 | .00032 | .00034
.01 | .00034 | .00035 | .00040
.05 | .00060 | .00075 | .00686
.1 | .00109 | .00108 | .01912
.5 | .00479 | .01076 | .26687
1 | .01199 | .01534 | 1.15766
(c) $n=14$ $\theta$ | $\bar{\eta}_{\hat{\theta}_{\text{new}}}$ | $\bar{\eta}_{\hat{\theta}_{\text{vMF}}}$
---|---|---
1e-04 | .00083 | .00084
5e-04 | .00056 | .00058
.001 | .00032 | .00034
.005 | .00025 | .00025
.01 | .00030 | .00039
.05 | .00058 | .00776
.1 | .00248 | .02001
.5 | .00635 | .27603
1 | .01248 | 1.17938
(d) $n=16$

$\theta$ | $\bar{\eta}_{\hat{\theta}_{\text{new}}}$ | $\bar{\eta}_{\hat{\theta}_{\text{vMF}}}$
---|---|---
1e-04 | .00051 | .00051
5e-04 | .00030 | .00028
.001 | .00022 | .00019
.005 | .00028 | .00018
.01 | .00041 | .00063
.05 | .00265 | .00842
.1 | .00650 | .02111
.5 | .03650 | .28727
1 | .03655 | 1.21505
(e) $n=20$ $\theta$ | $\bar{\eta}_{\hat{\theta}_{\text{new}}}$ | $\bar{\eta}_{\hat{\theta}_{\text{vMF}}}$
---|---|---
1e-04 | .00042 | .00040
5e-04 | .00023 | .00019
.001 | .00019 | .00013
.005 | .00036 | .00023
.01 | .00039 | .00097
.05 | .00439 | .00878
.1 | .00938 | .02168
.5 | .05099 | .29718
1 | .11892 | 1.24676
(f) $n=24$

Table 11: Average absolute error $\bar{\eta}_{\hat{\theta}}$ for the
estimation of the concentration parameter associated to inferential procedures
based on alternative approximations of the Spearman distance ditribution.
Inference was performed on 100 samples simulated from the MMS with varying
values of $\theta$ (column 1) and different values of $n\leq 14$ (upper panel)
and $n>14$ (lower panel).

## C. Effect of the approximation on the inference

Here we study the effect of the approximation on the inference on $\theta$ by
means of simulated data. We simulate 100 datasets of $N=1000$ rankings from
the MMS with $n\in\\{6,10,14,16,20,24\\}$ and concentration
$\theta\in\\{0.0001,0.0005,0.001,0.005,0.01,0.05,0.1,0.5,1\\}$ and then
produce the following three estimates of $\theta$:

  1. 1.

$\hat{\theta}_{\text{OEIS}}$, obtained from the EM procedure, when $Z(\theta)$
and $\text{E}_{\theta}[D]$ have the exact form (this is available only for
$n\leq 14$);

  2. 2.

$\hat{\theta}_{\text{new}}$, obtained from the EM procedure, when $Z(\theta)$
and $\text{E}_{\theta}[D]$ are estimated by using our new approximation
$\hat{N}_{d_{n}}$;

  3. 3.

$\hat{\theta}_{\text{vMF}}$, obtained from the procedure described in Xu et
al. (2018), that is, a Newton-based recursive approximation solution of (7)
when $Z(\theta)$ is approximated with $\hat{Z}_{\text{vMF}}(\theta)$, see Xu
et al. (2018), page 116.

In Table 11, we report the average absolute error
$\bar{\eta}_{\hat{\theta}}=\sum_{s=1}^{100}|\hat{\theta}_{s}-\theta|/100$ over
the 100 samples to evaluate the performance of the three estimates for $n\leq
14$ (upper panel) and $n>14$ (lower panel). From Table 11, we observe that
$\bar{\eta}_{\hat{\theta}_{\text{new}}}$ is smaller than
$\bar{\eta}_{\hat{\theta}_{\text{vMF}}}$ for all cases $n\leq 14$ considered,
and the difference among the two errors is quite small for low values of
$\theta$, but it grows fast for increasing values of the concentration
parameter. This result, in line with the findings described in Section 3.2,
holds generally true also for larger values of $n$, as one can notice in the
lower panel of Table 11. However, in this case, we obtained that for very
small values of $\theta$, $\bar{\eta}_{\hat{\theta}_{\text{vMF}}}$ is
sometimes smaller than $\bar{\eta}_{\hat{\theta}_{\text{new}}}$ The latter
happens for instance when $n=20,24$, and $\theta\leq 0.005$ (see Table 11 (e)
and (f)).