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. Here I analyse a dataset on NZ university rankings, taken from the Times Higher Educational Supplement. First load the data:

uni_table <- read.table("universities.txt",header=TRUE)
universities_table <- uni_table
head(universities_table,12)

Each row corresponds to a ranking in the THE ranking system. Thus the first row refers to The University of Auckland 2012, and we see that they scored 28.3 for their teaching, 25.0 for their research, and so on; rows 3 and 4 show that AUT and Lincoln were not given rankings that year. Rows 1-8 give results for 2012, rows 9-16 give 2013, and so on to 2020.

We can coerce each year to a hyper2 likelihood function by converting to an order table with order(), thence to a support function:

f <- function(a,y){
  a <- subset(a,a$year==y)
  o <- function(v){  # 'v' is a character vector
    good <- v != "DNP"
    v[good] <- rank(-as.numeric(v[good]),ties.method="last")
    return(v)
  }
  out <- data.frame(o(a[,3]),o(a[,4]),o(a[,5]),o(a[,6]),o(a[,7]))
  colnames(out) <- colnames(a)[3:7]
  rownames(out) <- a[,1]
  return(out)
}
f(universities_table,2012)

Now:

universities <- hyper2(pnames=universities_table[1:8,1])
for(i in 2013:2020){universities %<>% `+`(ordertable2supp(f(universities_table,i)))}
universities

Object universities is simply the support function for the universities' strengths for all years' observations (assumed to be independent) 2013-2020:

universities_maxp <- maxp(universities)
universities_maxp
pie(universities_maxp)

Note the dominance of Auckland (UoA). We observe that AUT has increased its strength since 2016 and we can test whether this increase is significant:

H1 <- hyper2(pnames=c("AUT_pre",universities_table[c(1:2,4:8),1]))
for(i in 2012:2016){
  b <- f(universities_table,i)
  rownames(b)[3] <- "AUT_pre"
  H1 <- H1 + ordertable2supp(b)
} # column loop closes


H2 <- hyper2(pnames=c("AUT_post",universities_table[c(1:2,4:8),1]))
for(i in 2017:2020){
  jj <- f(universities_table,i)
  rownames(jj)[3] <- "AUT_post"
  H2 <- H2 + ordertable2supp(jj)
}
maxp(H2)

H <- H1 + H2
maxp(H)

So from the above we see that the estimated post-2016 strength of AUT is larger than the pre-2016 strength. Is this change significant?

samep.test(H,c("AUT_pre","AUT_post"))

Package dataset {-}

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

save(universities_table,universities_maxp,universities,file="universities.rda")

RobinHankin/hyper2 documentation built on April 21, 2024, 11:38 a.m.