inst/extdata/scripts/equivalence_class_enumeration.R

In robertness/bninfo: Queries and information theoritic operations on Bayesian networks.

library(bnlearn)
node_names <- c("IGFR", "IRS", "PI3K", "GRB", "RAS", "MEK", "ERK", "AKT", "MTOR")
net <- empty.graph(nodes = node_names)
arcs(net) <- matrix(c("IGFR", "IRS",
                      "IRS", "PI3K",
                      "IRS", "GRB",
                      "GRB", "RAS",
                      "PI3K", "AKT",
                      "AKT", "MTOR",
                      "AKT", "RAS",
                      "RAS", "MEK",
                      "MEK", "ERK"), ncol = 2, byrow = T)
library(bninfo)
equiv_class <- cpdag(net, moral = FALSE)
undirected_arcs <- arcs(equiv_class) %>%
  apply(1, function(r) paste(sort(r), collapse = "--")) %>%
  {.[which(duplicated(.))]}
net_df <- arcs(net) %>%
  data.frame %>%
  mutate(edge_names = arcs2names(., directed_edges = F),
         undirected = ifelse(edge_names %in% undirected_arcs, TRUE, FALSE))

m <- expand.grid(rep(list(0:1), sum(net_df$undirected)))
arc_list <- list()
for(i in 1:nrow(m)){
  binary <- m[i, ]
  new_arcs <- net_df %>%
    mutate(change = 0) %>%
    {.$change[.$undirected] <- binary; .} %>%
    apply(1, function(r){
      if(r$change == 1) return(c(r$to, r$from))
      return(c(r$from, r$to))
    }) %>%
    t
  arc_list <- c(arc_list, list(new_arcs))
}
# Though uneccessary in this case, here we could put a filter for DAGs
net_boot <- lapply(arc_list, function(arc_set){
  net <- empty.graph(nodes = node_names)
  arcs(net) <- arc_set
  net
}) %>%
  custom.strength(nodes = node_names, cpdag = FALSE)