R/adj_intLik_p.R

In BriceonWiley/IntegratedLikelihood.R:

adj_intLik_p <- function(
  n00 = 49, n01 = 1, n10 = 1, n11 = 2, X = 14, Y = 433, k = 0.5, alpha = 0.05,
  closed_form = FALSE
) {
  # browser()


  # data preparation
  if (k <= 0 | k >= 1) stop("k must be bigger than 0 and smaller than 1.")
  data <- c(n00, n01, n10, n11, X, Y)
  n <- sum(data)
  crit <- n^k * (1 - exp(-qchisq(1 - alpha, 1) / n^k))
  pretty_k <- str_pad(as.character(k), 4, "right", pad = "0")
  CI_names <- c(paste0("Lower_k=", pretty_k), paste0("Upper_k=", pretty_k))


  # check if CI is possible
  if (any(is.na(mWald_p(n00, n01, n10, n11, X, Y, alpha)))) {
    interval <- tibble(lower = NaN, upper = NaN)
    names(interval) <- CI_names
    return(interval)
  }


  # find MLE
  if (closed_form) {
    phat <- optimize(
      likelihood_p_closed, c(0, 1), data = data, maximum = TRUE
    )$maximum
    num <- likelihood_p_closed(phat, data)
  } else {
    phat <- optimize(
      likelihood_p, c(0, 1), data = data, maximum = TRUE
    )$maximum
    num <- likelihood_p(phat, data)
  }


  # build CI
  if (phat != 0) {
    lower <- tryCatch({
      uniroot(
        likelihood_ratio_p,
        lower = 0,
        upper = phat,
        num = num,
        data = data,
        crit = crit,
        closed_form = closed_form,
        maxiter = 5e4
      )$root
    },
    error = function(e) {
      message(
        paste0(
          "There was a problem with the lower bound for the sample:\n",
          "n00 = ", n00, ", n01 = ", n01, ", n10 = ", n10, ", n11 = ", n11,
          ", X = ", X, " and Y = ", Y, "."
        )
      )
      0
    }
    )
  } else {
    lower <- 0
  }
  if (phat != 1) {
    upper <- tryCatch({
      uniroot(
        likelihood_ratio_p,
        lower = phat,
        upper = 1,
        num = num,
        data = data,
        crit = crit,
        closed_form = closed_form,
        maxiter = 5e4
      )$root
    },
    error = function(e) {
      message(
        paste0(
          "There was a problem with the upper bound for the sample:\n",
          "n00 = ", n00, ", n01 = ", n01, ", n10 = ", n10, ", n11 = ", n11,
          ", X = ", X, " and Y = ", Y, "."
        )
      )
      1
    }
    )
  } else {
    upper <- 1
  }


  # return CI
  interval <- tibble(lower, upper)
  names(interval) <- CI_names
  interval
}