test-pltree.R
In PlackettLuce: Plackett-Luce Models for Rankings

coef_tol <- 1e-06
loglik_tol <- 1e-07

if (require(psychotree) & require(sandwich)){
    data("Topmodel2007", package = "psychotree")
    preference <- Topmodel2007$preference

    # plfit vs btmodel
    btmod <- btmodel(preference, ref = "Anja")

    G <- as.grouped_rankings(preference)
    plmod <- plfit(G, npseudo = 0, ref = "Anja", estfun = TRUE, object = TRUE)

    # some conditions with psychotree 0.15.1 and psychotools 0.4.2
    test_that("plfit and bt model agree [Topmodel2007]",
              {
                  # coef
                  expect_equal(coef(btmod), coef(plmod)[-6])
                  # estfun - only agree if last level taken as ref in btmodel
                  E1 <- psychotools::estfun.btmodel(btmod)
                  E2 <- plmod$estfun
                  expect_equal(unname(E1), unname(E2))
                  # worth and vcov
                  # - vcov only agrees if first level taken as ref in btmodel
                  btmod <- btmodel(preference, ref = "Barbara")
                  w1 <- worth(btmod, log = FALSE, ref = 1)
                  w2 <- worth(plmod$object, log = FALSE, ref = 1)
                  expect_equal(coef(w1), coef(w2))
                  expect_equal(attr(w1, "vcov"), attr(w2, "vcov"),
                               tolerance = coef_tol)
              })
    # pltree vs bttree
    bt_tree <- bttree(preference ~ ., data = Topmodel2007, minsize = 5,
                      ref = "Anja")
    pl_tree <- pltree(G ~ ., npseudo = 0, data = Topmodel2007[, -1],
                      minsize = 5, ref = "Anja")

    test_that("pltree and bttree agree [Topmodel2007]",
              {
                 expect_equal(itempar(bt_tree), itempar(pl_tree))
              })
    # predict pltree
    # work round some bugs in psychotree 0.15.1
    newdata <- Topmodel2007[1:3,]
    test_that('predict.pltree works for type = "itempar" [Topmodel2007]',
              {
                  # probabilities
                  expect_equal(itempar(bt_tree)[
                      as.character(predict(bt_tree, type = "node")),],
                      predict(pl_tree, type = "itempar"),
                      ignore_attr = TRUE)
                  expect_equal(itempar(bt_tree)[
                      as.character(predict(bt_tree, newdata = newdata,
                                           type = "node")),],
                      predict(pl_tree,  newdata = newdata, type = "itempar"),
                      ignore_attr = TRUE)
                  # log-abilities
                  expect_equal(itempar(bt_tree, log = TRUE)[
                      as.character(predict(bt_tree, type = "node")),],
                      predict(pl_tree, type = "itempar", log = TRUE),
                      ignore_attr = TRUE)
                  expect_equal(itempar(bt_tree, log = TRUE)[
                      as.character(predict(bt_tree, newdata = newdata,
                                           type = "node")),],
                      predict(pl_tree, newdata = newdata, type = "itempar",
                              log = TRUE),
                      ignore_attr = TRUE)
              })
    test_that('predict.pltree works for type = "rank" [Topmodel2007]',
              {
                  rank <- t(apply(-itempar(bt_tree), 1, rank))
                  expect_equal(rank[
                      as.character(predict(bt_tree, newdata = newdata,
                                           type = "node")),],
                      predict(pl_tree, newdata = newdata, type = "rank"),
                      ignore_attr = TRUE)
              })
    test_that('predict.pltree works for type = "best" [Topmodel2007]',
              {
                  expect_equal(names(predict(bt_tree, type = "best")),
                               names(predict(pl_tree, type = "best")))
                  expect_equal(as.character(predict(bt_tree, type = "best")),
                               as.character(predict(pl_tree, type = "best")))
                  expect_equal(names(predict(bt_tree, newdata = newdata,
                                             type = "best")),
                               names(predict(pl_tree,  newdata = newdata,
                                             type = "best")))
                  expect_equal(as.character(predict(bt_tree, newdata = newdata,
                                                    type = "best")),
                               as.character(predict(pl_tree,  newdata = newdata,
                                                    type = "best")))
              })
    test_that('predict.pltree works for type = "node" [Topmodel2007]',
              {
                  tmp <- predict(pl_tree, type = "node")
                  mode(tmp) <- "numeric"
                  expect_equal(predict(bt_tree, type = "node"), tmp)
                  tmp <- predict(pl_tree,  newdata = newdata, type = "node")
                  mode(tmp) <- "numeric"
                  expect_equal(predict(bt_tree, newdata = newdata,
                                       type = "node"), tmp)
              })
    test_that('AIC.pltree works [Topmodel2007]',
              {
                  Topmodel2007$G <- G
                  aic1 <- AIC(pl_tree)
                  aic2 <- AIC(pl_tree, newdata = Topmodel2007)
                  expect_equal(aic1, aic2)

                  node <- predict(pl_tree, type = "node")
                  aic1 <- AIC(pl_tree, newdata = Topmodel2007[node == "3", -1])
                  aic2 <- -2*as.numeric(logLik(pl_tree[[3]])) +
                      2*attr(logLik(pl_tree), "df")
                  expect_equal(aic1, aic2)
              })
    # coef pltree
    test_that('coef pltree works with log = FALSE [Topmodel2007]',
              {
                  expect_equal(coef(pl_tree, log = FALSE),
                               itempar(pl_tree, log = FALSE, vcov = FALSE))
                  expect_equal(coef(pl_tree)[, -6],
                               coef(bt_tree))
              })
}

# example with weights

data(beans, package = "PlackettLuce")

# Fill in the missing ranking
beans$middle <- complete(beans[c("best", "worst")],
                         items = c("A", "B", "C"))

# Use these names to decode the orderings of order 3
order3 <- decode(beans[c("best", "middle", "worst")],
                 items = beans[c("variety_a", "variety_b", "variety_c")],
                 code = c("A", "B", "C"))

# Convert these results to a vector and get the corresponding trial variety
outcome <- unlist(beans[c("var_a", "var_b", "var_c")])
trial_variety <- unlist(beans[c("variety_a", "variety_b", "variety_c")])

# Create a data frame of the implied orderings of order 2
order2 <- data.frame(Winner = ifelse(outcome == "Worse",
                                     "Local", trial_variety),
                     Loser = ifelse(outcome == "Worse",
                                    trial_variety, "Local"),
                     stringsAsFactors = FALSE, row.names = NULL)

# Finally combine the rankings of order 2 and order 3
R <- rbind(as.rankings(order3, input = "ordering"),
           as.rankings(order2, input = "ordering"))

# Group the rankings by the corresponding farm
G <- group(R, rep(seq_len(nrow(beans)), 4))

weights <- c(rep(0.3, 177), rep(2, 481), rep(0.3, 184))

pl_tree <- pltree(G ~ season,
                  data = beans, alpha = 0.05, weights = weights )

test_that('grouped_rankings work w/ weights [beans]', {
    mod1 <- PlackettLuce(G, weights = weights)
    mod2 <- PlackettLuce(R, weights =rep(weights, 4))
    # iter can be different on some platform due to small difference in rowsums
    nm <- setdiff(names(mod1), c("call", "ranker", "iter"))
    expect_equal(mod1[nm], mod2[nm])
})

test_that('itempar.pltree works w/ weights [beans]',
          {
              # same results with newdata as original data
              pred1 <- predict(pl_tree, newdata = beans)
              pred2 <- predict(pl_tree)
              expect_equal(pred1, pred2)

              itempar1 <- unique(pred2)
              itempar2 <- itempar(pl_tree)
              # make sure nodes ordered the same way around
              expect_equal(itempar1[order(itempar1[,1]),],
                           itempar2[order(itempar2[,1]),],
                           ignore_attr = TRUE)
          })

test_that('AIC.pltree works w/ weights [beans]',
          {
              beans$G <- G
              aic1 <- AIC(pl_tree)
              aic2 <- AIC(pl_tree, newdata = beans, weights = weights)
              expect_equal(aic1, aic2)

              node <- predict(pl_tree, type = "node")
              id <- node == "3"
              aic1 <- AIC(pl_tree, newdata = beans[id, -1],
                          weights = weights[id])
              aic2 <- -2*as.numeric(logLik(pl_tree[[3]])) +
                  2*attr(logLik(pl_tree), "df")
              expect_equal(aic1, aic2)
          })

pl_tree1 <- pltree(G ~ maxTN, data = beans, alpha = 0)

test_that('itempar.pltree works w/ single node [beans]',  {
    # same results with newdata as original data
    pred1 <- predict(pl_tree1, newdata = beans)
    pred2 <- predict(pl_tree1)
    expect_equal(pred1, pred2)

    itempar1 <- unique(pred2)
    itempar2 <- itempar(pl_tree1)
    expect_equal(itempar1, itempar2, ignore_attr = TRUE)
})

test_that('AIC.pltree works w/ single node [beans]',
          {
              beans$G <- G
              aic1 <- AIC(pl_tree1)
              aic2 <- AIC(pl_tree1, newdata = beans)
              expect_equal(aic1, aic2)
          })

test_that('pltree works w/ estimated adherence [beans]', {
    pl_tree <- pltree(G ~ season,
                      data = beans, alpha = 0.05,
                      gamma = list(shape = 10, rate = 10))
    expect_snapshot_output(print(pl_tree))
})

test_that('pltree fails w/ fixed adherence [beans]', {
    expect_error(pltree(G ~ season,
                        data = beans, alpha = 0.05,
                        adherence = seq(0.75, 1.25, length.out = length(G))))
})