Simulating from Bradley-Terry models and comparing variance estimators

robust_QS_test.R

simulate <- function(means, ..., sampler = rpois) {
  # Simulate from a dataset and return
  # a sample that is in the same format
  sample <- means
  sample[] <- sampler(length(means), means, ...)
  return(sample)
}

make_QS <- function(n, mean = 10) {
  D <- diag(sample.int(n))
  S <- matrix(0, n, n)
  S[lower.tri(S)] <- rpois(n * (n - 1) / 2, mean)
  S <- S + t(S)
  D %*% S
}

simulate_BT <- function(pi, N) {
  n <- length(pi)
  P <- matrix(pi, n, n)
  P <- P / (P + t(P))
  out <- N
  counts <- rbinom(n * (n-1) / 2, N[lower.tri(N)], P[lower.tri(P)])
  complement <- N[lower.tri(N)] - counts
  out[lower.tri(out)] <- complement
  out <- t(out)
  out[lower.tri(out)] <- counts
  return(out)
}

BT4 <- BradleyTerry2::citations
pi4 <- scrooge::BTscores(BT4)
N4 <- round((BT4 + t(BT4)) / 2)
simulate_BT(pi4, N4)

pi10 <- runif(10)
pi10 <- pi10 / sum(pi10)
N10 <- simulate(matrix(10, 10, 10))
N10 <- round((N10 + t(N10)) / 2)
simulate_BT(pi10, N10)

n <- 10
A <- (diag(n) - matrix(1/n, n, n))[, -1]
sum2zero <- function(x) A %*% x %*% t(A)

local({
  #set.seed(5)
  op <- par(mfrow = c(4, 4))
  for (scale in c(1, (2:16) * 9)) {
    cite <- simulate_BT(pi10, N10 * scale)
    BT <- BTmodel$new(cite)
    MLE <- vcov(BT$model)
    Robust <- sandwich::sandwich(BT$model)
    BiasReduced <- BT$bias_reduced()
    x <- sum2zero(MLE)
    y <- sum2zero(Robust)
    y <- sum2zero(BiasReduced)
    x <- x[lower.tri(x, T)]
    y <- y[lower.tri(y, T)]
    plot(x, y, asp = 1,
         main = paste('Covariances x', scale),
         xlab = 'MLE',
         ylab = 'Bias Reduced')
    abline(0, 1, col = 'grey', lty = 2)
    abline(lm(y ~ x), col = 'red', lty = 3)
    abline(v = 0, col = 'lightgrey')
  }
  par(op)
})