set.seed(0) knitr::opts_chunk$set(echo = TRUE) library("hyper2") library("magrittr") options("digits" = 5)
This analysis is essentially without value, see the very end for why
(takes about an hour to process without cache)
![](`r system.file("help/figures/hyper2.png", package = "hyper2")`){width=10%}
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')
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)
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.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.