R/calcu_adjZ.R

In DPComb: Discrete p-Value Combination Tests

# This file contains functions for calculating the adjusted Z statistics and their mean and variance.  
# Created by ZWu 2025-02-10 

#' @title Calculate Adjusted Z Statistic and Its Moments
#' @description Calculates the adjusted Z statistic and its mean and variance corresponding to an observed p-value and the distribution of the p-value (which is fully characterized by its support vector).
#' @param p An observed p-value.
#' @param p_support A valid p-value support vector containing increasingly sorted possible p-values.
#' @param method The p-value combination method, one of "fisher_mean" (default), "fisher_median", "pearson", "edgington", "stouffer", or "george".
#' @return A list with the following elements:
#'   \item{Z}{The adjusted Z statistic corresponding to the observed p-value.}
#'   \item{Zmean}{The mean of the adjusted Z statistic.}
#'   \item{Zvar}{The variance of the adjusted Z statistic.}
#' @details
#'    The input `p` should be in `p_support`. If not, the `p_support` element closest to `p` will be used to calculate `Z`, with a warning message. `Zmean` and `Zvar` are calculated based on `p_support`. The adjustment is made based on a specific p-value combination method.
#' @examples
#' # Example usage:
#' methods <- c("fisher_mean", "fisher_median", "pearson", "george", "stouffer", "edgington")
#' p_support <- seq(0.01, 1, length.out = 100)
#' sapply(methods, function(m) adjZ_moments(p_support[10], p_support, m))
#' sapply(methods, function(m) adjZ_moments(0.105, p_support, m)) # Warning if p is not in p_support.
#' @export
adjZ_moments <- function(p, p_support, method="fisher_mean") {
    # Find the index in p_support corresponding to the observed p-value
    idx <- which.min(abs(p_support - p)) # Find the index of the closest match
    if (abs(p_support[idx] - p) > 1e-8) {
        warning("The closest p_support value differs from p by more than 1e-8.")
    }

  # Find the previous p-value less than the given p-value
  p_prev <- if (idx == 1) 0 else p_support[idx-1]

  Z = computeZ(f=p, f_prev=p_prev, method=method) 
  # Compute the mean and variance of the adjusted Z statistic
  mts = computeZmoment(p_support=p_support, method=method)

  return(list(Z = Z, Zmean = mts$Zmean, Zvar = mts$Zvar))
}


#' @title Compute Adjusted Z Statistic
#'
#' @description This core function computes the combination method-related adjusted Z statistic for a given observed discrete p-value.
#' @param f A given legitimate discrete p-value. It should be > 0 and <= 1.
#' @param f_prev The next smaller p-value (i.e., the previous element in the p-value support vector). It is zero if \code{f} is the smallest possible p-value in its support.
#' @param method The combination method that the adjusted Z statistic is related to. One of "fisher_mean", "fisher_median", "pearson", "george", "stouffer", or "edgington".
#' @return The adjusted Z statistic.
#' @details
#' The following are the formulas for the adjusted Z statistic when the p-value \eqn{P=F_i}. 
#' Notations: \eqn{\overline{F_i} \equiv 1-F_i}; \eqn{K_i \equiv (2\pi)^{-1/2}\exp [-\Phi^{-1}(F_i)^2/2]}, where \eqn{\Phi} is the cumulative distribution function of the standard normal distribution.
#' 
#' \tabular{lccc}{
#' \strong{Method} \tab \strong{Statistic} \tab \strong{Value when \eqn{P=F_i}} \cr
#' Fisher \tab \eqn{Z_F} \tab \eqn{2-2(F_i-F_{i-1})^{-1}(F_i\log F_i -F_{i-1}\log F_{i-1})} \cr
#' Pearson \tab \eqn{Z_P} \tab \eqn{2-2(F_i-F_{i-1})^{-1}(\overline{F_{i-1}}\log \overline{F_{i-1}}-\overline{F_i}\log \overline{F_i})} \cr
#' George \tab \eqn{Z_G} \tab \eqn{(Z_P-Z_F)/2} \cr
#' Stouffer \tab \eqn{Z_S} \tab \eqn{(F_i-F_{i-1})^{-1}\left[K_{i-1}-K_i\right]} \cr
#' Edgington \tab \eqn{Z_E} \tab \eqn{(F_i+F_{i-1})/2} \cr
#' }
#'
#' @examples
#' methods = c("fisher_mean", "fisher_median", "pearson", "george", "stouffer", "edgington")
#' sapply(methods, function(m) computeZ(0.1, 0.05, m))
#' sapply(methods, function(m) computeZ(0.1, 0, m)) # f_prev = 0
#' sapply(methods, function(m) computeZ(1, 0.9, m)) # f = 1
#' sapply(methods, function(m) computeZ(0.1, 0.1 - 1e-10, m)) # f - f_prev is small
#' @export
computeZ <- function(f, f_prev, method) {
    #stop if f or f_prev is > 1
    if (f > 1 || f_prev > 1) stop("Detected p-value > 1, which is invalid.")
        
    eps <- 1e-10 # small value for numerical stability
    d <- f - f_prev #the difference is also the probability for p-value = f

    if (method == "fisher_mean") {
        if (f_prev == 0) {
            return(2 * (1 - log(f))) #Simplify formula to avoid log(f_prev) if f_prev = 0
        } 
        if (d < eps) {
            return(-2 * log((f + f_prev) / 2)) #use midP to avoid division by d.
        } else {
            return(2 * (1 - (f * log(f) - f_prev * log(f_prev)) / d))
        }

    } else if (method == "fisher_median") {
        return(-2 * log((f + f_prev) / 2))

    } else if (method == "pearson") {
        if (f > 1 - eps) { #f == 1 case, but include eps for numerical stability
            return(2 * (1 - log(1-f_prev))) #Simplify formula to avoid log(1-f) if f = 1
        }
        if (d < eps) {
            return(-2 * log(((1 - f) + (1 - f_prev)) / 2))  #use 1-midP if to avoid division by d.
        } else {
            return(2 - 2 * ((1 - f_prev) * log(1 - f_prev) - (1 - f) * log(1 - f)) / d)
        }

    } else if (method == "george") {
        return( (computeZ(f, f_prev, "pearson") - computeZ(f, f_prev, "fisher_mean") ) / 2)

    } else if (method == "edgington") {
        return((f + f_prev)/2)

    } else if (method == "stouffer") { 
        if (d < eps) {
            return( pnorm((qnorm(f)+qnorm(f_prev))/2)/sqrt(2*pi) ) #use midP to avoid division by d.  
        } else {
            return( (-dnorm(qnorm(f)) + dnorm(qnorm(f_prev)))/(f - f_prev) ) #This formula allows f_prev = 0. It is same as (-exp(-qnorm(f)^2/2)+exp(-qnorm(f_prev)^2/2))/(f-f_prev)/sqrt(2*pi) in Gonzalo's code.
        }
    } else {
        stop ("Invalid method. Choose one of 'fisher_mean', 'fisher_median', 'pearson', 'george', 'stouffer', or 'edgington'.")
    }
}


#' @title Compute Moments of Adjusted Z Statistic
#' @description This function calculates the moments (mean and variance) of the adjusted Z statistic based on a given p-value support vector and a specified combination method.
#' @param p_support A numeric vector of p-values. Its elements must be nonnegative, nondecreasing, and not larger than 1.
#' @param method The combination method that the adjusted Z statistic is related to. One of "fisher_mean" (default), "fisher_median", "pearson", "george", "stouffer", or "edgington".
#' @return A list containing the mean (`Zmean`) and variance (`Zvar`) of the adjusted Z statistic.
#' @details
#' This function relies on the `computeZ` function to calculate the adjusted Z statistic for each element in `p_support`. The mean and variance are then computed based on these Z statistic values and their corresponding probabilities.
#' @examples
#' methods <- c("fisher_mean", "fisher_median", "pearson", "george", "stouffer", "edgington")
#' 
#' # Toy example
#' p_support <- c(0.1, 0.2, 0.5, 1)
#' sapply(methods, function(m) computeZmoment(p_support, m))
#' 
#' # Example for p_support containing many 0's and 1's due to numerical limitations.
#' p_support <- pbinom(0:100000, size = 100000, prob = 0.7)
#' sapply(methods, function(m) computeZmoment(p_support, m))
#' @export
computeZmoment  <- function(p_support, method="fisher_mean") {
    # Validate all elements in p_support are strictly larger than 0, increasingly ordered, and contain 1 (adjusted for numerical stability)
    if (any(p_support < -1e-8) || any(diff(p_support) < -1e-8) || any(p_support > 1 + 1e-8)) {
        stop("Invalid p_support vector. Its elements must be nonnegative, nondecreasing, and not larger than 1.")
    }
    
    # Compute mean and variance of the adjusted Z statistic
    p_support_prev <- c(0, p_support[1:(length(p_support) - 1)]) # Previous p-value for each p-value in p_support
    Zvals <- mapply(computeZ, f = p_support, f_prev = p_support_prev, method = method)
    probs <- p_support - p_support_prev #Probability of each p-value in p_support

    #Remove 0 probabilities or Inf Z values, for numerical stability (0*Inf = NaN in R)
    validIdx = which(probs > 0 & Zvals < Inf)
    Zvals <- Zvals[validIdx]
    probs <- probs[validIdx]
    Zmean <- sum(probs * Zvals)
    Zvar <- sum(probs * (Zvals - Zmean)^2)
    
    return(list(Zmean = Zmean, Zvar = Zvar))
}


#' @title Type I Error Accuracy Metrics for Method Selection
#' @description Computes the two metrics recommended for comparing and selecting discrete p-value combination methods: the variance ratio \eqn{\mathrm{Var}(Z)/\mathrm{Var}(Y)} and the normalized (scaled) Wasserstein distance \eqn{W_2(Z,\tilde{Y})/\mathrm{SD}(Y)}. Here \eqn{Z} is the adjusted discrete statistic, \eqn{Y} is the continuous reference null, and \eqn{\tilde{Y}} is the moment-matched continuous surrogate used in the testing procedure. A higher variance ratio (closer to 1) and a smaller normalized distance both indicate more accurate finite-sample Type I error control.
#' @param p_support A valid p-value support vector: nonnegative, nondecreasing, and ending at 1, characterizing the null distribution of the discrete p-value.
#' @param method The combination method, one of "fisher_mean" (default), "fisher_median", "pearson", "george", "stouffer", or "edgington".
#' @return A list with the following elements:
#'   \item{var_Z}{The variance of the adjusted statistic \eqn{Z}.}
#'   \item{var_Y}{The variance of the continuous reference null \eqn{Y}, a known constant per method: 4 for Fisher and Pearson, \eqn{\pi^2/3} for George, 1 for Stouffer, and 1/12 for Edgington.}
#'   \item{var_ratio}{The variance ratio \eqn{\mathrm{Var}(Z)/\mathrm{Var}(Y)}, which lies in \eqn{[0,1]}.}
#'   \item{W2}{The Wasserstein-2 distance \eqn{W_2(Z,\tilde{Y})} between \eqn{Z} and the moment-matched surrogate \eqn{\tilde{Y}}.}
#'   \item{norm_dist}{The normalized distance \eqn{W_2(Z,\tilde{Y})/\mathrm{SD}(Y)}.}
#' @details
#' The surrogate \eqn{\tilde{Y}} matches the first two moments of \eqn{Z}: a gamma distribution for Fisher's and Pearson's methods, and a normal distribution for George's, Stouffer's, and Edgington's methods. The Wasserstein distance is computed under the optimal (quantile) coupling,
#' \deqn{W_2^2(Z,\tilde{Y}) = \sum_i \int_{\tilde{G}^{-1}(P(Z<z_i))}^{\tilde{G}^{-1}(P(Z\le z_i))} (z_i-y)^2 \tilde{g}(y)\,dy,}
#' where \eqn{\tilde{G}} and \eqn{\tilde{g}} are the CDF and density of \eqn{\tilde{Y}}. The integral over each atom is evaluated numerically. For the non-i.i.d. case, the average variance ratio \eqn{\sum_j \mathrm{Var}(Z_j)/(n\,\mathrm{Var}(Y))} can be obtained by averaging \code{var_ratio} across the component supports.
#' @references
#' 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
#' # Distribution P_L (large mass at small p-values), cf. Table 4 of Contador and Wu (2026).
#' p_support <- c(0.40, (41:100) / 100)
#' methods <- c("fisher_mean", "pearson", "stouffer", "edgington", "george")
#' round(sapply(methods, function(m) unlist(accuracy_metrics(p_support, m))), 3)
#' @importFrom stats integrate qgamma dgamma qnorm dnorm
#' @export
accuracy_metrics <- function(p_support, method = "fisher_mean") {
  # Validate p_support (same checks as computeZmoment)
  if (any(p_support < -1e-8) || any(diff(p_support) < -1e-8) || any(p_support > 1 + 1e-8)) {
    stop("Invalid p_support vector. Its elements must be nonnegative, nondecreasing, and not larger than 1.")
  }

  # Variance of the continuous reference null Y = G^{-1}(U)
  var_Y <- switch(method,
                  fisher_mean   = 4,
                  fisher_median = 4,
                  pearson       = 4,
                  george        = pi^2 / 3,
                  stouffer      = 1,
                  edgington     = 1 / 12,
                  stop("Invalid method. Choose one of 'fisher_mean', 'fisher_median', 'pearson', 'george', 'stouffer', or 'edgington'."))

  # Mean and variance of the adjusted statistic Z (canonical, from computeZmoment)
  mts   <- computeZmoment(p_support, method = method)
  meanZ <- mts$Zmean
  varZ  <- mts$Zvar

  # Degenerate case: Z is (essentially) a point mass, so the surrogate matches it exactly.
  if (varZ <= .Machine$double.eps) {
    return(list(var_Z = varZ, var_Y = var_Y, var_ratio = varZ / var_Y,
                W2 = 0, norm_dist = 0))
  }

  # Reconstruct the atoms of Z and their probability masses for the coupling
  p_prev <- c(0, p_support[-length(p_support)])
  z      <- mapply(computeZ, f = p_support, f_prev = p_prev, method = method)
  probs  <- p_support - p_prev
  ok     <- which(probs > 0 & is.finite(z))
  z      <- z[ok]
  probs  <- probs[ok] / sum(probs[ok])

  # Moment-matched continuous surrogate Ytilde (mean = meanZ, var = varZ):
  # gamma for Fisher/Pearson, normal for George/Stouffer/Edgington.
  if (method %in% c("fisher_mean", "fisher_median", "pearson")) {
    shape <- meanZ^2 / varZ
    scale <- varZ / meanZ
    Qfun  <- function(u) qgamma(u, shape = shape, scale = scale)
    dfun  <- function(y) dgamma(y, shape = shape, scale = scale)
  } else {
    sdZ  <- sqrt(varZ)
    Qfun <- function(u) qnorm(u, mean = meanZ, sd = sdZ)
    dfun <- function(y) dnorm(y, mean = meanZ, sd = sdZ)
  }

  # W2^2(Z, Ytilde) under the optimal (quantile) coupling, integrated atom by atom.
  ord     <- order(z)
  zs      <- z[ord]
  cumv    <- cumsum(probs[ord])
  cumprev <- c(0, cumv[-length(cumv)])
  cumv[length(cumv)] <- 1
  W2sq <- 0
  for (i in seq_along(zs)) {
    lo <- Qfun(cumprev[i])
    hi <- Qfun(cumv[i])
    if (!is.finite(lo)) lo <- -Inf
    if (!is.finite(hi)) hi <-  Inf
    if (hi <= lo) next
    W2sq <- W2sq + integrate(function(y) (zs[i] - y)^2 * dfun(y),
                             lower = lo, upper = hi,
                             rel.tol = 1e-8, subdivisions = 1000L)$value
  }
  W2 <- sqrt(max(W2sq, 0))

  list(var_Z     = varZ,
       var_Y     = var_Y,
       var_ratio = varZ / var_Y,
       W2        = W2,
       norm_dist = W2 / sqrt(var_Y))
}