R/Pval.R

In MHTdiscrete: Multiple Hypotheses Testing for Discrete Data

#' Calculating p-values for discrete data
#'
#' The function for calculating the original available p-values and all attaianble p-values for the corresponding hypothesis.
#'
#' @usage
#' getPval(raw.data, test.type, alternative)
#' @param raw.data original data set with count number for treatment group and study group. The data set type could be \code{\link[base]{matrix}} or \code{\link[base]{data.frame}}.
#' @param test.type there are two discrete test available now, must be one of \code{"FET"} for Fisher's Exact Test and \code{"BET"} for Binomial Exact Test.
#' @param alternative indicates the alternative hypothesis and must be one of \code{"two.sided"}, \code{"greater"} or \code{"less"}.
#' @return
#' A numeric vector of the adjusted p-values (of the same length as \eqn{p}).
#' @author Yalin Zhu
#' @references
#' Zhu, Y., & Guo, W. (2017).
#' Familywise error rate controlling procedures for discrete data
#' \emph{arXiv preprint} arXiv:1711.08147.
#'
#' Clopper, C. J. & Pearson, E. S. (1934).
#' The use of confidence or fiducial limits illustrated in the case of the binomial.
#' \emph{Biometrika}, \strong{26}: 404-413.
#'
#' Fisher, R. A. (1922).
#' On the Interpretation of \eqn{\chi^2} from Contingency Tables, and the Calculation of P.
#' \emph{Journal of the Royal Statistical Society}, \strong{85}: 87-94.
#'
#' @examples
#'  ## Using Fisher's Exact Test to get the avaiable and attainablep-values
#'  # import raw data set as data.frame type
#'  df <-  data.frame(X1=c(4, 2, 2, 13, 6, 8, 4, 0, 1), N1 = rep(148, 9),
#'  	X2 = c(0, 0, 1, 3, 2, 1, 2, 2, 2), N2 = rep(132, 9))
#'  # obtain the avaiable p-values and attainable p-values using two-sided Fisher's Exact Test
#'  getPval(raw.data=df, test.type = "FET",alternative = "two.sided")
#'  # store the avaiable p-values
#' p <- getPval(raw.data=df, test.type = "FET",alternative = "two.sided")[[1]]; p
#'  # store the attainable p-values
#' p.set <- getPval(raw.data=df, test.type = "FET",alternative = "two.sided")[[2]]; p.set
#' @export

getPval <- function(raw.data, test.type = c("FET", "BET"), alternative = c("two.sided","greater","less")){
  if (class(raw.data)!="data.frame" & class(raw.data)!="matrix"){
    stop("The raw data type has to be 'data.frame' or 'matrix'")
  }
  else{
    test.type <- match.arg(test.type)
    alternative <- match.arg(alternative)
    len <- nrow(raw.data)
    p <- numeric(len)
    p.set <- vector("list", len)
    if (test.type=="FET"){
      ssum <- raw.data[,1]+raw.data[,3]
      slim <- pmin(ssum, raw.data[,2])

      for (i in 1:len){
        p[i] <- stats::fisher.test(matrix(c(raw.data[i,1],
                                            raw.data[i,2]-raw.data[i,1],
                                            raw.data[i,3],
                                            raw.data[i,4]-raw.data[i,3]),2,2),
                                   alternative = alternative)[[1]]
        s<-numeric(slim[i] + 1)
        for (j in 0:slim[i]){
          s[j+1] <- stats::fisher.test(matrix(c(j,
                                                raw.data[i,2]-j,
                                                ssum[i]-j,
                                                raw.data[i,4]-(ssum[i]-j)),2,2), alternative = alternative)[[1]]
        }
        p.set[[i]] <- sort(s)
      }
    }
    if (test.type=="BET"){
      ssum <- raw.data[,1]+raw.data[,2]
      for (i in 1:len){
        p[i] <- stats::binom.test(c(raw.data[i,1],raw.data[i,2]), p = 0.5, alternative = alternative)[[3]]
        s <- numeric(ssum[i] + 1)
        for (j in 0:ssum[i]){
          s[j+1] <- stats::binom.test(c(j, ssum[i]-j), p = 0.5, alternative = alternative)[[3]]
        }
        p.set[[i]] <- sort(s)
      }
    }
    return(list(available.pvalue=p, attainable.pvalue=p.set))
  }
}