Person with advantage functionality for hyper3 objects: pwa3()
Robin K. S. Hankin
(takes about five minutes to run without cache)
pwa3## function (H3, ...)
## {
## l <- list(...)
## if (length(l) == 1) {
## x <- l[[1]]
## if (is.vector(x) && (!is.list(x)) && is.numeric(x) &&
## (!is.null(names(x)))) {
## l <- as.list(l[[1]])
## }
## }
## for (i in seq_along(l)) {
## H3 <- pwa3_single(H3, names(l[i]), l[[i]])
## }
## return(H3)
## }
## <bytecode: 0x55fd645e2350>
## <environment: namespace:hyper2>pwa3_single## function (H3, pwa, lambda)
## {
## stopifnot(pwa %in% pnames(H3))
## stopifnot(length(pwa) == 1)
## H3 <- as.hyper3(H3)
## B <- elements(brackets(H3))
## W <- elements(weights(H3))
## for (i in seq_along(B)) {
## W[[i]][B[[i]] %in% pwa] <- lambda
## }
## return(hyper3_bw(B, W, powers(H3), pnames(H3)))
## }
## <bytecode: 0x55fd645e11c0>
## <environment: namespace:hyper2>To cite the hyper2 package in publications, please use
Hankin (2017). Here we analyse some simple Plackett-Luce
observations, using hyper3 likelihood functions.
Our observations are a total of six order statistics for four
competitors a,b,c,d. These competitors race three times under
Plackett-Luce probabilities. They then race three times, but with
competitor b having some advantage (steroids?), parametrized by his
Plackett-Luce strength being \(\lambda b\).
Our object of inference is \(\lambda\); the Plackett-Luce strengths
a,b,c,d are nuisance parameters.
lam <- function(lambda){
race_nohelp <-
race(c("b","a","d","c")) +
race(c("c","b","a","d")) +
race(c("b","d","a","c"))
race_help <-
race(c("a","c","d","b")) |> pwa3(b = lambda) +
race(c("c","d","b","a")) |> pwa3(b = lambda) +
race(c("d","c","a","b")) |> pwa3(b = lambda)
maxp(as.hyper3(race_nohelp) + race_help, give=TRUE)$value
}Above we see that in the unaided case (the first three races), b
does quite well, coming first twice and second once. In the aided
case, b does much worse, coming either third or fourth. Evidently,
the “help” actually hinders b and we suspect that \(\lambda<1\).
l <- exp(seq(from=log(0.0002), to=log(2), len=17))
like <- sapply(l,lam)optimal_sol <- optimize(lam, c(0.01, 3), maximum=TRUE)optimal_sol## $maximum
## [1] 0.024679
##
## $objective
## [1] -15.486f <- function(l){lam(l) - optimal_sol$objective + 2}
lower <- uniroot(f, interval = c(exp(-8), exp(-7)))
upper <- uniroot(f, interval = c(exp(-2), exp( 0)))plot(log(l), like - optimal_sol$objective,
type='b', pch=16,
xlab = expression(log(lambda)),
ylab = "log-likelihood")
abline(h = c(0, -2), col='gray')
abline(v = log(optimal_sol$maximum))
abline(v = 0, lty=2)
segments(x0 = log(c(lower$root, upper$root)), y0 = -2.25, y1 = -1.75, col='gray')
Figure 1: Log-likelihood function for lambda
Figure 1 shows a log-likelihood function for \(\log\lambda\). We see that the evaluate is about \(\lambda=\exp(-3.7)\simeq 0.025\). Further, the support interval runs from about \(\exp(-7.56)\simeq 5.2\times 10^{-4}\) to about \(\exp(-0.87)\simeq 0.42\) thereby excluding the null of \(\lambda=1\) which has a support of about -3.54.
Another example
In this case, competitor a receives help \(\lambda\) or \(\lambda^n\) if
he gets \(n\) sets of help. Below, in function lam2() we see that, in
general, the more help a gets, the better he does. See how a
negative power corresponds to some sort of hindrance or penalty, and a
zero power (of course) corresponds to no help at all.
lam2 <- function(lam){
race_help <-
race(c("d","c","a","b")) |> pwa3(a = lam^0 ) +
race(c("c","b","d","a")) |> pwa3(a = lam^-3) +
race(c("c","b","d","a")) |> pwa3(a = lam^-2) +
race(c("c","b","a","d")) |> pwa3(a = lam^-1) +
race(c("a","c","b","d")) |> pwa3(a = lam ) +
race(c("a","c","b","d")) |> pwa3(a = lam^3 ) +
race(c("a","b","d","c")) |> pwa3(a = lam ) +
race(c("c","a","d","b")) |> pwa3(a = lam^2 )
maxp(race_help, give=TRUE)$value
}l <- exp(seq(from=log(0.8), to=log(19.1), len=11))
like2 <- sapply(l,lam2)plot(log(l),like2-max(like2),type="b",pch=16)
abline(h=c(0,-2), col='gray')
abline(v=0, lty=2)
Another other example
In this example we see help being given to two competitors at the same time.
myr <- function(a1,a2,a3,a4,a5){race(letters[c(a1,a2,a3,a4,a5)])}
lam3 <- function(l){
race_help <-
myr(3,1,2,4,5) |> pwa3(a=l, b=l) +
myr(2,3,1,4,5) |> pwa3(a=l, b=l) +
myr(3,1,5,4,2) |> pwa3(a=1, b=1) + # no help!
myr(4,2,5,3,1) |> pwa3(a=1, b=1) + # no help!
myr(3,1,2,4,5) |> pwa3(a=l, b=l)
maxp(race_help, give=TRUE)$value
}l <- exp(seq(from=log(0.5), to=log(84.1), len=13))
like3 <- sapply(l,lam3)plot(log(l),like3-max(like3), type="b", pch=16)
abline(h=c(0,-2), col='gray')
abline(v=0, lty=2)
Yet another example
Here we consider help in turn to each of five competitors. From the function below, if someone gets help they generally do quite well.
myr <- function(a1,a2,a3,a4,a5){race(letters[c(a1,a2,a3,a4,a5)])}
lam4 <- function(x, givelike=FALSE){
race_help <-
myr(3,1,2,4,5) |> pwa3(a=x) +
myr(2,3,1,4,5) |> pwa3(b=x) +
myr(3,1,5,4,2) |> pwa3(c=x) +
myr(4,2,5,3,1) |> pwa3(d=x) +
myr(3,5,2,4,1) |> pwa3(e=x)
if(givelike){
return(race_help)
} else {
return(maxp(race_help, give=TRUE)$value)
}
}l <- exp(seq(from= -0.5, to=log(250), len=15))
like4 <- sapply(l,lam4)optimal_sol4 <- optimize(lam4, exp(2:3), maximum=TRUE)It makes sense to consider our likelihood function at the estimate for \(\lambda\):
plot(log(l),like4-max(like4), type="b", pch=16)
abline(h=c(0,-2), col='gray')
abline(v=0, lty=2)
abline(v = log(optimal_sol4$maximum))
H <- lam4(optimal_sol4$maximum, givelike=TRUE)
H## log( (a=1)^2 * (a=1, b=1, c=1, d=1, e=16.705)^-1 * (a=1, b=1, c=1, d=16.705, e=1)^-1 * (a=1, b=1, c=1,
## e=1)^-1 * (a=1, b=1, c=16.705, d=1, e=1)^-1 * (a=1, b=1, d=1)^-1 * (a=1, b=1, d=1, e=1)^-1 * (a=1, b=1,
## d=1, e=16.705)^-1 * (a=1, b=16.705, c=1, d=1, e=1)^-1 * (a=1, c=1)^-1 * (a=1, c=1, d=1, e=1)^-1 * (a=1,
## c=1, e=1)^-1 * (a=1, d=1)^-1 * (a=1, d=1, e=1)^-1 * (a=16.705)^1 * (a=16.705, b=1, c=1, d=1, e=1)^-1 *
## (a=16.705, b=1, d=1, e=1)^-1 * (b=1)^3 * (b=1, d=1)^-1 * (b=1, d=1, e=1)^-2 * (b=16.705)^1 * (c=1)^4 *
## (c=16.705)^1 * (d=1)^4 * (d=1, e=1)^-2 * (d=16.705)^1 * (e=1)^2 * (e=16.705)^1)summary(H)## A hyper3 object of size 5.
## pnames: a b c d e
## Number of brackets: 28
## Sum of powers: 0
##
## Table of bracket lengths:
## 1 2 3 4 5
## 10 4 4 5 5
##
## Table of powers:
## -2 -1 1 2 3 4
## 2 16 5 2 1 2
##
## Table of weights:
## 1 16.7049918243144
## 63 12maxp(H)## a b c d e
## 0.073862 0.171647 0.589727 0.106656 0.058108Above we see that the evaluate is very high for c having a BT
strength of about 0.6. This is consistent with our observations for
twice we see c winning [that is, coming first] without any help.
Is there statistical evidence to support this assertion?
The simplest approach would be to consider the likelihood function at \(\hat{\lambda}\simeq 16.705\):
f <- function(o){samep.test(H, c(o, "c"))$statistic}
stat4 <- c(a=f("a"), b=f("b"), d=f("d"), e=f("e"))stat4## a b d e
## 2.39028 0.74467 1.51336 2.64778Above we see that there is no strong evidence for c being the
strongest, on the grounds that we are unable to reject \(H_0\colon
p_b=p_c\) [the support difference is about 0.74, short of the two-units
of support criterion required for rejection], at least with
\(\lambda=16.705\).
Another further idea
Now we have a model for the starting grid in formula 1. In the following function we specify the finishing order and also the starting grid order.
tolam <- function(lam,o){setNames(lam^(seq_along(o)-1),o)}
lam5 <- function(lam, givelike=FALSE){
race_help <-
race(c("a","b","d","c","e")) |> pwa3(tolam(lam, c("a","b","c","d","e"))) +
race(c("a","c","b","d","e")) |> pwa3(tolam(lam, c("c","b","a","d","e"))) +
race(c("a","d","c","b","e")) |> pwa3(tolam(lam, c("a","c","b","d","e"))) +
race(c("e","c","b","a","d")) |> pwa3(tolam(lam, c("a","c","b","e","d"))) +
race(c("b","c","e","a","d")) |> pwa3(tolam(lam, c("c","b","e","d","a"))) +
race(c("a","d","c","b","e")) |> pwa3(tolam(lam, c("a","b","c","d","e")))
if(givelike){
return(race_help)
} else {
return(maxp(race_help, give=TRUE)$value)
}
}l <- exp(seq(from= -2, to=log(1.5), len=15))
like5 <- sapply(l,lam5)
optimal_sol5 <- optimize(lam5, exp(c(-1,0)), maximum=TRUE)f <- function(l){lam5(l) - optimal_sol5$objective + 2}
lower <- uniroot(f, interval = c(exp(-1.5), exp(-1.0)))
upper <- uniroot(f, interval = c(exp(-0.5), exp(+0.5)))plot(log(l), like5 - max(like5), type='b', pch=16)
abline(h = c(0, -2), col='gray')
abline(v = log(optimal_sol5$maximum))
abline(v = 0, lty=2)
segments(x0 = log(c(lower$root, upper$root)), y0 = -2.5, y1 = -1.5, col='red')
Figure 2: Log-likelihood function for lambda, the pole position measure
H <- lam5(optimal_sol5$maximum, givelike=TRUE)
mP <- maxp(H)mP## a b c d e
## 0.30881 0.16141 0.16307 0.24192 0.12478References
hyper2 Package: Likelihood Functions for Generalized Bradley-Terry Models.” The R Journal 9 (2): 429–39.