SEset2.R

library("SEset")
library("qgraph")
library("pcalg")

# For true DAG:
A <- matrix(c(
  0,0.25,0.25,
  0,0,0,
  0,0,0
),3,3,byrow=TRUE)
S <- diag(3)

# True DAG:
graph_dag <- qgraph(t(A), layout = "circle", theme = "colorblind", title = "True DAG", asize = 10,
       labels = c("A","B","C"), vsize = 20, edge.label.cex = 2, edge.labels = TRUE)

# True precision:
I <- diag(3)
omega <- t(I - A) %*% (I - A)

# True network:
pcor <- wi2net(omega)
graph_ggm <- qgraph(pcor, layout = "circle", theme = "colorblind", title = "True GGM",
       labels = c("A","B","C"), edge.labels = TRUE, vsize = 20, edge.label.cex = 2)

# pc algorithm:
pc <- pc(suffStat = list(C = solve(omega), n = 99999), indepTest = gaussCItest,
         labels = c("A","B","C"), alpha = 0.05)
graph_pc <- qgraph(pc, layout = "circle", theme = "colorblind", title = "pcalg result", asize = 10,
                   labels = c("A","B","C"), vsize = 20, edge.label.cex = 2, esize = 5,
                   unCol = "black")

# SE set:
SE <- precision_to_SEset(omega, digits = 2, rm_duplicates = TRUE)

# Directed edge frequency
propd <- propcal(SE, rm_duplicate = FALSE, directed = TRUE)
graph_se <- qgraph(propd, layout = "circle", theme = "colorblind", title = "SE frequency", asize = 10,
       parallelEdge = TRUE, edge.labels = TRUE,
       labels = c("A","B","C"), posCol = "black", vsize = 20, edge.label.cex = 2)

# Make output file:
pdf("SEset2.pdf",width = 8, height = 8)
layout(matrix(1:4,2,2,byrow=TRUE))
plot(graph_dag)
box("figure")
plot(graph_ggm)
box("figure")
plot(graph_pc)
box("figure")
plot(graph_se)
box("figure")
dev.off()