R/CI_calc.R

In koohyun-kwon/OptACI: Methods to Construct Optimal Average Coverage Confidence Intervals

#' Optimal Average Coverage Confidence Interval
#'
#' Calculates the estimated optimal average coverage CIs, and
#' tuning parameters associated with those CIs.
#'
#' @param x_vec a vector of observed values.
#' @param sig_vec a vector of individual standard deviations.
#' @param alpha a desired level of non-coverage.
#' @param Jn the order of series approximation;
#' do not specify its value when using the simple estimator.
#' @param lam_init the candidate vector of minimum shrinkage factors to be used;
#' do not specify its value when using the simple estimator.
#' @inheritParams lam_adj
#'
#' @return a list with four components; \code{lam_opt} gives the value of the
#' optimal shrinkage factor, \code{chi_opt} gives the value of the normalized half-length
#' associated with the optimal shrinkage factor, and \code{CI_l} and \code{CI_u} gives the vector of
#' lower and upper ends of the confidence intervals, respectively.
#' @export
#'
#' @examples th_vec <- rnorm(50) + 1
#' sig_vec <- rnorm(50)^2 + 1
#' x_vec <- rnorm(50, th_vec, sig_vec)
#' lam_init <- c(0.3, 0.35, 0.4, 0.45)
#' ACI(x_vec, sig_vec, 0.05)
#' ACI(x_vec, sig_vec, 0.05, 2, lam_init, 0.025)
ACI <- function(x_vec, sig_vec, alpha, Jn, lam_init, tol, chilen = 50){

  x_vec_norm <- (x_vec - mean(x_vec))/sig_vec

  if(missing(Jn)){

    optres <- stats::optimize(f = chilam, interval = c(1/2, 1),
                              xvec = x_vec_norm, alpha = alpha)

  }else{

    lammin_new <- lam_adj(lam_init, x_vec_norm, Jn, alpha, chilen, tol)
    optres <- stats::optimize(f = chilam, interval = c(lammin_new, 1),
                              xvec = x_vec_norm, alpha = alpha, Jn = Jn)
  }

  lam_opt <- optres$minimum
  chi_opt <- optres$objective
  CI_l <- lam_opt * x_vec + (1 - lam_opt) * mean(x_vec) - chi_opt * sig_vec
  CI_u <- lam_opt * x_vec + (1 - lam_opt) * mean(x_vec) + chi_opt * sig_vec

  res <- list(lam_opt = lam_opt, chi_opt = chi_opt, CI_l = CI_l, CI_u = CI_u)
  return(res)
}

#' Oracle Average Coverage Confidence Interval
#'
#' Calculates the true optimal average coverage CIs, and
#' tuning parameters associated with those CIs.
#'
#' This function was created for the purpose of simulations.
#'
#' @param th_vec the true mean vector
#' @param min_lam the minimum value of the shrinkage factor to be used in optimization;
#' default is 0.01.
#' @param x_vec a vector of observed values; it can be unspecified if only the shrinkage factor
#' and the half-length are calculated.
#' @inheritParams ACI
#'
#' @return a list with four components; \code{lam_opt} gives the value of the
#' optimal shrinkage factor, \code{chi_opt} gives the value of the normalized half-length
#' associated with the optimal shrinkage factor, and \code{CI_l} and \code{CI_u} gives the vector of
#' lower and upper ends of the confidence intervals, respectively.
#' @export
#'
#' @examples th_vec <- rnorm(50) + 1
#' sig_vec <- rnorm(50)^2 + 1
#' x_vec <- rnorm(50, th_vec, sig_vec)
#' ACI_orc(th_vec, sig_vec, 0.05)
#' ACI_orc(th_vec, sig_vec, 0.05, x_vec)
ACI_orc <- function(th_vec, sig_vec, alpha, x_vec, min_lam = 0.01){

  th_vec_norm <- (th_vec - mean(th_vec))/sig_vec
  optres <- stats::optimize(f = chilam_orc, interval = c(min_lam, 1),
                            th_vec = th_vec_norm, alpha = alpha)

  lam_opt <- optres$minimum
  chi_opt <- optres$objective
  if(missing(x_vec)){

    CI_l <- NULL
    CI_u <- NULL

  }else{

    CI_l <- lam_opt * x_vec + (1 - lam_opt) * mean(x_vec) - chi_opt * sig_vec
    CI_u <- lam_opt * x_vec + (1 - lam_opt) * mean(x_vec) + chi_opt * sig_vec
  }


  res <- list(lam_opt = lam_opt, chi_opt = chi_opt, CI_l = CI_l, CI_u = CI_u)
  return(res)

}

#' Optimal Average Coverage Confidence Interval under Normality Assumption
#'
#' Calculates the estimated optimal average coverage CIs, and
#' tuning parameters associated with those CIs, using the function \code{gnB_norm}.
#'
#' @inheritParams ACI
#' @inheritParams chilam_norm
#'
#' @return a list with four components; \code{lam_opt} gives the value of the
#' optimal shrinkage factor, \code{chi_opt} gives the value of the normalized half-length
#' associated with the optimal shrinkage factor, and \code{CI_l} and \code{CI_u} gives the vector of
#' lower and upper ends of the confidence intervals, respectively.
#' @export
#'
#' @examples th_vec <- rnorm(50) + 1
#' sig_vec <- rnorm(50)^2 + 1
#' x_vec <- rnorm(50, th_vec, sig_vec)
#' lam_init <- c(0.3, 0.35, 0.4, 0.45)
#' ACI_aux(x_vec, sig_vec, 0.05, 2, lam_init, 4)
ACI_aux <- function(x_vec, sig_vec, alpha, Jn, lam_init, Jn2){

  x_vec_norm <- (x_vec - mean(x_vec))/sig_vec

  lammin_new <- min(lam_init) # No adjustment of lambda


  optres <- stats::optimize(f = chilam_norm, interval = c(lammin_new, 1),
                            xvec = x_vec_norm, alpha = alpha, Jn = Jn, Jn2 = Jn2)

  lam_opt <- optres$minimum
  chi_opt <- optres$objective
  CI_l <- lam_opt * x_vec + (1 - lam_opt) * mean(x_vec) - chi_opt * sig_vec
  CI_u <- lam_opt * x_vec + (1 - lam_opt) * mean(x_vec) + chi_opt * sig_vec

  res <- list(lam_opt = lam_opt, chi_opt = chi_opt, CI_l = CI_l, CI_u = CI_u)
  return(res)
}


#' Optimal Average Coverage Confidence Interval Using AKP Procedure
#'
#' Calculates the estimated optimal average coverage CIs, and
#' tuning parameters associated with those CIs, using \code{ebci} package.
#'
#' @inheritParams ACI
#' @param min_mu2 the cut-off value.
#'
#' @return a list with four components; \code{lam_opt} gives the value of the
#' optimal shrinkage factor, \code{chi_opt} gives the value of the normalized half-length
#' associated with the optimal shrinkage factor, and \code{CI_l} and \code{CI_u} gives the vector of
#' lower and upper ends of the confidence intervals, respectively.
#' @export
#'
#' @examples th_vec <- rnorm(50) + 1
#' sig_vec <- rnorm(50)^2 + 1
#' x_vec <- rnorm(50, th_vec, sig_vec)
#' ACI_AKP(x_vec, sig_vec, 0.05, 0.04)
ACI_AKP <- function(x_vec, sig_vec, alpha, min_mu2){

  x_vec_norm <- (x_vec - mean(x_vec))/sig_vec

  eps = x_vec_norm - mean(x_vec_norm)
  mu2 = max(min_mu2, mean(eps^2 - 1))

  w_opt_res = ebci::w_opt(mu2, kappa = Inf, alpha)
  chi_opt = w_opt_res$length
  lam_opt = w_opt_res$w

  CI_l <- lam_opt * x_vec + (1 - lam_opt) * mean(x_vec) - chi_opt * sig_vec
  CI_u <- lam_opt * x_vec + (1 - lam_opt) * mean(x_vec) + chi_opt * sig_vec

  res <- list(lam_opt = lam_opt, chi_opt = chi_opt, CI_l = CI_l, CI_u = CI_u)
  return(res)

}


#' Hybrid ACI Procedure
#'
#' @inheritParams ACI_aux
#' @inheritParams ACI_AKP
#'
#' @return a list with four components; \code{lam_opt} gives the value of the
#' optimal shrinkage factor, \code{chi_opt} gives the value of the normalized half-length
#' associated with the optimal shrinkage factor, and \code{CI_l} and \code{CI_u} gives the vector of
#' lower and upper ends of the confidence intervals, respectively.
#' @export
#'
#' @examples th_vec <- rnorm(50) + 1
#' sig_vec <- rnorm(50)^2 + 1
#' x_vec <- rnorm(50, th_vec, sig_vec)
#' lam_init <- c(0.3, 0.35, 0.4, 0.45)
#' ACI_hybr(x_vec, sig_vec, 0.05, 2, lam_init, 4, 0.1)
ACI_hybr <- function(x_vec, sig_vec, alpha, Jn, lam_init, Jn2, min_mu2){

  norm_res <- ACI_aux(x_vec, sig_vec, alpha, Jn, lam_init, Jn2)
  AKP_res <- ACI_AKP(x_vec, sig_vec, alpha, min_mu2)

  if(norm_res$chi_opt <= AKP_res$chi_opt | AKP_res$lam_opt >= min(lam_init)){

    res <- norm_res

  }else{

    res<- AKP_res
  }

  return(res)
}