library(IsingSampler)
library(IsingFit)
library(bootnet)
### Simulate binary datasets under null hypothesis:
### underlying network structures are similar
# 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
# Simulate:
data1 <- IsingSampler(nSample, Graph, Thresh)
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).
# we don't run these two examples by default as they take too long for the R CMD check
# but they are still interesting.
if (FALSE) {
## 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)
}
Run the code above in your browser using DataLab