R/intLik_phi.R

In BriceonWiley/IntegratedLikelihood.R:

intLik_phi <- function(
  z1 = 40, z2 = 77, m01 = 15, m02 = 12, y01 = 6, y02 = 1,
  N0 = 0.7, N = 27.89, alpha = 0.05, closed_form = TRUE
) {
  # browser()


  # data preparation and find MLE
  data <- c(z1, z2, m01, m02, y01, y02, N0, N)
  crit <- qchisq(1 - alpha, 1)
  init <- calculate_mles_pois(data)["phihat"]
  if(closed_form) opt <-
      optimize(likelihood_phi_closed, c(1e-4, init + 3), data = data, maximum = T)
  else opt <- optimize(
    likelihood_phi, c(1e-4, init + 3), data = data, init = init, maximum = T
  )
  phihat <- opt$maximum; H_1 <- opt$objective


  # build CI
  lower <- tryCatch(
    uniroot(
      likelihood_ratio_phi, c(1e-6, phihat),
      H_1 = H_1, data = data, crit = crit, closed = closed_form
    )$root,
    error = function(e) NaN
  )
  upper <- tryCatch(
    uniroot(
      likelihood_ratio_phi, c(phihat, 7),
      H_1 = H_1,
      data = data, crit = crit, closed = closed_form
    )$root,
    error = function(e) {
      tryCatch(
        uniroot(
          likelihood_ratio_phi, c(phihat, 100),
          H_1 = H_1,
          data = data, crit = crit, closed = closed_form
        )$root,
        error = function(e) {
          tryCatch(
            uniroot(
              likelihood_ratio_phi, c(phihat, 1000),
              H_1 = H_1,
              data = data, crit = crit, closed = closed_form
            )$root,
            error = function(e) NaN
          )
        }
      )
    }
  )


  # return CI
  tibble(ILR_Lower = lower, ILR_Upper = upper)
}