permChacko issue #13

issue-13.R

# Simulations to address https://github.com/ocbe-uio/permChacko/issues/13

# Setup =================================================================
remotes::install_github("ocbe-uio/permChacko@issue-13")
runs <- 1000L
results <- matrix(
  NA, nrow = runs, ncol = 6L,
  dimnames = list(
    "run" = seq_len(runs),
    "p_values" = c(
      "analytic", "numeric", "numeric_alt", "tabular", "bias_num", "bias_alt"
    )
  )
)

# Simulations ===========================================================

# Samples
parameters <- data.frame(
  "k" = sample(3:10, runs, replace = TRUE),
  "lambda" = sample(100:500, runs, replace = TRUE)
)
vecs <- list()
for (r in seq_len(runs)) {
  vecs[[r]] <- rpois(parameters[r, "k"], parameters[r, "lambda"])
}

# Paralellized calculations
library(doParallel)
registerDoParallel()
message("Starting ", runs, " simulations on ", getDoParWorkers(), " cores")
pvals <- foreach (r = seq_len(runs)) %dopar% {
  require(permChacko)
  permChacko(vecs[[r]])
}

# Calculating results
for (i in seq_len(runs)) {
  results[i, 1:4] <- pvals[[i]]$p_values
  if (!is.na(results[i, "analytic"])) {
    # If analytic p-value is present, use it as the true p-value
    results[i, "bias_num"] <- results[i, "numeric"] - results[i, "analytic"]
    results[i, "bias_alt"] <- results[i, "numeric_alt"] - results[i, "analytic"]
  } else {
    # If analytic p-value is missing, use tabular p-value as the true p-value
    results[i, "bias_num"] <- results[i, "numeric"] - results[i, "tabular"]
    results[i, "bias_alt"] <- results[i, "numeric_alt"] - results[i, "tabular"]
  }
}

# Presenting results ====================================================
message("General results")
print(summary(results))

message("Results for cases where the analytic p-value is missing")
print(summary(results[is.na(results[, "analytic"]), ]))

message("Results for cases where the analytic p-value is present")
print(summary(results[!is.na(results[, "analytic"]), ]))