Bayesian nonparametric Plackett-Luce models for the analysis of preferences for college degree programmes

Introduction

In this paper we consider partial ranking data consisting of ordered lists of the top-$m$ items among a set of objects. Data in the form of partial rankings arise in many contexts. For example, in this paper we shall consider data pertaining to the top ten preferences of Irish secondary school graduates who are applying to undergraduate degree programmes offered in Irish third level institutions. The third level institutions consist of universities, institutes of technologies and private colleges. This application is described in detail in Section 2.

The Plackett–Luce model [Luce (1959); Plackett (1975)] is a popular model for modeling such partial rankings of a finite collection of $M$ items. It has found many applications, including choice modeling [Luce (1977); Chapman and Staelin (1982)], sport ranking [Hunter (2004)] and voting [Gormley and Murphy (2008)]. Diaconis (1988), Chapter 9, provides detailed discussions on the statistical foundations of this model.

In the Plackett–Luce model, each item $k\in[M]=\{1,\ldots,M\}$ is assigned a positive rating parameter $w_{k}$, which represents the desirability or rating of a product in the case of choice modeling, or the skill of a player in sport rankings. The Plackett–Luce model assumes the following generative story for a top-$m$ list $\rho=(\rho_{1},\ldots,\rho_{m})$ of items $\rho_{i}\in[M]$: at each stage $i=1,\ldots,m$, an item is chosen to be the $i$th item in the list from among the items that have not yet been chosen, with the probability that $\rho_{i}$ is selected being proportional to its desirability $w_{\rho_{i}}$. The overall probability of a given partial ranking $\rho$ is then

with the denominator in (1) being the sum over all items not yet selected at stage $i$.

In many situations the collection of available items can be very large and/or potentially unknown. In this case a nonparametric approach can be sensible, where the pool of items is assumed to be infinite and the model allows for the possibility of items not observed in previous top-$m$ lists to appear in future ones. A naïve approach, building upon recent work on Bayesian inference for the (finite) Plackett–Luce model and its extensions [Gormley and Murphy (2009); Guiver and Snelson (2009); Caron and Doucet (2012)], is to first derive a Markov chain Monte Carlo sampler for the finite model, then to “take the infinite limit” of the sampler, where the number of available items becomes infinite, but such that all unobserved items are grouped together for computational tractability.

Such an approach, outlined in Section 3, is reminiscent of a number of previous approaches deriving the (Gibbs sampler for the) Dirichlet process mixture model as the infinite limit of (a Gibbs sampler for) finite mixture models [Neal (1992); Rasmussen (2000); Ishwaran and Zarepour (2002)]. Although intuitively appealing, this is not a satisfying approach since it is not clear what the underlying nonparametric model actually is, as it is actually the algorithm whose infinite limit was taken. It also does not directly lead to more general and flexible nonparametric models with no obvious finite counterpart, nor does it lead to alternative perspectives and characterisations of the same model, or resultant alternative inference algorithms. Orbanz (2009) further investigates the approach of constructing nonparametric Bayesian models from finite-dimensional parametric Bayesian models.

Caron and Teh (2012) recently proposed a Bayesian nonparametricPlackett–Luce model based on a natural representation of items along with their ratings as an atomic measure. Specifically, the model assumes the existence of an infinite pool of items $\{X_{k}\}_{k=1}^{\infty}$, each with its own rating parameter, $\{w_{k}\}_{k=1}^{\infty}$. The atomic measure then consists of an atom located at each $X_{k}$ with a mass of $w_{k}$:

The probability of a top-$m$ list of items, say, $(X_{\rho_{1}},\ldots,X_{\rho_{m}})$, is then a direct extension of the finite case (1):

Using this representation, note that the top item $X_{\rho_{1}}$ in the list is simply a draw from the probability measure obtained by normalising $G$, while subsequent items in the top-$m$ list are draws from probability measures obtained by first removing from $G$ the atoms corresponding to previously picked items and normalising. Described this way, it is clear that the Plackett–Luce model is none other than a partial size-biased permutation of the atoms in $G$ [Patil and Taillie (1977)], and the existing machinery of random measures and exchangeable random partitions [Pitman (2006); Lijoi and Prünster (2010)] can be brought to bear on our problem.

For example, we may use a variety of existing stochastic processes to specify a prior over the atomic measure $G$. Caron and Teh (2012) considered the case, described in Section 4, where $G$ is a gamma process. This is a completely random measure [Kingman (1967); Lijoi and Prünster (2010)] with gamma marginals, such that the corresponding normalised probability measure is a Dirichlet process [Ferguson (1973)]. They showed that with the introduction of a suitable set of auxiliary variables, it is possible to characterise the posterior law of $G$ given observations of top-$m$ lists distributed according to (3). A simple Gibbs sampler can then be derived to simulate from the posterior distribution which corresponds to the infinite limit of the Gibbs sampler for finite models. In the Appendix, we show that the construction can be extended from gamma processes to general completely random measures, and we discuss extensions of the Gibbs sampler to this more general case.

In Section 5 we describe a Dirichlet process mixture model [Ferguson (1973); Lo (1984)] for heterogeneous partial ranking data, where each mixture component is a gamma process nonparametric Plackett–Luce model. As shown in Section 2, such a model is relevant for capturing heterogeneity in preferences for college degree programmes. As we will see, in this model it is important to allow the same atoms to appear across the different random measures of the mixture components, otherwise the model becomes degenerate with all observed items that ever appeared together in some partial ranking being assigned to the same mixture component. To allow for this, we use a tree-structured extension of the time-varying model of Caron and Teh (2012). In Section 6 we apply this mixture model to the Irish college degree programme preferences data, showing that the model is able to recover clusters of students with similar and interpretable preferences.

Finally, we conclude in Section 7 with a discussion of the important contributions of this paper and proposals for future work.

Irish college degree programmes

Applications to college degree programmes in Ireland are handled by a centralised applications system called the College Application Office (CAO) (\surlwww.cao.ie); a degree programme involves studying a specific subject (broad or focussed) in a particular third level institution. The CAO handles applications for 35 different third level institutions including universities, institutes of technologies and private colleges. In the autumn of each year, a list of all degree programmes for the subsequent year is made available to applicants. Quite often new degree programmes are added to the list of potential choices after the initial list has been published, thus meaning that the potential list of degree programme choices is evolving and not always completely known. Applications are completed early in the year in which the students plan to enter their college degree programme. The list of available degree programmes changes from year to year but has been generally growing in size year on year. Many degree programmes have a specific subject area, for example, Mathematics, History or Computer Science, but others are more general, for example, Science, Commerce or Arts. In the year 2000, which we are examining herein, there were 533 degree programmes available to be selected by the applicants. When students apply for degree programmes they rank up to ten degree programmes, in order of preference, from the list of all degree programmes that are being offered. Two examples of such applications for two different applicants are shown in Table 1.

Places in these degree programmes are allocated on the basis of the applicants’ performance in the Irish Leaving Certificate examination. Students typically take between seven and nine subjects in the Leaving Certificate examination. Points between zero and one hundred are awarded for each applicant’s best six subjects in the Leaving Certificate examination and the points are totalled to give an overall points score. The allocation of applicants to most degree programmes is solely on the basis of the applicant’s points score and applicants with a high points score are more likely to get their high preference choices. The minimum points score of all applicants accepted into a degree programme is publicly available and is called the points requirement. It is worth mentioning that even though degree programmes may have required Leaving Certificate subjects and grades as part of the minimum entry requirements, the subjects used in the applicant’s points score calculation can be any six Leaving Certificate subjects.

The college applications system in Ireland is much debated in the educational sector and it receives much attention in the Irish media. The debate has two main parts: one part of the debate is whether the current system of allocating points to students on the basis of a single Leaving Certificate examination is a fair method, especially when the points can be gained from any Leaving Certificate subjects; the other part of the debate explores the choice behaviour of the applicants and whether students are choosing degree programmes in a coherent manner. We focus on the applicant’s choices which are core to the second part of the debate.

Many people feel that students do not necessarily pick degree programmes on the basis of the courses offered but that they choose on other grounds, like the perceived prestige of the degree programme. However, other factors like geographical location of the third level institution may also have an impact on the applicant’s choice behaviour. The two example applications in Table 1 illustrate that a number of factors influence applicants choices. The first applicant has selected degree programmes in medicine and other health sciences, so their choices appear to be largely based on the course material. However, the second application includes a wide variety of different degree programmes; the applicant’s first choice degree programme leads to a career in Primary Teaching, whereas the other degree programmes are in different areas. However, the institutions that have been chosen are geographically close (within 100 km).

In the year 1997, the Department of Education and Science commissioned a review of the Irish college applications system. A report [Hyland (1999)] reviewed the current system and made some recommendations concerning the future of the system. In addition, four research reports were published, one of which [Tuohy (1998)] examined the applicant’s choices. Tuohy (1998) used a number of exploratory data analysis techniques to investigate the degree programmes selected, but without reference to the preference ordering, and he found that subject matter was an important factor in applicant choices. More recently, Gormley and Murphy (2006) used a finite mixture of Plackett–Luce models to find clusters of applications with similar choice profiles. They fitted their model using maximum likelihood and chose the number of mixture components using the Bayesian Information Criterion (BIC). Their results also indicated that subject matter and geographical location were strong determinants of student choices. However, the model fitting paradigm used in their analysis could not find small clusters of applicants because of the manner that BIC penalises each additional mixture component. Further, McNicholas (2007) used association rule mining to further explore college applicant choices, but he restricted his attention to degree programme choice combinations that were selected by at least 0.5% of the applicants; thus, that analysis emphasised only high frequency choice behaviour.

O’Connell, Clancy and McCoy (2006) conducted a survey of new college entrants (as opposed to applicants) in 2004 and found that the choice of college where they commenced their degree programme was influenced primarily by reputation and geographical location of the third level institution, and that the choice of degree programme was influenced by intrinsic interest in the subject matter and, to a lesser extent, future career prospects. Whilst that study only looks at students who entered college and the degree programme that they ultimately studied, it provides a further insight into the factors that influence choice of degree programme.

We investigate the complete degree programme choice data for the year 2000 cohort of applications to the College Application Office; these data correspond to top-10 rankings of college degree programmes for 53,757 applicants. The model proposed herein has a number of appealing properties because it can account for choosing from the large number of degree programmes on offer, it allows for small differences in preference between degree programmes, it facilitates discovering large and small clusters of applicants with similar preferences, and the fitting in the Bayesian paradigm facilitates a deep exploration of the clustering and co-clustering of applicants.

An extension of the Plackett–Luce model to countably infinite choice sets

We start this section with a review of a Bayesian approach to inference in finite Plackett–Luce models [Gormley and Murphy (2009); Guiver and Snelson (2009); Caron and Doucet (2012)] and take the infinite limit to arrive at a nonparametric model. This will give good intuitions for how the model operates, before we rederive the same nonparametric model more formally in the next section using gamma processes.

being exponentially distributed. The M step of the EM algorithm can be easily derived as well. The resulting algorithm was first proposed by Hunter (2004) as an instance of the MM (majorisation–maximisation) algorithm [Lange, Hunter and Yang (2000)] and its reinterpretation as an EM algorithm was recently given by Caron and Doucet (2012).

Taking a further step, we note that the joint probability (3.1) is conjugate to a factorised gamma prior over the parameters, say, $w_{k}\sim\operatorname{Gamma}(\frac{\alpha}{M},\tau)$ with hyperparameters $\alpha,\tau>0$. Now Bayesian inference can be carried out, for example, using a variational Bayesian EM algorithm or a Gibbs sampler. In this paper we shall consider only Gibbs sampling algorithms. By regrouping the terms in the exponential in (3.1), the parameter updates are derived to be [Caron and Doucet (2012)]:

where $n_{k}$ is the number of occurrences of item $k$ among the observed partial rankings and

2 Taking the infinite limit

A Gibbs sampler for a nonparametricPlackett–Luce model can now be easily derived by taking the limit as the number of choice items $M\rightarrow\infty$. If item $k$ has appeared among the observed partial rankings, the limiting conditional distribution (5) is well defined since $n_{k}>0$. For items that did not appear in the observations, (5) becomes degenerate at 0. Instead we can define $w_{*}=\sum_{k\dvtx n_{k}=0}w_{k}$ to be the total desirability among all the infinitely many unobserved items. Making use of the fact that sums of independent gammas with the same scale parameter is a gamma with shape parameter given by the sum of the shape parameters,

This nonparametric model allows us to estimate the probability of seeing new items appearing in future partial rankings in a coherent manner. While intuitive, the derivation is ad hoc, in the sense that it arises as the infinite limit of the Gibbs sampler for finite Plackett–Luce models, and is unsatisfying, as it did not directly capture the structure of the underlying infinite-dimensional object, which we will show in the next section to be a gamma process.

A Bayesian nonparametric Plackett–Luce model based on the gamma process

where $\alpha>0$ is the concentration parameter and $\tau>0$ the inverse scale. We write this as $G\sim\Gamma(\alpha,\tau,H)$. Under this parametrisation, we have that $G(A)\sim\operatorname{Gamma}(\alpha H(A),\tau)$. $\lambda(w)h(x)$ is known as the Lévy intensity of the homogeneous CRM $G$. The jump part $\lambda(w)$ of the Lévy intensity verifies the necessary condition

and plays a significant role in characterising the properties of the gamma process.

We shall interpret each atom $X_{k}$ as a choice item, with its mass $w_{k}>0$ corresponding to the desirability parameter. The Thurstonian view described in the finite model can be easily extended to the nonparametric one, where a partial ranking $(X_{\rho_{1}},\ldots,X_{\rho_{m}})$ can be generated as the first $m$ items to arrive in a race. In particular, for each atom $X_{k}$ let $z_{k}\sim\operatorname{Exp}(w_{k})$ be the time of arrival of $X_{k}$ and $X_{\rho_{i}}$ the $i$th item to arrive. The first $m$ items to arrive $(X_{\rho_{1}},\ldots,X_{\rho_{m}})$ then constitute our partial ranking, with probability as given in (3). This construction is depicted in Figure 1. The top row of Figure 2 visualises some top-5 rankings generated from the model, with $\tau=1$ and different values of $\alpha$. Figure 3 shows the mean number of items appearing in $L$ top-$m$ rankings. For $m=1$, one recovers the well-known result on the number of clusters for a Dirichlet process model.

Again reparametrising using inter-arrival durations, let $Z_{i}=z_{\rho_{i}}-z_{\rho_{i-1}}$ for $i=1,2,\ldots$ (with $z_{\rho_{0}}=0$). The joint probability of an observed partial ranking of length $m$ along with the $m$ associated latent variables can be derived to be

Marginalising out $(Z_{1},\ldots,Z_{m})$ gives the probability of $(X_{\rho_{1}},\ldots,X_{\rho_{m}})$ as in (3). Further, conditional on $\rho=(\rho_{i})_{i=1}^{m}$, it is seen that the inter-arrival durations $Z_{1},\ldots,Z_{m}$ are mutually independent, with

In the next section we shall characterise the posterior distribution over $G$ given observed partial rankings and their associated latent variables. We end this subsection with two observations.

First, note that the jump part $\lambda(w)$ of the Lévy intensity of the gamma process satisfies the following property:

The joint probability of the item lists and auxiliary variables is then [cf. (4)]

Then the joint probability under the nonparametric Plackett–Luce model is

Taking expectation of (4.1) with respect to $G$ gives the following:

The marginal probability of the $L$ partial rankings and latent variables is

where $\psi(z)$ is the Laplace transform of $\lambda(w)$,

and $\kappa(n,z)$ is the $n$th moment of the exponentially tilted intensity $\lambda(w)e^{-zw}$:

The proof, using the Poisson process characterisation of completely random measures and the Palm formula [James, Lijoi and Prünster (2009)], is given in the Appendix.

Another application of the Palm formula [James, Lijoi and Prünster (2009)] now allows us to derive a posterior characterisation of $G$. The posterior CRM can be decomposed as the sum of a CRM with fixed atoms and a CRM whose jump part of the Lévy intensity is updated to $\lambda^{\ast}(w)$ in a conjugate fashion, similar to deriving a conjugate posterior for a parametric distribution.

where $G^{*}$ and $w^{*}_{1},\ldots,w^{*}_{K}$ are mutually independent. The law of $G^{*}$ is still a gamma process,

The denominator is as given in Theorem 1, while the numerator is obtained using the same Palm formula technique as Theorem 1, with the inclusion of the term $e^{-\int f(x)G(dx)}$. Some algebra then shows that the resulting characteristic functional of the posterior $G$ coincides with that of (13). The proof details are given in the Appendix.

2 Gibbs sampling

Note that the latent variables are conditionally independent given the masses and vice versa. Hyperparameters of the gamma process can be simply derived from the joint distribution in Theorem 1. Since the marginal probability of the partial rankings is invariant to rescaling of the masses, it is sufficient to keep $\tau$ fixed at 1. As for $\alpha$, if a $\operatorname{Gamma}(a,b)$ prior is placed on it, its conditional distribution is still gamma:

Note that this update was derived with $w^{*}_{*}$ marginalised out, so after an update to $\alpha$ it is necessary to immediately update $w^{*}_{*}$ via (4.2) before proceeding to update other variables.

In the Appendix C, we show that the construction can be extended from gamma processes to general completely random measures, and we discuss extensions of the Gibbs sampler to this more general case. In particular, we show that a simple Gibbs sampler can still be derived for the generalised gamma class of completely random measures.

Mixtures of nonparametric Plackett–Luce components

In this section we propose a mixture model for heterogeneous ranking data consisting of nonparametric Plackett–Luce components. Using the same data augmentation scheme, we show that an efficient Gibbs sampler can be derived and apply the model to a data set of preferences for Irish degree programmes by high school graduates.

where $\operatorname{GEM}(\gamma)$ denotes the Griffiths–Engen–McCloskey (GEM) distribution [Pitman (2006)] with concentration parameter $\gamma$ (also known as the stick-breaking construction) and $\operatorname{PL}(G)$ denotes the nonparametric Plackett–Luce model parameterised by the atomic measure $G$ described in Section 4. The $j$th cluster in the mixture model is parameterised by an atomic measure $G_{j}$ and has mixing proportion $\pi_{j}$.

To complete the model, we have to specify the prior on the component atomic measures $G_{j}$. An obvious choice would be to use independent draws from a gamma process $\Gamma(\alpha,\tau,H)$ for each $G_{j}$. This unfortunately does not work. The reason is because if $H$ is smooth, then different atomic measures will never share the same atoms. On the other hand, notice that all items appearing in some observed partial ranking have to come from the same Plackett–Luce model, and thus have to appear as atoms in the corresponding atomic measure. Putting these two observations together, the result is that any observed pair of partial rankings that share a common item will have to be assigned to the same component, and the mixture model will degenerate to using a few much larger components only. In consequence, the model will not capture the fine-scale preference structure that may be present in the partial rankings. This is a similar problem that motivated the hierarchical DP [Teh et al. (2006)], and the solution there, as in here, is to allow different atomic measures to share the same set of atoms, but to allow different atom masses.

Our solution, which is different from Teh et al. (2006), is to make use of the Pitt–Walker [Pitt and Walker (2005)] dependence model for gamma processes. Consider a tree-structured model where there is a single root $G_{0}$ and each component atomic measure $G_{j}$ is a leaf which connects directly to $G_{0}$. The Pitt–Walker model allows us to construct the dependence structure between the root $G_{0}$ and the leaves $(G_{j})$ such that each $G_{j}$ marginally follows a gamma process $\Gamma(\alpha,\tau,H)$. At the root, $G_{0}$ is first given a gamma process prior:

Since $G_{0}$ is atomic, we can write it in the form

Now for each $j$, define a random measure $U_{j}$ with conditional law:

where $\phi>0$ is a parameter which, as we shall see, governs the strength of dependence between $G_{0}$ and each $G_{j}$. Note that since $G_{0}$ has finite total mass, $U_{j}$ consists only of a finite number of atoms with positive masses; the other atoms all have masses equal to zero. Using the same Palm formula method as Section 4.1, we can show the following proposition:

Suppose the prior law of $G_{0}$ is $\Gamma(\alpha,\tau,H)$ and $U_{j}$ has conditional law given by (5.1). The posterior law of $G_{0}$ given $U_{j}$ is then

where $G_{0}^{*}$ and $(w^{*}_{0k})_{k=1}^{\infty}$ are all mutually independent. The law of $G_{0}^{*}$ is given by a gamma process while the masses are conditionally gamma,

Note that if $u_{jk}=0$, we define $w^{*}_{0k}$ to be degenerate at 0, thus, the posterior of $G_{0}$ consists of a finite number of atoms in common with $U_{j}$, along with an infinite number of atoms (those in $G^{*}_{0}$) not in common. The total mass of $G_{0}^{*}$ has distribution $\operatorname{Gamma}(\alpha,\tau+\phi)$.

The idea, inspired by Pitt and Walker (2005), is to define the conditional law of $G_{j}$ given $G_{0}$ and $U_{j}$ to be independent of $G_{0}$ and to coincide with the conditional law of $G_{0}$ given $U_{j}$ as in Proposition 3. In other words, define

where $G_{j}^{*}\sim\Gamma(\alpha,\tau+\phi,H)$ and $w_{jk}^{*}\sim\operatorname{Gamma}(u_{jk},\tau+\phi)$ are mutually independent. Note that if $u_{jk}=0$, the conditional distribution of $w_{jk}^{*}$ will be degenerate at 0. Hence, $G_{j}$ has an atom at $X_{k}$ if and only if $U_{j}$ has an atom at $X_{k}$, that is, if $u_{jk}>0$. In addition, it also has an infinite number of atoms (those in $G_{j}^{*}$) which are in neither $U_{j}$ nor $G_{0}$.

Since the conditional laws of $G_{j}$ and $G_{0}$ given $U_{j}$ coincide, and $G_{0}$ has prior $\Gamma(\alpha,\tau,H)$, it can be seen that $G_{j}$ will marginally follow the same law $\Gamma(\alpha,\tau,H)$ as well. More compactly, we can write the dependence model as

As a final observation, the parameter $\phi$ can be interpreted as controlling the strength of dependence between $G_{0}$ and each $G_{j}$. Indeed, it can be shown that

so that larger $\phi$ corresponds to each $G_{j}$ being more similar to $G_{0}$. Larger $\phi$ may also favour a larger number of clusters, as similar partial rankings are more likely to be clustered in different groups.

Our construction to inducing sharing of atoms has a number of qualitative differences from that of the hierarchical DP [Teh et al. (2006)]. First, the marginal law of each $G_{j}$ is known: it is marginally a gamma process. For the hierarchical DP the marginal laws of the individual random measures are not of simple analytical forms. Since normalising a gamma process gives a DP, our construction can be used as an alternative method to induce sharing of atoms across multiple random measures, each of which still has marginal DP law. Second, in our construction only a finite number of atoms will be shared across random measures (though the number shared can be controlled by the dependence parameter $\phi$), while in the hierarchical DP all infinitely many atoms are shared. In Caron and Teh (2012) we used the Pitt–Walker construction for a different purpose: we constructed a dynamical nonparametric Plackett–Luce model, where at each time $t$, $G_{t}$ is a gamma process, with the Pitt–Walker construction used to define a Markov dependence structure for the sequence of random measures $(G_{t})$.

The structure of (5.1), with a DP mixture with each component specified by a random atomic measure, is reminiscent of the nested DP of Rodríguez, Dunson and Gelfand (2008) as well, though our model has an additional hierarchical structure allowing the sharing of atoms among different component measures. In this respect, it also shares similarities with the hierarchical Dirichlet process model of Müller, Quintana and Rosner (2004).

We focused here on a DP mixture for its simplicity, with a single parameter $\gamma$ tuning the clustering structure. The model can be generalised to more flexible random measures, such as Pitman–Yor processes [Pitman (1995)] or normalised random measures [Regazzini, Lijoi and Prünster (2003); Lijoi, Mena and Prünster (2007)].

2 Posterior characterisation and Gibbs sampling

The overall graphical model is described in Figure 4.

where $G_{j}^{*}$ and $w_{j1},\ldots,w_{jK}$ are mutually independent. The law of $G_{j}^{*}$ is a gamma process,

Note that if $n_{jk}+u_{jk}=0$, then $w_{jk}=0$ and $G_{j}$ will not have a fixed atom at $X^{*}_{k}$. To complete the posterior characterisation, note that, conditioned on $G_{0}$ and $G_{j}$, the variables $u_{j1},\ldots,u_{jK}$ and $u_{j*}$ are independent, with $u_{jk}$ dependent only on $w_{0k}$ and $w_{jk}$ and similarly for $u_{j*}$. The conditional probabilities are

where $1_{a,b}=1$ if $a=b$, otherwise, and $\mathcal{I}$ is the modified Bessel function of the first kind. It is therefore possible to sample exactly from the discrete distributions (23) and (24) using standard retrospective sampling for discrete distributions; see, for example, Papaspiliopoulos and Roberts (2008). Alternatively, we describe in the Appendix a Metropolis–Hastings procedure that worked well in the applications.

Armed with the posterior characterisation, a Gibbs sampler can now be derived. Each iteration of the Gibbs sampler proceeds in the following order (details are in the Appendix): {longlist}[(2)]

The individual atom masses $(w_{jk},w_{j*})$ are scaled along with the update to the total masses. Then the Poisson masses $(u_{jk})$, $(u_{j*})$ are updated using (23) and (24).

The concentration parameter $\alpha$ and the masses $w_{0*}$, $(w_{j*})$ and $(u_{j*})$ associated with other unobserved items are updated efficiently using a forward–backward recursion detailed in the Appendix.

The masses $(w_{0k})$ and $w_{0*}$ of the atoms in $G_{0}$ are updated via an extension of Proposition 3. In particular, for each item $k=1,\ldots,K$, the masses are conditionally independent with distributions

while the total mass of the remaining atoms have conditional distribution

Finally, the scale parameter $\gamma$ of the Dirichlet process is updated using West (1992) and the dependence parameter $\phi$ is updated by a Metropolis–Hastings step using (5.2) and (LABEL:eqcondw2) with the latent $(u_{jk})$ and $(u_{j*})$ marginalised out. The resulting algorithm is a valid partially collapsed Gibbs sampler [van Dyk and Park (2008)]. Note, however, that permutations of the above steps could result in an invalid sampler. The computational cost scales as $O(K\times J\times m\times L)$, where $J$ is the average number of clusters. However, it is possible to parallelise over the different items in the algorithm to obtain an algorithm that scales as $O(J\times m\times L)$.

Application: Irish college degree programmes

We now consider the application of the proposed model to study the choices made by the 53,757 degree programme applicants to the College Application Office (CAO) in the year 2000.

The following flat priors are used for the hyperparameters

We run the Gibbs sampler with $N={}$20,000 iterations. In order to obtain a point estimate of the partition from the posterior distribution, we use the approach proposed by Dahl (2006). Let $c^{(i)},i=1,\ldots,N$ be the Monte Carlo samples. The point estimate $\hat{c}$ is obtained by

where the co-clustering matrix $\zeta$ is obtained with

2 Results

An examination of the Plackett–Luce parameter for each cluster reveals that the subject matter of the degree programme is a strong determinant of the clustering of students (Table 2). For example, clusters 6, 18 and 23 are characterised as construction, medicine and law, respectively. Besides the type of degree, geographical location is a strong determinant of degree programme choice. Clusters 12, 21 and 22 are, respectively, concerned with applications to college degree programmes in Cork, Galway and Limerick. There is a lot of heterogeneity in the subject area of the college degree programmes for these clusters, as can be seen, for example, for the Cork cluster 12 in Table 3. A number of clusters are also defined by a combination of both subject area and location, for example, for clusters 7 and 8 in Tables 4 and 5, which correspond to computer science, respectively, outside and inside Dublin.

As mentioned in Section 2, there is a common perception in the Irish society and media that students pick degree programme based on prestige rather than subject area. Another perception is that the points requirement for a degree programme is a measure of prestige; in fact, the points requirement is determined by a number of factors including the number of available places, the number of applicants who list the degree programme in their top-10 preferences and the quality of the applicants who apply for the degree programme. Such a selection-by-prestige phenomenon should be evidenced by a cluster of students picking degree programmes in medicine and law, both of which have very high points requirements, but no such cluster was found. In fact, medicine and law applicants are clustered separately into clusters 18 and 23, respectively. Therefore, the clustering suggests that students are primarily picking degree programmes on the basis of subject area and geographical considerations; this finding is in agreement with the results found in Gormley and Murphy (2006); McNicholas (2007).

It is also of interest to look at the variability of the student choices within each cluster. This can be quantified by the normalised entropy, which takes its values between 0 and 1, and defined for each cluster $j$ by

where $\hat{w}_{jk}$ are the averaged normalised weights of item $k$ in cluster $j$ obtained from the second MCMC run; the normalised entropy values for each cluster are reported in Table 2. A low value indicates low variability in the choices within a cluster, whereas a large value indicates a lot of variability. Interestingly, cluster 15 has very low normalised entropy, where 56% of the students in that cluster are likely to take one of the three most popular degree programmes of that cluster (Drumcondra, Froebel or Marina) as their first choice; these degree programmes are the main primary teacher education degree programmes in Dublin and, thus, many members of this cluster have a strong interest in teacher education as a degree choice. Further, there is much more variability in cluster 7, where students choices are spread across various computing degree programmes, and only 23% of the students are likely to take one of the three most popular degree programmes as their first choice.

The co-clustering matrix reveals some interesting connections between clusters, which have not been explored in previous analyses of the CAO data. For example, the plot reveals that a number of applicants have high probability of belonging to clusters 4 and 19 which are both in the arts. Cluster 4 is characterised by arts degrees which do not require the applicants to select their major in advance, whereas cluster 19 is characterised by arts degrees where the student needs to specify their major in advance. It is worth observing that the clusters are fairly well separated, and very few clusters exhibit the phenomenon of sharing applicants, which is further evidence that the applicants are only selecting degree programmes of a particular type (as described by the cluster names in Table 2).

Marginal Posterior distributions of the hyperparameters $\alpha$, $\gamma$ and $\phi$ are, respectively, in the ranges $$,$$ and $$. Correlation parameter$\phi$is rather high. This is due to the fact that some degree programmes, such as Arts in University College Dublin or Cork, often appear in the top-ten list of applicants, whatever their main subject matter is. Parameter$\gamma$is associated to the number of clusters, which is around 35. Parameter$\alpha$ relates to the variability of the weights within clusters (and thus to the entropy of the clusters).

Discussion

We have proposed a Bayesian nonparametric Plackett–Luce model for ranked data. Our approach is based on the theory of completely random measures, where we showed that the Plackett–Luce generative model corresponds exactly to a size-biased permutation of the atoms in the random measure. We characterised the posterior distribution and derived a simple MCMC sampling algorithm for posterior simulation. Our approach can be seen as a multi-stage generalisation of posterior inference in normalised random measures [Regazzini, Lijoi and Prünster (2003); James, Lijoi and Prünster (2009); Griffin and Walker (2011); Favaro and Teh (2013)].

We also developed a nonparametric mixture model consisting of nonparametric Plackett–Luce components to model heterogeneity in partial ranking data. In order to allow atoms to be shared across components, we made use of the Pitt–Walker construction, which was previously only used to define Markov dynamical models. Applying our model to a data set of preferences for Irish college degree programmes, we find interesting clustering structure supporting the observation that students were choosing programmes mainly based on subject area and geographical considerations.

It is worthwhile comparing our mixture model to another nonparametric mixture model, DPM-GM, where each component is a generalised Mallows model [Busse, Orbanz and Buhmann (2007); Meilă and Bao (2008); Meilă and Chen (2010)]. In the generalised Mallows model the component distributions are characterised by a (discrete) permutation parameter, whereas in the Plackett–Luce model the component distributions are characterised by a continuous rating parameter. Thus, the Plackett–Luce model offers greater modelling flexibility to capture the strength of preferences for each item. On the other hand, the scale parameters in the generalised Mallows model can accommodate varying precision in the ranking. Additionally, inference for the generalised Mallows models can be difficult.

The mixture model established the existence of clusters of applicants with similar degree programme preferences and characterises these clusters and their coherence in terms of choices. The results support the previous hypotheses that subject matter and geographical location are the primary drivers of degree programme choice [Gormley and Murphy (2006); McNicholas (2007)]. These factors are important because they reflect the intrinsic interest in the subject matter of the degree programmes and the economic and practical aspects of choosing a third level institution for study. The geographical location influence is further supported by results on acceptances to degree programmes [O’Connell, Clancy and McCoy (2006)] and studies on how students fund their education which found that 45% of Irish university students live in their family home [Clancy and Kehoe (1999)] and thus attend an institution that is geographically close by.

An interesting extension of the proposed model would be to consider inhomogeneous completely random measures, where the preferences would depend on a set of covariates (e.g., location).

Appendix A Proof of Theorem 1

The marginal probability (12) is obtained by taking the expectation of (4.1) with respect to $G$. Note however that (4.1) is a density, so to be totally precise here we need to work with the probability of infinitesimal neighborhoods around the observations instead, which introduces significant notational complexity. To keep the notation simple, we will work with densities, leaving it to the careful reader to verify that the calculations indeed carry over to the case of probabilities.

We now recall the Palm formula [see e.g., Bertoin (2006), Lemma 2.3].

Palm Formula. Let $N$ be a Poisson process on $S$ with mean measure $\nu$. Let $S_{p}$ denote the set of point measures on $S$, $f\dvtx S\rightarrow{}[0,+\infty{}[$ and $\mathcal{G}\dvtx S\times S_{p}\rightarrow{}[0,+\infty{}[$ be some measurable functional. Then we have the so-called Palm formula

where the expectation is with respect to $N$.

Applying the Palm formula for Poisson processes to pull the $k=1$ term out of the expectation,

Now iteratively pull out terms $k=2,\ldots,K$ using the same idea, and we get:

Appendix B Proof of Theorem 2

The proof is essentially obtained by calculating the numerator and denominator of (14). The denominator is already given in Theorem 1. The numerator is obtained using the same technique with the inclusion of the term $e^{\int f(x)G(dx)}$, which gives

By the Lévy–Khintchine theorem (using the fact that $G$ has a Poisson process representation $N$),

Dividing the numerator (B) by the denominator (A), the characteristic functional of the posterior $G$ is

Appendix C Generalisation to completely random measures

Both Theorems 1 and 2 generalise naturally to homogeneous CRMs. In fact the statements and the proofs in the appendix still hold with the more general Lévy intensity, along with its Laplace transform $\psi(z)$ and moment function $\kappa(n,z)$:

The marginal probability of the $L$ partial rankings and latent variables is

where $\psi(z)$ is the Laplace transform of $\lambda(w)$,

and $\kappa(n,z)$ is the $n$th moment of the exponentially tilted Lévy intensity $\lambda(w)e^{-zw}$:

where $G^{*}$ and $w^{*}_{1},\ldots,w^{*}_{K}$ are mutually independent. The law of $G^{*}$ is a homogeneous CRM with an exponentially tilted Lévy intensity:

Examples of CRMs that have been explored in the literature for Bayesian nonparametric modelling include the stable process [Kingman (1975)], the inverse Gaussian process [Lijoi, Mena and Prünster (2005)], the generalised gamma process [Brix (1999)], and the beta process [Hjort (1990)]. The generalised gamma process forms the largest known simple and tractable family of CRMs, with the gamma, stable and inverse Gaussian processes included as subfamilies. It has a Lévy intensity of the form

where the concentration parameter is $\alpha>0$, the inverse scale is $\tau\geq 0$, and the index is $0\leq\sigma<1$. The gamma process is recovered when $\sigma=0$, the stable when $\tau=0$, and the inverse Gaussian when $\sigma=1/2$. The Laplace transform and the moment function of the generalised gamma process are

where $f(w)$ is the density (assumed to exist) of the total mass $w^{*}_{*}$ under a CRM with the prior Lévy intensity $\lambda(w)$. Note that integrating out the parameters $(w^{*}_{k}),w^{*}_{*}$ from (C) gives the marginal probability in Theorem 1′. From the joint probability (C), the Gibbs sampler can now be derived:

Appendix D Gibbs sampler for the mixture of nonparametric Plackett–Luce components

Let $J$ be the number of different values taken by $c$ (number of clusters). Please note that the number of clusters is not set in advance and its value may change at each iteration. The Gibbs sampler proceeds with each of the following updates in turn:

For $j=1,\ldots,J$, update $(u_{j},u_{j*})$ given ($w_{0},w_{0*},w_{j},w_{j*},\phi,\alpha,c)$.

Update $w_{0*}$ given $(Z,\phi,c,\alpha)$.

For $j=1,\ldots,J$, update $u_{j*}$ given $(Z,\phi,c,\alpha,w_{0*})$.

For $j=1,\ldots,J$, update $w_{j*}$ given $(Z,\alpha,u_{j*},\phi,c)$.

Update $(w_{0k}),w_{0*}$ given $(U_{1\dvtx J},\alpha)$.

For $j=1,\ldots,J$, update $(w_{j},w_{j*})$ given $(Z,\alpha,u_{j},u_{j*},\phi,c)$.

Update $\phi$ given $w_{0},w_{0*},w_{1\dvtx J},w_{1\dvtx J*},\alpha,\phi$.

For $j=1,\ldots,J$, update $(u_{j},u_{j*})$ given ($w_{0},w_{0*},w_{j},w_{j*},\phi,\alpha,c)$

Consider first the sampling of $u_{j}$. We have, for $j=1,\ldots,J$ and $k=1,\ldots,K$

Hence we can have the following MH update. If $w_{jk}>0$, then we necessarily have $u_{jk}>0$. We sample $u_{jk}^{\ast}\sim$zPoisson$(\phi w_{0k})$ where zPoisson$(\phi w_{0k})$ denotes the zero-truncated Poisson distribution and accept $u_{jk}^{\ast}$ with probability

If $w_{jk}=0$, we only have two possible moves: $u_{jk}=0$ or $u_{jk}=1$, given by the following probabilities

Note that the above Markov chain is not irreducible, as the probability is zero to go from a state $(u_{jk}>0,w_{jk}>0)$ to a state $(u_{jk}=0,w_{jk}=0)$, even though the posterior probability of this event is nonzero in the case item $k$ does not appear in cluster $j$. We can add such moves by jointly sampling $(u_{jk},w_{jk})$. For each $k$ that does not appear in cluster $j$, sample $u_{jk}^{\ast}\sim\operatorname{Poisson}(\phi w_{0k})$ then set $w_{jk}^{\ast}=0$ if $u_{jk}^{\ast}=0$ otherwise sample $w_{jk}^{\ast}\sim\operatorname{Gamma}(u_{jk},\tau+\phi)$. Accept $(u_{jk}^{\ast},w_{jk}^{\ast})$ with probability

We now consider sampling of $u_{j*}$, $j=1,\ldots,J$. We can use a MH step. Sample $w_{j*}^{\ast}\sim\operatorname{Poisson}(\phi w_{0*})$ and accept with probability

We can sample from the full conditional which is given by

Update $w_{0*}$ given $(Z,\phi,c,\alpha)$

We can sample from the full conditional which is given by

where $x_{0}$ is defined above. {longlist}

For $j=1,\ldots,J$, update $u_{j*}$ given $(Z,\phi,c,\alpha,w_{0*})$

We can sample from the full conditional which is given, for $j=1,\ldots,J$ by

where $\widetilde{Z}_{j}$ is defined above. {longlist}

For $j=1,\ldots,J$, update $w_{j*}$ given $(Z,\alpha,u_{j*},\phi,c)$

We can sample from the full conditional which is given, for $j=1,\ldots,J$ by

where $\widetilde{Z}_{j}$ is defined above. {longlist}

Update $(w_{0k}),w_{0*}$ given $(U_{1\dvtx J},\alpha)$

For $j=1,\ldots,J$, update $(w_{jk}),w_{j*}$ given $(Z,\alpha,u_{j},u_{j*},\phi,c)$

if $u_{jk}+n_{jk}>0$, otherwise, set $w_{jk}=0$.

The allocation variables $(c_{1},\ldots,c_{L})$ are updated using the slice sampling technique described in [Walker (2007); Kalli, Griffin and Walker (2011); Fall and Barat (2012)]. It builds on the Introduction of additional latent slice variables, and does not require to set any truncation. For completeness, we briefly recall here the details of the sampler. From equation (5.1), we have

where the $\pi_{k}$ admit the following stick-breaking representation

Sample $(\pi_{1},\ldots,\pi_{J},\pi_{*})\sim\operatorname{Dirichlet}(\mu_{1},\ldots,\mu_{J},\gamma)$.

Set $k=J$. While $\sum_{j=1}^{k}\pi_{k}<(1-\min(\omega_{1},\ldots,\omega_{L}))$.

Sample $v_{k}\sim\operatorname{Beta}(1,\gamma)$.

Set $\pi_{k}=\pi_{*}v_{k}\prod_{j=J+1}^{k-1}(1-v_{j})$.

Sample $G_{k}$ given $G_{0}$ using equation (5.1).

The scale parameter $\gamma$ of the Dirichlet process is updated using the data augmentation technique of West (1992). {longlist}

Update $\phi$ given $w_{0},w_{0*},w_{1\dvtx J},w_{1\dvtx J*},\alpha,\phi$

We sample $\phi$ using a MH step. Propose $\phi^{\ast}=\phi\exp(\sigma\varepsilon)$ where $\sigma>0$ and $\varepsilon\sim\mathcal{N}(0,1)$. And accept it with probability

Acknowledgments

The authors thank Igor Prünster for very helpful feedback on an earlier version of this work. François Caron acknowledges the support of the European Commission under the Marie Curie Intra-European Fellowship Programme.The contents reflect only the authors views and not the views of the European Commission.

References