In RobinHankin/hyper2: The Hyperdirichlet Distribution, Mark 2

set.seed(0)
knitr::opts_chunk$set(echo = TRUE)
library("hyper2")
library("magrittr")
options("digits" = 5)

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

To cite the hyper2 package in publications, please use @hankin2017_rmd. Here I analyse a karate dataset taken from wikipedia. Short story: the dataset has too many degrees of freedom to be easily analysed with hyper2 and Zermelo's iteration scheme as implemented by the BradleyTerry package exploits the special pairwise structure of the likelihood function very effectively.

karate_table  <- read.table("karate.txt",header=TRUE)
head(karate_table)

Converting it to a likelihood function is straightforward:

karate <- hyper2()
for(i in seq_len(nrow(karate_table))){

  p1 <- karate_table$karateka1[i]
  p2 <- karate_table$karateka2[i]

  w1 <- karate_table$wins1[i]
  w2 <- karate_table$wins2[i]

  karate[p1] %<>% inc(w1)
  karate[p2] %<>% inc(w2)
  karate[c(p1,p2)] %<>% dec(w1+w2)
}
summary(karate)

In the summary above (in the summary of bracket lengths), we see how the likelihood function comprises only singletons [which appear on the numerator] and pairs [on the denominator].

karate_maxp <- maxp(karate)

dotchart(log10(karate_maxp))

MM analysis

We will define a matrix of win-loss statistics:

pk <- pnames(karate)
M <- matrix(0,length(pk),length(pk))
rownames(M) <- pk
colnames(M) <- pk
o <- karate_table # saves typing
for(i in seq_len(nrow(o))){
  M[o[i,1],o[i,3]] %<>% `+`(o[i,2])
  M[o[i,3],o[i,1]] %<>% `+`(o[i,4])
}

pM <- pairwise(M)
pM == karate  # consistency check; should be TRUE

Now use Zermelo's iteration method and compare with hyper2's maxp():

karate_zermelo <- zermelo(M)
loglik(karate_zermelo,karate)
loglik(karate_maxp,karate)

Above we see that Zermelo gives a much more likely point in parameter space than maxp().

par(pty='s')
plot(log10(karate_zermelo),log10(karate_maxp),asp=1,xlim=c(-4,0),ylim=c(-4,0))
abline(0,1)
grid()

Package dataset {-}

Following lines create karate.rda, residing in the data/ directory of the package.

save(karate_table,karate,karate_maxp,karate_zermelo,file="karate.rda")

References {-}

RobinHankin/hyper2 documentation built on May 6, 2024, 4:25 p.m.