R/BallDivergence.R

In DataSimilarity: Quantifying Similarity of Datasets and Multivariate Two- And k-Sample Testing

################################################################################
##                              BALL DIVERGENCE                               ##
##                                                                            ##
################################################################################
BallDivergence <- function (X1, X2,..., n.perm = 0, seed = NULL, num.threads = 0, 
                            kbd.type = "sum", 
                            weight = c("constant", "variance"), 
                            args.bd.test = NULL) {
  if(!requireNamespace("Ball", quietly = TRUE)) {
    stop("Package \"Ball\" required for using method BallDivergence().")
  }
  if(!is.null(seed)) {
    set.seed(seed)
  }
  genv <- globalenv()
  # Make sure to leave '.Random.seed' as-is on exit
  old_seed <- genv$.Random.seed
  on.exit(suspendInterrupts({
    if (is.null(old_seed)) {
      rm(".Random.seed", envir = genv, inherits = FALSE)
    } else {
      assign(".Random.seed", value = old_seed, envir = genv, inherits = FALSE)
    }
  }))
  data.list <- c(list(X1, X2), list(...))
  K <- length(data.list)
  if(any(!sapply(data.list, function(x) inherits(x, "matrix") | inherits(x, "data.frame")))) {
    stop("All datasets must be provided as data.frames or matrices.")
  }
  p <- sapply(data.list, ncol)
  if(length(unique(p)) > 1) { 
    stop("All datasets must have the same number of variables")
  }
  if(K == 2) {
    dname <- c(deparse1(substitute(X1)), deparse1(substitute(X2)))
  } else {
    mc <- as.list(match.call())
    mc <- mc[!names(mc) %in% c("n.perm", "seed", "num.threads", "kbd.type", 
                               "weight", "args.bd.test")]
    dname <- sapply(mc[-1], deparse)
  }
  NAs <- sapply(data.list, anyNA)
  if(sum(NAs) > 1) {
    warning(paste0("Missing values in ", paste0(dname[NAs], collapse = ", "), 
                   " were omitted from the calculation."))
    data.list <- lapply(data.list, na.omit)
  }
  n.vec <- sapply(data.list, nrow)
  data.list <- lapply(data.list, function(X) {
    colnames(X) <- paste0("X", 1:p[1])
    X
  })
  res <- do.call(Ball::bd.test, c(list(x = data.list, num.permutations = max(n.perm, 1),
                                       method = ifelse(n.perm > 0, "permutation", "limit"), 
                                       distance = FALSE, size = NULL, seed = seed, 
                                       num.threads = num.threads, kbd.type = kbd.type, 
                                       weight = weight), args.bd.test))#
  
  dname <- paste0(dname, collapse = ifelse(length(data.list) > 2, ", ", " and "))
  if(K > 2 & n.perm <= 0) {
    res <- list(statistic = res[paste0("kbd.", kbd.type, ".", weight[1])],
                p.value = NULL, replicates = n.perm, size = n.vec, 
                alternative = NA, method = paste0(K, "-sample Ball Divergence Test (", 
                                                  ifelse(n.perm > 0, "Permutation", "Limit"), 
                                                  " Distribution)"), 
                data.name = paste0(dname, 
                                   "\nnumber of observations = ", sum(n.vec), 
                                   ", group sizes: ", paste0(n.vec, collapse = " "), 
                                   "\nreplicates = ", n.perm, ", weight: ", 
                                   match.arg(weight)))
  } else {
    res$data.name <- paste0(dname, res$data.name)
  }
  
  names(res)[3] <- "n.perm"
  res <- res[names(res) != "complete.info"]
  res$alternative <- ifelse(K > 2, "At least one pair of distributions are unequal.", 
                            paste0("The distributions of ", dname, " are unequal."))
  class(res) <- "htest"
  return(res)
}