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"))

(takes about an hour to run without cache)

This analysis is essentially without value, see the very end for why

To cite the hyper2 package in publications, please use @hankin2017_rmd. The global liveability ranking shows the nicest places to live.

(glr <- readLines("global_liveability_ranking.txt"))

First we will turn it into a hyper2 object:

(o <- glr  %>% strsplit(" ") %>% lapply(function(l){l[-1]}))

Will remove Helsinki, Brisbane, and Hamburg on account of their having very low rank scores (if we leave them in, the Hessian becomes non-positive-definite at the evaluate).

o %<>%  lapply(function(x){x[!(x %in% c("Hamburg","Helsinki","Brisbane"))]})
o

H <- hyper2()
for(i in seq_along(o)){
  H <- H  + race(o[[i]])
}

And do some tests:

(mH <- maxp(H))
pie(mH)
equalp.test(H)

jjm <- maxp(H,n=1,give=TRUE,hessian=TRUE)

jjm$hessian
diag(jjm$hessian)
eigen(jjm$hessian)$values

So the Hessian is positive-definite [all the eigenvalues are positive]. Also, the evaluate is robust, as indicated by consistency():

consistency(H)

Now we tether some of the cities together:

elastic <- list(
Austria     = c("Vienna"),
Denmark     = c("Copenhagen"),
Switzerland = c("Zurich","Geneva"),
Canada      = c("Calgary","Vancouver","Toronto"),
Germany     = c("Frankfurt","Hamburg"),
Netherlands = c("Amsterdam"),
Japan       = c("Osaka","Tokyo"),
Australia   = c("Melbourne","Adelaide","Perth","Brisbane","Sydney"),
NZ          = c("Auckland","Wellington"),
Finland     = c("Helsinki")
)

Now create a hyper3 object:

ML <- function(h){
   H3 <- hyper3()
   for(r in o){ H3 <- H3 + cheering3(r,e=list2nv(elastic),help=rep(h,10))}
   return(H3)
}
f <- function(h){maxp(ML(h),give=TRUE,n=1)$likes}

And plot a profile likelihood

l <- seq(from=1,to=5,len=10)
L <- sapply(l,f)

and plot it:

plot(log(l),L-max(L),type='b')

A coarser division

We will try dividing the countries into continents:

elastic_continent <- list(
Europe      = c("Vienna","Copenhagen","Zurich","Geneva","Frankfurt","Hamburg","Amsterdam","Helsinki"),
America     = c("Calgary","Vancouver","Toronto"),
Asia        = c("Osaka","Tokyo"),
Australasia = c("Melbourne","Adelaide","Perth","Brisbane","Sydney","Auckland","Wellington")
)

ec <- NULL
for(i in seq_along(elastic_continent)){
  jj <- rep(i,length(elastic_continent[[i]]))
  names(jj) <- elastic_continent[[i]]
  ec <- c(ec,jj)
}
ec 
MLc <- function(h){
   H3 <- hyper3()
   for(r in o){ H3 <- H3 + cheering3(r,e=ec,help=rep(h,4))}
   return(H3)
}
fc <- function(h){maxp(ML(h),give=TRUE,n=1)$likes}

l <- seq(from=0.8,to=3,len=10)
L <- sapply(l,fc)

plot(l,L-max(L),type="b")
abline(v=1)

Country-level analysis

We need to replace each city with its country:

elastic
(jj <- list2nv(elastic,FALSE))
o_countries <- lapply(o, function(x){jj[names(jj) %in% x]})
(o_countries <- lapply(o_countries,as.vector))

We see that NZ and Australia are dominated so we should remove them:

(o_countries <- lapply(o_countries,function(x){x[!(x %in% c("NZ","Australia","Japan"))]}))

But hang on, if we continue to remove dominated countries we eventually remove every country that is not Austria. So that was a big fat waste of time.