R/intLik_p.R

In BriceonWiley/IntegratedLikelihood.R:

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


  # data preparation
  data <- c(n00, n01, n10, n11, X, Y)
  crit <- qchisq(1 - alpha, 1)


  # check if CI is possible
  if (any(is.na(mWald_p(n00, n01, n10, n11, X, Y, alpha)))) {
    return(tibble("ILR_Lower" = NaN, "ILR_Upper" = NaN))
  }


  # 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
  tibble("ILR_Lower" = lower, "ILR_Upper" = upper)
}