(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
hyper2 Package: Likelihood Functions for Generalized Bradley-Terry Models.” The R Journal 9 (2): 429–39.