R/simulate_gamma.R

In BriceonWiley/IntegratedLikelihood.R:

simulate_gamma <- function(
  N1, p1 = 0.1, phi1 = 0.1, theta1 = 0.1, N2 = N1, p2 = 0.1, phi2 = 0.1,
  theta2 = 0.1, alpha = 0.05, f = 0.1, len = 1e4, check = NULL, summary = FALSE
) {
  # browser()

  print_n <- str_pad(as.character(N1), 3, "left", "0")
  print_p <- str_pad(as.character(p1), 4, "right", "0")

  message(paste0("N1 is ", N1, ", and p1 is ", print_p, "."))


  # producing and checking samples
  sample_1 <- produce_samples(N1, p1, phi1, theta1, len, f)
  sample_2 <- produce_samples(N2, p2, phi2, theta2, len, f)
  if(!is.null(check)) {
    sample_1 <- check_zeros(sample_1, check)
    sample_2 <- check_zeros(sample_2, check)
  }
  samples <- suppressMessages(bind_cols(sample_1, sample_2)) %>%
    `colnames<-`(c(
      "n100", "n101", "n110", "n111", "X1", "Y1",
      "n200", "n201", "n210", "n211", "X2", "Y2"
      ))

  # producing intervals
  Walds <-
    pmap_dfr(
      samples,
      ~ Wald_gam(
        c(..1, ..2, ..3, ..4, ..5, ..6),
        c(..7, ..8, ..9, ..10, ..11, ..12),
        alpha
      )
    )
  intLiks <-
    pmap_dfr(
      samples,
      ~ intLik_gam(
        c(..1, ..2, ..3, ..4, ..5, ..6),
        c(..7, ..8, ..9, ..10, ..11, ..12),
        alpha
        )
      )
  counts <-
    c(
      "N1" = N1, "M1" = N1 * 0.9, "n1" = N1 * 0.1,
      "N2" = N2, "M2" = N2 * 0.9, "n2" = N2 * 0.1,
      "B" = len, "NaN" = sum(is.nan(Walds$W_Lower))
      )
  proportions <- c(
    "p1" = p1, "phi1" = phi1, "theta1" = theta1,
    "p2" = p2, "phi2" = phi2, "theta2" = theta2,
    "log odds" = log((p1 * (1 - p2)) / (p2 * (1 - p1))),
    "odds" = (p1 * (1 - p2)) / (p2 * (1 - p1))
    )
  intervals <- bind_cols(Walds, intLiks)

  results <- list(
    "counts" = counts,
    "proportions" = proportions,
    "sample 1" = sample_1,
    "sample 2" = sample_2,
    "intervals" = intervals
  )

  # producing final product
  if (summary) {
    coverages <- coverage_gam(results, ignore_NaN = TRUE)
    widths <- width_gam(results, ignore.NaN = TRUE)
    summarized <-
      bind_rows(coverages, widths) %>%
      select(W, ILR) %>%
      t() %>%
      `colnames<-`(c("Coverage Probability", "Average Width"))
    message(
      paste0(
        "Removed ", results$counts['NaN'],
        " bad intervals in summary calculations."
      )
    )
    list(
      "counts" = results$counts,
      "proportions" = results$proportions,
      "summary" = round(summarized, 4)
    )
  } else results
}