R/DPComb_Functions.R

In DPComb: Discrete p-Value Combination Tests

# This file contains user's functions for the DPCombine package
# Created by ZWu, Feb 7, 2025. 


#' @title Discrete P-value Combination Tests 
#'
#' @description This function combines evidence of significance to test a global null hypothesis that the given discrete null distributions are true. The tests can be based on discrete p-values or the corresponding discrete X statistics, and their null distributions.  
#'
#' @param ps Optional. A vector of observed discrete p-values to be combined.  
#' @param p_supports Optional. A vector or list of p-value supports, characterizing the null distribution of discrete p-values. A vector can be used when all p-values are from the same distribution under the null, i.e., the i.i.d. case. A list allows for non-identical distributions, containing p-value support vectors matching the distributions of the elements in `ps`. 
#' @param method The combination testing method. One of "fisher_mean" (default), "fisher_median", "pearson", "edgington", "stouffer", or "george". See the details. 
#' @param xs Optional. A vector of observed discrete X statistics, from which the p-values in `ps` are obtained. It is not used if `ps` is given, but required if `ps` is NULL.
#' @param side sidedness of p-values to be combined. One of "two" (default), "right" or "left". 
#' @param x_support_probs Optional. A list containing `x_support` and `x_prob` vectors to characterize the null distributions of X statistics under the null distribution. For the i.i.d. case, where each element in `xs` follows the same distribution, it is a list of two vectors `x_support` and `x_prob`. For non-identical distributions, `x_support_probs` is required to be a list of lists, each sublist containing the `x_support` and `x_prob` vectors to define the null distributions of X, corresponding to the elements of `xs`. The number of sublists should equal the length of `xs` or `ps`. 
#' @param x_distn_params Optional. A list containing the description of null distributions for the discrete X statistics. See the details. For the i.i.d. case, it is a simple list. For non-identical distributions, `x_distn_params` is required to be a list of lists, each sublist defining the null distributions of one X statistic, consistent with `xs` or `ps`. The number of sublists should equal the length of `xs`. 
#' @return A list with elements:
#'   \item{Sn}{The combination statistic}
#'   \item{pval}{The p-value of the combination statistic}
#' @details
#' This function calculates the following types of p-value combination statistics `Sn` and their testing p-values `pval`. 
#' \tabular{ll}{
#' \strong{Method} \tab \strong{Statistic} \cr
#' Fisher \tab \eqn{T_F = -2\sum_{j=1}^n \log P_j} \cr
#' Pearson \tab \eqn{T_P = -2\sum_{j=1}^n \log (1 - P_j)} \cr
#' George \tab \eqn{T_G = \frac{T_P - T_F}{2} = \sum_{j=1}^n \log \frac{P_j}{1 - P_j}} \cr
#' Stouffer \tab \eqn{T_S = \sum_{j=1}^n \Phi^{-1}(P_j)} \cr
#' Edgington \tab \eqn{T_E = \sum_{j=1}^n P_j} \cr
#' }
#' Smaller p-values represent higher significance against the null hypothesis. By its formula, a larger `Sn` in Fisher's combination method, or a smaller `Sn` in other combination methods, indicates a higher significance level against the null hypothesis. 
#' These statistics are adjusted by an "adjusted Z statistic" obtained by a Wasserstein distance optimization process. For Fisher's combination, the two methods "fisher_mean" and "fisher_median" correspond to Lancaster's mean-value chi-squared method and median-value chi-squared method, respectively.
#' 
#' The supported distribution descriptions and 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 = ).
#' 
#' @references
#' Lancaster, HO (1949). The combination of probabilities arising from data in discrete distributions. Biometrika, 36(3/4), 370-382.
#' 
#' Contador, Gonzalo and Wu, Zheyang (2025). A minimum Wasserstein distance approach to Fisher's combination of independent, discrete p-values. Scandinavian Journal of Statistics, 52(3), 1281-1300.
#'
#' Contador, Gonzalo and Wu, Zheyang (2026). Optimal Adjustment and Combination of Independent Discrete p-Values. Under revision at the Journal of Computational and Graphical Statistics.
#' @examples
#' # Example 1: p-values from the same distribution
#' p_supports <- seq(0.01, 1, length.out = 100)
#' methods <- c("fisher_mean", "fisher_median", "pearson", "george", "stouffer", "edgington")
#' 
#' ps <- c(0.1, 0.2, 0.21, 0.35)
#' sapply(methods, function(m) DPComb_tests(ps=ps, p_supports=p_supports, method=m))
#' 
#' # Example 2: p-values from different distributions
#' p_supports1 <- seq(0.01, 1, length.out = 10)
#' p_supports2 <- seq(0.3, 1, length.out = 5)
#' p_supports3 <- seq(0.01, 1, length.out = 100)
#' p_supports <- list(p_supports1, p_supports2, p_supports3)
#' ps <- c(0.12, 0.475, 0.21)
#' sapply(methods, function(m) DPComb_tests(ps=ps, p_supports=p_supports, method=m))
#' 
#' # Example 3: input xs and x_support_probs from the same distribution
#' xs <- c(0, 1, 2)
#' x_support_probs <- list(x_support = 0:5, x_prob = dbinom(0:5, size = 5, prob = 0.1))
#' sapply(methods, function(m) DPComb_tests(xs=xs, x_support_probs=x_support_probs, 
#'                                          side="two", method=m))
#' 
#' # Example 4: input xs and x_support_probs from different distributions
#' xs <- c(0, 1, 2)
#' x_supports <- list(0:5, 0:5, 0:5)
#' x_probs <- list(dbinom(0:5, size = 5, prob = 0.1), dbinom(0:5, size = 5, prob = 0.2), 
#'                 dbinom(0:5, size = 5, prob = 0.3))
#' x_support_probs <- lapply(1:length(x_supports), 
#'                          function(i) list(x_support = x_supports[[i]], 
#'                                           x_prob = x_probs[[i]]))
#' sapply(methods, function(m) DPComb_tests(xs=xs, x_support_probs=x_support_probs, 
#'                                          side="two", method=m))
#' sapply(methods, function(m) DPComb_tests(xs=xs, x_support_probs=x_support_probs, 
#'                                          side="right", method=m))
#' 
#' # Example 5: input xs and x_distn_params from the same distribution
#' xs <- c(0, 1, 2)
#' x_distn_params <- list(distn = "binom", size = 5, prob = 0.1)
#' sapply(methods, function(m) DPComb_tests(xs=xs, x_distn_params=x_distn_params, 
#'                                          side="two", method=m))
#' 
#' # Example 6: input xs and x_distn_params from different distributions
#' xs <- c(0, 1, 2)
#' x_distn_params <- list(list(distn = "binom", size = 5, prob = 0.1), 
#'                             list(distn = "binom", size = 5, prob = 0.2), 
#'                             list(distn = "binom", size = 5, prob = 0.3))
#' sapply(methods, function(m) DPComb_tests(xs=xs, x_distn_params=x_distn_params, 
#'                                          side="two", method=m))
#' sapply(methods, function(m) DPComb_tests(xs=xs, x_distn_params=x_distn_params, 
#'                                          side="right", method=m))
#' 
#' # Example 7: input xs and x_distn_params from different types of distributions: 
#' # binomial, poisson, and hypergeometric
#' xs <- c(0, 1, 2)
#' x_distn_params <- list(list(distn = "binom", size = 5, prob = 0.1), 
#'                        list(distn = "pois", lambda = 5), 
#'                        list(distn = "hyper", m = 10, n = 5, k = 3))
#' sapply(methods, function(m) DPComb_tests(xs=xs, x_distn_params=x_distn_params, 
#'                                          side="two", method=m))
#' # Example 8: input ps and x_support_probs. 
#' # Require ps and x_support_probs are consistent, i.e., ps are in the the p_supports 
#' # generated from x_support_probs. Avoid using this if unsure. 
#' ps <- c(0.00001, 0.00032, 0.00243)
#' x_support_probs <- list(list(x_support = 0:5, x_prob = dbinom(0:5, size = 5, prob = 0.1)),
#'             list(x_support = 0:5, x_prob = dbinom(0:5, size = 5, prob = 0.2)),
#'             list(x_support = 0:5, x_prob = dbinom(0:5, size = 5, prob = 0.3)))
#' sapply(methods, function(m) DPComb_tests(ps=ps, x_support_probs=x_support_probs, method=m))
#' sapply(methods, function(m) DPComb_tests(ps=ps, x_support_probs=x_support_probs, 
#'                                          side="right", method=m))
#' # Example 9: input ps and x_distn_params
#' # Require ps and x_distn_params are consistent, i.e., ps are in the the p_supports 
#' # generated from x_distn_params. Avoid using this if unsure. 
#' ps <- c(0.00001, 0.00032, 0.00243)
#' x_distn_params <- list(list(distn = "binom", size = 5, prob = 0.1), 
#'                        list(distn = "binom", size = 5, prob = 0.2), 
#'                        list(distn = "binom", size = 5, prob = 0.3))
#' sapply(methods, function(m) DPComb_tests(ps=ps, x_distn_params=x_distn_params, method=m))
#' 
#' @importFrom stats pgamma dbinom dgeom dhyper dnbinom dnorm dpois pnorm qgeom qnbinom qnorm qpois
#' @export
#' 
DPComb_tests <- function(ps = NULL, p_supports = NULL, 
                           xs = NULL, side = "two", 
                           x_support_probs = NULL,  
                           x_distn_params = NULL, 
                           method = "fisher_mean") {
  
  # Case 1: If ps is provided
  if (!is.null(ps)) {
    n <- length(ps)
    
    # If p_supports is not provided, use x_support_probs or x_distn_params to generate it
    if (is.null(p_supports)) {
      if (is.null(x_support_probs) && is.null(x_distn_params)) {
        stop("Either provide p_supports or specify x_support_probs or x_distn_params.")
      }
      
      p_supports <- vector("list", n) #initialize p_supports
      
      if (!is.null(x_support_probs)) {
        # Validate structure of x_support_probs: should be a list of length n or a single list to be replicated.
        if (!is.list(x_support_probs[[1]])) {
          # Single support/prob pair provided; replicate it for all ps
          x_support_probs <- rep(list(x_support_probs), n)
        } else if (length(x_support_probs) != n) {
          stop("Length of x_support_probs must equal length of ps.")
        }
      } else {
        # Generate support/probability list using distn_to_x_support_probs
        if (!is.list(x_distn_params[[1]])) {
          # Single distribution parameter list provided; replicate it for all ps
          x_support_probs <- distn_to_x_support_probs(x_distn_params, repN = n)
        } else if (length(x_distn_params) != n) {
          stop("Length of x_distn_params must equal length of ps.")
        } else {
          # Generate a list of lists, each containing x_support and x_prob
          x_support_probs <- lapply(x_distn_params, distn_to_x_support_probs)
        }
      }
      
      # Compute p_supports from x_support_probs using convert_x_to_p
      for (i in seq_len(n)) {
        p_supports[[i]] <- convert_x_to_p(x_support = x_support_probs[[i]]$x_support, 
                                          x_prob = x_support_probs[[i]]$x_prob, 
                                          side = side)$p_support
      }
    } else if (is.vector(p_supports) && !is.list(p_supports)) {
      p_supports <- rep(list(p_supports), n)
    } else if (length(p_supports) != n) {
      stop("Length of p_supports as a list must match length of ps.")
    }
    
  } else {
    # Case 2: If ps is not provided: require xs
    if (is.null(xs)) stop("xs must be provided if ps is NULL.")
    
    n <- length(xs)
    ps <- numeric(n)                #initialize ps
    p_supports <- vector("list", n) #initialize p_supports
    
    # If x_support_probs is not provided, use x_distn_params to generate it
    if (is.null(x_support_probs)) {
      if (is.null(x_distn_params))
        stop("Either provide x_support_probs or specify x_distn_params.")
      # Generate support/probability list using distn_to_x_support_probs
      if (!is.list(x_distn_params[[1]])) {
        # Single distribution parameter list provided; replicate it for all xs
        x_support_probs <- distn_to_x_support_probs(x_distn_params, repN = n)
      } else if (length(x_distn_params) != n) {
        stop("Length of x_distn_params must equal length of xs.")
      } else {
        # Generate a list of lists, each containing x_support and x_prob
        x_support_probs <- lapply(x_distn_params, distn_to_x_support_probs)
      }
    } else {
      # Validate structure of x_support_probs: should be a list of length n or a single list to be replicated.
      if (!is.list(x_support_probs[[1]])) {
        # Single support/prob pair provided; replicate it for all xs
        x_support_probs <- rep(list(x_support_probs), n)
      } else if (length(x_support_probs) != n) {
        stop("Length of x_support_probs must equal length of xs.")
      }
    }
    
    # Compute ps and p_supports from xs using convert_x_to_p
    for (i in seq_len(n)) {
      out <- convert_x_to_p(x = xs[i], 
                            x_support = x_support_probs[[i]]$x_support, 
                            x_prob = x_support_probs[[i]]$x_prob, 
                            side = side)
      ps[i] <- out$p
      p_supports[[i]] <- out$p_support
    }
  }

  # Compute adjusted Z statistics and their moments for each p-value.
  Zs <- lapply(1:length(ps), function(i) adjZ_moments(ps[i], p_supports[[i]], method))
  Sn <- sum(sapply(Zs, function(z) z$Z))
  Sn_mean <- sum(sapply(Zs, function(z) z$Zmean))
  Sn_var <- sum(sapply(Zs, function(z) z$Zvar))
  
  # Compute combined p-value based on selected method.
  # Larger Sn in fisher_mean and fisher_median is more significant.
  # Smaller Sn in pearson, stouffer, george, and edgington is more significant.
  if(method %in% c("fisher_mean", "fisher_median")) {
    pval <- pgamma(Sn, shape = Sn_mean^2/Sn_var, scale = Sn_var/Sn_mean, lower.tail = FALSE)
  } else if(method == "pearson") {
    pval <- pgamma(Sn, shape = Sn_mean^2/Sn_var, scale = Sn_var/Sn_mean)
  } else if(method %in% c("stouffer", "george")) {
    pval <- pnorm(Sn, mean = 0, sd = sqrt(Sn_var))
  } else if(method == "edgington") {
    pval <- pnorm(Sn, mean = length(ps)/2, sd = sqrt(Sn_var))
  } else {
    stop("Invalid method. Choose one of 'fisher_mean', 'fisher_median', 'pearson', 'edgington', 'stouffer', or 'george'.")
  }
  
  return(list(Sn = Sn, pval = pval))
}

#' Combine Fisher's Exact Tests in Case-Control Analysis
#'
#' This function performs a combination test on case-control data with binary covariates 
#' by first obtaining p-values from Fisher's exact test for each covariate and then 
#' combining them using the \code{DPComb_tests} function.
#'
#' @param Data A dataframe containing case-control data. It must include a column
#'   with 1 for cases and 0 for controls, and one or more
#'   columns representing binary covariates (with values 1 or 0).
#' @param response A character string specifying the name of the response variable indicating cases and controls.
#' @param covariates A vector of character strings specifying the names of the covariate columns to analyze.
#' @param method A character string specifying the combination method to be used.
#'   Default is \code{"fisher_mean"}. Supported methods include \code{"fisher_mean"},
#'   \code{"fisher_median"}, \code{"pearson"}, \code{"george"}, \code{"stouffer"}, and \code{"edgington"}.
#' @param alternative A character string indicating the tail for Fisher's exact test. One of "two.sided" (default), "greater", or "less", consistent with the \code{fisher.test} function.
#'
#' @return A list with the following elements:
#' \describe{
#'   \item{Sn}{The testing statistic combining statistical significance from Fisher's exact tests.}
#'   \item{pval}{The testing p-value for the combination test.}
#' }
#'
#' @details
#' For each covariate, Fisher's exact test is performed to compute a p-value, based on the hypergeometric distribution under the null hypothesis. The parameters are derived from the total number of cases, controls, and the total count of 1's in each covariate. The \code{DPComb_tests} function is then applied to compute the test statistic and associated p-value for the combination of Fisher's exact tests.
#'
#' @examples
#'   # Load case-control data from DPComb
#'   data(case_control, package = "DPComb")
#'   covariates <- c("marker1", "marker2", "marker3", "marker4", "marker5")
#'   test_case_control_fisher(Data = case_control, response = "disease_status", 
#'                            covariates = covariates, 
#'                            method = "fisher_mean", alternative = "two.sided")
#' @importFrom stats fisher.test
#' @export
test_case_control_fisher <- function(Data, response, covariates, method = "fisher_mean", alternative = "two.sided") {
  # The case-control response 
  status <- Data[[response]]
  
  # Compute p-values for each covariate using Fisher's exact test
  ps <- sapply(covariates, function(m) {
    tbl <- table(status, Data[[m]])
    fisher.test(tbl, alternative = alternative)$p.value
  })
  ps[ps > 1] <- 1  # Ensure p-values do not exceed 1
  
  # Derive hypergeometric distribution parameters for each covariate
  total_cases <- sum(status == 1)
  total_controls <- sum(status == 0)
  x_distn_params <- lapply(covariates, function(m) {
    k <- sum(Data[[m]])
    list(distn = "hyper", m = total_cases, n = total_controls, k = k)
  })
  
  # Run the combination test using DPComb_tests
  # Determine appropriate side for DPComb_tests
  side <- switch(alternative,
                 "two.sided" = "two",
                 "greater" = "right",
                 "less" = "left",
                 stop("Invalid alternative. Use 'two.sided', 'greater', or 'less'."))
                      
  res <- DPComb_tests(ps = ps, side = side, x_distn_params = x_distn_params, method = method)
  
  return(res)
}