Modelling Preference Data with the Wallenius Distribution
Clara Grazian Fabrizio Leisen Brunero Liseo








 Modelling Preference Data with the Wallenius Distribution
    
 
Clara Grazian1 
Fabrizio Leisen2 
Brunero Liseo3
        
 
        

        

        1
        University of Oxford, U.K.   
         2 
        University of Kent, U.K. 
         

        3
         Sapienza Università di Roma, Italy 
         
 

     
    
=========================================================================================================================================================================================================================================================================







 
The Wallenius distribution is a generalisation of the Hypergeometric
distribution where weights are assigned to balls of different colours.
This naturally defines a model for ranking categories which can be
used for classification purposes. Since, in general,
the resulting likelihood is not analytically available, we adopt an approximate Bayesian computational (ABC) approach 
for estimating the importance of the categories. We illustrate the performance of the estimation procedure on simulated datasets. Finally, we use the new model for analysing two datasets concerning movies ratings and Italian academic statisticians' journal preferences. The latter
is a novel dataset collected by the authors.

Keywords: Approximate Bayesian Computation, Biased Urn, Movies ratings, Scientific Journals Preferences. 







§ INTRODUCTION AND MOTIVATIONS



Human beings naturally tend, in everyday life, to compare and rank concepts and objects such as food, shops, singers and football teams, according to their preferences. In general, to rank a set of objects means to arrange them in order with respect to some characteristic. Ranked data are often employed in contexts where objective and precise measurements are difficult, unreliable, or even impossible to obtain and the observer is bound to collect ordinal information about preferences, judgments, relative or absolute ranking among competitors, called items. 

Modern web technologies have made available a huge amount of ranked data, which can provide information about social and psychological behaviour, marketing strategies and political preferences. The codification of this information has been of interest to the statisticians since the beginning of the 20th century. The Thurstone model (TM) assumes that each
item i is associated with a score W_i on which the comparative judgment is based; examples of unidimensional scores are the unrecorded finishing times of players in a race or any possible preference/attitude measure towards items.
Item i is preferred to item j if W_i is greater than W_j, see <cit.>. From the modelling point of view, this corresponds to assigning a probability p_ij=(W_i > W_j). 
The Bradley-Terry model (BT) is a particular case of the TM model with p_ij=p_i(p_i+p_j)^-1 where p_i,p_j≥ 0  are the item parameters reflecting the rate of each item, see <cit.>. Paired comparison models are always applicable to rankings after converting the latter in a suitable set of pairwise preferences. Conversely, paired comparisons of K items do not necessarily correspond to a ranking, due to the potential presence of circularities. A popular extension of the BT model is the Plackett-Luce model (PL). 
Given a set of L items and a vector of probabilities (p_1,…, p_L), such that ∑_i=1^L p_i=1, the PL model assigns a probability distribution on all the set of possible rankings of these objects which is a function of the (p_1,…, p_K), see <cit.> and <cit.>. TM, BT and PL are not the only proposals in the field, and modelling ranking is an active area of research, see <cit.> and <cit.>. 


There is no wide consensus about the use of choice or ranking data for better representing preferences and, very often, the best solution is problem specific. In this paper, we consider a sort of hybrid situation; in fact, we assume that choices related to single items can be further classified into categories of different relevance, and the ranking of categories is the main goal of the statistical analysis. 
Our approach makes use of an extension of the Hypergeometric distribution, namely the Wallenius distribution <cit.> and can be used in the cases where data are available in the form of rankings, votes, preferences of items but the interest is in defining the importance of the categories in which the items can be 
clustered.

The Wallenius distribution arises quite naturally in situations where sampling is performed without replacement and units in the population have different probabilities to be drawn. To be more specific, consider a urn with balls of c different colours: for i = 1, …, c there are m_i balls of colour i. In addition, colour i has a priority ω_i>0 which specifies its relative importance with respect to the other colours. A sample of n balls, with n < ∑_i=1^c m_i, is drawn sequentially without replacement. The Wallenius distribution describes the probability distribution for all possible strings of balls of length n drawn from this urn. 
This experimental situation arises in very different contexts. For example, in auditing problems, transactions are examined by randomly selecting a single euro (or pound, or dollar) among the total amount, so larger transactions are more likely to be drawn and checked. 

The Wallenius distribution was introduced by <cit.> and it is also known as the noncentral Hypergeometric distribution; this alternative name is justified by the fact that, when all the priorities ω_i's are equal, one gets back to the classical Hypergeometric distribution. However this name should be avoided because, as extensively discussed by <cit.>, this is also the name of another distribution, proposed by <cit.>.
Although the Wallenius distribution is a very natural statistical model for the aforementioned situations, its popularity in applied settings has been prevented by the lack of a closed form expression of the probability mass function: see Section <ref> for details.

The gist of this paper is the use of the priorities vector ω=(ω_1, …, ω_c) of the Wallenius distribution as a measure of importance for different values of a categorical variable. 

In particular, we analyse two datasets, where we aim at ranking the categories rather than the items.
The first dataset considers data downloaded from the MovieLens website, which consists of 105,339 ratings across 10,329 movies performed by 668 users.
In this framework, it is of interest to classify the different genres in terms of  satisfaction, in order to provide some useful feedback to users and/or providers.

The second dataset considers data we collected between October and November 2016 among Italian academic statisticians.
They indicated their journal preferences from the 2015 ISI “Statistics and Probability” list of Journals. In this context, we are interested in ranking the journal categories in order to provide a description of the research interests of the Italian Statistical community.

We adopt a Bayesian methodology which allows us to overcome the computational problems related to the lack of a closed form expression of the probability mass function of the Wallenius distribution.
We propose a novel approximate Bayesian computational approach <cit.>, where the vector of summary statistics is represented by the relative frequencies of the different categories and the acceptance mechanism is based on the distance in variation <cit.>

The paper is organized as follows: in Section <ref> we introduce the Wallenius distribution; in Section <ref> our approximated inferential strategy is described, based on an ABC algorithm. The performance of the algorithm has been tested in several examples, first in an extensive simulation study (Section <ref>) and then on two real datasets (Section <ref>). A discussion concludes the paper. 






§ THE WALLENIUS DISTRIBUTION



Consider an urn with N balls of c different colours.  There are m_i balls of the i-th colour, so that ∑_i^c m_i=N. In this situation, the multivariate Hypergeometric distribution is the discrete probability distribution which describes the sampling without replacement of n balls. In this framework, the probability of drawing a ball of a certain colour is proportional to the number of balls of the same colour. It is possible to generalise the experiment with a biased sampling of balls. For instance, each colour may have a different priority or importance, say ω_i>0, i=1, …, c. Suppose we have drawn n balls without replacement from the urn and let X_n= (X_1n, X_2n, …, X_cn) denote the frequencies of balls of different colours in the sample. Let  Z_n be the colour of the ball drawn at time n. 
In this setting, the probability that the next ball is of colour i also depends on its priority and is defined as


    P( Z_n+1 = i |X_n  ) =  (m_i-X_in) ω_i/∑_j=1^c ( m_j - X_jn) ω_j.
  

<cit.> provided the above expression and the probability mass function of X_n for the case c=2. <cit.>
derived the following general expression. For a given integer n, and parameters
m= (m_1, …, m_c) and ω= (ω_1, …, ω_c),
the probability of observing a vector of colour frequencies x=(x_1, … , x_c) is 

    P(x; n, m, ω) = ∏_j=1^c m_jx_j∫_0^1 
    ∏_j=1^c  ( 1- t^ω_j/d )^x_j dt,


where ∑_i=1^c x_i=n and d= ∑_j=1^c ω_j (m_j - x_j). When ω_i = ω, for every i=1,…, c, the Wallenius distribution reduces to the multivariate Hypergeometric distribution. This can be easily shown by considering, without loss of generality, ω=1 and c=2. In this particular case, the probability mass function simplifies to  

    P(x; n, m )= mxN-mn-x∫_0^1  ( 1- t^1/d )^ndt.

The change of variable z=t^1/d leads to

    P(x; n, m )   = mxN-mn-x
    d ∫_0^1  ( 1- z  )^n   z^d-1  dz 
       =mxN-mn-xΓ(d+1) Γ(n+1)/Γ(n+d+1).


Since d=N-n, the probability mass function reduces to

    P(x; n, m )= mxN-mn-x /Nn,

which is the probability mass function of the Hypergeometric distribution when two colours are considered. 


The Wallenius distribution has been underemployed in the statistical literature mainly 
because the integral appearing in (<ref>) cannot be solved in a closed form and numerical approximations are necessary. <cit.> has made a substantial contributions in this direction, providing approximations based either on asymptotic expansions or numerical integration. To our knowledge, the Wallenius distribution has only been used in a limited number of applications, mainly devoted to auditing problems <cit.>, ecology <cit.>, vaccine efficacy <cit.> and modeling of RNA sequences <cit.>. 

In this work, we propose a novel look at the Wallenius distribution and we use it as statistical model, with the goal of ranking the values of a categorical random variable, based on preference data.

This is motivated by the sampling nature of the Wallenius distribution where an importance ω_j is associated with category j. The highest ω_j's represent the most popular categories. This naturally defines a new model which allows us to rank preferences.

Notice that we are implicitly assuming that all balls of the same colour have the same importance; this may not be the case in some applications: we will discuss this aspect in the final section.

 Recently, the development of social networks and the competitive pressure to provide customized services has motivated many new ranking problems involving hundreds or thousands of objects. 

Recommendations on products such as movies, books and songs are typical examples in which the number of objects is extraordinarily large. In recent years, many researchers in statistics and computer science have developed models to handle such big data. For instance, in Section <ref> we consider the problem of ranking customer movie choices in terms of genres  such as Comedy, Drama and Science Fiction. We consider data downloaded from the MovieLens website  which consists of 105,339 online ratings of 10,329 movies by 668 raters on a scale of 1-5. We rank the categories by estimating the priority parameters of the Wallenius distribution by using an approximate Bayesian approach.

In particular, in the next section, we introduce a simple ABC algorithm which allows us to avoid the direct computation of the integral in equation (<ref>). 






§ BAYESIAN INFERENCE FOR THE WALLENIUS MODEL



Let x_h=(x_h1, …, x_hc) be a draw of n_h balls from the Wallenius urn described in equation (<ref>), 
where h=1, …, k and ∑_j=1^c x_hj=n_h.
In this paper we adopt a Bayesian approach, where the parameter vector ω
is considered random. For a given prior distribution π(ω), the resulting posterior 
is

    π(ω|x_1, …, x_k ) ∝π(ω)
    ∏_h=1^k 
     [ ∫_0^1 ∏_j=1^c  ( 1- t_h^ω_j/d_h )^x_hj dt_h
     ],


with d_h = ∑_j=1^c ω_j (m_j - x_hj).
Here k represents the sample size, that is, the number of different and conditionally
independent preference lists provided by the interviewees, while n_h (h=1, …, k) is the number of items 
selected by the h-th interviewee. 
The above posterior distribution depends on 
k different integrals which cannot be reduced to a closed form. 
This makes the implementation of standard Markov Chain Monte Carlo (MCMC) methods for estimating ω rather complex.
Indeed, most MCMC methods rely on the direct evaluation of the unnormalized posterior distribution (<ref>).
Although there are many available routines, in different software packages, to evaluate univariate integrals, we noticed that they lack accuracy especially for large values 
of the n_h's and m.
We believe that this problem has had a strong negative impact   
on the popularization of the Wallenius distribution despite a need for 
interpretable models in the applied setting. For instance, the Wallenius distribution arises naturally in genetics as an alternative to the Fisher exact test, see <cit.> and the references therein. 

In this section, we propose an algorithm which allows to sample from the posterior distribution introduced in (<ref>). The algorithm belongs to the class of approximate Bayesian computational (ABC) methods. This approach is philosophically different from the standard MCMC methods since the implementation only requires to draw samples from the generating model for a given parameter value. In the case of the Wallenius distribution, the task of generating draws is not hard, making the use of ABC particularly straightforward. <cit.> provided methods and algorithms to sample from the Wallenius distribution. He also made available a reliable R package, called , which has been used extensively in this work.  

The ABC methodology can be considered as a (class of) popular algorithms that achieves posterior simulation by avoiding the computation of the likelihood function: see <cit.>, <cit.> and <cit.> for recent surveys. As remarked by <cit.>, the first genuine ABC algorithm was
introduced by <cit.> in a population genetics setting.
Explicitly, we consider a parametric model {f(·|θ), θ∈Θ} and suppose that a dataset
y∈𝒟⊂ℝ^n is observed. Let
ε>0 be a tolerance level, η a summary statistic
(which is often not sufficient) defined on 𝒟 and ρ a
distance or metric acting on the η space. 
Let π be a prior distribution for θ; the ABC algorithm is described in Algorithm 1.



The basic idea behind the ABC is that, for a small (enough)
ε and a representative summary statistic, we can obtain a
reasonable approximation of the posterior distribution. The practical implementation of an ABC algorithm requires the selection of a suitable summary statistic, a distance and a tolerance level. In our specific case we summarized the data by using the arithmetic mean of the 
observed and simulated frequency vectors, i.e., at the ℓ-th iteration of pseudo data generation, we have  


    η(x^(ℓ)) =p^(ℓ) = 1/k∑_h=1^k p^(ℓ)_h,
 
with 

    p^(ℓ)_h =  (x_h1^(ℓ)/n_h, …, x_hc^(ℓ)/n_h )

to be compared with the relative frequencies observed in the sample

    η(x^(t))=p^(t)= 1/k∑_h=1^k p^(t)_h.

with

    p^(t)_h = ( x_h1/n_h, …, x_hc/n_h ).

Since the frequencies p^(ℓ)=( p_1^(ℓ),…, p_c^(ℓ)) and p^(t)= ( p_1,…, p_c) can be interpreted as discrete probability distributions, it is natural to compare them through the “distance in variation” <cit.> metrics

    ρ(p^(ℓ),p^(t))=1/2∑_j=1^c |p_j^(ℓ)-p_j|


Regarding the setting of tolerance level we refer to the Section <ref> where the algorithm will be tested on simulated data.



  §.§.§ The prior distribution


The vector of parameters ω= (ω_1, …, ω_c) assumes values in ℝ_+^c and different priors can be considered. However, one must take into account that the priority parameters ω_j must be interpreted in a relative way. 
In fact, the quantity d in the  p.m.f. of the  Wallenius distribution (defined in equation (<ref>)) depends on the priority parameters ω. In particular,

    d= ∑_j=1^c ω_j (m_j - x_j).

If we consider two different vectors ω' and ω such that ω' = κω for κ > 0, we have that 

    ω_j'/d'=κω_j/∑_j=1^c κω_j (m_j - x_j)=ω_j/∑_j=1^c ω_j (m_j - x_j)=ω_j/d

where d' and d are computed respectively with ω' and ω. Equation (<ref>) implies that the p.m.f. of the Wallenius distribution does not change if we consider the vector of priorities ω' instead of ω. This induces an identifiability issue, which can be resolved by a normalization step. 
From this perspective, the most natural way to follow is to assume that ∑_j=1^cω_j=1, and to assume a Dirichlet prior on the normalized vector. Hereafter we will assume that the Dirichlet prior we adopt in the simulations and the real data examples are symmetric (i.e., all the hyperparameters are equal). Our default choice will be to set them all equal to 1, making the prior uniform on its support.
An alternative default choice, especially useful when c is large, is given by α = 1/c, as explained in <cit.>.
 
 


  §.§.§ Alternative computational approaches


The R package  allows the approximate numerical evaluation of the probability mass function of the Wallenius distribution. In a classical setting, this makes feasible the computation of the MLE. 
In a Bayesian setting this enables the implementation of standard MCMC algorithms, such as the Metropolis-Hastings sampler. 
Nonetheless, we deem more appropriate to use the ABC approach illustrated in this section for several reasons. 
First, the output of the Bayesian approach is far richer than the one available in a classical setting. For instance, in Section <ref> we are able to easily compute important summaries of the posterior distribution, i.e. the probability p_ij=(ω_i>ω_j). 
Second, standard MCMC methods require repeated evaluations of the likelihood function. This could lead to an unsustainable computational burden compared to ABC. 
Last but not least, we have performed a simulation study regarding the behaviour of the maximum likelihood estimator of the vector 
ω and we noticed that it typically tends to produce unreliable and unstable estimates when the “true” ω is close to the boundary of the simplex and/or when the number of categories is large.


 







 





§ SIMULATION STUDY



In order to test Algorithm §<ref> with the summary statistics shown in Section <ref>, we have conducted an extensive simulation study, with different scenarios. We performed 20 repeated simulations of k draws from the Wallenius distribution where each draw consists of a number n_h (h=1, ⋯, k) of balls. We use the prior distribution defined in Section <ref>, i.e. a Dirichlet prior 𝒟ir(1,…,1). As already stated in Section <ref>, we use the summary statistics and the distance in variation defined in equations (<ref>) and (<ref>). The tolerance level ε has been chosen with a pilot simulation where 10^5 values have been simulated by fixing the tolerance level to a very large value. Then, the distribution of the distances from the true values has been studied. The tolerance level is fixed as a small quantile of this distribution (it is common practice to fix it as the quantile of level 0.05). The complete procedure will be described in the following.



The simulated experiments have been performed for different values of c, ranging between 2 and 20, and using three configurations for both m and ω, as explained below: 


  * same number of balls for each colour, i.e. m_j=m, j=1, …,c;  uniform importance weights, i.e. ω_j=ω, j=1, …,c;  

  * increasing values for m_j's (all the integers between 1 and c) and ω's (all the integers between 1 and c, normalized to sum to one), j=1, …, c;

  * increasing values for m_j's (all the integers between 1 and c) and decreasing values for the ω's (all the integers between c and 1, normalized to sum to one), j=1, …, c;

Finally, we have used three different sample sizes, namely k=5, k=50 and k=1000. 
The value of n_h's has been taken to be half the total number of balls in the urn. The results are available in Tables <ref>, <ref> and <ref>. 

Surprisingly, as the sample size k increases, the root mean squared error (RMSE) remains relatively stable. 

Results are less accurate for those configurations where both ω and m are uniform, while they are more accurate for configurations where ω and m follow an opposite ordering. This may be explained by observing that data are carrying more information on ω in this situation. 

The RMSE is decreasing almost everywhere as the value of c increases: the only case where this is not true is the case of both ω and m uniform. 

This may suggest that the Wallenius distribution does not perform well when the “true” model is the simpler classical multivariate Hypergeometric model, especially when the number of categories c is large.

Table <ref>, <ref> and <ref> also show the average acceptance rates of the ABC algorithm used in the simulation experiments. The acceptance rate depends on the value of the tolerance level ϵ chosen in the experiment: we have followed the strategy described in <cit.>, where a pilot run is done to study the distribution of the distance between the summary statistics computed on the observed data and on the simulated data. Then, ε is chosen to be a quantile of the empirical distribution of this distance. We have chosen to consider the quantile of level 0.05. With this automatic choice of ε we obtain an acceptance rate of about 0.01-0.02 on average. We obtained lower acceptance rates in the case of a small number of colours. These rates are compatible with the average tolerance level. It could be possible to reduce the RMSE by reducing the tolerance level ε, however there is a balance between the goodness of the approximation and the computational cost. In an applied context, it is always advisable to compare several tolerance levels. We will propose this comparison in Section <ref>. In this context, we use only one threshold ε (in the automatic way above described) to focus the analysis on a Monte Carlo comparison by varying the sample size and the number of colours in the urn.

As a conclusive remark of the section, we have performed a sensitivity analysis regarding the common hyperparameter of the Dirichlet prior. For values ranging from 1/c (the choice suggested in <cit.>)  and 1 (the uniform prior), we have always obtained similar results in terms of RMSE, showing a sort of robustness of the model, at least with respect to this particular aspect.













§ REAL DATA APPLICATIONS



We now apply the proposed approach to two real datasets, in order to assess the applicability and the performance of the algorithm. 
In both cases, we obtain the ratings of a group of individuals about specific elements from a list. Each individual may choose the number of elements to rate. The elements are then grouped in categories and the goal is to provide a ranking of the categories. By using the urn terminology of Section <ref>, the categories are the colours and each element from the list is a ball; the aim of the analysis is to perform inference on the importance weights of each colour.  



 §.§ Movies dataset



This dataset describes 5-star (with half-star increments) rating from MovieLens, a movie recommendation service (http://grouplens.org/datasets/movielens/). The dataset may change over time. We consider the dataset which contains 105,339 ratings across 10,329 movies. These data were created by 668 users between April 03, 1996 and January 09, 2016. This dataset was generated on January 11, 2016. Users were randomly selected by MovieLens, with no demographic information, and each of them has rated at least 20 movies. The movies in the dataset were described by genre, following the IMDb information (https://www.themoviedb.org/); nineteen genres were considered in the dataset, including a “no genre” category; we have decided to eliminate the empty category from the analysis. In this case, we consider a movie to be "good" if its rating is at least 3.5 stars. Therefore, the vector X_n represents the frequencies of "good movies" in each category. Each film may be described by more than one genre. In this case we have proceeded as follows: we have ordered the genres in terms of their generality and then assigned to the movie the least general genre with which it was described. We have decided the following order (from the less general to the most general): Animation → Children → Musical → Documentary → Horror → Sci-Fi → Film Noir → Crime → Fantasy → War → Western → Mistery → Action → Thriller → Adventure → Romance → Comedy → Drama. Of course, this is an experimental choice, which may affect the results. Since the movies can be cross-classified, an interesting (and more realistic) development would be considering a model which can take into account this feature; this is left for further research. We have then replicated the same prior choice and the same choices of distance and vector of summary statistics described in Section <ref>. The tolerance level ε has been chosen with a pilot simulation in order to produce a sample of size 10^5, as described in Section <ref>. In this particular case, we have used ε=0.5. Table <ref> displays the posterior mean estimates of the vector of importance weights ω. The importance weights seem to be very close, with small differences among them. This suggests that there is not a category which is particularly popular. Nonetheless, we can observe a slightly preference for the Action and Sci-Fi genres and less interest in the Fantasy, War and Drama genres. We believe that this similarity in the importance weights is due to an excessive number of categories in the movies dataset. In this setting the graphical comparison of the marginal posterior distributions can provide a better insight on the customer preferences. Figure <ref> shows that there is more variability in the users preferences to choose a particular movie genre, such as Action or Romance.







 §.§ Statistical Journals dataset



The scientific areas (or “settori scientifici disciplinari”, S.S.D.) are a characterization used in the academic Italian system to classify knowledge in higher education. 
The sectors are determined by the Italian Ministry of Education. In particular, there are 367 S.S.D., divided into 14 macro-areas and each member of the academic staff pertains to a single sector. We have performed a survey on the preferences of the researchers in Statistics (Sector SECS-S/01) of Italian universities about the available scientific journals. It should be noted that researchers in Probability and Mathematical Statistics, Medical, Economic and Social Statistics are not included in this survey, because they pertain to different sectors. We have considered only staff with both teaching and research contracts. Postdoctoral fellows and PhD students have been excluded. In this survey we have used the 2015 “Statistics and Probability” list of journals of the Institute for Scientific Information (ISI). We have asked to SESC-S/01 researchers to indicate their preferences in this list, between a minimum of ten and a maximum of twenty. One difference from the Movies example of Section <ref> is that the participants do not have to indicate the level of their preference, only a list of journals which each of the participants considers either


	
  * prestigious and/or
	
  * likely for a potential submission and/or
	
  * professionally significant (in terms of frequency of readings).


The survey was conducted between 25th October 2016 and 4th November 2016. We have collected 174 responses, distributed, in terms of role, as follows: 49 Full professors (Professori Ordinari), 72 Associate Professors (Professori Associati) and 53 Assistant Professors, both fixed-term and tenure-track (Ricercatori a tempo indeterminato e a tempo determinato). We have then grouped the journals by category, considering five main classes of interest: Methodology, Probability, Applied Statistics, Computational Statistics and Econometrics and Finance. The list of journals and relative category is available in the Appendix. Among the 124 journals available in the “Statistics and Probability” ISI list, we have classified 23 journals in Probability, 45 in Methodology, 34 in Applied Statistics, 9 in Computational Statistics and 13 in Econometrics and Finance. We assume the Wallenius distribution for modelling the dataset, where c represents the number of the categories. The preferences of each respondent are summarized in a vector where the position of each entry represents the number of journals falling in the corresponding category. We consider that this vector is a realization of the Wallenius distribution. 



The results are available in Figure <ref>, Figure <ref> and Table <ref>, which show that there seems to be a preference for the research in Methodological and Applied Statistics among the researchers in Statistics and less interest in journals of Probability. As already stated, this should highlight the fact that researchers in Mathematical Statistics and Probability do not pertain to the investigated sector. 
These results also show that the effect of a decrease of the tolerance level seems to be a concentration of the posterior distributions of the importance weights ω, except for the weight relative to the Computational journals, for which there is a shift.  As a possible explanation of this fact, one should consider that this category is under-represented in the list (at least, according our classification) with respect to the others. Table <ref> shows the estimated pair comparison probabilities for the journal categories.  











§ DISCUSSION



In this paper we have considered the problem of ranking categories of items.  
We have proposed a novel model based on the Wallenius distribution. In terms of an urn scheme, it generalizes the Hypergeometric distribution with an additional vector of parameters ω, which represents the importance of the different types of balls in the urn. 

A referee noticed that “the model assumes that the balls of the same colours (eg. the journals in the same category) are equally likely to be drawn.”  This assumption may not be justified, since, in the Journal example, journals in the same category may have different standing. This is exactly the reason why we propose the Wallenius model for ranking categories rather than single items; the weight ω refers to the entire categories and they do not discriminate within categories. However, it is certainly of scientific interest to pursue the above issue and to conceive a nested model 
where items might be further ranked within categories; see, for example, <cit.>.  In a Bayesian nonparametric setting, this approach could be further generalized by using nested non-exchangeable species sampling sequences, see <cit.> and <cit.>.  

So far the Wallenius model has been definitely under-employed, due to the analytical intractability of the probability mass function. In this work we proposed an approximate Bayesian computational algorithm which provides a fast and reliable approach to the estimation of the vector of priorities ω. Our method is easy to implement and it might be very useful in several statistical applications where balls are drawn from the urn in a biased fashion. Paradigmatic examples of the importance of the Wallenius model especially appear in auditing where transactions are randomly checked with probability proportional to their monetary value. We analysed two datasets concerning movies ratings and Italian academic statisticians' journal preferences. The ABC algorithm allows us to estimate the importance of movies categories or journal preferences under the assumption of a Wallenius generating model. Future work will focus on the use of the Wallenius distribution to other areas of application and on the estimation of the category multiplicities m given the knowledge of the importance weights ω.



§ ACKNOWLEDGEMENTS

The authors are very grateful to Martin Ridout for his valuable comments on a first draft of the paper. This project has been funded by the Royal Society International Exchanges Grant “Empirical and Bootstrap Likelihood Procedures for Approximate Bayesian Inference”. F Applied Statistics, 9 in Computational Statistics and 13 in Econometrics and Finance. We assume the Wallenius distribution for modelling the dataset, where c represents the number of the categories. The preferences of each respondent are summarized in a vector where the position of each entry represents the number of journals falling in the corresponding category. We consider that this vector is a realization of the Wallenius distribution. 
abrizio Leisen was supported by the European Community's Seventh Framework Programme [FP7/2007-2013] under grant agreement no: 630677.








20
natexlab#1#1

[Airoldi  et al.(2014)]Airoldi
Airoldi, E., Costa, T., Bassetti, F., Leisen, F.  and Guindani, M. (2014).
Generalized species sampling priors with latent beta reinforcements.
Journal of the American Statistical Association 109, 1466–1480.

[Allingham et al.(2009)]allingham2009bayesian
Allingham, D., King, R.A.R., and Mengersen, K.L. (2010).
Bayesian estimation of quantile distributions. Applied Statistics, 9 in Computational Statistics and 13 in Econometrics and Finance. 
Statistics and Computing
  19(2), 189–201.

[Alvo and Yu(2014)]alvo
Alvo, M. and Yu, P. L. H. (2014).
Statistical Methods for Ranking Data.
Springer, New York.

[Bassetti, Crimaldi and Leisen (2010)]Bassetti
Bassetti, F., Crimaldi, I. and Leisen, F. (2010).
Conditionally identically distributed species sampling sequences.
Advances in Applied Probability 42, 433–459.

[Beaumont(2010)]bea:10
Beaumont, M. (2010).
Approximate Bayesian computation in evolution and Applied Statistics, 9 in Computational Statistics and 13 in Econometrics and Finance..
Annual Review of Ecology, Evolution, and Systematics
  41, 379–406.

[Berger et al.(2015)Berger, Bernardo and Sun]berger2015
Berger, J. O., Bernardo, J. M. and Sun, D. (2015).
Overall objective priors.
Bayesian Analysis 10, 189–221.

[Bradley and Terry(1952)]bt52
Bradley, R. A. and Terry, M. E. (1952).
Rank analysis of incomplete block designs. I: The method of paired
  comparisons.
Biometrika 39, 324–345.

[Bremaud(1998)]bremaud
Bremaud, P. (1998).
Markov Chains: Gibbs Fields, Monte Carlo Simulation
  and Queues.
Springer-Verlag: New York.

[Chesson(1976)]chesson
Chesson, J. (1976).
A non-central multivariate Hypergeometric distribution arising from
  biased sampling with application to selective predation.
Journal of Applied Probability 13, 795–797.

[Fisher(1935)]fisher1
Fisher, R. (1935).
The logic of inductive inference.
Journal of the Royal Statistical Society 98, 39–82.

[Fog(2008a)]fog1
Fog, A. (2008a).
Calculation methods for Wallenius' Noncentral Hypergeometric
  Distribution.
Communications in Statistics - Simulation and Computation
  37, 258–273.

[Fog(2008b)]fog2
Fog, A. (2008b).
Sampling methods for Wallenius' and Fisher's Noncentral
  Hypergeometric Distributions.
Communications in Statistics - Simulation and Computation
  37, 241–257.

[Gao et al.(2011)Gao, Fang, Zhang, Zhi and Cui]gao11
Gao, L., Fang, Z., Zhang, K., Zhi, D. and
  Cui, X. (2011).
Length bias correction for RNA-seq data in gene set analyses.
Bioinformatics 27, 662–669.

[Gillett(2000)]gillett
Gillett, P. R. (2000).
Monetary unit sampling: a belief-function implementation for audit
  and accounting applications.
International Journal of Approximate Reasoning 25,
  43–70.

[Hernández-Suárez and Castillo-Chavez(2000)]castillo2000urn
Hernández-Suárez, C. M. and Castillo-Chavez, C. (2000).
Urn models and vaccine efficacy.
Statistics in Medicine 19, 827–835.


[Inskip et al.(2013)]inskip
Inskip, C., Ridout, M., Fahad, Z., Tully, R., 
	Barlow, A., Greenwood Barlow C., Islam, M.A., Roberts, T., 
	MacMillan, D. (2013).
Human–Tiger Conflict in Context: Risks to Lives
and Livelihoods in the Bangladesh Sundarbans.
Human Ecology  41, 169–186.

[Karabatsos and Leisen (2018)]KL2018
Karabatsos, G. and Leisen, F. (2018).
An approximate likelihood perspective on ABC methods.
Statistics Surveys 12, 66–104.

[Luce(1959)]luce
Luce, R. D. (1959).
Individual Choice Behavior: A Theoretical Analysis.
John Wiley and Sons Inc., New York.

[Manly(1974)]manly
Manly, B. J. (1974).
A model for certain types of selection experiments.
Biometrics 30(2), 281–294.

[Marden(1995)]marden
Marden, J. (1995).
Analyzing and Modeling Rank Data.
Chapman and Hall, London.

[Marin et al.(2012)Marin, Robert and Pudlo]marin-abc
Marin, J. M., Robert, C. P. and Pudlo, P. (2012).
Approximate Bayesian computational methods.
Statistics and Computing 22, 1167–1180.

[Plackett(1975)]plack
Plackett, R. L. (1975).
The analysis of permutations.
Journal of the Royal Statistical Society  Series C 24,
  193–202.

[Pritchard et al.(1999)Pritchard, Seielstad, Perez-Lezaun and
  Feldman]Pritchard
Pritchard, J., Seielstad, M., Perez-Lezaun, A. and
  Feldman, M. (1999).
Population growth of human Y chromosomes: a study of Y chromosome micro-satellites.
Molecular Biology and Evolution 16, 1791–1798.

[Thurstone(1927)]thurs
Thurstone, L. L. (1927).
A law of comparative judgment.
Psychological review 34, 273–286.

[Wallenius(1963)]walle
Wallenius, K. T. (1963).
Biased Sampling: The Non-Central Hypergeometric Probability
  Distribution - Department of Statistics - Stanford University.
Ph.D. thesis, Department of Statistics - Stanford University.





















§ APPENDIX













