# In this code we used the same data set but permute the T nodes
# to generate multiple permuted data sets and apply the methods
# to check the variation in inferred graphs.
CompareMethodsNodeOrdering <- function (N, signal,model,n_data, n_nodeordering) {
# Parameters
p <- 0.45
b0.1 <- 0
b1.1 <- signal
b1.2 <- signal
b1.3 <- signal
sd.1 <- 1
switch(model,
truth1 = {
# Truth1 model (V1-->T1-->T2-->T3)
tarmat_s1 <- matrix(0,nrow=4,ncol = 4)
colnames(tarmat_s1) <- c("V1","T1","T2","T3")
rownames(tarmat_s1) <- colnames(tarmat_s1)
# Adjacency matrix from the true graph
tarmat_s1[1,2] <- 1
tarmat_s1[2,3] <- 1
tarmat_s1[3,4] <- 1
# The number of genetic variants in the graph.
GV <- 1
# The number of nodes
n.nodes <- ncol (tarmat_s1)
# Lists for the adjacency matrix output by each method.
# MRPC
Adjlist_MRPC <- vector (mode = 'list', length = n_data)
# pc
Adjlist_PC <- vector (mode = 'list', length = n_data)
# pc.stable
Adjlist_pc.stable <- vector (mode = 'list', length = n_data)
# mmpc
Adjlist_mmpc <- vector (mode = 'list', length = n_data)
# mmhc
Adjlist_mmhc <- vector (mode = 'list', length = n_data)
# hc
Adjlist_hc <- vector (mode = 'list', length = n_data)
# Matrix for the number of unique graphs for each data set (across n permutations).
countMatrix <- matrix (nrow = n_data, ncol = 6)
# name the columns according to the method.
colnames (countMatrix) <- c('MRPC', 'pc', 'pc.stable', 'mmpc', 'mmhc', 'hc')
# Loop through the independent data sets.
for (e in 1:n_data) {
cat ("Dataset:", e, "\n")
# Simulate data for the true1.
V1 <- c(sample(c(0, 1, 2),size = N,replace = TRUE,prob = c((1 - p)^2,2*p*(1 - p),p^2)))
T1 <- SimulateData1P(N=N,P1=V1,b0.1=b0.1,b1.1=b1.1,sd.1=sd.1)
T2 <- SimulateData1P(N=N,P1=T1,b0.1=b0.1,b1.1=b1.1,sd.1=sd.1)
T3 <- SimulateData1P(N=N,P1=T2,b0.1=b0.1,b1.1=b1.1,sd.1=sd.1)
Data1 <- cbind(V1,T1,T2,T3)
# MRPC
Adjlist_MRPC[[e]] <- vector(mode = 'list', length = n_nodeordering)
# pc
Adjlist_PC[[e]] <- vector(mode = 'list', length = n_nodeordering)
# pc.stable
Adjlist_pc.stable[[e]] <- vector(mode = 'list', length = n_nodeordering)
# mmpc
Adjlist_mmpc[[e]] <- vector(mode = 'list', length = n_nodeordering)
# mmhc
Adjlist_mmhc[[e]] <- vector(mode = 'list', length = n_nodeordering)
# hc
Adjlist_hc[[e]] <- vector(mode = 'list', length = n_nodeordering)
# Loop through all node permutations for one data set.
for (j in 1:n_nodeordering) {
cat ("Node ordering:", j, "\n")
# Create a new ordering for the T nodes.
temp.order <- c (GV, sample((GV + 1):n.nodes))
# New data matrix with permuted nodes.
Data2 <- Data1[, temp.order]
n <- nrow (Data2) # Sample size
V <- colnames (Data2) # Node labels
# Calculate Pearson correlation
suffStat <- list (C = cor (Data2), n = n)
# Infer the graph by MRPC
MRPC_Inferred <- MRPC (Data2, suffStat, GV = GV, FDR = 0.05, FDRcontrol = 'LOND', indepTest = 'gaussCItest', labels = V, verbose = FALSE)
# Adjacency matrix from the graph by MRPC
G_MRPC <- as (MRPC_Inferred@graph, "matrix")
# Save adjacency matrix
Adjlist_MRPC[[e]][[j]] <- AdjustMatrix (tarmat_s1, G_MRPC)
# Infer the graph by PC
PC_Inferred <- pc (suffStat, alpha = 0.05, indepTest = gaussCItest, labels = V, verbose = FALSE)
# Adjacency matrix from the graph by pc
G_PC <- as (PC_Inferred@graph, "matrix")
# Save adjacency matrix
Adjlist_PC[[e]][[j]] <- AdjustMatrix (tarmat_s1, G_PC)
# arcs not to be included from gene expression to genotype
to <- rep (colnames (Data2)[1:GV], each = (ncol (Data2) - GV))
from <- rep (colnames (Data2)[(GV + 1):ncol (Data2)], GV)
bl <- cbind (from, to)
# Infer the graph by pc.stable
pc.stable_Inferred <- pc.stable (data.frame (Data2), blacklist = bl, alpha = 0.05, B = NULL, max.sx = NULL, debug = FALSE, undirected = FALSE)
# Inferred graph object by pc.stable
G_pc.stable <- amat (pc.stable_Inferred)
# Save adjacency matrix
Adjlist_pc.stable[[e]][[j]] <- AdjustMatrix (tarmat_s1, G_pc.stable)
# Infer the graph by mmpc
mmpc_Inferred <- mmpc (data.frame (Data2), blacklist = bl, alpha = 0.05, B = NULL, max.sx = NULL, debug = FALSE, undirected = FALSE)
# Inferred graph object by mmpc
G_mmpc <- amat (mmpc_Inferred)
# Save adjacency matrix
Adjlist_mmpc[[e]][[j]] <- AdjustMatrix (tarmat_s1, G_mmpc)
# Infer the graph by mmhc
mmhc_Inferred <- mmhc (data.frame (Data2), blacklist = bl, debug = FALSE)
# Inferred graph object by mmhc
G_mmhc <- amat (mmhc_Inferred)
# Save adjacency matrix
Adjlist_mmhc[[e]][[j]] <- AdjustMatrix (tarmat_s1, G_mmhc)
# Infer the graph by hc
hc_Inferred <- hc (data.frame (Data2), blacklist = bl, debug = FALSE)
# Inferred graph object by mmhc
G_hc <- amat (hc_Inferred)
# Save adjacency matrix
Adjlist_hc[[e]][[j]] <- AdjustMatrix (tarmat_s1, G_hc)
}
# Calculate the number of unique graphs inferred by each method across all permutations.
countMatrix[e, ] <- c(length (unique (Adjlist_MRPC[[e]])), length (unique (Adjlist_PC[[e]])), length (unique (Adjlist_pc.stable[[e]])), length (unique (Adjlist_mmpc[[e]])), length( unique (Adjlist_mmhc[[e]])), length( unique (Adjlist_hc[[e]])))
}
return(countMatrix)
},
truth2 = {
# Truth2 model (V1-->T1<--T2-->T3)
tarmat_s2 <- matrix(0,nrow=4,ncol = 4)
colnames(tarmat_s2) <- c("V1","T1","T2","T3")
rownames(tarmat_s2) <- colnames(tarmat_s2)
# Adjacency matrix from the true graph
tarmat_s2[1,2] <- 1
tarmat_s2[3,2] <- 1
tarmat_s2[3,4] <- 1
# The number of genetic variants in the graph.
GV <- 1
# The number of nodes
n.nodes <- ncol (tarmat_s2)
# Lists for the adjacency matrix output by each method.
# MRPC
Adjlist_MRPC <- vector (mode = 'list', length = n_data)
# pc
Adjlist_PC <- vector (mode = 'list', length = n_data)
# pc.stable
Adjlist_pc.stable <- vector (mode = 'list', length = n_data)
# mmpc
Adjlist_mmpc <- vector (mode = 'list', length = n_data)
# mmhc
Adjlist_mmhc <- vector (mode = 'list', length = n_data)
# hc
Adjlist_hc <- vector (mode = 'list', length = n_data)
# Matrix for the number of unique graphs for each data set (across n permutations).
countMatrix <- matrix (nrow = n_data, ncol = 6)
# name the columns according to the method.
colnames (countMatrix) <- c('MRPC', 'pc', 'pc.stable', 'mmpc', 'mmhc', 'hc')
# Loop through the independent data sets.
for (e in 1:n_data) {
cat ("Dataset:", e, "\n")
# Simulate data for the true2.
V1 <- c(sample(c(0, 1, 2),size = N,replace = TRUE,prob = c((1 - p)^2,2*p*(1 - p),p^2)))
T2 <- SimulateDataNP(N=N,b0.1=b0.1,sd.1=sd.1)
T1 <- SimulateData2P(N=N,P1=V1,P2=T2,b0.1=b0.1,b1.1=b1.1,b1.2=b1.2,sd.1=sd.1)
T3 <- SimulateData1P(N=N,P1=T2,b0.1=b0.1,b1.1=b1.1,sd.1=sd.1)
# Combined
Data1 <- cbind(V1,T1,T2,T3)
# MRPC
Adjlist_MRPC[[e]] <- vector(mode = 'list', length = n_nodeordering)
# pc
Adjlist_PC[[e]] <- vector(mode = 'list', length = n_nodeordering)
# pc.stable
Adjlist_pc.stable[[e]] <- vector(mode = 'list', length = n_nodeordering)
# mmpc
Adjlist_mmpc[[e]] <- vector(mode = 'list', length = n_nodeordering)
# mmhc
Adjlist_mmhc[[e]] <- vector(mode = 'list', length = n_nodeordering)
# hc
Adjlist_hc[[e]] <- vector(mode = 'list', length = n_nodeordering)
# Loop through all node permutations for one data set.
for (j in 1:n_nodeordering) {
cat ("Node ordering:", j, "\n")
# Create a new ordering for the T nodes.
temp.order <- c (GV, sample((GV + 1):n.nodes))
# New data matrix with permuted nodes.
Data2 <- Data1[, temp.order]
n <- nrow (Data2) # Sample size
V <- colnames (Data2) # Node labels
# Calculate Pearson correlation
suffStat <- list (C = cor (Data2), n = n)
# Infer the graph by MRPC
MRPC_Inferred <- MRPC (Data2, suffStat, GV = GV, FDR = 0.05, FDRcontrol = 'LOND', indepTest = 'gaussCItest', labels = V, verbose = FALSE)
# Adjacency matrix from the graph by MRPC
G_MRPC <- as (MRPC_Inferred@graph, "matrix")
# Save adjacency matrix
Adjlist_MRPC[[e]][[j]] <- AdjustMatrix (tarmat_s2, G_MRPC)
# Infer the graph by PC
PC_Inferred <- pc (suffStat, alpha = 0.05, indepTest = gaussCItest, labels = V, verbose = F)
# Adjacency matrix from the graph by pc
G_PC <- as (PC_Inferred@graph, "matrix")
# Save adjacency matrix
Adjlist_PC[[e]][[j]] <- AdjustMatrix (tarmat_s2, G_PC)
# arcs not to be included from gene expression to genotype
to <- rep (colnames (Data2)[1:GV], each = (ncol (Data2) - GV))
from <- rep (colnames (Data2)[(GV + 1):ncol (Data2)], GV)
bl <- cbind (from, to)
# Infer the graph by pc.stable
pc.stable_Inferred <- pc.stable (data.frame (Data2), blacklist = bl, alpha = 0.05, B = NULL, max.sx = NULL, debug = FALSE, undirected = FALSE)
# Inferred graph object by pc.stable
G_pc.stable <- amat (pc.stable_Inferred)
# Save adjacency matrix
Adjlist_pc.stable[[e]][[j]] <- AdjustMatrix (tarmat_s2, G_pc.stable)
# Infer the graph by mmpc
mmpc_Inferred <- mmpc (data.frame (Data2), blacklist = bl, alpha = 0.05, B = NULL, max.sx = NULL, debug = FALSE, undirected = FALSE)
# Inferred graph object by mmpc
G_mmpc <- amat (mmpc_Inferred)
# Save adjacency matrix
Adjlist_mmpc[[e]][[j]] <- AdjustMatrix (tarmat_s2, G_mmpc)
# Infer the graph by mmhc
mmhc_Inferred <- mmhc (data.frame (Data2), blacklist = bl, debug = FALSE)
# Inferred graph object by mmhc
G_mmhc <- amat (mmhc_Inferred)
# Save adjacency matrix
Adjlist_mmhc[[e]][[j]] <- AdjustMatrix (tarmat_s2, G_mmhc)
# Infer the graph by hc
hc_Inferred <- hc (data.frame (Data2), blacklist = bl, debug = FALSE)
# Inferred graph object by hc
G_hc <- amat (hc_Inferred)
# Save adjacency matrix
Adjlist_hc[[e]][[j]] <- AdjustMatrix (tarmat_s2, G_hc)
}
# Calculate the number of unique graphs inferred by each method across all permutations.
countMatrix[e, ] <- c(length (unique (Adjlist_MRPC[[e]])), length (unique (Adjlist_PC[[e]])), length (unique (Adjlist_pc.stable[[e]])), length (unique (Adjlist_mmpc[[e]])), length( unique (Adjlist_mmhc[[e]])), length( unique (Adjlist_hc[[e]])))
}
return(countMatrix)
},
stop("Model not included or missing"))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.