R/cpPermuteEntropy.R

In chr1swallace/FasterCliquePercolation: Clique Percolation for Networks

#' Confidence Intervals Of Entropy Based On Random Permutations Of Network
#' 
#' Function for determining confidence intervals of entropy values calculated for
#' community partition from clique percolation based on randomly permuted networks
#' of original network.
#' 
#' @param W A qgraph object; see also \link[qgraph]{qgraph}
#' @param cpThreshold.object A cpThreshold object; see also \link{cpThreshold}
#' @param n number of permutations (default is 100)
#' @param interval requested confidence interval (larger than zero and smaller 1;
#'    default is 0.95)
#' @param CFinder logical indicating whether clique percolation for weighted networks
#'    should be performed as in CFinder ; see also \link{cpAlgorithm}
#'    
#' @return
#'   A list object with the following elements:
#'   \describe{
#'   \item{Confidence.Interval}{a data frame with lower and upper bound of
#'      confidence interval for each \code{k}}
#'   \item{Extracted.Rows}{rows extracted from \code{cpThreshold.object} that are
#'      larger than the upper bound of the specified confidence interval for
#'      each \code{k}}
#' }
#' 
#' @details
#'   The function generates \code{n} random permutations of the network
#'   specified in \code{W}. For each randomly permuted network, it runs \code{cpThreshold}
#'   (see \link{cpThreshold} for more information) with \code{k} and \code{I} values
#'   extracted from the cpThreshold object specified in \code{cpThreshold.object}.
#'   Across permutations, the confidence intervals of the entropy values are determined
#'   for each \code{k} separately.
#' 
#'   The confidence interval of the entropy values is determined separately for each \code{k}.
#'   This is because larger \code{k} have to produce less communities on average,
#'   which will decrease entropy. Comparing confidence intervals of smaller \code{k} to
#'   those of larger \code{k} would therefore be disadvantageous for larger \code{k}.
#' 
#'   In the output, one can check the confidence intervals of each \code{k}. Moreover,
#'   a data frame is produced that takes the cpThreshold object that was specified in
#'   \code{cpThreshold.object} and removes all rows that do not exceed the upper bound of the
#'   confidence interval of the respective \code{k}.
#' 
#' @examples
#' # create qgraph object
#' W <- matrix(c(0,1,1,1,0,0,0,0,
#'               0,0,1,1,0,0,0,0,
#'               0,0,0,0,0,0,0,0,
#'               0,0,0,0,1,1,1,0,
#'               0,0,0,0,0,1,1,0,
#'               0,0,0,0,0,0,1,0,
#'               0,0,0,0,0,0,0,1,
#'               0,0,0,0,0,0,0,0), nrow = 8, ncol = 8, byrow = TRUE)
#' W <- Matrix::forceSymmetric(W)
#' W <- qgraph::qgraph(W)
#' 
#' # create cpThreshold object
#' cpThreshold.object <- cpThreshold(W = W, method = "unweighted", k.range = c(3,4),
#'                                   threshold = "entropy")
#' 
#' # run cpPermuteEntropy with 100 permutations and 95% confidence interval
#' \donttest{
#' set.seed(4186)
#' results <- cpPermuteEntropy(W = W, cpThreshold.object = cpThreshold.object,
#'                             n = 100, interval = 0.95)
#' 
#' # check results
#' results$Confidence.Interval
#' results$Extracted.Rows
#' }
#' 
#' @author Jens Lange, \email{lange.jens@@outlook.com} 
#' 
#' @export cpPermuteEntropy

cpPermuteEntropy <- function(W, cpThreshold.object, n = 100, interval = 0.95,
                             CFinder = FALSE) {
  
  #error message if n is not larger than zero
  if (n <= 0) {
    stop("n must be larger than zero.")
  }
  
  #error message if interval is not larger than zero and smaller than 1
  if (!(interval < 1 & interval > 0)) {
    stop("interval must be larger than zero and smaller than 1.")
  }
  
  #for n random permutations of graph, the function runs cpTthreshold for k.range and I.range
  #random permutations generated by permuting the upper triangle of the weights matrix and forcing it symmetric
  #for each random permutation, the maximum entropy is extracted for each k and saved in data frame
  #information about parameters in original cpThreshold function is extracted from cpThreshold.object
  if (CFinder == FALSE) {
    if ("Intensity"%in%names(cpThreshold.object)) {method <- "weighted"} else {method <- "unweighted"}
  }
  if (CFinder == TRUE) {
    method <- "weighted.CFinder"
  }
  k.range <- sort(unique(cpThreshold.object$k))
  I.range <- sort(unique(cpThreshold.object$Intensity), decreasing = TRUE)
  max_ent <- data.frame(k = numeric(),
                        max = numeric())
  counter <- 1
  
  progress_bar_permute <- utils::txtProgressBar(min = 0, max = 100,
                                                style = 3) #progress bar definition
  progress_bar_permute_counter <- 0 #counter for progress bar
  
  for (i in 1:n) {
    Wmat <- qgraph::getWmat(W)
    Wmat[upper.tri(Wmat)] <- sample(Wmat[upper.tri(Wmat)])
    Wmat <- Matrix::forceSymmetric(Wmat, "U")
    graph <- qgraph::qgraph(Wmat, DoNotPlot = TRUE)
    
    utils::capture.output(
      results_cpThreshold <- cpThreshold(graph, method = method, k.range = k.range,
                                          I.range = I.range, threshold = "entropy")
    ) #capture.output suppresses display of progress bar from cpThreshold
    
    for (j in k.range) {
      max_ent[counter, "k"] <- j
      max_ent[counter, "max"] <- max(results_cpThreshold[results_cpThreshold$k == j, ]$Entropy.Threshold)
      counter <- counter + 1
    }
    
    #progress bar update
    progress_bar_permute_counter <- progress_bar_permute_counter + 1
    utils::setTxtProgressBar(progress_bar_permute, progress_bar_permute_counter)
  }
  
  #determine confidence interval based on normal distribution of extracted entropies for each k
  #extract rows of cpThreshold.object with entropy larger than upper confidence for each k
  results <- data.frame(k = numeric(),
                        l = numeric(),
                        u = numeric())
  extracted_rows <- list()
  counter2 <- 1
  for (k in k.range) {
    m <- mean(max_ent[max_ent$k == k, "max"])
    std <- stats::sd(max_ent[max_ent$k == k, "max"])
    lower_ci <- m - stats::qnorm((1 + interval)/2) * std/sqrt(n)
    upper_ci <- m + stats::qnorm((1 + interval)/2) * std/sqrt(n)
    results[counter2, "k"] <- k
    results[counter2, "l"] <- lower_ci
    results[counter2, "u"] <- upper_ci
    
    rows <- cpThreshold.object[cpThreshold.object$k == k & cpThreshold.object$Entropy.Threshold > upper_ci, ]
    extracted_rows[[counter2]] <- rows
    
    counter2 <- counter2 + 1
  }
  names(results) <- c("k",
                      paste(interval*100, "% CI lower", sep = ""),
                      paste(interval*100, "% CI upper", sep = ""))
  extracted_rows <- Reduce(rbind, extracted_rows)
  
  return(list(Confidence.Interval = results, Extracted.Rows = extracted_rows))
}