R/est.prop.R

Defines functions est.prop

Documented in est.prop

#' @title Proportion Estimation
#' @description Estimates false null hypothesis Proportion from multiple p-values using higher criticism test estimator.
#'
#' @param p.value A sequence of p-values from test data, not including p-values from covariates.
#' @param cn A value of bounding sequence generated by HCTR::bounding.seq().
#' @param adj A boolean algebra to decide whether to use adjusted Higher Criticism test statistic, the default value is TRUE.
#'
#' @return An estimated proportion of false null hypothesis.
#'
#' @export
#' @importFrom Rdpack reprompt
#'
#' @examples
#' set.seed(10)
#' X <- matrix(runif(n = 10000, min = 0, max = 1), nrow = 100)
#' result <- bounding.seq(p.value = X)
#' Y <- matrix(runif(n = 100, min = 0, max = 1), nrow = 100)
#' test <- est.prop(p.value = Y, cn = result)
#'
#' @references
#' \insertRef{meinshausen2006estimating}{HCTR}
est.prop <- function(p.value, cn, adj = TRUE) {
  n <- length(p.value)
  p_order <- sort(p.value, index.return = TRUE)$x
  pihat <- rep(0, n)
  ind <- 1 : n
  if(adj == T){
    pihat <- (ind / n - p_order - cn * sqrt(p_order * (1 - p_order))) / (1 - p_order)
  }
  else{
    pihat <- (ind / n - p_order - sqrt(2 * log(log(n))/n) * sqrt(p_order * (1 - p_order))) / (1 - p_order)
  }
  est <- max(pihat[2 : floor(n / 2)], 0)
  return(est)
}

Try the HCTR package in your browser

Any scripts or data that you put into this service are public.

HCTR documentation built on Dec. 1, 2019, 1:21 a.m.