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. This short document discusses a small dataset first presented by @hankin2010. There, results from tennis games among four players was considered:

       {p1,p2} vs {p3,p4}   9-2  
       {p1,p3} vs {p2,p4}   4-4  
       {p1,p4} vs {p2,p3}   6-7  
          {p1} vs {p3}     10-14 
          {p2} vs {p3}     12-14 
          {p1} vs {p4}     10-14 
          {p2} vs {p4}     11-10 
          {p3} vs {p4}     13-13 

It is suspected that $p_1$ and $p_2$ have some form of team cohesion and play better when paired than when either solo or with other players. As the scores show, each player and, apart from the first line, each doubles partnership is of approximately the same strength.

We will create two likelihood functions: a regular one that assumes the stength of a doubles partnership is the sum of the strength of the players; and one that accounts for team cohesion.

H <- hyper2()

## 1&3 vs 2&4, scoreline 4-4:
H[c("p1","p3")] %<>% inc(4)
H[c("p2","p4")] %<>% inc(4)
H[c("p1","p2","p3","p4")] %<>% dec(8)

## 1&4 vs 2&3, scoreline 6-7:
H[c("p1","p4")] %<>% inc(6)
H[c("p2","p3")] %<>% inc(7)
H[c("p1","p2","p3","p4")] %<>% dec(13)

## 1 vs 3, scoreline 10-14:
H["p1"] %<>% inc(10)
H["p3"] %<>% inc(14)
H[c("p1","p3")] %<>% dec(24)

## 2 vs 3, scoreline 12-14:
H["p2"] %<>% inc(12)
H["p3"] %<>% inc(14)
H[c("p2","p3")] %<>% dec(26)

## 1 vs 4, scoreline 10-14:
H["p1"] %<>% inc(10)
H["p4"] %<>% inc(14)
H[c("p1","p4")] %<>% dec(24)

## 2 vs 4, scoreline 11-10:
H["p2"] %<>% inc(11)
H["p4"] %<>% inc(10)
H[c("p2","p4")] %<>% dec(21)

## 3 vs 4, scoreline 13-13:
H["p3"] %<>% inc(13)
H["p4"] %<>% inc(13)
H[c("p3","p4")] %<>% dec(26)

Thus H is a hyper2 object corresponding to results except the first line. We include this in two ways. Firstly, treating the first scoreline like the others:

tennis <- H

## 1&2 vs 3&4, scoreline 9-2:
tennis[c("p1","p2")] %<>% inc(9)
tennis[c("p3","p4")] %<>% inc(2)
tennis[c("p1","p2","p3","p4")] %<>% dec(11)

tennis

And secondly, accounting for any possible team cohesion by way of a ghost player G:

tennis_ghost <- H

## 1&2 vs 3&4 (NB: includes ghost!), scoreline 9-2 (again):
tennis_ghost[c("p1","p2","G")] %<>% inc(9)
tennis_ghost[c("p3","p4")] %<>% inc(2)
tennis_ghost[c("p1","p2","p3","p4","G")] %<>% dec(11)

tennis_ghost

standard evaluates:

(tennis_maxp <- maxp(tennis))
(tennis_ghost_maxp <- maxp(tennis_ghost))

First test the null of equal strength:

equalp.test(tennis)

not significant. Is the ghost real?

specificp.test(tennis_ghost,"G",0.0001)  # nonzero value needed to avoid triggering known R bug

so we can infer that the ghost is indeed real as we reject $H_0\colon p_G=0$.

We can use the Hessian matrix:

(M <- hessian(tennis,indep(tennis_maxp),border=TRUE))

See how the fillup value p4 takes up the unit sum constraint [it is defined as $p_4=1-p_1-p_2-p_3$]. Matrix M should be negative-semidefinite:

eigen(M,TRUE,TRUE)$values # one value is zero to numerical precision; the rest should be <0

See how one eigenvalue is zero (to numerical precision) and the rest are strictly negative, corresponding to the support function having a well-defined maximum.

Analysis with hyper3 {-}

The likelihood for the 9-2 scoreline between $\left\lbrace p_1,p_2\right\rbrace$ and $\left\lbrace p_3,p_4\right\rbrace$ would have a likelihood function of

$$ \frac{ \left(\lambda p_1 +\lambda p_2\right)^{9} \left(p_2 + p_3\right)^{2} }{ \left(\lambda p_1 +\lambda p_2+p_2 + p_3\right)^{11} } $$

H3 <- as.hyper3(H)
H3

## 1&2 vs 3&4;  scoreline 9-2 (again), with lambda=1.2

f <- function(lambda){
  H3[c(p1=lambda,p2=lambda)] %<>% inc(9)
  H3[c(p3=1,p2=1)] %<>% dec(2)
  H3[c(p1=lambda,p2=lambda,p3=1,p2=1)] %<>% dec(11)
  return(H3)
}
f(1.1443345)
f(1.888888)

maximize:

tennis_H3_maxp <- maxp(H3)
tennis_H3_maxp

Compare the hyper2 evaluates:

tennis_maxp
tennis_ghost_maxp

Package dataset {-}

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

save(tennis,tennis_maxp,tennis_ghost, tennis_ghost_maxp,file="tennis.rda")

References {-}

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