From b588ea3a5e50e0167b9d284431b33e8f683de537 Mon Sep 17 00:00:00 2001 From: Michael Pinus Date: Sun, 3 Apr 2022 10:00:32 +0300 Subject: [PATCH 1/5] initial support for multicore processing --- DESCRIPTION | 9 +- NAMESPACE | 1 + R/NCT.R | 366 +++++++++++++++++++++++++---------- man/networkcomparisontest.rd | 5 +- 4 files changed, 277 insertions(+), 104 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 08f8f42..fe61c9d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -10,6 +10,8 @@ Authors@R: c( person("Alex", "Millner", role = c("ctb"))) Maintainer: Claudia van Borkulo Imports: + snow, + foreach, qgraph, IsingFit, IsingSampler, @@ -18,7 +20,10 @@ Imports: graphics, utils, Matrix, - methods + methods, + bigstatsr, + doParallel, + parallelly Suggests: networktools Description: This permutation based hypothesis test, suited for Gaussian and binary data, @@ -30,4 +35,4 @@ Description: This permutation based hypothesis test, suited for Gaussian and bin one group which is measured twice. See van Borkulo et al. (2017) . License: GPL-2 -RoxygenNote: 6.1.1 +RoxygenNote: 7.1.2 diff --git a/NAMESPACE b/NAMESPACE index b61f291..e33bcbb 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,6 +7,7 @@ import(qgraph) import(IsingFit) import(IsingSampler) import(reshape2) +import(foreach) importFrom("graphics", "hist", "points") importFrom("stats", "cor", "p.adjust") importFrom("utils", "setTxtProgressBar", "txtProgressBar") diff --git a/R/NCT.R b/R/NCT.R index 6383a3c..14cbfc3 100644 --- a/R/NCT.R +++ b/R/NCT.R @@ -9,7 +9,8 @@ NCT <- function(data1, data2, centrality=c("strength","expectedInfluence"),nodes="all", communities=NULL,useCommunities="all", estimator, estimatorArgs = list(), - verbose = TRUE){ + verbose = TRUE, + nCores = 1){ # store function call, including default arguments not explicitly set - credit to Neal Fultz match.call.defaults <- function(...) { @@ -148,7 +149,7 @@ NCT <- function(data1, data2, } - if (progressbar==TRUE) pb <- txtProgressBar(max=it, style = 3) + if (progressbar==TRUE) {pb <- txtProgressBar(max=it, style = 3)} x1 <- data1 x2 <- data2 nobs1 <- nrow(x1) @@ -171,11 +172,13 @@ NCT <- function(data1, data2, } } - glstrinv.perm <- glstrinv.real <- nwinv.real <- nwinv.perm <- c() - diffedges.perm <- matrix(0,it,nedges) - einv.perm.all <- array(NA,dim=c(nvars, nvars, it)) - corrpvals.all <- matrix(NA,nvars,nvars) - edges.pvalmattemp <- matrix(0,nvars,nvars) + glstrinv.real <- nwinv.real <- c() + glstrinv.perm <- nwinv.perm <- bigstatsr::FBM(nrow = it,ncol = 1) + diffedges.perm <- bigstatsr::as_FBM(matrix(0,it,nedges)) + einv.perm.all <- bigstatsr::as_FBM(matrix(NA,nvars^2, it)) + einv.perm.all_temp <- bigstatsr::as_FBM(matrix(NA,nvars,nvars)) + corrpvals.all <- matrix(NA,nvars,nvars) + edges.pvalmattemp <- matrix(0,nvars,nvars) validCentrality <- c("closeness", "betweenness", @@ -188,7 +191,7 @@ NCT <- function(data1, data2, }else { centrality } - diffcen.perm <- matrix(NA, it, nnodes*length(centrality)) + diffcen.perm <- bigstatsr::as_FBM(matrix(NA, it, nnodes*length(centrality))) ##################################### ### procedure for all data ### @@ -255,124 +258,287 @@ NCT <- function(data1, data2, ##################################### ##### Start permutations ##### ##################################### - - for (i in 1:it) - { - diffedges.permtemp <- matrix(0, nvars, nvars) - - # If not paired data - if(paired==FALSE) + if (nCores == 1) { + for (i in 1:it) { - # s <- sample(1:(nobs1+nobs2),nobs1,replace=FALSE) - # x1perm <- dataall[s,] - # x2perm <- dataall[b[-s], ] - - # Include variance check - okay <- FALSE - counter <- 0 + diffedges.permtemp <- matrix(0, nvars, nvars) - if(binary.data) { # if binary data we need to resample the permutation to ensure mininum required variance for glmnet - while(okay == FALSE) { - - # Permute + # If not paired data + if(paired==FALSE) + { + # s <- sample(1:(nobs1+nobs2),nobs1,replace=FALSE) + # x1perm <- dataall[s,] + # x2perm <- dataall[b[-s], ] + + # Include variance check + okay <- FALSE + counter <- 0 + + if(binary.data) { # if binary data we need to resample the permutation to ensure mininum required variance for glmnet + while(okay == FALSE) { + + # Permute + s <- sample(1:(nobs1+nobs2),nobs1,replace=FALSE) + x1perm <- dataall[s,] + x2perm <- dataall[b[-s], ] + + # check glmnet requirement: at least two instances of each category + ind <- all(apply(x1perm, 2, function(x) min(c(sum(x==0), sum(x==1)))) > 1) & all(apply(x2perm, 2, function(x) min(c(sum(x==0), sum(x==1)))) > 1) + if(ind) okay <- TRUE else counter <- counter + 1 + + } # end: while + } else{ s <- sample(1:(nobs1+nobs2),nobs1,replace=FALSE) x1perm <- dataall[s,] x2perm <- dataall[b[-s], ] - - # check glmnet requirement: at least two instances of each category - ind <- all(apply(x1perm, 2, function(x) min(c(sum(x==0), sum(x==1)))) > 1) & all(apply(x2perm, 2, function(x) min(c(sum(x==0), sum(x==1)))) > 1) - if(ind) okay <- TRUE else counter <- counter + 1 - - } # end: while - } else{ - s <- sample(1:(nobs1+nobs2),nobs1,replace=FALSE) - x1perm <- dataall[s,] - x2perm <- dataall[b[-s], ] + } + + + # Estimate the networks: + r1perm <- do.call(estimator,c(list(x1perm),estimatorArgs)) + if (is.list(r1perm)) r1perm <- r1perm$graph + + r2perm <- do.call(estimator,c(list(x2perm),estimatorArgs)) + if (is.list(r2perm)) r2perm <- r2perm$graph + + if(weighted==FALSE){ + r1perm=(r1perm!=0)*1 + r2perm=(r2perm!=0)*1 + } } - - # Estimate the networks: - r1perm <- do.call(estimator,c(list(x1perm),estimatorArgs)) - if (is.list(r1perm)) r1perm <- r1perm$graph - - r2perm <- do.call(estimator,c(list(x2perm),estimatorArgs)) - if (is.list(r2perm)) r2perm <- r2perm$graph - - if(weighted==FALSE){ - r1perm=(r1perm!=0)*1 - r2perm=(r2perm!=0)*1 + # If paired data + if(paired==TRUE) + { + if (verbose) message("Note: NCT for dependent data has not been validated.") + s <- sample(c(1,2),nobs1,replace=TRUE) + x1perm <- x1[s==1,] + x1perm <- rbind(x1perm,x2[s==2,]) + x2perm <- x2[s==1,] + x2perm <- rbind(x2perm,x1[s==2,]) + + # To do: add variance check for paired data + + + # Estimate the networks: + r1perm <- do.call(estimator,c(list(x1perm),estimatorArgs)) + if (is.list(r1perm)) r1perm <- r1perm$graph + + r2perm <- do.call(estimator,c(list(x2perm),estimatorArgs)) + if (is.list(r2perm)) r2perm <- r2perm$graph + + if(weighted==FALSE){ + r1perm=(r1perm!=0)*1 + r2perm=(r2perm!=0)*1 + } } - } - - # If paired data - if(paired==TRUE) - { - if (verbose) message("Note: NCT for dependent data has not been validated.") - s <- sample(c(1,2),nobs1,replace=TRUE) - x1perm <- x1[s==1,] - x1perm <- rbind(x1perm,x2[s==2,]) - x2perm <- x2[s==1,] - x2perm <- rbind(x2perm,x1[s==2,]) - - # To do: add variance check for paired data + ## Invariance measures for permuted data + if(abs){ + glstrinv.perm[i] <- abs(sum(abs(r1perm[upper.tri(r1perm)]))-sum(abs(r2perm[upper.tri(r2perm)]))) + } else { + glstrinv.perm[i] <- abs(sum(r1perm[upper.tri(r1perm)])-sum(r2perm[upper.tri(r2perm)])) + } - # Estimate the networks: - r1perm <- do.call(estimator,c(list(x1perm),estimatorArgs)) - if (is.list(r1perm)) r1perm <- r1perm$graph + diffedges.perm[i,] <- abs(r1perm-r2perm)[upper.tri(abs(r1perm-r2perm))] + diffedges.permtemp[upper.tri(diffedges.permtemp, diag=FALSE)] <- diffedges.perm[i,] + diffedges.permtemp <- diffedges.permtemp + t(diffedges.permtemp) + einv.perm.all[,,i] <- diffedges.permtemp + nwinv.perm[i] <- max(diffedges.perm[i,]) - r2perm <- do.call(estimator,c(list(x2perm),estimatorArgs)) - if (is.list(r2perm)) r2perm <- r2perm$graph - if(weighted==FALSE){ - r1perm=(r1perm!=0)*1 - r2perm=(r2perm!=0)*1 + if(test.centrality==TRUE){ + cen1permtemp <- centrality_auto(r1perm)$node.centrality + cen2permtemp <- centrality_auto(r2perm)$node.centrality + names(cen1permtemp) <- names(cen2permtemp) <- c("betweenness","closeness","strength","expectedInfluence") + if(TRUE %in% (bridgecen %in% centrality)){ + b1permtemp <- networktools::bridge(r1perm, communities=communities, useCommunities=useCommunities) + b2permtemp <- networktools::bridge(r2perm, communities=communities, useCommunities=useCommunities) + names(b1permtemp) <- names(b2permtemp) <- c(bridgecen, "bridgeExpectedInfluence2step", + "communities") + b1permtemp$communities <- b2permtemp$communities <- NULL + cen1permtemp <- data.frame(c(cen1permtemp,b1permtemp)) + cen2permtemp <- data.frame(c(cen2permtemp,b2permtemp)) + } + diffcen.permtemp <- as.matrix(cen1permtemp) - as.matrix(cen2permtemp) + if(nodes[1]=="all"){ + diffcen.perm[i,] <- reshape2::melt(diffcen.permtemp[,centrality])$value + } else { + diffcen.perm[i,] <- reshape2::melt(diffcen.permtemp[which(nodes%in%colnames(data1)),centrality])$value + } } + + + if (progressbar==TRUE) setTxtProgressBar(pb, i) } + } else if (nCores > 1) { - ## Invariance measures for permuted data - if(abs){ - glstrinv.perm[i] <- abs(sum(abs(r1perm[upper.tri(r1perm)]))-sum(abs(r2perm[upper.tri(r2perm)]))) - } else { - glstrinv.perm[i] <- abs(sum(r1perm[upper.tri(r1perm)])-sum(r2perm[upper.tri(r2perm)])) - } + ##################################### + ##### Start Parallel permutations ### + ##################################### + + nClust <- nCores - 1 + cl_temp <- parallelly::makeClusterPSOCK(workers = nClust,verbose = T) + doParallel::registerDoParallel(cl_temp) - diffedges.perm[i,] <- abs(r1perm-r2perm)[upper.tri(abs(r1perm-r2perm))] - diffedges.permtemp[upper.tri(diffedges.permtemp, diag=FALSE)] <- diffedges.perm[i,] - diffedges.permtemp <- diffedges.permtemp + t(diffedges.permtemp) - einv.perm.all[,,i] <- diffedges.permtemp - nwinv.perm[i] <- max(diffedges.perm[i,]) + excl <- c() # placeholder for arguments that should NOT be passed to the cluster + # pb <- txtProgressBar(0, it, style = 2) - if(test.centrality==TRUE){ - cen1permtemp <- centrality_auto(r1perm)$node.centrality - cen2permtemp <- centrality_auto(r2perm)$node.centrality - names(cen1permtemp) <- names(cen2permtemp) <- c("betweenness","closeness","strength","expectedInfluence") - if(TRUE %in% (bridgecen %in% centrality)){ - b1permtemp <- networktools::bridge(r1perm, communities=communities, useCommunities=useCommunities) - b2permtemp <- networktools::bridge(r2perm, communities=communities, useCommunities=useCommunities) - names(b1permtemp) <- names(b2permtemp) <- c(bridgecen, "bridgeExpectedInfluence2step", - "communities") - b1permtemp$communities <- b2permtemp$communities <- NULL - cen1permtemp <- data.frame(c(cen1permtemp,b1permtemp)) - cen2permtemp <- data.frame(c(cen2permtemp,b2permtemp)) - } - diffcen.permtemp <- as.matrix(cen1permtemp) - as.matrix(cen2permtemp) - if(nodes[1]=="all"){ - diffcen.perm[i,] <- reshape2::melt(diffcen.permtemp[,centrality])$value - } else { - diffcen.perm[i,] <- reshape2::melt(diffcen.permtemp[which(nodes%in%colnames(data1)),centrality])$value - } + comb <- function(x, ...) { + lapply(seq_along(x), + function(i) c(x[[i]], lapply(list(...), function(y) y[[i]]))) } + permResults <- foreach (i = 1:it, .combine = "comb", .multicombine=TRUE, + .init=list(list(), list()), .verbose = F, + # .export = ls()[!ls()%in%c(excl,"cl_temp")], + .errorhandling = "stop" + # , .packages = c("networktools","bootnet","reshape2") + ) %dopar% { + + diffedges.permtemp <- matrix(0, nvars, nvars) + + # If not paired data + if(paired==FALSE) + { + # s <- sample(1:(nobs1+nobs2),nobs1,replace=FALSE) + # x1perm <- dataall[s,] + # x2perm <- dataall[b[-s], ] + + # Include variance check + okay <- FALSE + counter <- 0 + + if(binary.data) { # if binary data we need to resample the permutation to ensure mininum required variance for glmnet + while(okay == FALSE) { + + # Permute + s <- sample(1:(nobs1+nobs2),nobs1,replace=FALSE) + x1perm <- dataall[s,] + x2perm <- dataall[b[-s], ] + + # check glmnet requirement: at least two instances of each category + ind <- all(apply(x1perm, 2, function(x) min(c(sum(x==0), sum(x==1)))) > 1) & all(apply(x2perm, 2, function(x) min(c(sum(x==0), sum(x==1)))) > 1) + if(ind) okay <- TRUE else counter <- counter + 1 + + } # end: while + } else{ + s <- sample(1:(nobs1+nobs2),nobs1,replace=FALSE) + x1perm <- dataall[s,] + x2perm <- dataall[b[-s], ] + } + + + # Estimate the networks: + r1perm <- do.call(estimator,c(list(x1perm),estimatorArgs)) + if (is.list(r1perm)) r1perm <- r1perm$graph + + r2perm <- do.call(estimator,c(list(x2perm),estimatorArgs)) + if (is.list(r2perm)) r2perm <- r2perm$graph + + if(weighted==FALSE){ + r1perm=(r1perm!=0)*1 + r2perm=(r2perm!=0)*1 + } + } + + # If paired data + if(paired==TRUE) + { + if (verbose) message("Note: NCT for dependent data has not been validated.") + s <- sample(c(1,2),nobs1,replace=TRUE) + x1perm <- x1[s==1,] + x1perm <- rbind(x1perm,x2[s==2,]) + x2perm <- x2[s==1,] + x2perm <- rbind(x2perm,x1[s==2,]) + + # To do: add variance check for paired data + + + # Estimate the networks: + r1perm <- do.call(estimator,c(list(x1perm),estimatorArgs)) + if (is.list(r1perm)) r1perm <- r1perm$graph + + r2perm <- do.call(estimator,c(list(x2perm),estimatorArgs)) + if (is.list(r2perm)) r2perm <- r2perm$graph + + if(weighted==FALSE){ + r1perm=(r1perm!=0)*1 + r2perm=(r2perm!=0)*1 + } + } + + ## Invariance measures for permuted data + if(abs){ + glstrinv.perm[i] <- abs(sum(abs(r1perm[upper.tri(r1perm)]))-sum(abs(r2perm[upper.tri(r2perm)]))) + } else { + glstrinv.perm[i] <- abs(sum(r1perm[upper.tri(r1perm)])-sum(r2perm[upper.tri(r2perm)])) + } + + diffedges.perm[i,] <- abs(r1perm-r2perm)[upper.tri(abs(r1perm-r2perm))] + diffedges.permtemp[upper.tri(diffedges.permtemp, diag=FALSE)] <- diffedges.perm[i,] + diffedges.permtemp <- diffedges.permtemp + t(diffedges.permtemp) + einv.perm.all[,i] <- c(diffedges.permtemp) + nwinv.perm[i] <- max(diffedges.perm[i,]) + for(ji in 1:nvars){ + for(ii in 1:nvars){ + einv.perm.all_temp[ii,ji] = diffedges.permtemp[ii,ji] + } + } + + + if(test.centrality==TRUE){ + cen1permtemp <- centrality_auto(r1perm)$node.centrality + cen2permtemp <- centrality_auto(r2perm)$node.centrality + names(cen1permtemp) <- names(cen2permtemp) <- c("betweenness","closeness","strength","expectedInfluence") + if(TRUE %in% (bridgecen %in% centrality)){ + b1permtemp <- networktools::bridge(r1perm, communities=communities, useCommunities=useCommunities) + b2permtemp <- networktools::bridge(r2perm, communities=communities, useCommunities=useCommunities) + names(b1permtemp) <- names(b2permtemp) <- c(bridgecen, "bridgeExpectedInfluence2step", + "communities") + b1permtemp$communities <- b2permtemp$communities <- NULL + cen1permtemp <- data.frame(c(cen1permtemp,b1permtemp)) + cen2permtemp <- data.frame(c(cen2permtemp,b2permtemp)) + } + diffcen.permtemp <- as.matrix(cen1permtemp) - as.matrix(cen2permtemp) + if(nodes[1]=="all"){ + diffcen.perm[i,] <- reshape2::melt(diffcen.permtemp[,centrality])$value + } else { + diffcen.perm[i,] <- reshape2::melt(diffcen.permtemp[which(nodes%in%colnames(data1)),centrality])$value + } + } else { + diffcen.perm[i,] = NA + } + + # res <- list( + # glstrinv.perm = glstrinv.perm_iter, + # diffedges.perm = diffedges.perm_iter, + # einv.perm.all = einv.perm.all_iter, + # nwinv.perm = nwinv.perm_iter, + # diffcen.perm = diffcen.perm_iter + # ) + if (progressbar==TRUE) setTxtProgressBar(pb, i) + return(NULL) + } + parallelly::autoStopCluster(cl_temp) + gc() - if (progressbar==TRUE) setTxtProgressBar(pb, i) } + + glstrinv.perm = as.numeric(glstrinv.perm[]) + diffedges.perm = diffedges.perm[] + einv.perm.all2 <- aperm(`dim<-`(t(einv.perm.all[]), c(nvars, nvars, it)), c(2, 1, 3)) + nwinv.perm = nwinv.perm[] + diffcen.perm = diffcen.perm[] + # einv.perm.all_temp + ##################################### ##### End permutations ##### ##################################### - # browser() + browser() + ##################################### ##### Calculate p-values ##### diff --git a/man/networkcomparisontest.rd b/man/networkcomparisontest.rd index c7592d3..3c7d2f7 100644 --- a/man/networkcomparisontest.rd +++ b/man/networkcomparisontest.rd @@ -19,7 +19,8 @@ NCT(data1, data2, centrality=c("strength","expectedInfluence"),nodes="all", communities=NULL,useCommunities="all", estimator, estimatorArgs = list(), - verbose = TRUE) + verbose = TRUE, + nCores = 1) } \arguments{ @@ -89,7 +90,7 @@ A function that takes data as input and returns a network structure. This can be Arguments to the \code{estimator} function } \item{verbose}{Logical: Should some warnings and notes be printed?} - +\item{nCores}{Number of cores to use in computing results. Set to 1 to not use parallel computing.} } From 702188f75151acc9bb0d628094630e92eea57103 Mon Sep 17 00:00:00 2001 From: Michael Pinus Date: Sun, 3 Apr 2022 13:16:14 +0300 Subject: [PATCH 2/5] result objects' init by nCores --- DESCRIPTION | 2 +- R/NCT.R | 72 ++++++++++++++---------------------- man/networkcomparisontest.rd | 2 +- 3 files changed, 30 insertions(+), 46 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index fe61c9d..6eba206 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: NetworkComparisonTest Type: Package Title: Statistical Comparison of Two Networks Based on Three Invariance Measures -Version: 2.2.1 +Version: 2.2.2 Authors@R: c( person("Claudia", "van Borkulo", email = "cvborkulo@gmail.com",role = c("aut","cre")), person("Sacha", "Epskamp", role = c("aut")), diff --git a/R/NCT.R b/R/NCT.R index 14cbfc3..eb34cfb 100644 --- a/R/NCT.R +++ b/R/NCT.R @@ -50,6 +50,12 @@ NCT <- function(data1, data2, # Fix for networktools example: if (missing(edges)) edges <- "all" + #progressbar isn't supported when using parallel computing + if (progressbar & nCores > 1){ + message("Note: Progress bar isn't supported when using parallel computing.") + } + + # Small test to warn people of arguments that are ignored if the input is a bootnet object: # Test for bootnet objects: if (is(data1,"bootnetResult") || is(data2,"bootnetResult")){ @@ -173,14 +179,9 @@ NCT <- function(data1, data2, } glstrinv.real <- nwinv.real <- c() - glstrinv.perm <- nwinv.perm <- bigstatsr::FBM(nrow = it,ncol = 1) - diffedges.perm <- bigstatsr::as_FBM(matrix(0,it,nedges)) - einv.perm.all <- bigstatsr::as_FBM(matrix(NA,nvars^2, it)) - einv.perm.all_temp <- bigstatsr::as_FBM(matrix(NA,nvars,nvars)) corrpvals.all <- matrix(NA,nvars,nvars) edges.pvalmattemp <- matrix(0,nvars,nvars) - - + validCentrality <- c("closeness", "betweenness", "strength", "expectedInfluence", "bridgeStrength", "bridgeCloseness", "bridgeBetweenness", "bridgeExpectedInfluence") @@ -191,7 +192,6 @@ NCT <- function(data1, data2, }else { centrality } - diffcen.perm <- bigstatsr::as_FBM(matrix(NA, it, nnodes*length(centrality))) ##################################### ### procedure for all data ### @@ -259,6 +259,12 @@ NCT <- function(data1, data2, ##### Start permutations ##### ##################################### if (nCores == 1) { + + glstrinv.perm <- nwinv.perm <- c() + diffedges.perm <- matrix(0,it,nedges) + einv.perm.all <- array(NA,dim=c(nvars, nvars, it)) + diffcen.perm <- matrix(NA, it, nnodes*length(centrality)) + for (i in 1:it) { diffedges.permtemp <- matrix(0, nvars, nvars) @@ -376,25 +382,18 @@ NCT <- function(data1, data2, ##################################### ##### Start Parallel permutations ### ##################################### - + + glstrinv.perm <- nwinv.perm <- bigstatsr::FBM(nrow = it,ncol = 1) + diffedges.perm <- bigstatsr::as_FBM(matrix(0,it,nedges)) + einv.perm.all <- bigstatsr::as_FBM(matrix(NA,nvars^2, it)) + diffcen.perm <- bigstatsr::as_FBM(matrix(NA, it, nnodes*length(centrality))) + nClust <- nCores - 1 cl_temp <- parallelly::makeClusterPSOCK(workers = nClust,verbose = T) doParallel::registerDoParallel(cl_temp) - excl <- c() # placeholder for arguments that should NOT be passed to the cluster - - # pb <- txtProgressBar(0, it, style = 2) - - comb <- function(x, ...) { - lapply(seq_along(x), - function(i) c(x[[i]], lapply(list(...), function(y) y[[i]]))) - } - - permResults <- foreach (i = 1:it, .combine = "comb", .multicombine=TRUE, - .init=list(list(), list()), .verbose = F, - # .export = ls()[!ls()%in%c(excl,"cl_temp")], + permResults <- foreach (i = 1:it, .errorhandling = "stop" - # , .packages = c("networktools","bootnet","reshape2") ) %dopar% { diffedges.permtemp <- matrix(0, nvars, nvars) @@ -481,12 +480,6 @@ NCT <- function(data1, data2, diffedges.permtemp <- diffedges.permtemp + t(diffedges.permtemp) einv.perm.all[,i] <- c(diffedges.permtemp) nwinv.perm[i] <- max(diffedges.perm[i,]) - for(ji in 1:nvars){ - for(ii in 1:nvars){ - einv.perm.all_temp[ii,ji] = diffedges.permtemp[ii,ji] - } - } - if(test.centrality==TRUE){ cen1permtemp <- centrality_auto(r1perm)$node.centrality @@ -511,34 +504,25 @@ NCT <- function(data1, data2, diffcen.perm[i,] = NA } - # res <- list( - # glstrinv.perm = glstrinv.perm_iter, - # diffedges.perm = diffedges.perm_iter, - # einv.perm.all = einv.perm.all_iter, - # nwinv.perm = nwinv.perm_iter, - # diffcen.perm = diffcen.perm_iter - # ) - if (progressbar==TRUE) setTxtProgressBar(pb, i) return(NULL) } parallelly::autoStopCluster(cl_temp) gc() + glstrinv.perm = as.numeric(glstrinv.perm[]) + diffedges.perm = diffedges.perm[] + einv.perm.all <- `dim<-`(einv.perm.all[], c(nvars, nvars, it)) + nwinv.perm = nwinv.perm[] + diffcen.perm = diffcen.perm[] } - glstrinv.perm = as.numeric(glstrinv.perm[]) - diffedges.perm = diffedges.perm[] - einv.perm.all2 <- aperm(`dim<-`(t(einv.perm.all[]), c(nvars, nvars, it)), c(2, 1, 3)) - nwinv.perm = nwinv.perm[] - diffcen.perm = diffcen.perm[] - # einv.perm.all_temp ##################################### ##### End permutations ##### ##################################### - - browser() - + + # browser() + ##################################### ##### Calculate p-values ##### diff --git a/man/networkcomparisontest.rd b/man/networkcomparisontest.rd index 3c7d2f7..2dd2777 100644 --- a/man/networkcomparisontest.rd +++ b/man/networkcomparisontest.rd @@ -58,7 +58,7 @@ Logical. Can be TRUE of FALSE to indicate whether or not differences in individu Character or list. When 'all', differences between all individual edges are tested. When provided a list with one or more pairs of indices referring to variables, the provided edges are tested. } \item{progressbar}{ -Logical. Should the pbar be plotted in order to see the progress of the estimation procedure? Defaults to TRUE. +Logical. Should the pbar be plotted in order to see the progress of the estimation procedure? Defaults to TRUE. Not supported when using parallel computing (nCores > 1). } \item{make.positive.definite}{ If \code{make.positive.definite = TRUE}, the covariance matrices used for the glasso are projected to the nearest positive definite matrices, if they are not yet positive definite. This is useful for small n, for which it is very likely that at least one of the bootstrap comparisons involves a covariance matrix that is not positive definite. From 236c34b9c9420f93964dde4b906eef0d3b673555 Mon Sep 17 00:00:00 2001 From: Michael Pinus Date: Sun, 3 Apr 2022 13:48:25 +0300 Subject: [PATCH 3/5] remove snow dependece --- DESCRIPTION | 1 - 1 file changed, 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 6eba206..45809a1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -10,7 +10,6 @@ Authors@R: c( person("Alex", "Millner", role = c("ctb"))) Maintainer: Claudia van Borkulo Imports: - snow, foreach, qgraph, IsingFit, From 6d1c9aefdf234048f45747ae912bebf8735ffc03 Mon Sep 17 00:00:00 2001 From: Michael Pinus Date: Wed, 22 Jun 2022 00:34:58 +0300 Subject: [PATCH 4/5] Update NetworkComparisonTest.Rproj --- NetworkComparisonTest.Rproj | 1 + 1 file changed, 1 insertion(+) diff --git a/NetworkComparisonTest.Rproj b/NetworkComparisonTest.Rproj index 21a4da0..176148b 100644 --- a/NetworkComparisonTest.Rproj +++ b/NetworkComparisonTest.Rproj @@ -14,4 +14,5 @@ LaTeX: pdfLaTeX BuildType: Package PackageUseDevtools: Yes +PackageCleanBeforeInstall: Yes PackageInstallArgs: --no-multiarch --with-keep.source From a56712e7d8eb91904ca192ee5b55b4b1e56fc9c5 Mon Sep 17 00:00:00 2001 From: Michael Pinus Date: Wed, 22 Jun 2022 01:00:13 +0300 Subject: [PATCH 5/5] Revert "Merge remote-tracking branch 'upstream/master'" This reverts commit 0b4795319d5170e2c80d901e44025920e8b06d88, reversing changes made to 6d1c9aefdf234048f45747ae912bebf8735ffc03. --- DESCRIPTION | 20 +++--- NAMESPACE | 5 +- R/NCT.R | 62 +++++++--------- R/s3Methods.R | 11 ++- inst/CITATION | 11 ++- man/networkcomparisontest-package.Rd | 10 ++- man/networkcomparisontest.rd | 102 +++++++++------------------ 7 files changed, 83 insertions(+), 138 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 2727b9b..45809a1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: NetworkComparisonTest Type: Package -Title: Statistical Comparison of Two Networks Based on Several Invariance Measures +Title: Statistical Comparison of Two Networks Based on Three Invariance Measures Version: 2.2.2 Authors@R: c( person("Claudia", "van Borkulo", email = "cvborkulo@gmail.com",role = c("aut","cre")), @@ -10,7 +10,6 @@ Authors@R: c( person("Alex", "Millner", role = c("ctb"))) Maintainer: Claudia van Borkulo Imports: - bootnet, foreach, qgraph, IsingFit, @@ -26,14 +25,13 @@ Imports: parallelly Suggests: networktools -Description: This permutation based hypothesis test, suited for several types of data - supported by the estimateNetwork function of the bootnet package (Epskamp & Fried, 2018), - assesses the difference between two networks based on several invariance measures (network - structure invariance, global strength invariance, edge invariance, several centrality - measures, etc.). Network structures are estimated with l1-regularization. The Network - Comparison Test is suited for comparison of independent (e.g., two different groups) and - dependent samples (e.g., one group that is measured twice). See van Borkulo et al. (2021, - in press; the final article will be available, upon publication, via its DOI: - 10.1037/met0000476). +Description: This permutation based hypothesis test, suited for Gaussian and binary data, + assesses the difference between two networks based on several invariance measures + (e.g., network structure invariance, global strength invariance, edge invariance). + Network structures are estimated with l1-regularized partial correlations (Gaussian data) + or with l1-regularized logistic regression (eLasso, binary data). Suited for comparison + of independent and dependent samples. For dependent samples, only supported for data of + one group which is measured twice. See van Borkulo et al. (2017) + . License: GPL-2 RoxygenNote: 7.1.2 diff --git a/NAMESPACE b/NAMESPACE index 0123008..e33bcbb 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,7 +7,6 @@ import(qgraph) import(IsingFit) import(IsingSampler) import(reshape2) -import(bootnet) import(foreach) importFrom("graphics", "hist", "points") importFrom("stats", "cor", "p.adjust") @@ -15,4 +14,6 @@ importFrom("utils", "setTxtProgressBar", "txtProgressBar") importFrom("Matrix", "nearPD") importFrom("reshape2","melt") importFrom("stats", "na.omit") -importFrom("methods", "is") + + +importFrom("methods", "is") \ No newline at end of file diff --git a/R/NCT.R b/R/NCT.R index cc24084..eb34cfb 100644 --- a/R/NCT.R +++ b/R/NCT.R @@ -195,6 +195,7 @@ NCT <- function(data1, data2, ##################################### ### procedure for all data ### + ##################################### # Estimate the networks: @@ -257,24 +258,16 @@ NCT <- function(data1, data2, ##################################### ##### Start permutations ##### ##################################### -==== BASE ==== - - for (i in 1:it) - { - diffedges.permtemp <- matrix(0, nvars, nvars) -==== BASE ==== + if (nCores == 1) { + + glstrinv.perm <- nwinv.perm <- c() + diffedges.perm <- matrix(0,it,nedges) + einv.perm.all <- array(NA,dim=c(nvars, nvars, it)) + diffcen.perm <- matrix(NA, it, nnodes*length(centrality)) for (i in 1:it) { -==== BASE ==== - # s <- sample(1:(nobs1+nobs2),nobs1,replace=FALSE) - # x1perm <- dataall[s,] - # x2perm <- dataall[b[-s], ] - - # Include variance check - okay <- FALSE - counter <- 0 -==== BASE ==== + diffedges.permtemp <- matrix(0, nvars, nvars) # If not paired data if(paired==FALSE) @@ -345,24 +338,13 @@ NCT <- function(data1, data2, r2perm=(r2perm!=0)*1 } } -==== BASE ==== - } - - # If paired data - if(paired==TRUE) - { - if (verbose) message("Note: NCT for dependent data has not been validated.") - s <- sample(c(1,2),nobs1,replace=TRUE) - x1perm <- x1[s==1,] - x1perm <- rbind(x1perm,x2[s==2,]) - x2perm <- x2[s==1,] - x2perm <- rbind(x2perm,x1[s==2,]) - # To do: add variance check for paired data -==== BASE ==== - -==== BASE ==== -==== BASE ==== + ## Invariance measures for permuted data + if(abs){ + glstrinv.perm[i] <- abs(sum(abs(r1perm[upper.tri(r1perm)]))-sum(abs(r2perm[upper.tri(r2perm)]))) + } else { + glstrinv.perm[i] <- abs(sum(r1perm[upper.tri(r1perm)])-sum(r2perm[upper.tri(r2perm)])) + } diffedges.perm[i,] <- abs(r1perm-r2perm)[upper.tri(abs(r1perm-r2perm))] diffedges.permtemp[upper.tri(diffedges.permtemp, diag=FALSE)] <- diffedges.perm[i,] @@ -538,6 +520,9 @@ NCT <- function(data1, data2, ##################################### ##### End permutations ##### ##################################### + + # browser() + ##################################### ##### Calculate p-values ##### @@ -558,9 +543,10 @@ NCT <- function(data1, data2, rownames(corrpvals.all) <- colnames(corrpvals.all) <- colnames(x1) einv.pvals <- melt(corrpvals.all, na.rm=TRUE, value.name = 'p-value') einv.perm <- einv.perm.all - einv.real <- diffedges.realoutput - einv.pvals <- cbind(einv.pvals, round(einv.real[upper.tri(einv.real)],8)) - colnames(einv.pvals) <- c('Var1', 'Var2', 'p-value', "Test statistic E") + einv.real <- diffedges.realoutput + + edges.tested <- "all" + } ## If a selection of edges should be tested @@ -588,10 +574,11 @@ NCT <- function(data1, data2, corrpvals_mat[,3] <- corrpvals corrpvals_mat[,1:2] <- pairs einv.pvals <- as.data.frame(corrpvals_mat) - einv.pvals <- cbind(einv.pvals, einv.real) - colnames(einv.pvals) <- c('Var1', 'Var2', 'p-value', "Test statistic E") + colnames(einv.pvals) <- c('Var1', 'Var2', 'p-value') } + edges.tested <- colnames(einv.perm) + res <- list(glstrinv.real = glstrinv.real, glstrinv.sep = glstrinv.sep, glstrinv.pval = (sum(glstrinv.perm >= glstrinv.real) + 1) / (it + 1), @@ -599,6 +586,7 @@ NCT <- function(data1, data2, nwinv.real = nwinv.real, nwinv.pval = (sum(nwinv.perm >= nwinv.real) + 1) / (it + 1), nwinv.perm = nwinv.perm, + edges.tested = edges.tested, einv.real = einv.real, einv.pvals = einv.pvals, einv.perm = einv.perm, diff --git a/R/s3Methods.R b/R/s3Methods.R index 7499ebe..ba6525f 100644 --- a/R/s3Methods.R +++ b/R/s3Methods.R @@ -14,7 +14,7 @@ print.NCT <- function(x,...){ print(x$einv.pvals) } if(x$info$call$test.centrality){ - cat("\n CENTRALITY INVARIANCE TEST p-value\n") + cat("\n CENTRALITY INVARIANCE TEST \n") print(x$diffcen.pval) } } @@ -46,11 +46,10 @@ summary.NCT <- function(object,...){ "\n p-value", object$nwinv.pval, global_stat_message, object$glstrinv.sep, "\n Test statistic S: ", object$glstrinv.real, - "\n p-value", object$glstrinv.pval) - if(object$info$call$test.edges){ - cat("\n\n EDGE INVARIANCE TEST \n") - print(object$einv.pvals) - } + "\n p-value", object$glstrinv.pval, + "\n\n EDGE INVARIANCE TEST \n Edges tested:", object$edges.tested, + "\n Test statistic E:", object$einv.real, + "\n p-value:", object$einv.pvals$`p-value`) if(object$info$call$test.centrality){ cenTest <- c(reshape2::melt(object$diffcen.real)$value, reshape2::melt(object$diffcen.pval)$value) diff --git a/inst/CITATION b/inst/CITATION index 65204dd..5cba3fa 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -3,19 +3,16 @@ citHeader("To cite NCT in publications use:") citEntry(entry = "Article", title = "Comparing network structures on three aspects: A permutation test", author = personList(as.person("Claudia D. van Borkulo"), - as.person("Riet van Bork"), as.person("Lynn Boschloo"), as.person("Jolanda J. Kossakowski"), as.person("Pia Tio"), as.person("Robert A. Schoevers"), as.person("Denny Borsboom"), as.person("Lourens J. Waldorp")), - year = "2021", - journal = "Psychological Methods", - doi = "DOI: 10.1037/met0000476", + year = "2017", + doi = "10.13140/RG.2.2.29455.38569", textVersion = - paste("Claudia D. van Borkulo, Riet van Bork, Lynn Boschloo, Jolanda J. Kossakowski, Pia Tio, Robert A. Schoevers, Denny Borsboom, Lourens J. Waldorp (2021).", + paste("Claudia D. van Borkulo, Lynn Boschloo, Jolanda J. Kossakowski, Pia Tio, Robert A. Schoevers, Denny Borsboom, Lourens J. Waldorp (2017).", "Comparing network structures on three aspects: A permutation test.", - "Psychological Methods", - "DOI: 10.1037/met0000476.") + "doi:10.13140/RG.2.2.29455.38569.") ) diff --git a/man/networkcomparisontest-package.Rd b/man/networkcomparisontest-package.Rd index 4a876b9..6dba699 100644 --- a/man/networkcomparisontest-package.Rd +++ b/man/networkcomparisontest-package.Rd @@ -2,16 +2,16 @@ \alias{NetworkComparisonTest-package} \docType{package} \title{ -Statistical Comparison of Two Networks Based on Several Invariance Measures +Statistical Comparison of Two Networks Based on Three Invariance Measures } \description{ -This permutation based hypothesis test, suited for several types of data supported by the estimateNetwork function of the bootnet package (Epskamp & Fried, 2018), assesses the difference between two networks based on several invariance measures (network structure invariance, global strength invariance, edge invariance, several centrality measures, etc.). Network structures are estimated with l1-regularization. The Network Comparison Test is suited for comparison of independent (e.g., two different groups) and dependent samples (e.g., one group that is measured twice). +This permutation based hypothesis test, suited for gaussian and binary data, assesses the difference between two networks based on several invariance measures (network structure invariance, global strength invariance, edge invariance). Network structures are estimated with l1-regularized partial correlations (gaussian data) or with l1-regularized logistic regression (eLasso, binary data). Suited for comparison of independent and dependent samples. For dependent samples, only supported for data of one group which is measured twice. } \details{ \tabular{ll}{ Package: \tab NetworkComparisonTest\cr Type: \tab Package\cr -Version: \tab 2.2.2\cr +Version: \tab 2.2.1\cr License: \tab GPL-2\cr } } @@ -21,15 +21,13 @@ Claudia D. van Borkulo, with contributions from Jonas Haslbeck, Sacha Epskamp, P Maintainer: Claudia D. van Borkulo } \references{ -Epskamp, S., & Fried, E. I. (2018). A tutorial on regularized partial correlation networks. Psychological Methods, 23(4), 617-634. - Ernst, M.D. Permutation methods: A basis for exact inference. Stat Sci. 2004;19(4):676-685. Good, P.I. Permutation, parametric and bootstrap tests of hypotheses. Vol. 3. New York:: Springer, 2005. van Borkulo, C. D., Boschloo, L., Borsboom, D., Penninx, B. W. J. H., Waldorp, L. J., & Schoevers, R.A. (2015). Association of symptom network structure with the course of depression. JAMA Psychiatry. 2015;72(12). doi:10.1001/jamapsychiatry.2015.2079 -van Borkulo, C. D., van Bork, R., Boschloo, Kossakowski, J., Tio, P., L., Schoevers, R.A., Borsboom, D., & , Waldorp, L. J. (2021). Comparing network structures on three aspects: A permutation test. Psychological Methods, in press. +van Borkulo, C. D., Boschloo, Kossakowski, J., Tio, P., L., Schoevers, R.A., Borsboom, D., & , Waldorp, L. J. (2016). Comparing network structures on three aspects: A permutation test. doi:10.13140/RG.2.2.29455.38569 } % ~~ Optionally other standard keywords, one per line, from ~~ % ~~ file KEYWORDS in the R documentation directory ~~ diff --git a/man/networkcomparisontest.rd b/man/networkcomparisontest.rd index d110cbf..2dd2777 100644 --- a/man/networkcomparisontest.rd +++ b/man/networkcomparisontest.rd @@ -2,10 +2,10 @@ \alias{NetworkComparisonTest} \alias{NCT} \title{ -Statistical Comparison of Two Networks Based on Several Invariance Measures +Statistical Comparison of Two Networks Based on Three Invariance Measures } \description{ -This permutation based hypothesis test, suited for several types of data supported by the estimateNetwork function of the bootnet package (Epskamp & Fried, 2018), assesses the difference between two networks based on several invariance measures (network structure invariance, global strength invariance, edge invariance, several centrality measures, etc.). Network structures are estimated with l1-regularization. The Network Comparison Test is suited for comparison of independent (e.g., two different groups) and dependent samples (e.g., one group that is measured twice). +This permutation based hypothesis test, suited for gaussian and binary data, assesses the difference between two networks based on several invariance measures (network structure invariance, global strength invariance, edge invariance). Network structures are estimated with l1-regularized partial correlations (gaussian data) or with l1-regularized logistic regression (eLasso, binary data). Suited for comparison of independent and dependent samples. For dependent samples, only supported for data of one group which is measured twice. } \usage{ NCT(data1, data2, @@ -37,7 +37,7 @@ A single value between 0 and 1. When not entered, gamma is set to 0.25 for binar The number of iterations (permutations). } \item{binary.data}{ -Logical. Can be TRUE or FALSE to indicate whether the data is binary or not. If binary.data is FALSE, the data is regarded gaussian. This argument is ignored when using estimateNetwork() output as input for NCT. +Logical. Can be TRUE or FALSE to indicate whether the data is binary or not. If binary.data is FALSE, the data is regarded gaussian. } \item{paired}{ Logical. Can be TRUE of FALSE to indicate whether the samples are dependent or not. If paired is TRUE, relabeling is performed within each pair of observations. If paired is FALSE, relabeling is not restricted to pairs of observations. Note that, currently, dependent data is assumed to entail one group measured twice. @@ -105,6 +105,7 @@ NCT returns a 'NCT' object that contains the following items: \item{nwinv.perm}{The values of the maximum difference in edge weights of the permuted networks} \item{nwinv.pval }{The p value resulting from the permutation test concerning the maximum difference in edge weights.} \item{einv.pvals}{p-values (corrected for multiple testing or not according to 'p.adjust.methods') per edge from the permutation test concerning differences in edges weights. Only returned if test.edges = TRUE.} +\item{edges.tested}{The pairs of variables between which the edges are called to be tested. Only if test.edges = TRUE.} \item{einv.real}{The value of the difference in edge weight of the observed networks (multiple values if more edges are called to test). Only if test.edges = TRUE.} \item{einv.perm}{The values of the difference in edge weight of the permuted networks. Only if test.edges = TRUE.} \item{diffcen.real}{The values of the difference in centralities of the observed networks. Only if test.centrality = TRUE.} @@ -119,7 +120,7 @@ Good, P.I. Permutation, parametric and bootstrap tests of hypotheses. Vol. 3. Ne van Borkulo, C. D., Boschloo, L., Borsboom, D., Penninx, B. W. J. H., Waldorp, L. J., & Schoevers, R.A. (2015). Association of symptom network structure with the course of depression. JAMA Psychiatry. 2015;72(12). doi:10.1001/jamapsychiatry.2015.2079 -van Borkulo, C. D., van Bork, R., Boschloo, Kossakowski, J., Tio, P., L., Schoevers, R.A., Borsboom, D., & , Waldorp, L. J. (2021). Comparing network structures on three aspects: A permutation test. DOI: 10.1037/met0000476 +van Borkulo, C. D., Boschloo, Kossakowski, J., Tio, P., L., Schoevers, R.A., Borsboom, D., & , Waldorp, L. J. (2016). Comparing network structures on three aspects: A permutation test. doi:10.13140/RG.2.2.29455.38569 } \author{ Claudia D. van Borkulo, with contributions from Jonas Haslbeck, Sacha Epskamp, Payton Jones and Alex Millner @@ -133,24 +134,16 @@ See also my website: http://cvborkulo.com \examples{ library("IsingSampler") library("IsingFit") -library("bootnet") ### Simulate binary datasets under null hypothesis: -### underlying network structures are similar +### underlying network structures have the same strength # Input: N <- 6 # Number of nodes nSample <- 500 # Number of samples # Ising parameters: -set.seed(123) Graph <- matrix(sample(0:1,N^2,TRUE,prob = c(0.8, 0.2)),N,N) * runif(N^2,0.5,2) Graph <- pmax(Graph,t(Graph)) -Graph[4,1] <- Graph[4,1]*-1 -Graph[1,4] <- Graph[1,4]*-1 -Graph[5,1] <- Graph[5,1]*-1 -Graph[1,5] <- Graph[1,5]*-1 -Graph[6,1] <- Graph[6,1]*-1 -Graph[1,6] <- Graph[1,6]*-1 diag(Graph) <- 0 Thresh <- -rowSums(Graph) / 2 @@ -160,62 +153,33 @@ data2 <- IsingSampler(nSample, Graph, Thresh) colnames(data1) <- colnames(data2) <- c('V1', 'V2', 'V3', 'V4', 'V5', 'V6') ### Compare networks of data sets using NCT ### -## Networks can be compared by either (1) feeding the data directly into NCT (whereby -## you need to specify arguments such as "gamma" and "binary.data") or (2) by using -## estimateNetwork() (bootnet package) and feeding that output into NCT. For the latter -## option, we refer to the help file of estimateNetwork() for its usage. Below, both -## options are illustrated. We recommend using estimateNetwork(), since this function -## has implemented many network estimation methods. - -## gamma = 0 (in estimateNetwork this hyperparameter is called "tuning"; to illustrate -# how to specify a different value than the default) -## iterations (it) set to 10 to save time -## Note: Low number of iterations can give unreliable results; should be 1000 at least - -## Testing whether there are differences in the three aspects that are validated -# (network invariance, global strength, edge weight) -## 2 edges are tested here: between variable 1 and 2, and between 3 and 6 (can be -# "list(c(2,1),c(6,3))" as well) - -## (1) Feeding data directly into NCT -set.seed(123) -NCT_a <- NCT(data1, data2, gamma=0, it=10, binary.data = TRUE, - test.edges=TRUE, edges=list(c(1,2),c(3,6))) -summary(NCT_a) -## Plot results of global strength invariance test (not reliable with only 10 -# permutations!) -plot(NCT_a, what="strength") - -## (2) Feeding the estimateNetwork() output into NCT -est_1 <- estimateNetwork(data1, default = "IsingFit", tuning = 0) -est_2 <- estimateNetwork(data2, default = "IsingFit", tuning = 0) -## When using estimateNetwork() output, there is no need to specify gamma and binary.data -## This yields similar output as NCT_a -set.seed(123) -NCT_b <- NCT(est_1, est_2, it=10, test.edges=TRUE, -edges=list(c(1,2),c(3,6))) -summary(NCT_b) - -## Next, an example of testing whether there are differences in node strength -# when data is paired (e.g., a group which is measured pre- and post-treatement). -# Also, here you can see how to specify that you want to take the sign of node strength -# into account (by default, the absolute value is taken and, therefore, the sign is -# ignored). - -## abs = FALSE -set.seed(123) -NCT_c = NCT(est_1, est_2, paired = TRUE, abs = FALSE, test.edges = TRUE, -edges = list(c(1,2),c(3,6)), test.centrality = TRUE, -centrality = c("strength"), nodes = "all", it=10) -summary(NCT_c) - -## Finally, an example how to test for differences in centrality (e.g., expectedInfluence) - -set.seed(123) -NCT_d = NCT(est_1, est_2, paired = TRUE, abs = FALSE, test.edges = TRUE, -edges = list(c(1,2),c(3,6)), test.centrality = TRUE, -centrality = c("expectedInfluence"), nodes = "all", it=10) -summary(NCT_d) +# with gamma = 0. +# Iterations (it) set to 10 to save time. +# Low number of iterations can give unreliable results. Should be 1000 at least. + +# Testing the three aspects that are validated (network invariance, global strength, edge weight) +# 2 edges are tested here: between variable 1 and 2, +# and between 3 and 6 (can be list(c(2,1),c(6,3)) as well) +Res_1 <- NCT(data1, data2, gamma=0, it=10, binary.data = TRUE, +test.edges=TRUE, edges=list(c(1,2),c(3,6))) + +## Plotting of NCT results +## See the help file of plot.NCT for more information about the plotting function and its arguments + +# Plot results of the network structure invariance test (not reliable with only 10 permutations!): +plot(Res_1, what="network") + +# Plot results of global strength invariance test (not reliable with only 10 permutations!): +plot(Res_1, what="strength") + +# Plot results of the edge invariance test (not reliable with only 10 permutations!): +# Note that two distributions are plotted +plot(Res_1, what="edge") + +# Without testing for (an) individual edge(s) +# The arguments 'test.edges' and 'edges' don't need to be specified +# Not run +# Res_0 <- NCT(data1, data2, gamma=0, it=10, binary.data = TRUE) } % Add one or more standard keywords, see file 'KEYWORDS' in the % R documentation directory.