Assesses the goodness of fit of competing statistical models
# S3 method for class 'PlackettLuce'
anova(object, ...)
likelihood_ratio(x, split, ...)library("PlackettLuce")
example("beans", package = "PlackettLuce")
#>
#> beans> # Consider the best and worst rankings. These give the variety the
#> beans> # farmer thought was best or worst, coded as A, B or C for the
#> beans> # first, second or third variety assigned to the farmer
#> beans> # respectively.
#> beans> data(beans)
#>
#> beans> head(beans[c("best", "worst")], 2)
#> best worst
#> 1 C A
#> 2 B A
#>
#> beans> # Fill in the missing item
#> beans> beans$middle <- complete(beans[c("best", "worst")],
#> beans+ items = c("A", "B", "C"))
#>
#> beans> head(beans[c("best", "middle", "worst")], 2)
#> best middle worst
#> 1 C B A
#> 2 B C A
#>
#> beans> # This gives an ordering of the three varieties the farmer was
#> beans> # given. The names of these varieties are stored in separate
#> beans> # columns
#> beans> varieties <- beans[c("variety_a", "variety_b", "variety_c")]
#>
#> beans> 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
#>
#> beans> # Use these names to decode the orderings of order 3
#> beans> order3 <- decode(beans[c("best", "middle", "worst")],
#> beans+ items = beans[c("variety_a", "variety_b", "variety_c")],
#> beans+ code = c("A", "B", "C"))
#>
#> beans> # Now consider the paired comparisons agains the local variety
#> beans> head(beans[c("var_a", "var_b", "var_c")], 2)
#> var_a var_b var_c
#> 1 Worse Worse Better
#> 2 Worse Better Better
#>
#> beans> # Convert these results to a vector and get the corresponding trial variety
#> beans> outcome <- unlist(beans[c("var_a", "var_b", "var_c")])
#>
#> beans> trial_variety <- unlist(beans[c("variety_a", "variety_b", "variety_c")])
#>
#> beans> # Create a data frame of the implied orderings of order 2
#> beans> order2 <- data.frame(Winner = ifelse(outcome == "Worse",
#> beans+ "Local", trial_variety),
#> beans+ Loser = ifelse(outcome == "Worse",
#> beans+ trial_variety, "Local"),
#> beans+ stringsAsFactors = FALSE, row.names = NULL)
#>
#> beans> head(order2, 2)
#> Winner Loser
#> 1 Local BRT 103-182
#> 2 Local INTA Rojo
#>
#> beans> # Finally combine the rankings of order 2 and order 3
#> beans> R <- rbind(as.rankings(order3, input = "orderings"),
#> beans+ as.rankings(order2, input = "orderings"))
#>
#> beans> 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"
#>
#> beans> tail(R)
#> [1] "INTA Sequia > Local" "INTA Sequia > Local" "BRT 103-182 > Local"
#> [4] "Local > INTA Matagalpa" "Local > INTA Rojo" "Local > SJC 730-79"
G = group(R, rep(seq_len(nrow(beans)), 4))
d = cbind(G, beans)
split = ifelse(d$maxTN < 18.7175, TRUE, FALSE)
likelihood_ratio(G, split)
#> deviance DF_delta Pr(>Chisq)
#> <dbl> <dbl> <dbl> <chr>
#> 1: 42.5998 10.0000 0.0000 ***
mod = PlackettLuce(G)
anova(mod)
#> logLik DF Statistic Pr(>Chisq)
#> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1: NULL -3259.5513 5052.0000 NA NA
#> 2: Model -3192.5834 5042.0000 133.9360 0.0000 ***