R/simulate_phi.R

In BriceonWiley/IntegratedLikelihood.R:

simulate_phi <- function(
  N0, N, theta1, theta2 = theta1, lam1 = 20, lam2 = 15, len = 10000,
  alpha = 0.05, adjustment = 0, closed, summary = FALSE
) {
  # browser()

  # messaging
  # print_n0 <- str_pad(as.character(N0), 2, "left", "0")
  # print_n <- str_pad(as.character(N), 2, "left", "0")
  print_theta <- str_pad(as.character(theta1), 4, "right", "0")

  message(
    paste0("N0 is ", N0, ", N is ", N, " and theta is ", print_theta, ".")
  )

  if (adjustment < 0 | adjustment > 1)
    stop("The adjustment must be between 0 and 1.")


  # producing and checking samples
  samples <- rerun(len, {
    m01 <- rpois(1, N0 * lam1)
    m02 <- rpois(1, N0 * lam2)
    y01 <- rbinom(1, m01, theta1)
    y02 <- rbinom(1, m02, theta2)
    z1 <- rpois(1, N * (lam1 * (1 - theta1) + lam2 * theta2))
    z2 <- rpois(1, N * (lam2 * (1 - theta2) + lam1 * theta1))
    if (m01 == 0) m01 <- adjustment
    if (m02 == 0) m02 <- adjustment
    if (y01 == 0) y01 <- adjustment
    if (y02 == 0) y02 <- adjustment
    tibble(z1 = z1, z2 = z2, m01 = m01, m02 = m02, y01 = y01, y02 = y02)
  }) %>%
    bind_rows()


  # producing intervals
  walds <- pmap_dfr(
    samples, ~ Wald_phi(..1, ..2, ..3, ..4, ..5, ..6, N0, N, alpha)
  )
  scores <- pmap_dfr(
    samples, ~ Score_phi(..1, ..2, ..3, ..4, ..5, ..6, N0, N, alpha)
  )
  ILR <- pmap_dfr(
    samples, ~ intLik_phi(
      ..1, ..2, ..3, ..4, ..5, ..6, N0, N, alpha, closed
    )
  )
  intervals <- bind_cols(walds, scores, ILR)
  NaNs <- intervals %>%
    pmap_lgl(~ any(is.nan(c(..1, ..2, ..3, ..4, ..5, ..6)))) %>%
    sum()
  values <- c(
    "lambda1" = lam1, "lambda2" = lam2, "phi" = round(lam1 / lam2, 4),
    "theta1" = theta1, "theta2" = theta2, "N0" = N0, "N" = N,
    "B" = len, "NaN" = NaNs
  )
  results <- list(
    "values" = values,
    "samples" = samples,
    "intervals" = intervals
  )

  # producing final product
  if (summary) {
    coverages <- coverage_phi(results, ignore_NaN = TRUE)
    widths <- width_phi(results, ignore_NaN = TRUE)
    summarized <-
      bind_rows(coverages, widths) %>%
      select(-N0, -N, -phi, -theta1, -theta2) %>%
      t() %>%
      `colnames<-`(c("Coverage Probability", "Average Width"))
    message(
      paste0(
        "Removed ", results$counts['NaN'],
        " bad intervals in summary calculations."
      )
    )
    list(
      "values" = results$values,
      "summary" = round(summarized, 4)
    )
  } else results
}