(takes a “long time” to run without cache; the time is spent in chunk calculateprofilelike)

To cite the hyper2 package in publications, please use Hankin (2017). The file eurovision2022_jury_pre_null.txt, used below, refers to the jury vote for the 2022 Eurovision Song contest, pre-nullification; data supplied by Riley Uttley. First we specify the matrix as appearing in the Wikipedia page:

eurovision_table <- as.matrix(read.table("eurovision2022_jury_pre_null.txt"))
excluded <- c(2, 8, 11, 14,15,16)

contestants <- c(
"AU", # Australia        1
"AZ", # Azerbaijan       2
"BE", # Belgijm          3
"CY", # Cyprus           4
"CZ", # Czech republic   5
"EE", # Estonia          6
"FI", # Finland          7
"GE", # Georgia          8 
"IE", # Ireland          9 
"IL", # Israel          10 
"ME", # Montenegro      11
"MK", # North Macedonia 12
"MT", # Malta           13 
"PL", # Poland          14 
"RO", # Romania         15
"SM", # San Marino      16
"RS", # Serbia          17
"SE"  # Sweden          18
)
removed <- contestants[excluded]
non_competitors <- c("DE","ES","UK") # Germany, Spain, UK
juries <- c(contestants, non_competitors)
dimnames(eurovision_table) <- list(contestants=contestants, juries=juries)
eurovision_table
##            juries
## contestants AU AZ BE CY CZ EE FI GE IE IL ME MK MT PL RO SM RS SE DE ES UK
##          AU NA  4  8 10  8 10 10  1  5 10 10  0  5  1  5  4  2 12  6 10  3
##          AZ  3 NA  3  6  4  6  7 12  0  0  4  8  1 10  7  6  8  5 10 12  8
##          BE  8  0 NA  3  6  5  0  0  2  7  0  0 10  0  0  0  6  8  5  8  2
##          CY  0  0  0 NA  0  0  3  0  4  0  0  1  0  0  0  0  0  2  0  0  0
##          CZ  2  0  4  4 NA  7  6  0 10  6  7  0  0  0  0  0  4  6  0  0 10
##          EE  5  0  0  7  7 NA  4  3  8  3  8  0  2  3  4  1 10 10  4  0  4
##          FI  0  5  5  2  2  8 NA  2  0  5  3  2  4  5  0  0  5  4  1  2  0
##          GE  1 10  0  0  5  3  0 NA  0  0  1 12  0  4  6 10  0  0  3  0  0
##          IE  6  0  1  1  0  4  0  4 NA  0  0  0  0  0  1  0  0  0  0  0  0
##          IL 10  0  7  0  0  0  0  5  3 NA  2  0  0  0  0  2  0  3  2  3  1
##          ME  0  0  0  0  1  1  0  0  6  2 NA  0  6  0  0  0  0  7  0  0  0
##          MK  0  7  0  0  0  0  0 10  0  0  0 NA  3  7  8  7  7  0  0  1  0
##          MT  7  1  0  0 10  0  2  0  1  0  5  3 NA  0  3  0  1  0 12  0  7
##          PL  0 12  6  8  3  0  1  6  0  8  6  5  8 NA 10  8  3  1  0  7  6
##          RO  0  8  0  0  0  2  0  8  0  4  0  7  0  8 NA 12  0  0  8  4  0
##          SM  0  6 12  0  0  0  5  7  0  0  0  6  0 12 12 NA  0  0  0  0  0
##          RS  4  2  2  5  0  0  8  0  7  1  0 10 12  2  2  5 NA  0  0  5  5
##          SE 12  3 10 12 12 12 12  0 12 12 12  4  7  6  0  3 12 NA  7  6 12
(points <- sort(unique(c(eurovision_table)),decreasing=TRUE))
##  [1] 12 10  8  7  6  5  4  3  2  1  0
preference <- eurovision_table*0  
for(i in seq_along(points)){ preference[eurovision_table == points[i]] <- i }
preference
##            juries
## contestants AU AZ BE CY CZ EE FI GE IE IL ME MK MT PL RO SM RS SE DE ES UK
##          AU NA  7  3  2  3  2  2 10  6  2  2 11  6 10  6  7  9  1  5  2  8
##          AZ  8 NA  8  5  7  5  4  1 11 11  7  3 10  2  4  5  3  6  2  1  3
##          BE  3 11 NA  8  5  6 11 11  9  4 11 11  2 11 11 11  5  3  6  3  9
##          CY 11 11 11 NA 11 11  8 11  7 11 11 10 11 11 11 11 11  9 11 11 11
##          CZ  9 11  7  7 NA  4  5 11  2  5  4 11 11 11 11 11  7  5 11 11  2
##          EE  6 11 11  4  4 NA  7  8  3  8  3 11  9  8  7 10  2  2  7 11  7
##          FI 11  6  6  9  9  3 NA  9 11  6  8  9  7  6 11 11  6  7 10  9 11
##          GE 10  2 11 11  6  8 11 NA 11 11 10  1 11  7  5  2 11 11  8 11 11
##          IE  5 11 10 10 11  7 11  7 NA 11 11 11 11 11 10 11 11 11 11 11 11
##          IL  2 11  4 11 11 11 11  6  8 NA  9 11 11 11 11  9 11  8  9  8 10
##          ME 11 11 11 11 10 10 11 11  5  9 NA 11  5 11 11 11 11  4 11 11 11
##          MK 11  4 11 11 11 11 11  2 11 11 11 NA  8  4  3  4  4 11 11 10 11
##          MT  4 10 11 11  2 11  9 11 10 11  6  8 NA 11  8 11 10 11  1 11  4
##          PL 11  1  5  3  8 11 10  5 11  3  5  6  3 NA  2  3  8 10 11  4  5
##          RO 11  3 11 11 11  9 11  3 11  7 11  4 11  3 NA  1 11 11  3  7 11
##          SM 11  5  1 11 11 11  6  4 11 11 11  5 11  1  1 NA 11 11 11 11 11
##          RS  7  9  9  6 11 11  3 11  4 10 11  2  1  9  9  6 NA 11 11  6  6
##          SE  1  8  2  1  1  1  1 11  1  1  1  7  4  5 11  8  1 NA  4  5  1
H_excluded <- hyper2()
H_included <- hyper2()
text <- "remove this for R >= 4.6-0"
print(text)
## [1] "remove this for R >= 4.6-0"
`%notin%` <- Negate(`%in%`)
for(i in excluded){
    d <- preference[,i,drop=TRUE]
    d[d==11] <- 0
    newH <- suppfun(d[!is.na(d)])
    H_excluded <- H_excluded + newH
} # i loop closes

for(i in seq_along(ncol(preference))){
    if(i %notin% excluded){
        d <- preference[,i,drop=TRUE]
        d[d==11] <- 0
        newH <- suppfun(d[!is.na(d)])
        H_excluded <- H_excluded + newH
    }
} # i loop closes
H_excluded <- H_excluded |>
    pwa("AZ", "a_az") |> pwa("GE", "a_GE") |>
    pwa("ME", "a_ME") |> pwa("PL", "a_PL") |>
    pwa("RO", "a_RO") |> pwa("SM", "a_SM")
H <- H_excluded + H_included
pnames(H_excluded)
##  [1] "a_az" "a_GE" "a_ME" "a_PL" "a_RO" "a_SM" "AU"   "AZ"   "BE"   "CY"  
## [11] "CZ"   "EE"   "FI"   "GE"   "IE"   "IL"   "ME"   "MK"   "MT"   "PL"  
## [21] "RO"   "RS"   "SE"   "SM"
S_likelihood <- function(S, give=FALSE){

    make_strengths <- function(x){  # x is independent
        stopifnot(length(x) == 17)
        jj <- c(rep(S,6), x)
        setNames(c(jj, 1-sum(jj)),  # sum = 1
                 pnames(H_excluded))
        
    }
    
    objective <- function(x){ return(loglik(make_strengths(x), H)) }
    
    x_start <- rep(0.01, 17)
    
    n <- 17
    SMALL <- 1e-5
    UI <- rbind(diag(nrow = n), -1)
    CI <- c(rep(SMALL, n), -1 + SMALL)
    CI[18] <- 6*S - 1
    
    jj <- constrOptim(theta = x_start, f = objective, grad = NULL, ui = UI, ci = CI)

    if(give){
        return(jj)
    } else {
        return(jj$value)
    }
}    
jj6 <- S_likelihood(0.06)
jj7 <- S_likelihood(0.07)
jj6
## [1] -467.24
jj7
## [1] -447.6
pchisq(2*(jj7-jj6),df=1, lower.tail=FALSE)
## [1] 3.672e-10
Hankin, R. K. S. 2017. “Partial Rank Data with the hyper2 Package: Likelihood Functions for Generalized Bradley-Terry Models.” The R Journal 9 (2): 429–39.