Nothing
#' 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 or a symmetric matrix; 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}
#' @param ncores Numeric.
#' Number of cores to use in computing results.
#' Defaults to \code{parallel::detectCores() / 2} or half of your
#' computer's processing power.
#' Set to \code{1} to not use parallel computing
#' @param seed Numeric.
#' Set seed for reproducible results.
#' Defaults to \code{NULL}
#'
#' @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}}
#' \item{Settings}{user-specified settings}
#' }
#'
#' @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
#' ## Example with fictitious data
#'
#' # 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{
#' results <- cpPermuteEntropy(W = W, cpThreshold.object = cpThreshold.object,
#' n = 100, interval = 0.95, ncores = 1, seed = 4186)
#'
#' # check results
#' results
#' }
#'
#' ## Example with Obama data set (see ?Obama)
#'
#' # get data
#' data(Obama)
#'
#' # estimate network
#' net <- qgraph::EBICglasso(qgraph::cor_auto(Obama), n = nrow(Obama))
#'
#' # create cpThreshold object
#' \donttest{
#' threshold <- cpThreshold(net, method = "weighted",
#' k.range = 3:4,
#' I.range = seq(0.1, 0.5, 0.01),
#' threshold = "entropy")
#' }
#'
#' # run cpPermuteEntropy with 50 permutations and 99% confidence interval
#' \donttest{
#' permute <- cpPermuteEntropy(net, cpThreshold.object = threshold,
#' interval = 0.99, n = 50, ncores = 1, seed = 4186)
#'
#' # check results
#' permute
#' }
#'
#' @author Jens Lange, \email{lange.jens@@outlook.com}
#'
#' @export cpPermuteEntropy
cpPermuteEntropy <- function(W, cpThreshold.object, n = 100, interval = 0.95,
CFinder = FALSE, ncores, seed = NULL) {
###check whether W is a qgraph object
###if not check whether matrix is symmetric and convert to qgraph object
if (!isTRUE(methods::is(W, "qgraph"))) {
#error if W is a matrix but not symmetric
if (isSymmetric(W) == FALSE) {
stop("If W is a matrix, it must be symmetric.")
}
#converts matrix to qgraph if input is not a qgraph object
W <- qgraph::qgraph(W, DoNotPlot = TRUE)
#### OLD ERROR BEGIN ####
#stop("W (network object) must be a qgraph object.")
#### OLD ERROR END ####
}
#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.")
}
#set default for missing processing cores
if(missing(ncores)){
ncores <- ceiling(parallel::detectCores() / 2)
}
#for n random permutations of graph, the function runs cpThreshold 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)
#initialize iterations
iterations <- as.list(1:n)
#parallel processing
cl <- parallel::makeCluster(ncores)
#set seed for reproducibility
if(!is.null(seed)){
parallel::clusterSetRNGStream(cl = cl, iseed = seed)
}
#export variables
parallel::clusterExport(cl = cl,
varlist = c("W", "k.range", "I.range",
"method", "cpThreshold"),
envir=environment())
#let user know permutations has started
message("Permutating...\n", appendLF = FALSE)
#Permutate
perms <- pbapply::pblapply(
X = iterations, cl = cl,
FUN = function(X, W, k.range, I.range, method){
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
max_ent <- matrix(nrow = length(k.range), ncol = 2)
colnames(max_ent) <- c("k", "max")
for (j in 1:length(k.range)) {
max_ent[j,] <- c(k.range[j], max(results_cpThreshold[results_cpThreshold$k == k.range[j], ]$Entropy.Threshold))
}
return(max_ent)
},
W = W, k.range = k.range, I.range = I.range, method = method
)
parallel::stopCluster(cl)
max_ent <- as.data.frame(do.call(rbind, perms))
#### OLD CODE BEGIN ####
#max_ent <- data.frame(k = numeric(),
# max = numeric())
#counter <- 1
#progress_bar_permute <- utils::txtProgressBar(min = 0, max = n,
# 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)
#}
#close progress bar
#close(progress_bar_permute)
#### OLD CODE END ####
#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)
returned_object <- list(Confidence.Interval = results,
Extracted.Rows = extracted_rows,
Settings = list(n = n,
interval = interval,
CFinder = CFinder,
ncores = ncores,
seed = seed))
class(returned_object) <- "cpPermuteEntropy"
return(returned_object)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.