R/distn_process.R

In DPComb: Discrete p-Value Combination Tests

#This file contains supportive functions for processing discrete distributions.
#Create by ZWu, Feb 8, 2025.

#' @title Convert X Distribution to the p-Value Distribution
#' @description This function converts a discrete X statistic and its distribution (characterized by its support vector and probability mass vector) to the p-value and its distribution (characterized by the p-value support, with increasingly sorted unique possible elements).
#' @param x An observed X value. Default is NULL (for the purpose of getting the `p_support` only).
#' @param x_support Vector of support of X, i.e., all possible unique X values sorted increasingly.
#' @param x_prob Vector of probability masses, matching with elements of `x_support`.
#' @param side Side of p-value. One of "two" (default), "left", or "right", specifying the tail for the p-value computation.
#' @return A list with elements:
#'   \item{p}{The observed p-value corresponding to the input `x`. It is NULL if `x` is NULL.}
#'   \item{p_support}{Unique sorted vector of p-value support.}
#'   \item{p_by_x}{Vector of possible p-values matching the X support. This is not used in the main function, but can be useful for debugging.}
#' @details
#' Let X be a random discrete variable. The p-values are defined as:
#' 
#' For left-sided:
#'   \deqn{p(x) = P(X \le x)}
#' 
#' For right-sided:
#'   \deqn{p(x) = P(X \ge x)}
#' 
#' For two-sided:
#'   \deqn{p(x) = \sum_{y: P(X \le y) \le P(X \le x)} P(X = y)}
#' 
#' By definition, `p_support` does not contain 0 and always contains 1. However, the `p_support` may contain repeated 0's and 1's due to numerical limitations. The package has tried to account for this issue.
#' @examples
#' x_support <- 0:5
#' x_prob <- dbinom(x_support, size = 5, prob = 0.1)
#' convert_x_to_p(3, x_support, x_prob, side = "left")
#' convert_x_to_p(3, x_support, x_prob, side = "right")
#' convert_x_to_p(3, x_support, x_prob, side = "two")
#' convert_x_to_p(NULL, x_support, x_prob, side = "two") # x can be NULL: output p is NULL
#' 
#' # x must be in x_support if not NULL. Error message otherwise.
#' # convert_x_to_p(6, x_support, x_prob, side = "two") #This gives Error as expected.
#' 
#' # symmetric distribution leads to duplicated values in p_by_x and a shorter p_support.
#' x_prob <- dbinom(x_support, size = 5, prob = 0.5)
#' convert_x_to_p(3, x_support, x_prob, side = "two")
#' 
#' # For large distributions, the p_support is truncated for numerical stability.
#' # It is still valid for practical hypothesis testing purposes.
#' n = 100000
#' x_support = 0:n
#' x_prob = dbinom(x_support, size = n, prob = 0.7)
#' result = convert_x_to_p(3, x_support, x_prob, side = "left")
#' result$p
#' length(result$p_support) # truncated
#' @export
convert_x_to_p <- function(x=NULL, x_support, x_prob, side = "two") {
  #validations...
  if(length(x_support) != length(x_prob)) stop("x_support and x_prob must have equal length.")
  
  # Compute the p-values matching with all elements in x_support
  if(side == "left"){
    p_by_x <- cumsum(x_prob)
  } else if(side == "right"){
    p_by_x <- rev(cumsum(rev(x_prob)))
  } else if(side == "two"){
    # For two-sided, sum of probabilities <= p at each point 
    p_by_x <- sapply(1:length(x_prob), function(i) sum(x_prob[x_prob <= x_prob[i]]))
  } else {
    stop("Invalid side. Choose one of 'left', 'right', or 'two'.")
  }

  # Convert any elements > 1 (which is possible due to numerical additions) to be 1
  p_by_x[p_by_x > 1] <- 1
  
  # Compute unique sorted p-value support
  # This also removes many 0's and 1's in distributions with very small probabily mass at the tails
  p_support <- sort(unique(p_by_x))

  if (is.null(x)) {
    p <- NULL
  } else {
    if(!(x %in% x_support)) stop("Input x must be NULL or in the x_support vector.")
    idx <- which(x_support == x)
    p <- p_by_x[idx[1]]
  }
  
  return(list(p = p,
              p_support = p_support,
              p_by_x = p_by_x))
}



#' @title Get Support and Probability Vectors from a Distribution Description
#' @description Generates the support vector (i.e., increasingly sorted unique possible values) and the corresponding probability mass vector for a given discrete distribution description, including the name of the distribution and its parameters. Allows for replication.
#' @param distn_params A list containing: 1) `distn`: the distribution name. One of "binom", "pois", "hyper", "noncenhypergeom", "nbinom", or "geom". 2) Parameters for the specified distribution. They are consistent with the corresponding R density functions. See details.
#' @param repN The number of support/probability vectors to generate. Default is 1.
#' @return A list of `repN` replicated lists, each containing the `x_support` and `x_prob` vectors.
#' @details 
#' Supported distributions and their parameters:
#' - Binomial: list(distn = "binom", size = , prob = ).
#' - Poisson: list(distn = "pois", lambda = ).
#' - Hypergeometric: list(distn = "hyper", m = , n = , k = ).
#' - Noncentral Hypergeometric: list(distn = "noncenhypergeom", n1 = , n2 = , m1 = , psi = ). Requires the `MCMCpack` package.
#' - Negative Binomial: list(distn = "nbinom", size = , prob = ).
#' - Geometric: list(distn = "geom", prob = ).
#' @examples
#' # For binomial distribution
#' distn_params <- list(distn = "binom", size = 5, prob = 0.1)
#' distn_to_x_support_probs(distn_params)
#' distn_to_x_support_probs(distn_params, 2) # replicate 2 times
#' 
#' # For Poisson distribution
#' distn_params <- list(distn = "pois", lambda = 100)
#' distn_to_x_support_probs(distn_params)
#' 
#' # For hypergeometric distribution
#' distn_params <- list(distn = "hyper", m = 10, n = 5, k = 3)
#' distn_to_x_support_probs(distn_params)
#' 
#' # For noncentral hypergeometric distribution
#' distn_params <- list(distn = "noncenhypergeom", n1 = 10, n2 = 5, m1 = 3, psi = 1)
#' distn_to_x_support_probs(distn_params) # same as the hypergeometric distribution when psi = 1
#' distn_params <- list(distn = "noncenhypergeom", n1 = 10, n2 = 5, m1 = 3, psi = 0.5)
#' distn_to_x_support_probs(distn_params) # smaller psi makes higher prob for smaller values
#' distn_params <- list(distn = "noncenhypergeom", n1 = 10, n2 = 5, m1 = 3, psi = 2)
#' distn_to_x_support_probs(distn_params) # larger psi makes higher prob for larger values
#' 
#' # For negative binomial distribution
#' distn_params <- list(distn = "nbinom", size = 5, prob = 0.1)
#' distn_to_x_support_probs(distn_params)
#' 
#' # For geometric distribution
#' distn_params <- list(distn = "geom", prob = 0.1)
#' distn_to_x_support_probs(distn_params)
#' @export
distn_to_x_support_probs <- function(distn_params, repN=1) {
  
  # Validate that distn_params has a 'distn' element.
  if (is.null(distn_params$distn)) {
    stop("distn_params must include the element 'distn'.")
  }
  
  # Generate support and probability vector based on the specified distribution.
  if (distn_params$distn == "binom") {
    if (!all(c("size", "prob") %in% names(distn_params)))
      stop("For 'binom', distn_params must include 'size' and 'prob'.")
    support_i <- 0:distn_params$size
    prob_i <- dbinom(support_i, size = distn_params$size, prob = distn_params$prob)

  } else if (distn_params$distn == "pois") {
    if (!"lambda" %in% names(distn_params))
      stop("For 'pois', distn_params must include 'lambda'.")
    support_i <- 0:qpois(0.9999, distn_params$lambda)
    prob_i <- dpois(support_i, lambda = distn_params$lambda)
    
  } else if (distn_params$distn == "hyper") {
    if (!all(c("m", "n", "k") %in% names(distn_params)))
      stop("For 'hyper', distn_params must include 'm', 'n', and 'k'.")
    m=distn_params$m
    n=distn_params$n
    k=distn_params$k
    support_i <- seq(max(0, k - n), min(k, m))
    prob_i <- dhyper(support_i, m = m, n = n, k = k)
    
  } else if (distn_params$distn == "noncenhypergeom") {
    # Check for required package
    if (!requireNamespace("MCMCpack", quietly = TRUE)) {
      stop("Package 'MCMCpack' required for noncentral hypergeometric distribution. Please install it.")
    }
    # Mapping: use m1 as m, n1 as n, n2 as k, and psi as odds.
    if (!all(c("n1", "n2", "m1", "psi") %in% names(distn_params)))
      stop("For 'noncenhypergeom', distn_params must include 'n1', 'n2', 'm1', and 'psi'.")
    n1=distn_params$n1
    n2=distn_params$n2
    m1=distn_params$m1
    psi=distn_params$psi
    support_i <- seq(max(0, m1 - n2), min(m1, n1))
    prob_i <- unlist(lapply(support_i, function(s) MCMCpack::dnoncenhypergeom(s, n1=n1, n2=n2, m1=m1, psi=psi)))   
    
  } else if (distn_params$distn == "nbinom") {
    if (!all(c("size", "prob") %in% names(distn_params)))
      stop("For 'nbinom', distn_params must include 'size' and 'prob'.")
    support_i <- 0:qnbinom(0.9999, size = distn_params$size, prob = distn_params$prob)
    prob_i <- dnbinom(support_i, size = distn_params$size, prob = distn_params$prob)
    
  } else if (distn_params$distn == "geom") {
    if (!"prob" %in% names(distn_params))
      stop("For 'geom', distn_params must include 'prob'.")
    support_i <- 0:qgeom(0.9999, prob = distn_params$prob) +1 ### geometric is defined as number of successes in R
    prob_i <- dgeom(support_i-1, prob = distn_params$prob)
    
  } else if (distn_params$distn == "multinom") {
          stop("The 'multinom' distribution is not supported in this function.")
    
  } else {
    stop("Unsupported distribution. Choose one of 'binom', 'pois', 'hyper', 'noncenhypergeom', 'nbinom', or 'geom'.")
  }
  
  # Create a list of length repN, each component contains both support and probability.
  if(repN == 1) {
    return(list(x_support = support_i, x_prob = prob_i))
  } else if (repN > 1)  {
    res <- vector("list", repN)
    for(i in seq_len(repN)) {
      res[[i]] <- list(x_support = support_i, x_prob = prob_i)
    }
    return(res)
  } else {
    stop("repN must be a positive integer.")
  }

}