R/telescoped.gtest.R

In Upsilon: Another Test of Association for Count Data

# Author: Xuye Luo
# Date: December 11, 2025

#' @title Find Best G-test Result in Telescoped Data
#'
#' @description Performs G-tests (Likelihood Ratio Tests) on telescoped data at multiple 
#' resolution levels (or a specific focus level) to identify the strongest association.
#' The "best" result is selected based on either the most significant p-value or 
#' the largest effect size (Mutual Information), depending on the \code{esize} parameter.
#'
#' @param x A numeric vector without NA values, same length as y.
#' @param y A numeric vector without NA values, same length as x.
#' @param focus An integer specifying a specific resolution level to test. 
#'              If \code{NULL} (default), tests are performed across all available levels.
#' @param base Integer. The base for exponential scaling in telescoping. Defaults to 2.
#' @param log.p Logical. If \code{TRUE}, returns the natural logarithm of the p-value. Defaults to \code{TRUE}.
#' @param esize Logical. If \code{TRUE}, sorts results primarily by Effect Size (Mutual Information) in descending order. 
#'              If \code{FALSE} (default), sorts primarily by P-value (significance) in descending order (most significant first).
#'
#' @return A data frame containing the best result with columns:
#' \item{p.value}{The p-value (or log p-value) of the G-test.}
#' \item{estimate}{The effect size (Mutual Information).}
#' \item{focus}{The resolution level (focus) where this result was found.}
#'
#'
#' @examples
#' library("Upsilon")
#' set.seed(123)
#' n <- 100
#' x <- rnorm(n)
#' y <- rnorm(n)
#' 
#' # Best result by significance
#' telescoped.gtest(x, y)
#' 
#' @keywords internal
#' export
#' @noRd
telescoped.gtest <- function(x, 
                             y, 
                             focus = NULL, 
                             base = 2, 
                             log.p = TRUE, 
                             esize = FALSE) {
  
  # Get telescoped data
  telescope_list <- telescoped.data(x, y, focus = focus, base = base)
  
  # Search across multiple levels
  if (is.list(telescope_list) && is.null(focus)) {
    
    # Remove the first element (Raw Data) as it is continuous and not suitable for G-test directly
    if (length(telescope_list) > 1) {
      telescope_list <- telescope_list[-1]
    }
    
    # Iterate through each level
    gtest_result_list <- lapply(telescope_list, function(mat) {
      # Use matrix indexing instead of unlist for safety
      val_x <- mat[, 1]
      val_y <- mat[, 2]
      
      results <- fast.gtest(val_x, val_y, log.p = log.p)
      
      data.frame(
        p.value = results$p.value, 
        estimate = results$estimate
      )
    })
    
    gtest_result <- do.call(rbind, gtest_result_list)
    
    # Assign focus levels (0 to N-1 based on list index)
    gtest_result$focus <- seq_len(nrow(gtest_result)) - 1
    
    # Sorting Logic
    if (esize) {

      ord <- order(
        gtest_result$estimate, 
        -gtest_result$p.value, 
        -gtest_result$focus, 
        decreasing = TRUE
      )
    } else {

      ord <- order(
        -gtest_result$p.value, 
        -gtest_result$focus, 
        decreasing = TRUE
      )
    }
    
    gtest_result_sorted <- gtest_result[ord, , drop = FALSE]
    gtest_result_best <- gtest_result_sorted[1, , drop = FALSE]
    
  } else {
    # Specific focus level
    mat <- if (is.list(telescope_list)) telescope_list[[1]] else telescope_list
    
    results <- fast.gtest(mat[, 1], mat[, 2], log.p = log.p)
    
    gtest_result_best <- data.frame(
      p.value = results$p.value,
      estimate = results$estimate,
      focus = ifelse(is.null(focus), -1, focus)
    )
  }
  
  return(gtest_result_best)
}