R/Diagnostics.R

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

#' Average Coverage Probability Diagnostics
#'
#' Calculates the true average coverage probability of
#' calculated confidence intervals.
#'
#' This function was created for the purpose of simulations.
#'
#' @param th_vec the true mean vector
#' @param CI_l the lower end of confidence intervals, in the same length
#' of \code{th_vec}.
#' @param CI_u the upper end of confidence intervals, in the same length
#' of \code{th_vec}.
#'
#' @return the true average coverage probability
#' @export
#'
#' @examples th_vec <- rnorm(50) + 1
#' sig_vec <- rnorm(50)^2 + 1
#' x_vec <- rnorm(50, th_vec, sig_vec)
#' ACIres <- ACI(x_vec, sig_vec, 0.05)
#' CI_l <- ACIres$CI_l
#' CI_u <- ACIres$CI_u
#' ACP_check(th_vec, CI_l, CI_u)
ACP_check <- function(th_vec, CI_l, CI_u){

  cov_L <- as.numeric(th_vec >= CI_l)
  cov_U <- as.numeric(th_vec <= CI_u)
  cov_ind <- cov_L * cov_U
  cov_avg <- mean(cov_ind)
}


#' Mean Vector Generator
#'
#' Generates a mean vector from a specified design.
#'
#' This function was created for the purpose of simulations.
#'
#' @param n length of the vector
#' @param M dispersion parameter
#' @param method design to be used
#'
#' @return a vector of length \code{n}.
#' @export
#'
#' @examples th_gen(100, 1, "unif")
th_gen <- function(n, M, method){

  if(method == "unif"){

    th_vec <- seq(from = -M, to = M, length.out = n)

  }else if(method == "tnorm"){

    th_vec <- truncnorm::rtruncnorm(n = n, a = -M, b = M, mean = 0, sd = M/2)

  }else if(method == "norm"){

    th_vec <- stats::rnorm(n = n, mean = 0, sd = M/2)

  }else if(method == "t"){

    th_vec <- sgt::rsgt(n, mu = 0, sigma = M/2, lambda = 0, p = 2, q = 2)

  }else if(method == "sk-t"){

    th_vec <- sgt::rsgt(n, mu = 0, sigma = M/2, lambda = 0.9, p = 2, q = Inf)

  }

  return(th_vec)
}

#' Simulation Result Data Cleaning
#'
#' Generates cleaned data frames from simulation results.
#'
#' This function was created for the purpose of simulations.
#'
#' @param config_names names of configurations
#' @param M_range a vector of dispersion parameters
#' @param CI_res_all a simulation output matrix, with the number of rows equal to
#' \code{2 * length(M_range) * length(config_names)} and the number of columns
#' equal to number of simulations.
#' @param chi_orc a vector of optimal half-lengths, corresponding to
#' each value in \code{M_range}.
#' @param alpha a desired level of non-coverage
#'
#' @return a list of two dataframes, corresponding to coverage (\code{cov_data}) and
#' half-length (\code{chi_data}) results, respectively.
#' @export
#'
#' @examples config_names <- c("a", "b")
#' M_range <- c(1, 2)
#' nrow <- 2 * length(M_range) * length(config_names)
#' CI_res_all <- matrix(rnorm(nrow * 50), nrow, 50)
#' chi_orc <- c(1, 2)
#' res_format(config_names, M_range, CI_res_all, chi_orc, 0.05)
res_format <- function(config_names, M_range, CI_res_all, chi_orc, alpha){

  cov_names <- config_names
  chi_names <- c(config_names, "Orc. Chi")

  config_len <- length(config_names)

  covval <- numeric(0)
  chival <- numeric(0)

  grid_len <- length(M_range)

  for (i in 0:(config_len - 1)){

    ind1 <- 2 * i * grid_len + 1
    ind2 <- (2 * i + 1) * grid_len
    covres <- CI_res_all[ind1:ind2, ]
    chires <- CI_res_all[(ind1 + grid_len):(ind2 + grid_len), ]
    cov_mean <- rowMeans(covres)
    chi_mean <- rowMeans(chires)

    covval = c(covval, cov_mean)
    chival = c(chival, chi_mean)
  }

  chival = c(chival, chi_orc)

  cov_data <- data.frame(M = rep(M_range, length(cov_names)), covval = covval,
                         ind = rep(cov_names, each = grid_len))

  chi_data <- data.frame(M = rep(M_range, length(chi_names)),
                         chival = chival,
                         ind = rep(chi_names, each = grid_len))

  res <- list(cov_data = cov_data, chi_data = chi_data)
  return(res)
}

#' Specification Setup
#'
#' Sets up simulation settings.
#'
#' This function was created for the purpose of simulations.
#'
#' @param spec_num specification number to be used.
#' @param design_name the design for the true mean vector;
#' currently supports \code{c("unif", "tnorm", "norm", "t", "sk-t")}.
#' @param n the length of the true mean vector.
#'
#' @return a list of components to be used in the simulation
#' @export
#'
#' @examples set_spec(1, "unif", 100)
#' set_spec(2, "unif", 100)
set_spec <- function(spec_num, design_name, n){

  alpha <- 0.05

  Jn_exp <- 1/3
  Jn_til <- 2 * (log(n)/log(log(n)))^Jn_exp
  Jnval <- floor(Jn_til) + (1 - floor(Jn_til)%%2)  # To make Jn odd

  Jnval2 <- Jnval + 2 # Used for gnB_norm method

  # Default lambda values
  expnt_eps = 0.001
  lammin_cons <- 1
  lamlen <- 5

  expnt = 1/2 - expnt_eps
  lammin <- lammin_cons * 1 / (1 + Jn_til^(expnt))
  lamvec <- seq(from = lammin, to = 0.45, length.out = lamlen)

  # Setting theta designs
  grid_len <- 20
  Mmax <- 3
  M_range <- seq(from = 0.1, to = Mmax, length.out = grid_len)
  th_mat <- matrix(0, nrow = n, ncol = grid_len)
  set.seed(1)
  for(i in 1:grid_len){

    th_mat[ , i] <- th_gen(n = n, M = M_range[i], method = design_name)
  }

  # Tolerance level for bias
  btol <- alpha / 2

  # Standard deviations

  sig_vec <- rep(1, n)

  if(spec_num == 2){

    btol <- 1

  }else if(spec_num == 3){

    btol <- Inf
    lammin_cons <- 0.5
    lammin <- lammin_cons * 1 / (1 + Jn_til^(expnt))
    lamvec <- seq(from = lammin, to = 0.45, length.out = lamlen)

  }else if(spec_num == 4){

    btol <- Inf
    lammin_cons <- 0.25
    lammin <- lammin_cons * 1 / (1 + Jn_til^(expnt))
    lamvec <- seq(from = lammin, to = 0.45, length.out = lamlen)

    Jnval2 <- Jnval + 6
  }

  res <- list(Jnval = Jnval, lamvec = lamvec, M_range = M_range, th_mat = th_mat, btol = btol,
              alpha = alpha, sig_vec = sig_vec, Jnval2 = Jnval2)
  return(res)
}


#' CI Result Generation for Multiple Methods
#'
#' Generates half-length and coverage results for multiple methods.
#'
#' This function was created for the purpose of simulations.
#'
#' @param methods a vector of method names to use; currently support
#' \code{c("Df", "Sp", "Norm", "ebci", "hybr")}.
#' @param spec_list a list object created by \code{set_spec}.
#' @param min_mu2 the cut-off point for ebci method
#'
#' @return a vector of half-length and coverage results
#' @export
#'
#' @examples spec_list <- OptACI::set_spec(1, "unif", 100)
#' methods <- c("Df", "Sp", "Norm", "ebci", "hybr")
#' res_gen(methods, spec_list, 0.03)
res_gen <- function(methods, spec_list, min_mu2){

  alpha <- spec_list$alpha
  Jnval <- spec_list$Jnval
  lamvec <- spec_list$lamvec
  M_range <- spec_list$M_range
  th_mat <- spec_list$th_mat
  btol <- spec_list$btol
  sig_vec <- spec_list$sig_vec

  grid_len <- length(M_range)
  c_len <- length(methods)

  Jnval2 <- spec_list$Jnval2 # Used for gnB_norm method

  res <- numeric(grid_len * 2 * c_len) # Result matrix

  for(i in 1:grid_len){

    # x generation -------------------------------------------------
    th_vec <- th_mat[ , i]
    n <- length(th_vec)
    x_vec <- stats::rnorm(n, th_vec, sig_vec)

    # choice of method -------------------------------------------------
    ACI_res_gen <- function(c){

      if(c == "Df"){

        res <- ACI(x_vec, sig_vec, alpha, Jnval, lamvec, tol = btol)

      }else if(c == "Sp"){

        res <- ACI(x_vec, sig_vec, alpha)

      }else if(c == "Norm"){

        res <- ACI_aux(x_vec, sig_vec, alpha, Jnval, lamvec, Jnval2)

      }else if(c == "ebci"){

        res <- ACI_AKP(x_vec, sig_vec, alpha, min_mu2)

      }else if(c == "hybr"){

        res <- ACI_hybr(x_vec, sig_vec, alpha, Jnval, lamvec, Jnval2, min_mu2)
      }

      return(res)
    }

    # Storing results -------------------------------------------------

    for(c in 1:c_len){

      met_c = methods[c]
      ACIres <- ACI_res_gen(met_c)

      res[i + (2 * c - 2) * grid_len] <- ACP_check(th_vec, ACIres$CI_l, ACIres$CI_u)
      res[i + (2 * c - 1) * grid_len] <- ACIres$chi_opt
    }
  }

  return(res)
}

#' Oralce Half-legnth Generation
#'
#' @inheritParams res_gen
#'
#' @return a vector of oracle half-length values
#' @export
#'
#' @examples spec_list <- OptACI::set_spec(1, "unif", 100)
#' res_chi_orc(spec_list)
res_chi_orc <- function(spec_list){

  alpha <- spec_list$alpha
  M_range <- spec_list$M_range
  th_mat <- spec_list$th_mat
  sig_vec <- spec_list$sig_vec

  grid_len <- length(M_range)

  chi_orc <- numeric(grid_len)
  for(i in 1:grid_len){

    th_vec <- th_mat[ , i]
    ACI_orc_res <- ACI_orc(th_vec, sig_vec, alpha)
    chi_orc[i] <- ACI_orc_res$chi_opt

  }

  return(chi_orc)
}