R/conf_int_hybrid_et.R

In adviksh/winfer:

#' Equal-Tailed Hybrid Confidence Interval
#' @export
#' @inheritParams ci_docs
#'
ci_hybrid_unbiased <- function(estimates, st_errs, conf_level = 0.95,
                               beta = (1 - conf_level) / 10) {

  # Defense -----------------------------------------------------------------
  check_estimates(estimates = estimates)
  check_st_errs(st_errs = st_errs, estimates = estimates)
  check_conf_level(conf_level = conf_level)
  check_beta(beta = beta, conf_level = conf_level)


  # Fast Returns ------------------------------------------------------------
  if (beta == 0) {
    cis <- ci_conditional_unbiased(estimates  = estimates,
                                   st_errs    = st_errs,
                                   conf_level = conf_level)
    return(cis)
  }


  # Function Body -----------------------------------------------------------
  # see discussion before prop 6 in Andrews, Kitagawa, McCloskey (2019)
  adj_conf_level <- conf_level / (1 - beta)

  half_alpha <- (1 - adj_conf_level) / 2
  tail_tiles <- c(half_alpha, 1 - half_alpha)

  ordering <- order(estimates, decreasing = TRUE)

  ci_hybrid_unbiased_single(estimates_sorted = estimates[ordering],
                            st_errs_sorted   = st_errs[ordering],
                            tail_tiles       = tail_tiles,
                            beta             = beta)

}

ci_hybrid_unbiased_single <- function(estimates_sorted, st_errs_sorted,
                                      tail_tiles, beta) {

  if (is.unsorted(rev(estimates_sorted))) {
    rlang::abort("estimates_unsorted is not sorted in decreasing order")
  }

  tp <- conditional_truncation_points(estimates_sorted)

  purrr::map_dbl(tail_tiles, mu_alpha_hybrid,
                 beta  = beta,
                 n_est = length(estimates_sorted),
                 y     = estimates_sorted[1],
                 se    = st_errs_sorted[1],
                 trunc_lo = min(tp),
                 trunc_hi = max(tp))
}

mu_alpha_hybrid <- function(alpha, beta, n_est, y, se,
                            trunc_lo = -Inf, trunc_hi = Inf) {

  c_beta <-  q_max_absnorm(p    = 1 - beta,
                           size = n_est)

  ur_result <- stats::uniroot(f         = tcdf_minus_left_tail_hybrid,
                              interval  = c(y +  3 * se, y - 3 * se),
                              y         = y,
                              se        = se,
                              alpha     = alpha,
                              c_beta    = c_beta,
                              trunc_lo  = trunc_lo,
                              trunc_hi  = trunc_hi,
                              extendInt = "yes")
  ur_result$root
}

tcdf_minus_left_tail_hybrid <- function(mu, se, y, alpha, c_beta,
                                        trunc_lo = -Inf,
                                        trunc_hi = Inf) {

  conf_int_null_uncond <- c(mu - c_beta * se,
                            mu + c_beta * se)

  tp_hybrid <- intersect_intervals(c(trunc_lo, trunc_hi),
                                   conf_int_null_uncond)

  if (y <= min(tp_hybrid)) {

    p <- 0

  } else {

    p <- TruncatedNormal::ptmvnorm(q     = y,
                                   mu    = mu,
                                   sigma = se,
                                   lb    = min(tp_hybrid),
                                   ub    = max(tp_hybrid),
                                   log   = FALSE,
                                   type  = "mc",
                                   B     = 1e5)

  }

  if (is.na(p)) browser()

  p - (1 - alpha)
}