R/hypervolume_n_occupancy_test.R

In hypervolume: High Dimensional Geometry, Set Operations, Projection, and Inference Using Kernel Density Estimation, Support Vector Machines, and Convex Hulls

hypervolume_n_occupancy_test <- function(observed, path, alternative = "two_sided", significance = 0.05, cores = 1, p_adjust = "none", multi_comp_type = "pairwise") {
  
  # intitialize multi core calculations
  
  exists_cluster = TRUE
  if(cores > 1 & getDoParWorkers() == 1) {
    # If no cluster is registered, create a new one based on use input
    cl = makeCluster(cores)
    clusterEvalQ(cl, {
      library(hypervolume)
    })
    registerDoParallel(cl)
    exists_cluster = FALSE
  }
  
  on.exit({
    # If a cluster was created for this specific function call, close cluster and register sequential backend
    if(!exists_cluster) {
      stopCluster(cl)
      registerDoSEQ()
    }
  }, add = TRUE)
  
  
  # check if the first element of the first pairwise combination is called permutation1.rds
  # n <- list all the files in the first pa

  if(list.files(file.path(path, list.files(path)[1]))[1] != "permutation1.rds" ) {
    if(!exists_cluster) {
      stopCluster(cl)
      registerDoSEQ()
    }
    stop("Invalid input")
  }


  # name of the groups of the observed HypervolumeList
  group_names <- unlist(lapply(observed@HVList, function(x)x@Name))
  
  # combination of these names
  observed_combn <- combn(group_names, 2)
  
  # number of the pairwise combination inferred from the parent directory
  groups <- length(list.files(path))
  
  # check if the number of columns of observed_combn is the same as the one in groups
  if(ncol(observed_combn) != groups){
    stop("The number of pairwise combination does not match between observed and path.")
  }
  
  # store the names of the permuted files. Names are the same for each pairwise combination
  n <- list.files(file.path(path, list.files(path)[1]))
  
  # initialize an empty list for storing the pairwise differences between
  # ValueAtRandomPoints of the observed argument
  observed_res <- vector(ncol(observed_combn), mode = "list")
  
  # calculate pairwise differences
  for(i in 1:ncol(observed_combn)){
    group_1 <- which(group_names == observed_combn[1, i])
    group_2 <- which(group_names == observed_combn[2, i])
    result <- lapply(observed[[c(group_1, group_2)]]@HVList, function(x) x@ValueAtRandomPoints)
    result <- do.call(cbind, result)
    observed_res[[i]] <- result[ , 1] - result [ , 2]
  }
  
  

  # calculate the difference in the resampled hypervolumes
  # perform a nested loop with foreach and store the results
  # this loop will return a list of length equal to the number of pairwise combinations
  # each list contains a matrix, where the rose are the RandomPoints and columns represent
  # the difference between groups
  
  simu_res <- foreach(i = list.files(path)) %:% foreach(j = n, .combine = cbind) %dopar% {
    h1 = readRDS(file.path(path, i, j))
    result <- lapply(h1@HVList, function(x) x@ValueAtRandomPoints)
    result <- do.call(cbind, result)
    simulated_res <- result[ , 1] - result [ , 2]
  }
  
  
  # initialize a list to store the results
  result_list <- vector(groups, mode = "list")
  
  # number of permutations
  n <- length(list.files(file.path(path, list.files(path)[1])))
  
  # calculate the probability that observed differences are significantly different
  # from simulated differences
  # retain only ValueAtRandomPoints for which the difference is significant
  # se the other ValueAtRandomPoints to 0
  
  for(i in 1:groups){
    # simulated results for the i-th combination
    distribution <- simu_res[[i]]
    
    # observed result for the i-th combination
    obs <- observed_res[[i]]
    
    # calculates probabilities: can be more, less or two_sided
    # credits to vegan oecosimu
    if(alternative == "more"){
      result <- sweep(distribution, 1, obs, ">=")
      result <- apply(result, 1, sum)
      p <- pmin(1, (result  + 1)/(n + 1))
    }
    
    if(alternative == "less"){
      result <- sweep(distribution, 1, obs, "<=")
      result <- apply(result, 1, sum)
      p <- pmin(1, (result + 1)/(n + 1))
    }
    
    # if(alternative == "more_less"){
    #   result_less <- sweep(distribution, 1, obs, ">=")
    #   result_more <- sweep(distribution, 1, obs, "<=")
    #   result_less <- apply(result_less, 1, sum)
    #   result_more <- apply(result_more, 1, sum)
    #   p_less <- pmin(1, (result_less  + 1)/(n + 1))
    #   p_more <- pmin(1, (result_more  + 1)/(n + 1))
    #   p <- apply(data.frame(p_less, p_more), 1, max)
    # }
    # 
    
    if(alternative == "two_sided"){
      # result <- sweep(abs(distribution), 1, abs(obs), ">=")
      # result <- apply(result, 1, sum)
      result_less <- sweep(distribution, 1, obs, "<=")
      result_more <- sweep(distribution, 1, obs, ">=")
      result_less <- apply(result_less, 1, sum)
      result_more <- apply(result_more, 1, sum)
      result <- 2*pmin(result_less, result_more)
      p <- pmin(1, (result + 1)/(n + 1))
    }
    
    if(! identical(p_adjust, "none")){
      if(identical(multi_comp_type, "pairwise")){
        p <- p.adjust(p, method = p_adjust, n = length(p))
      }
      
      if(identical(multi_comp_type, "all")){
        p <- p.adjust(p, method = p_adjust, n = length(p) * groups)
      } 
      
    }
    
    # retain only those probabilities that are greater than significance
    p_obs <- p <= significance
    

    
    # retain only ValueAtRandomPoints for which the difference is significant
    if(sum(p_obs) == 0){
      result_list[[i]] <- rep(0, length(obs))
    } else {
      obs[! p_obs] <- 0
      result_list[[i]] <- obs
    }
  }
  

  # labels of the pairwise combinations
  label_res <- unlist(apply(observed_combn, 2, function(x) paste(x[1], x[2], sep = "__")))
  
  # cbind the results of pairwise combinations
  result_combn <- do.call(cbind, result_list)
  
  # assign colnames
  colnames(result_combn) <- label_res

  
  # create an empty list to store the hypervolume results 
  hv_list_res <- vector(ncol(observed_combn), mode = "list")
  
  # create an empty hypervolume as a model for storing results
  empty_hypervolume <- new("Hypervolume")
  
  # get mindensity from observed, it is the same ofor all the hypervolumes
  mindensity <- observed@HVList[[1]]@PointDensity
  
  
  vol_list <- unlist(lapply(observed@HVList, function(x) x@Volume))
  res <- lapply(observed@HVList, function(x) x@ValueAtRandomPoints)
  res <- do.call(cbind, res)
  res[res > 0] <- 1
  
  # volume of the significant fraction for th groups under comparison
  intersection_weights <- sweep(res, 1, apply(res, 1, sum), "/")
  vol_tot <- sum(apply(intersection_weights, 2, function(x) mean(x[x > 0], na.rm = TRUE)) * vol_list, na.rm = TRUE)

  
  # dimnames
  cn <- dimnames(observed@HVList[[1]]@RandomPoints)[[2]]
  dn <- list(NULL, cn)
  
  # create the hypervolumes to return to the user
  for(i in 1:ncol(observed_combn)){
    
    # position of the groups under comparison in observed_combn
    group_1 <- which(group_names == observed_combn[1, i])
    group_2 <- which(group_names == observed_combn[2, i])
    
    # vol_list <- c(observed@HVList[group_names == observed_combn[1, i]][[1]]@Volume, 
    #               observed@HVList[group_names == observed_combn[2, i]][[1]]@Volume)
    # 

    result <- lapply(observed[[c(group_1, group_2)]]@HVList, function(x) x@ValueAtRandomPoints)
    result <- do.call(cbind, result)
    # result <- result[colSums(result) != 0, ]
    result[result != 0] <- 1
    
    
    res_vol <- sum(result_combn[, i] != 0) / nrow(res) * vol_tot
    
    # get Data from observed, merge data for the two groups under comparison
    Data <- lapply(observed[[c(group_1, group_2)]]@HVList, function(x) x@Data)
    Data <- unique(do.call(rbind, Data))
    
    # volume of the significant fraction for th groups under comparison
    # res_vol <- sum(rowSums(result_combn[, i, drop = FALSE]) != 0)
    
    # ValueAtRandomPoints, obtained from result_combn
    empty_hypervolume@ValueAtRandomPoints <- result_combn[, i]
    
    # the name of the method is now n_occupancy_test
    empty_hypervolume@Method <- "n_occupancy_test"
    
    # label of the two groups under comparison
    empty_hypervolume@Name <- label_res[i]
    
    # assign merged Data
    empty_hypervolume@Data <- Data
    
    # number of dimensions
    empty_hypervolume@Dimensionality <- ncol(observed@HVList[[1]]@RandomPoints)
    
    # empty parameters
    empty_hypervolume@Parameters = list()
    
    # The coordinates of RandomPoints are the same across all the hypervolumes
    empty_hypervolume@RandomPoints = observed@HVList[[1]]@RandomPoints
    
    # assign the volume
    empty_hypervolume@Volume <- res_vol
    
    # meandensity is the same across all the hypervolumes 
    empty_hypervolume@PointDensity <- sum(result_combn[, i] != 0) / res_vol
    
    # set dimnames
    dimnames(empty_hypervolume@RandomPoints) = dn
    dimnames(empty_hypervolume@Data) = dn
    
    hv_list_res[[i]] <- empty_hypervolume
    
    
  }
  
  
  # return the results, if only 1 comparison is present return an object of class Hypervolume
  # otherwise returns a HypervolumeList
  if(ncol(observed_combn) == 1){
    result <- empty_hypervolume
  } else {
    hv_list_res <- new("HypervolumeList", HVList = hv_list_res)
    result <- hv_list_res
  }

  return(result)

}