In RobinHankin/hyper2: The Hyperdirichlet Distribution, Mark 2

knitr::opts_chunk$set(echo = TRUE)

knitr::include_graphics(system.file("help/figures/hyper2.png", package = "hyper2"))

To cite the hyper2 package in publications, please use @hankin2017_rmd. Suppose there are five people in a family. Each day for five days, one person cooks a meal. After the fifth day, each person considers the four meals cooked by other family members, and puts them in order of preference. This very short vignette shows how to apply generalized Bradley-Terry to this situation using the hyper2 package.

library(hyper2)

Now add some data:

H <- hyper2() + 
  race(c("p5","p2","p3","p4")) + # person 1 says p5 is the best, p4 the worst.
  race(c("p1","p5","p3","p4")) + # person 2 says p1 is the best, p4 the worst.
  race(c("p1","p4","p2","p5")) + # person 3 says p1 is the best, p5 the worst.
  race(c("p1","p2","p3","p5")) + # person 4 says p1 is the best, p5 the worst.
  race(c("p2","p1","p3","p4"))   # person 5 says p2 is the best, p4 the worst.
H

Note that noone is allowed to vote for themselves. We may find maximum likelihood estimate for the strengths, $\operatorname{argmax}\mathcal{L}\left(p_1,\ldots,p_5\right)$:

maxp(H)

Graphically:

dotchart(maxp(H),pch=16,main="probabilities")
dotchart(log(maxp(H)),pch=16,main="log probabilities")

We can assess the hypothesis that all players have the same strength:

equalp.test(H)

Thus we reject the hypothesis of equal strength. Now we can assess the hypothesis that person a does in fact have the highest strength of the five. We can follow the reasoning in the icons demonstration and test $H_0\colon p_1=p_2$:

samep.test(H,c("p1","p2"))

There is no strong evidence to support the assertion that person 1 is actually stronger than person 2 in the sense that $p_1>p_2$.

Suspect observation

Suppose that we subsequently observe order statistic p5,p4,p3,p2,p1, that is, person 5 is the best, 4 is the second best, and so on to 1 being the worst. Is this observation consistent with the previous dataset?

It is possible to perform a permutation test on this as follows. First, calculate the probability of each of the $5!=120$ possible observations (up to a constant):

library(partitions)
M <- perms(5)
M[] <- paste("p",M,sep="")
imH <- indep(maxp(H))
f <- function(o){loglik(imH,race(o),log=FALSE)}
LL <- apply(perms(5),2,f)
length(LL)
sum(LL)
head(LL)

So LL is the probability of observing each permutation of 5 objects. We then calculate the $p$-value as the probability of obtaining the observation or an observation more extreme; in this case we define a permutation to be "more extreme" if it has a smaller probability than the observed permutation of 5,4,3,2,1 (under the null). This is easily calculated:

obs <- c(5,4,3,2,1)
(pval <- sum(LL[LL <= f(obs)])/sum(LL))

Thus it is reasonable to reject the hypothesis that this particular observation was drawn from the same distribution as the others, and infer that it was suspect.

Changes in cooking strength

We now consider a sequence of observations where one person, say person 1, is suspected of an increase in strength (perhaps they attended a domestic cookery course). We make some observations before and after the training and seek evidence that it was effective. Unlike the previous cases, here we have a full order statistic without the complication of not being allowed to vote for oneself (perhaps a professional chef puts the dishes in order). First, before:

# observations before training:
Hbefore <- 
  race(c("p5","p2","p3","p4","p1")) + #  5 is the best, 4 the worst.
  race(c("p5","p3","p2","p4","p1")) + #  5 is the best, 1 the worst.
  race(c("p4","p5","p3","p2","p1")) + #  4 is the best, 1 the worst.
  race(c("p5","p3","p4","p1","p2")) + #  5 is the best, 2 the worst.
  race(c("p2","p5","p3","p1","p4"))   #  2 is the best, 4 the worst.
Hbefore

Now, after the training we have two finishing order observations: p1,p2,p3,p4,p5 and p3,p1,p5,p2,p4. Thus in the first observation, p1 comes first and in the second observation, p1 comes second second (after p3). The likelihood for the first observation would be

[ \frac{p_1+S}{p_1+S+p_2+p_3+p_4+p_5}\cdot \frac{p_2 }{ p_2+p_3+p_4+p_5}\cdot \frac{p_3 }{ p_3+p_4+p_5}\cdot \frac{p_4 }{ p_4+p_5}\cdot \frac{p_5 }{ p_5} ]

and the likelihood for the second observation would be

[ \frac{p_3 }{p_1+S+p_2+p_3+p_4+p_5}\cdot \frac{p_1+S}{p_1+S+p_2+ p_4+p_5}\cdot \frac{p_5 }{ p_2+ p_4+p_5}\cdot \frac{p_2 }{ p_2+ p_4 }\cdot \frac{p_4 }{ p_4 }\cdot ]

We can translate these into package idiom with function pwa() ("person with advantage"):

Hafter <- race(c("p1","p2","p3","p4","p5")) + race(c("p3","p1","p5","p2","p4"))
Hafter %<>%  pwa("p1")
Hafter

In the above, the additional strength conferred by the training is represented by S. We can estimate the effect of training:

maxp(Hbefore+Hafter)

[In the above, note the use of the overloaded "+", meaning to add two likelihood functions for independent observations]. The evaluate suggests that the training had a sizable effect; we may use \code{specficp.gt.test()} to test a null of zero strength:

specificp.gt.test(Hbefore+Hafter, "S", 0)

thus we have strong evidence for the training being effective.

Analysis with `hyper3` objects

We can take the "changes in cooking strength" augmented Plackett-Luce likelihood functions above, and apply hyper3 idiom to it. Instead of using a reified entity $S$ we have likelihood functions

[ \frac{\lambda p_1}{\lambda p_1+p_2+p_3+p_4+p_5}\cdot \frac{p_2 }{ p_2+p_3+p_4+p_5}\cdot \frac{p_3 }{ p_3+p_4+p_5}\cdot \frac{p_4 }{ p_4+p_5}\cdot \frac{p_5 }{ p_5} ]

[ \frac{p_3 }{\lambda p_1+p_2+p_3+p_4+p_5}\cdot \frac{\lambda p_1}{\lambda p_1+p_2+ p_4+p_5}\cdot \frac{p_5 }{ p_2+ p_4+p_5}\cdot \frac{p_2 }{ p_2+ p_4 }\cdot \frac{p_4 }{ p_4 }\cdot ]

and then test $H_0\colon\lambda=1$. The hyper3 likelihood function is mathematically the same as the reified hyper2 BT likelihood function provided $p_1+S=\lambda p_1$ (the likelihood functions differ only if more than one player has increased strength); note that $\lambda < 1$ is not possible with reified BT because it implies negative $S$. We will create a hyper3 object for the post-training observations:

Hpost1 <- function(lambda){
   out <- hyper3()
   out[c(p1=lambda)] %<>% inc
   out[c(p1=lambda,p2=1,p3=1,p4=1,p5=1)] %<>% dec

   out[c(p2=1)] %<>% inc
   out[c(p2=1,p3=1,p4=1,p5=1)] %<>% dec

   out[c(p3=1)] %<>% inc
   out[c(p3=1,p4=1,p5=1)] %<>% dec

   out[c(p4=1)] %<>% inc
   out[c(p4=1,p5=1)] %<>% dec
   return(out)
}

Hpost2 <- function(lambda){
   out <- hyper3()
   out[c(p3=1)] %<>% inc
   out[c(p1=lambda,p2=1,p3=1,p4=1,p5=1)] %<>% dec

   out[c(p1=lambda)] %<>% inc
   out[c(p1=lambda,p2=1,p4=1,p5=1)] %<>% dec

   out[c(p5=1)] %<>% inc
   out[c(p2=1,p4=1,p5=1)] %<>% dec

   out[c(p2=1)] %<>% inc
   out[c(p2=1,p4=1)] %<>% dec
   return(out)
}

Now The likelihood functions for the "before" observations simply add, so for example

as.hyper3(Hbefore) + Hpost1(1.8888) + Hpost2(1.8888)

Then we can calculate a profile likelihood curve:

H3 <- function(lambda){as.hyper3(Hbefore) + Hpost1(lambda) + Hpost2(lambda)}
supplam <- function(lambda){maxp(H3(lambda),n=1,give=TRUE)$likes}
lam <- exp(seq(from=log(0.8),to=log(100),len=50))
L <- sapply(lam,supplam)

plot(lam,L-max(L),type="l",lwd=2)
abline(h=c(0,-2),col='gray')
abline(v=1,col='red')
text(0,-0.2,"lambda=1",col="red",pos=4)

plot(log(lam),L-max(L),type="l",lwd=2)
abline(h=c(0,-2),col='gray')
abline(v=1,col='red')
text(0,-0.2,"lambda=1",col="red",pos=4)

We see that we can reject $H_0\colon\lambda=1$ but for this dataset there does not appear to be a clear upper bound on $\lambda$.

We can calculate the maximum likelihood estimate as follows:

o <- optimize(supplam,c(40,70),maximum=TRUE)
o

To calculate a formal likelihood ratio we need to calculate the support for $\lambda=1$ and compare it with the support for $\lambda=\hat{\lambda}$

supplam(o$maximum)
supplam(1)

and we get a $p$-value for $H_0\colon\lambda=1$ using Wilks, specifically that $2\log\left(\frac{\mathcal{L}{\lambda=\hat{\lambda}}}{\mathcal{L}{\lambda=1}}\right)\sim\chi^2_1$

pchisq(2*(supplam(o$maximum)-supplam(1)),lower.tail=FALSE,df=1)

Actually we should be able to analyse $\lim_{\lambda\longrightarrow\infty}H_\lambda$:

lambig <- 10^(0:6)
LL <- sapply(lambig,supplam)

plot(log(lambig),LL-max(LL),pch=16)
for(i in seq(from=0,to=8,by=2)){abline(i,-1,col='gray',lwd=0.5)}
abline(v=log(o$maximum))

Looking at the two parts of the likelihood function that includes lambda, and thinking about $\lim_{\lambda\longrightarrow\infty}$ we see that the first product is approximately constant and the second is $\propto\lambda^{-1}$. This analysis is a little bit naive because as $\lambda$ grows very large, we might expect $\hat{p_3}$ to grow too. But:

M <- rbind(
maxp(H3(10^0    ),n=1,give=TRUE)$par,
maxp(H3(10^1    ),n=1,give=TRUE)$par,
maxp(H3(10^2    ),n=1,give=TRUE)$par,
maxp(H3(10^3    ),n=1,give=TRUE)$par,
maxp(H3(10^4    ),n=1,give=TRUE)$par,
maxp(H3(10^5    ),n=1,give=TRUE)$par,
maxp(H3(10^6    ),n=1,give=TRUE)$par
)
M

matplot(M,type='b')

Above we see that the evaluate hardly changes when $\lambda > 10^4$, and indeed from the figure we see that the slope of the support graph is about $-1$ for large values of $\lambda$ [the points follow a line parallel to one of the slanted gray lines, which have a slope of exactly $-1$; note that the aspect ratio of the graph is not unity]