1.1 Comparison with other packages

Even though the Plackett-Luce model is a well-established method for analysing rankings, the software available to fit the model is limited. By considering each choice in the ranking as a multinomial observation, with one item observed out of a possible set, the “Poisson trick” (see, for example, Baker 1994) can be applied to express the model as a log-linear model, where the response is the count (one or zero) of each possible outcome within each choice. In principle, the model can then be fitted using standard software for generalized linear models. However there are a number of difficulties with this. Firstly, dummy variables must be set up to represent the presence or absence of each item in each choice and a factor created to identify each choice, which is a non-standard task. Secondly the factor identifying each choice will have many levels: greater than the number of rankings, for rankings of more than two objects. Thus there are many parameters to estimate and a standard function such as glm will be slow to fit the model, or may even fail as the corresponding model matrix will be too large to fit in memory. This issue can be circumvented by using the gnm function from gnm, which provides an eliminate argument to efficiently estimate the effects of such a factor. Even then, the model-fitting may be relatively slow, given the expansion in the number of observations when rankings are converted to counts. For example, the ranking {item 3 \(\prec\) item 1 \(\prec\) item 2} expands to two choices with five counts all together:

##      choice item 1 item 2 item 3 count
## [1,]      1      1      0      0     0
## [2,]      1      0      1      0     0
## [3,]      1      0      0      1     1
## [4,]      2      1      0      0     1
## [5,]      2      0      1      0     0

It is possible to aggregate observations of the same choice from the same set of alternatives, but the number of combinations increases quickly with the number of items.

Given the issues with applying general methods, custom algorithms and software have been developed. One approach is using Hunter’s (2004) minorization-maximization (MM) algorithm to maximize the likelihood, which is equivalent to an iterative scaling algorithm; this algorithm is used by the StatRank package. Alternatively the likelihood of the observed data under the PlackettLuce model can be maximised directly using a generic optimisation method such the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm, as is done by the pmr and hyper2 packages. Finally, Bayesian methods can be used to either maximize the posterior distribution via an Expectation Maximization (EM) algorithm or to simulate the posterior distribution using Markov-chain Monte-Carlo (MCMC) techniques, both of which are provided by PLMIX. PlackettLuce offers both iterative scaling and generic optimization using either BFGS or a limited memory variant (L-BFGS) via the lbfgs package.

Even some of these specialized implementations can scale poorly with the number of items and/or the number of rankings as shown by the example timings in Table 2. Specifically pmr::pl becomes impractical to use with a moderate number of items (~10), while the functions from hyper2 and StatRank take much longer to run with a large number (1000s) of unique rankings. PlackettLuce copes well with these moderately-sized data sets, though is not quite as fast as PLMIX when both the number of items and the number of unique rankings is large.

## Warning in kableExtra::kable_styling(kab): Please specify format in
## kable. kableExtra can customize either HTML or LaTeX outputs. See https://
## haozhu233.github.io/kableExtra/ for details.

(ref:timings) Timings in seconds for fitting the Plackett-Luce model to data sets summarised in Table 1 using different packages. See Appendix 5.1 for details and code.

## Warning in kableExtra::kable_styling(kab): Please specify format in
## kable. kableExtra can customize either HTML or LaTeX outputs. See https://
## haozhu233.github.io/kableExtra/ for details.

## Warning in kableExtra::add_header_above(kab, c(` ` = 1, `Time elapsed (s)` =
## 5)): Please specify format in kable. kableExtra can customize either HTML or
## LaTeX outputs. See https://haozhu233.github.io/kableExtra/ for details.

Note: ^a Function fails to complete.

As the number of items increases, it is typically more common to observe partial rankings than complete rankings. Partial rankings can be of two types: sub-rankings, where only a subset of items are ranked each time, and incomplete rankings, where the top \(n\) items are selected and the remaining items are unranked, but implicitly ranked lower than the top \(n\). Sub-rankings are accommodated by the standard Plackett-Luce model, while incomplete rankings require an extension of the Plackett-Luce model. PlackettLuce handles sub-rankings only, while PLMIX handles incomplete rankings only and hyper2 can handle both types. StatRank seems to support partial rankings, but the extent of this support is not clear. The timings in Table 3 for fitting the Plackett-Luce model on the NASCAR data from Hunter (2004) illustrate that PlackettLuce is more efficient than hyper2 for modelling sub-rankings of a relatively large number of items.

## Warning in kableExtra::kable_styling(kab): Please specify format in
## kable. kableExtra can customize either HTML or LaTeX outputs. See https://
## haozhu233.github.io/kableExtra/ for details.

## Warning in kableExtra::add_header_above(kab, c(`Features of NASCAR data` =
## 3, : Please specify format in kable. kableExtra can customize either HTML or
## LaTeX outputs. See https://haozhu233.github.io/kableExtra/ for details.

PlackettLuce is the only package out of those based on maximum likelihood estimation with the functionality to compute standard errors for the item parameters and thereby the facility to conduct inference about these parameters. PLMIX allows for inference based on the posterior distribution. In some cases, when the network of wins and losses is disconnected or weakly connected, the maximum likelihood estimate does not exist, or has infinite standard error; such issues are handled in PlackettLuce by utilising pseudo-rankings. This is similar to incorporating prior information as in the Bayesian approach.

PlackettLuce is also the only package that can accommodate tied rankings, through a novel extension of the Plackett-Luce model. On the other hand hyper2 is currently the only package that can handle rankings of combinations of items, for example team rankings in sports. PLMIX offers the facility to model heterogeneous populations of subjects that have different sets of worth parameters via mixture models. This is similar in spirit to the model-based partitioning offered by PlackettLuce, except here the sub-populations are defined by binary splits on subject attributes. A summary of the features of the various packages for Plackett-Luce models is given in Table 4.

## Warning in kableExtra::kable_styling(kab): Please specify format in
## kable. kableExtra can customize either HTML or LaTeX outputs. See https://
## haozhu233.github.io/kableExtra/ for details.

2.1 Extended Plackett-Luce model

The PlackettLuce package permits rankings of the form

\[R = \{C_1, C_2, \ldots, C_J\}\]

where the items in set \(C_1\) are ranked higher than (better than) the items in \(C_2\), and so on. If there are multiple objects in set \(C_j\) these items are tied in the ranking. For a set \(S\), let

\[f(S) = \delta_{|S|} \left(\prod_{i \in S} \alpha_i \right)^\frac{1}{|S|}\]

where \(|S|\) is the cardinality of the set \(S\), \(\delta_n\) is a parameter representing the prevalence of ties of order \(n\), and \(\alpha_i\) is a parameter representing the worth of item \(i\). Then under an extension of the Plackett-Luce model allowing ties up to order \(D\), the probability of the ranking \(R\) is given by

where \(D_j\) is the cardinality of \(C_j\), \(A_j\) is the set of alternatives from which \(C_j\) is chosen, and \(A_j \choose k\) is all the possible choices of \(k\) items from \(A_j\). The value of \(D\) can be set to the maximum number of tied items observed in the data, so that \(\delta_n = 0\) for \(n > D\).

When the worth parameters are constrained to sum to one, they represent the probability that the corresponding item comes first in a ranking of all items, given that first place is not tied.

The 2-way tie prevalence parameter \(\delta_2\) is interpretable via the probability that two given items of equal worth tie for first place, given that the first place is not a 3-way or higher tie. Specifically, that probability is \(\delta_2/(2 + \delta_2)\).

The 3-way and higher tie-prevalence parameters are interpretable similarly, in terms of tie probabilities among equal-worth items.

library(PlackettLuce)
head(pudding)

##   i j r_ij w_ij w_ji t_ij
## 1 1 2   57   19   22   16
## 2 1 3   47   16   19   12
## 3 2 3   48   19   19   10
## 4 1 4   54   18   23   13
## 5 2 4   51   23   19    9
## 6 3 4   54   19   20   15

nr <- 3*nrow(pudding)
R <- matrix(0, nrow = nr, ncol = 6,
            dimnames = list(NULL, seq_len(6)))
i <- rep(pudding$i, 3)
j <- rep(pudding$j, 3)
R[cbind(seq_len(nr), i)] <- rep(c(1, 2, 1), each = nrow(pudding))
R[cbind(seq_len(nr), j)] <- rep(c(2, 1, 1), each = nrow(pudding))
head(R, 3)

##      1 2 3 4 5 6
## [1,] 1 2 0 0 0 0
## [2,] 1 0 2 0 0 0
## [3,] 0 1 2 0 0 0

tail(R, 3)

##       1 2 3 4 5 6
## [43,] 0 0 1 0 0 1
## [44,] 0 0 0 1 0 1
## [45,] 0 0 0 0 1 1

R <- as.rankings(R)
head(R)

## [1] "1 > 2" "1 > 3" "2 > 3" "1 > 4" "2 > 4" "3 > 4"

tail(R)

## [1] "4 = 5" "1 = 6" "2 = 6" "3 = 6" "4 = 6" "5 = 6"

w <- unlist(pudding[c("w_ij", "w_ji", "t_ij")])

mod <- PlackettLuce(R, weights = w, npseudo = 0, maxit = 7)

## Warning in PlackettLuce(R, weights = w, npseudo = 0, maxit = 7): Iterations
## have not converged.

coef(mod, log = FALSE)

##         1         2         3         4         5         6      tie2 
## 0.1388005 0.1729985 0.1617420 0.1653930 0.1586805 0.2023855 0.7468147

2.2 Inference

A standard way to report model parameter estimates is to report them along with their corresponding standard error. This is an indication of the estimate’s precision, however implicitly this invites comparison with zero. Such comparison is made explicit in many summary methods for models in R, with the addition of partial t or Z tests testing the null hypothesis that the parameter is equal to zero, given the other parameters in the model. However this hypothesis is generally not of interest for the worth parameters in a Plackett-Luce model: we expect most items to have some worth, the question is whether the items differ in their worth. In addition, a Z test based on asymptotic normality of the maximum likelihood estimate will not be appropriate for worth parameters near zero or one, since it does not take account of the fact that the parameters cannot be outside of these limits.

On the log scale however, there are no bounds on the parameters and we can set a reference level to provide meaningful comparisons. By default, the summary method for "PlackettLuce" objects sets the first item (the first element of colnames(R)) as the reference:

summary(mod)

## Call: PlackettLuce(rankings = R, npseudo = 0, weights = w, maxit = 7)
## 
## Coefficients:
##      Estimate Std. Error z value Pr(>|z|)    
## 1      0.0000         NA      NA       NA    
## 2      0.2202     0.1872   1.176 0.239429    
## 3      0.1530     0.1935   0.790 0.429271    
## 4      0.1753     0.1882   0.931 0.351683    
## 5      0.1339     0.1927   0.695 0.487298    
## 6      0.3771     0.1924   1.960 0.049983 *  
## tie2  -0.2919     0.0825  -3.539 0.000402 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual deviance:  1619.4 on 1484 degrees of freedom
## AIC:  1631.4 
## Number of iterations: 7

None of the Z tests for the item parameters provides significant evidence against the null hypothesis of no difference from the worth of item 1, which is consistent with the test for equal preferences presented in Davidson (1970). The tie parameter is also shown on the log scale here, but it is an integral part of the model rather than a parameter of interest for inference, and its scale is not as relevant as that of the worth parameters.

The reference level for the item parameters can be changed via the ref argument, for example setting to NULL sets the mean worth as the reference:

summary(mod, ref = NULL)

## Call: PlackettLuce(rankings = R, npseudo = 0, weights = w, maxit = 7)
## 
## Coefficients:
##       Estimate Std. Error z value Pr(>|z|)    
## 1    -0.176581   0.121949  -1.448 0.147619    
## 2     0.043664   0.121818   0.358 0.720019    
## 3    -0.023617   0.126823  -0.186 0.852274    
## 4    -0.001295   0.122003  -0.011 0.991533    
## 5    -0.042726   0.127054  -0.336 0.736657    
## 6     0.200555   0.126594   1.584 0.113140    
## tie2 -0.291938   0.082499  -3.539 0.000402 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual deviance:  1619.4 on 1484 degrees of freedom
## AIC:  1631.4 
## Number of iterations: 7

As can be seen from the output above, the standard error of the item parameters changes with the reference level. Therefore in cases where there is not a natural reference (like for example in comparisons of own brand versus competitor’s brands), inference can depend on an arbitrary choice. This problem can be handled through the use of quasi standard errors that remain constant for a given item regardless of the reference. In addition quasi standard errors are defined for the reference item, so even in cases where there is a natural reference, the uncertainty around the worth of that item can still be represented.

Quasi standard errors for the item parameters are implemented via a method for the qvcalc function from the qvcalc package:

qv <- qvcalc(mod)
summary(qv)

## Model call:  PlackettLuce(rankings = R, npseudo = 0, weights = w, maxit = 7) 
##        estimate        SE   quasiSE   quasiVar
##     1 0.0000000 0.0000000 0.1328950 0.01766108
##     2 0.2202447 0.1872168 0.1327373 0.01761919
##     3 0.1529644 0.1935181 0.1395740 0.01948091
##     4 0.1752864 0.1882110 0.1330240 0.01769538
##     5 0.1338550 0.1927043 0.1399253 0.01957908
##     6 0.3771362 0.1924059 0.1392047 0.01937796
## Worst relative errors in SEs of simple contrasts (%):  -0.8 0.8 
## Worst relative errors over *all* contrasts (%):  -1.7 1.7

Again by default, the first item is taken as the reference, but this may be changed via a ref argument. The plot method for the returned object visualizes the item parameters (log-worth parameters) along with comparison intervals - item parameters for which the comparison intervals do not cross are significantly different:

plot(qv, xlab = "Brand of pudding", ylab = "Worth (log)", main = NULL)

Figure 1: Worth of brands of chocolate pudding
Intervals based on quasi-standard errors.

The quasi-variances allow comparisons that are approximately correct, for every possible contrast among the parameters. The routine error report in the last two lines printed above by summary(qv) tells us that, in this example, the approximation error has been very small: the approximation error for the standard error of any simple contrast among the parameters is less than 0.8%.

2.3 Disconnected networks

The wins and losses between items can be represented as a directed network. For example, consider the following set of paired comparisons:

R <- matrix(c(1, 2, 0, 0,
              2, 0, 1, 0,
              1, 0, 0, 2,
              2, 1, 0, 0,
              0, 1, 2, 0), byrow = TRUE, ncol = 4,
            dimnames = list(NULL, LETTERS[1:4]))
R <- as.rankings(R)

The adjacency function from PlackettLuce can be used to convert these to an adjacency matrix where element \((i, j)\) is the number of times item \(i\) is ranked higher than item \(j\):

A <- adjacency(R)
A

##   A B C D
## A 0 1 0 1
## B 1 0 1 0
## C 1 0 0 0
## D 0 0 0 0
## attr(,"class")
## [1] "adjacency" "matrix"

Using functions from igraph we can visualise the corresponding network:

library(igraph)
net <- graph_from_adjacency_matrix(A)
plot(net, edge.arrow.size = 0.5, vertex.size = 30)

Figure 2: Network representation of toy rankings

A sufficient condition for the worth parameters (on the log scale) to have finite maximum likelihood estimates (MLEs) and standard errors is that the network is strongly connected, i.e. there is a path of wins and a path of losses between each pair of items. In the example above, A, B and C are strongly connected. For example, C directly loses against B and although C never directly beats B, it does beat A and A in turn beats B, so C indirectly beats B. Similar paths of wins and losses can be found for all pairs of A, B and C. On the other hand D is only observed to lose, therefore the MLE of the log-worth would be \(-\infty\), with infinite standard error.

If one item always wins, the MLE of the log-worth would be \(+\infty\) with infinite standard error. Or if there are clusters of items that are strongly connected with each other, but disconnected or connected only by wins or only by loses (weakly connected) to other clusters, then the maximum likelihood estimates are undefined, because there is no information on the relative worth of the clusters or one cluster is infinitely worse than the other.

The connectivity of the network can be checked with the connectivity function from PlackettLuce

connectivity(A)

## $membership
## A B C D 
## 1 1 1 2 
## 
## $csize
## [1] 3 1
## 
## $no
## [1] 2

If the network is not strongly connected, information on the clusters within the network is returned. In this case a model could be estimated by excluding item D:

R2 <- R[, -4]
R2

## [1] "A > B" "C > A" "B > A" "B > C"

mod <- PlackettLuce(R2, npseudo = 0)
summary(mod)

## Call: PlackettLuce(rankings = R2, npseudo = 0)
## 
## Coefficients:
##   Estimate Std. Error z value Pr(>|z|)
## A   0.0000         NA      NA       NA
## B   0.8392     1.3596   0.617    0.537
## C   0.4196     1.5973   0.263    0.793
## 
## Residual deviance:  5.1356 on 2 degrees of freedom
## AIC:  9.1356 
## Number of iterations: 3

Note that since R is a rankings object, the rankings are automatically updated when items are dropped, so in this case the paired comparison with item D is dropped.

By default however PlackettLuce provides a way to handle disconnected/weakly connected networks, through the addition of pseudo-rankings. This works by adding a win and a loss between each item and a hypothetical or ghost item with fixed worth. This makes the network strongly connected so all the worth parameters are estimable. It also has an interpretation as a Bayesian prior, in particular an exchangeable prior in which all items have equal worth.

The npseudo argument defines the number of wins and loses with the ghost item that are added for each real item. Setting npseudo = 0 means that no pseudo-rankings are added, so PlackettLuce will return the standard MLE if the network is strongly connected and throw an error otherwise. The larger npseudo is, the stronger the influence of the prior, by default npseudo is set to 0.5, so each pseudo-ranking is weighted by 0.5. This is enough to connect the network, but is a weak prior. In this toy example, the item parameters change quite considerably:

mod2 <- PlackettLuce(R)
coef(mod2)

##          A          B          C          D 
##  0.0000000  0.5184185  0.1354707 -1.1537565

This is because there are only 5 rankings, so there is not much information in the data. In more realistic examples, the default prior will have a weak shrinkage effect, shrinking the items’ worth parameters towards \(1/N\), where \(N\) is the number of items.

For a practical example, we consider the NASCAR data from Hunter (2004). This collects the results of the 36 races in the 2002 NASCAR season in the United States. Each race involves 43 drivers out of a total of 87 drivers. The d_nascar data provided by PlackettLuce records the results as an ordering of the drivers in each race:

data(nascar)
nascar[1:2, 1:45]

##      rank1 rank2 rank3 rank4 rank5 rank6 rank7 rank8 rank9 rank10 rank11
## [1,]    83    18    20    48    53    51    67    72    32     42      2
## [2,]    52    72     4    82    60    31    32    66     3     44      2
##      rank12 rank13 rank14 rank15 rank16 rank17 rank18 rank19 rank20 rank21
## [1,]     31     62     13     37      6     60     66     33     77     56
## [2,]     48     83     67     41     77     33     61     45     38     51
##      rank22 rank23 rank24 rank25 rank26 rank27 rank28 rank29 rank30 rank31
## [1,]     63     55     70     14     43     71     35     12     44     79
## [2,]     14     42     62     35     12     25     37     34      6     18
##      rank32 rank33 rank34 rank35 rank36 rank37 rank38 rank39 rank40 rank41
## [1,]      3     52      4      9     45     41     61     34     39     49
## [2,]     79     39     59     43     55     49     56      9     53      7
##      rank42 rank43 rank44 rank45
## [1,]     15     82      0      0
## [2,]     13     71      0      0

For example, in the first race, driver 83 came first, followed by driver 18 and so on. Ranks 43 to 87 are zero for all races. We can convert these orderings to rankings using as.rankings with input = "ordering":

R <- as.rankings(nascar, input = "ordering")
R[1:2,]

## [1] "83 > 18 > 20 > 48 > 53 > 51 > 67 > 7 ..."
## [2] "52 > 72 > 4 > 82 > 60 > 31 > 32 > 66 ..."

The names corresponding to the driver IDs are available as an attribute of nascar; since the columns of R now correspond to the drivers, we can use these as the column names:

colnames(R) <- attr(nascar, "drivers")
R[1:3, 1:3, as.rankings = FALSE]

##      Austin Cameron Bill Elliott Bobby Hamilton
## [1,]              0           11             32
## [2,]              0           11              9
## [3,]              0            8             43

R[1:3]

## [1] "Ward Burton > Elliott Sadler > Geoff ..."
## [2] "Matt Kenseth > Sterling Marlin > Bob ..."
## [3] "Sterling Marlin > Jeremy Mayfield >  ..."

Maximum likelihood estimation cannot be used in this example, because four drivers placed last in each race they entered. So Hunter (2004) dropped these four drivers to fit the Plackett-Luce model, which we can reproduce as follows:

keep <- seq_len(83)
R2 <- R[, keep]
mod <- PlackettLuce(R2, npseudo = 0)

In order to demonstrate the correspondence with the results from Hunter (2004), we order the item parameters by the driver’s average rank:

avRank <- apply(R, 2, function(x) mean(x[x > 0]))
coefs <- round(coef(mod)[order(avRank[keep])], 2)
head(coefs, 3)

##     PJ Jones Scott Pruett  Mark Martin 
##         4.15         3.62         2.08

tail(coefs, 3)

##  Dave Marcis Dick Trickle    Joe Varde 
##         0.03        -0.31        -0.15

Now we fit the Plackett-Luce model to the full data, using the default pseudo-rankings method.

mod2 <- PlackettLuce(R)

For items that were in the previous model, we see that the log-worth parameters generally shrink towards zero:

coefs2 <- round(coef(mod2), 2)
coefs2[names(coefs)[1:3]]

##     PJ Jones Scott Pruett  Mark Martin 
##         3.20         2.77         1.91

coefs2[names(coefs)[81:83]]

##  Dave Marcis Dick Trickle    Joe Varde 
##         0.02        -0.38        -0.12

The new items have relative large negative log worth

coefs2[84:87]

## Andy Hillenburg  Gary Bradberry  Jason Hedlesky   Randy Renfrow 
##           -2.17           -1.74           -1.59           -1.77

Nonetheless, the estimates are finite and have finite standard errors:

coef(summary(mod2))[84:87,]

##                  Estimate Std. Error    z value  Pr(>|z|)
## Andy Hillenburg -2.171065   1.812994 -1.1975028 0.2311106
## Gary Bradberry  -1.744754   1.855365 -0.9403828 0.3470212
## Jason Hedlesky  -1.590764   1.881708 -0.8453828 0.3978972
## Randy Renfrow   -1.768629   1.904871 -0.9284767 0.3531604

Note that the reference here is simply the driver that comes first alphabetically: A. Cameron. We can plot the quasi-variances for a better comparison:

qv <- qvcalc(mod2)
qv$qvframe <- qv$qvframe[order(coef(mod2)),]
plot(qv, xlab = NULL, ylab = "Ability (log)", main = NULL, xaxt="n")
axis(1, at = seq_len(87), labels = rownames(qv$qvframe), las = 2, cex.axis = 0.7)

Figure 3: Ability of drivers based on NASCAR 2002 season
Intervals based on quasi-standard errors.

As in the previous example, we can use summary(qv) to see a report on the accuracy of the quasi-variance approximation. In this example the error of approximation, across the standard errors of all of the 3741 possible simple contrasts (contrasts between pairs of the 87 driver-specific parameters), ranges between -0.7% and +6.7% — which is still remarkably accurate, and which means that the plot of comparison intervals is a good visual guide to the uncertainty about drivers’ relative abilities. The results of using summary(qv) are not shown here, as they would occupy too much space.

Although the prior pseudo-rankings are only necessary to add when the network is incomplete, the default behaviour is always to use them (with a weight of 0.5) because the small shrinkage effect that the pseudo-data delivers has also a beneficial impact in terms of reducing the bias and the variance of the estimators.

5.1 Timings

Data for the package comparison in Table 2 was downloaded from PrefLib (Mattei and Walsh 2013) using the read.soc function provided in PlackettLuce to read in files with the “Strict Orders - Complete List” format.

library(PlackettLuce)
# read in example data sets
preflib <- "http://www.preflib.org/data/election/"
netflix <- read.soc(file.path(preflib, "netflix/ED-00004-00000101.soc"))
tshirt <- read.soc(file.path(preflib, "shirt/ED-00012-00000001.soc"))
sushi <- read.soc(file.path(preflib, "sushi/ED-00014-00000001.soc"))

A wrapper was defined for each function in the comparison, to prepare the rankings and run each function with reasonable defaults. The Plackett-Luce model was fitted to aggregated rankings where possible (for PlackettLuce, hyper2, and pmr). Arguments were set to obtain the maximum likelihood estimate, with the default convergence criteria. The default iterative scaling algorithm was used for PlackettLuce.

pl <- function(dat, ...){
    # convert ordered items to ranking
    R <- as.rankings(dat[,-1], "ordering")
    # fit without adding pseudo-rankings, weight rankings by count
    PlackettLuce(R, npseudo = 0, weights = dat$n)
}
hyper2 <- function(dat, ...){
    requireNamespace("hyper2")
    # create likelihood object based on ordered items and counts
    H <- hyper2::hyper2(pnames = paste0("p", seq_len(ncol(dat) - 1)))
    for (i in seq_len(nrow(dat))){
        x <-  dat[i, -1][dat[i, -1] > 0]
        H <- H + hyper2::order_likelihood(x, times = dat[i, 1])
    }
    # find parameters to maximise likelihood
    p <- hyper2::maxp(H)
    structure(p, loglik = hyper2::lhyper2(H, p[-length(p)]))
}
plmix <- function(dat, ...){
    requireNamespace("PLMIX")
    # disaggregate data (no functionality for weights or counts)
    r <- rep(seq_len(nrow(dat)), dat$n)
    # maximum a posteriori estimate, with non-informative prior
    # K items in each ranking, single component distribution
    # default starting values do not always work so specify as uniform
    K <- ncol(dat)
    PLMIX::mapPLMIX(as.matrix(dat[r, -1]), K = K, G = 1,
                    init = list(p = rep.int(1/K, K)), plot_objective = FALSE)
}
pmr <- function(dat, ...){
    requireNamespace("pmr")
    # convert ordered items to ranking
    R <- as.rankings(dat[,-1], "ordering")
    # create data frame with counts as required by pl
    X <- as.data.frame(unclass(R))
    X$n <- dat$n
    capture.output(res <- pmr::pl(X))
    res
}
statrank <- function(dat, iter){
    requireNamespace("StatRank")
    # disaggregate data (no functionality for weights or counts)
    r <- rep(seq_len(nrow(dat)), dat$n)
    capture.output(res <- StatRank::Estimation.PL.MLE(as.matrix(dat[r, -1]),
                                                      iter = iter))
    res
}

When recording timings, the number of iterations for StatRank was set so that the log-likelihood on exit was equal to the log-likelihood returned by the other functions with relative tolerance 1e-6.

timings <- function(dat, iter = NULL,
                    fun = c("pl", "hyper2", "plmix", "pmr", "statrank")){
    res <- list()
    for (nm in c("pl", "hyper2", "plmix", "pmr", "statrank")){
        if (nm %in% fun){
            res[[nm]] <- suppressWarnings(
                system.time(do.call(nm, list(dat, iter)))[["elapsed"]])
        } else res[[nm]] <- NA
    }
    res
}
netflix_timings <- timings(netflix, 6)
tshirt_timings <- timings(tshirt, 341,
                          fun = c("pl", "hyper2", "plmix", "statrank"))
sushi_timings <- timings(sushi, 5,
                         fun = c("pl", "hyper2", "plmix", "statrank"))

5.2 beans data preparation

First we handle the best and worst rankings. These give the variety the farmer thought was best or worst, coded as A, B or C for the first, second or third variety assigned to the farmer respectively.

data(beans)
head(beans[c("best", "worst")], 2)

## # A tibble: 2 x 2
##   best  worst
##   <chr> <chr>
## 1 C     A    
## 2 B     A

We convert these to numeric values, allowing us to impute the middle-ranked variety (a strict ranking is assumed, so the sum of each row must be 6)

beans <- within(beans, {
    best <- match(best, c("A", "B", "C"))
    worst <- match(worst, c("A", "B", "C"))
    middle <- 6 - best - worst
})
head(beans[c("best", "middle", "worst")], 2)

## # A tibble: 2 x 3
##    best middle worst
##   <int>  <dbl> <int>
## 1     3     2.     1
## 2     2     3.     1

This gives an ordering of the three varieties the farmer was given. The names of these varieties are stored in separate columns

varieties <- as.matrix(beans[c("variety_a", "variety_b", "variety_c")])
head(varieties, 2)

##      variety_a     variety_b         variety_c    
## [1,] "BRT 103-182" "SJC 730-79"      "PM2 Don Rey"
## [2,] "INTA Rojo"   "INTA Centro Sur" "INTA Sequia"

So we can convert the variety IDs to the variety names

n <- nrow(beans)
beans <- within(beans, {
    best <- varieties[cbind(seq_len(n), best)]
    worst <- varieties[cbind(seq_len(n), worst)]
    middle <- varieties[cbind(seq_len(n), middle)]
})
head(beans[c("best", "middle", "worst")], 2)

## # A tibble: 2 x 3
##   best            middle      worst      
##   <chr>           <chr>       <chr>      
## 1 PM2 Don Rey     SJC 730-79  BRT 103-182
## 2 INTA Centro Sur INTA Sequia INTA Rojo

Next we convert these orderings to sub-rankings of the full set of varieties, including the local variety as an additional item, so that we can add the paired comparisons shortly:

lab <- c("Local", sort(unique(as.vector(varieties))))
R <- as.rankings(beans[c("best", "middle", "worst")],
                 input = "ordering", labels = lab)
head(R)

## [1] "PM2 Don Rey > SJC 730-79 > BRT 103-182"  
## [2] "INTA Centro Sur > INTA Sequia > INTA ..."
## [3] "INTA Ferroso > INTA Matagalpa > BRT  ..."
## [4] "INTA Rojo > INTA Centro Sur > ALS 0532-6"
## [5] "PM2 Don Rey > INTA Sequia > SJC 730-79"  
## [6] "ALS 0532-6 > INTA Matagalpa > INTA Rojo"

The comparisons with the local variety are stored in another set of columns

head(beans[c("var_a", "var_b", "var_c")], 2)

## # A tibble: 2 x 3
##   var_a var_b  var_c 
##   <chr> <chr>  <chr> 
## 1 Worse Worse  Better
## 2 Worse Better Better

The following converts each of these columns to a matrix of ordered pairs:

paired <- list()
for (id in c("a", "b", "c")){
    ordering <- matrix("Local", nrow = n, ncol = 2)
    worse <- beans[[paste0("var_", id)]] == "Worse"
    ## put trial variety in column 1 when it is not worse than local
    ordering[!worse, 1] <- beans[[paste0("variety_", id)]][!worse]
    ## put trial variety in column 2 when it is worse than local
    ordering[worse, 2] <- beans[[paste0("variety_", id)]][worse]
    paired[[id]] <- ordering
}
head(paired[["a"]])

##      [,1]              [,2]         
## [1,] "Local"           "BRT 103-182"
## [2,] "Local"           "INTA Rojo"  
## [3,] "INTA Ferroso"    "Local"      
## [4,] "INTA Centro Sur" "Local"      
## [5,] "INTA Sequia"     "Local"      
## [6,] "ALS 0532-6"      "Local"

Again we convert these orderings to sub-rankings of the full set of varieties and combine them with the rankings of order three:

paired <- lapply(paired, as.rankings, input = "ordering", labels = lab)
R <- rbind(R, paired[["a"]], paired[["b"]], paired[["c"]])
head(R)

## [1] "PM2 Don Rey > SJC 730-79 > BRT 103-182"  
## [2] "INTA Centro Sur > INTA Sequia > INTA ..."
## [3] "INTA Ferroso > INTA Matagalpa > BRT  ..."
## [4] "INTA Rojo > INTA Centro Sur > ALS 0532-6"
## [5] "PM2 Don Rey > INTA Sequia > SJC 730-79"  
## [6] "ALS 0532-6 > INTA Matagalpa > INTA Rojo"

tail(R)

## [1] "INTA Sequia > Local"    "INTA Sequia > Local"   
## [3] "BRT 103-182 > Local"    "Local > INTA Matagalpa"
## [5] "Local > INTA Rojo"      "Local > SJC 730-79"