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

library(bninfo)

Bootstrap inference of Bayesian network structure provides enables more robust inference, and enables a edge-wise view of inference in terms of estimates of probability that an edge is present (presence probability), and conditional probability that an edge is oriented in a given direction conditional on its presence (orientation probability.). The approach uses resampled data, random starting networks, and sometimes stochastic simulation introduce randomness into the network seach, resuling in generation of a diverse set of high-scoring networks. If using a scoring metric such as Bayesian Dirichlet equivalent score, then this approach is equivalent to sampling from the posterior distribution of network structures.

Cooper GF, Yoo C (1999) introduced the modified Bayesian Dirichlet equivalent score, which adjusts the BDe calculation to account for targeted interventions, which fix nodes a specific values. This enables the learning of causal networks form data.

Problem: When sample size is high, search algorithms, particularly greedy search algorithms converge on very similar network structures, resulting in poor sampling from the posterior of network structures. This can be somewhat remedied by converting each search result to a partially-directed network which represents the entire equivalence class of networks. However, in the case interventions are used, converting to the equivalence class losses the causal information gained by interventions, because edges oriented by interventions are converted to undirected edges.

Solution: Create a function that converts directed networks to a TS-equivalence class. Two networks G1 and G2 are TS-equivalent if and only if they have the same structure ignoring arc directions, the same set of v-structures, and the same sets of parents for all manipulated variable. So given a DAG learned with an mBDe score, and a target node used in mBDe, then edges are assumed fixed and cannot be converted to undirected edges.

My approach to this is to create PDAG, attempt to coerce edges connected to targeted nodes back to what they were.

I'll start with a tree to avoid v-structures.

set.seed(10)
dag <- make_tree(10, 3) %>%
  igraph.to.graphNEL %>%
  as.bn
nodes(dag) <- paste0("V", nodes(dag))
graphviz.plot(dag)

Generate the PDAG

pdag <- cpdag(dag, moral = F)
graphviz.plot(pdag)

Which edges are undirected?

undirected.arcs(pdag)

In bnlearn, set.arc compells arc direction. Lets give that a shot for node 3.

dag_arcs <- arcs(dag)
incoming_arcs <- dag_arcs[dag_arcs[, 2] == "V3", , drop = FALSE]
outgoing_arcs <- dag_arcs[dag_arcs[, 1] == "V3", , drop = FALSE]
change_arcs <- rbind(incoming_arcs, outgoing_arcs)
tsdag <- pdag
for(i in 1:nrow(change_arcs)){
  tsdag <- set.arc(tsdag, change_arcs[i, 1], change_arcs[i, 2])
}
graphviz.plot(tsdag)

In the tsdag object, the incoming and outgoing arcs to V3 in dag have replaced the undirected arcs in pdag.

Testing that we can still convert this back into an instance of the equivalence class:

graphviz.plot(cextend(tsdag))

Test cases

Next, some test cases.

base <- empty.graph(LETTERS[1:5])
net1 <- `modelstring<-`(base, value = "[A|C][B|C][D|C][E|D][C]") 
net2 <- `modelstring<-`(base, value = "[A|C][B|C][C|D][E|D][D]") 
net3 <- `modelstring<-`(base, value = "[A|C][B|C][C|D][D|E][E]") 
net4 <- `modelstring<-`(base, value = "[A|C][C|B][D|C][E|D][B]") 
net5 <- `modelstring<-`(base, value = "[C|A][B|C][D|C][E|D][A]") 
graphviz.plot(net1, main = "net1")
graphviz.plot(net2, main = "net2")
graphviz.plot(net3, main = "net3")
graphviz.plot(net4, main = "net4")
graphviz.plot(net5, main = "net5")

These are each of the same equivalence class.

all.equal(cpdag(net1, F), cpdag(net2, F))
all.equal(cpdag(net2, F), cpdag(net3, F))
all.equal(cpdag(net3, F), cpdag(net4, F))
all.equal(cpdag(net4, F), cpdag(net5, F))

• An intervention on C in net1 should not change anything, because D -> E is coerced otherwise it would introduce a v-structure.
• Similar results are true for net4 and net5
• An intervention on C in net2 and net3 produce the same tsdag, with E-D undirected

tsdag_23_C <- base
arcs(tsdag_23_C) <- matrix(c("C", "A",
                             "C", "B",
                             "D", "C",
                             "E", "D",
                             "D", "E"), byrow = T, ncol = 2)
graphviz.plot(tsdag_23_C)

• An intervention in A in net5 should produce the same graph, again because of avoidance of v-structures.
• An intervention in A for the other nets should produce undirected edges in all but the incoming edge to A.

tsdag_1234_A <- base
arcs(tsdag_1234_A) <- matrix(c("C", "A",
                               "C", "B",
                               "B", "C",
                               "D", "C",
                               "C", "D",
                               "E", "D",
                               "D", "E"), byrow = T, ncol = 2)
graphviz.plot(tsdag_1234_A)

With these test cases, I implement a function.

ctsdag <- function(dag, interventions){
  dag_arcs <- arcs(dag)
  pdag <- cpdag(dag, moral = FALSE)
  for(int in interventions){
    incoming_arcs <- dag_arcs[dag_arcs[, 2] == int, , drop = FALSE]
    outgoing_arcs <- dag_arcs[dag_arcs[, 1] == int, , drop = FALSE]
    change_arcs <- rbind(incoming_arcs, outgoing_arcs)
    tsdag <- pdag
    for(i in 1:nrow(change_arcs)){
      tsdag <- set.arc(tsdag, change_arcs[i, 1], change_arcs[i, 2])
    }
  }
  tsdag
}

Applying the function to the test cases:

all.equal(ctsdag(net1, "C"), net1)

Failed! The reason is that by converting to the pdag first, some edges become undirected when they should stay directed because otherwise it would induce a v-structure.

So we need another way. Trying this algorithm:

• annotate each edges connected to targeted nodes as fixed.
• annotate edges involved in v-structures as fixed.
• For the remaining edges, if the from node is the to node in a fixed edge, reversing it would induce a v-structure, so annotate these nodes as fixed.

irreversable <- function(from_node, to_nodes){
  (from_node %in% to_nodes)
}
ctsdag <- function(dag, interventions){
  arcs_in_vstructs <- dag %>%
    vstructs %>%
    {rbind(.[, c(1, 2)], .[, c(3, 2)])} %>%
    arcs2names
  dag_arcs <- arcs(dag) %>%
    data.frame(stringsAsFactors = FALSE) %>%
    mutate(arcs = arcs2names(.),
           vstruct = ifelse(arcs %in% arcs_in_vstructs, TRUE, FALSE),
           target = ifelse(from %in% interventions, TRUE, FALSE),
           target = ifelse(to %in% interventions, TRUE, target),
           fixed = ifelse(vstruct + target > 0, TRUE, FALSE))
  to_nodes <- dag_arcs %>%
    filter(fixed) %$%
    to
  dag_arcs <- mutate(dag_arcs, fixed = ifelse(irreversable(from, to_nodes), TRUE, fixed),
                     variable = !fixed)
  tsdag <- dag
  for(i in 1:nrow(dag_arcs)){
    if(dag_arcs[i, "variable"]){
      dag_arc <- as.character(dag_arcs[i, c(1, 2)])
      tsdag <- set.edge(tsdag, dag_arc[1], dag_arc[2])
    }
  }
  tsdag
}

Trying again, applying the function to the test cases:

all.equal(ctsdag(net1, "C"), net1 )
all.equal(ctsdag(net4, "C"), net4)
all.equal(ctsdag(net5, "C"), net5)
all.equal(ctsdag(net2, "C"), tsdag_23_C)
all.equal(ctsdag(net3, "C"), tsdag_23_C)
all.equal(ctsdag(net1, "A"), tsdag_1234_A)
all.equal(ctsdag(net2, "A"), tsdag_1234_A)
all.equal(ctsdag(net3, "A"), tsdag_1234_A)
all.equal(ctsdag(net4, "A"), tsdag_1234_A)
all.equal(ctsdag(net5, "A"), net5) #!!!!

All but net5 works. This is because the ifelse approach I use doesn't account for the cascade of orientations. I need to think recursively.

Lets imagine a more extreme version of net5;

base <- empty.graph(LETTERS[1:8])
net_long <- `modelstring<-`(base, value = "[A][C|A][B|C][D|C][E|D][F|E][G|F][H|G]") 
graphviz.plot(net_long)

I need to think in terms of how orientation of an edge cascades down this change. A->C is fixed. If A->C is fixed, C->D can't be reversed, then D->E can't be reversed, etc. The stopping criteria of the recursion can be no more children or hitting another fixed edge.

How might this look in with propagation in Lucy? Lucy's propagation works by landing on a given edge and applying a callback. The callback depends on the states of determiner edges. In this case, I believe I'm just looking at the next upstream edge. If I am on edge F->G and edge E->F is updated and fixed, then I fix E->F and set its state to updated. If I am on edge F->G and edge E->F is updated but not fixed, then I set E->F as variable and set its state to updated.

get_upstream_edges <- function(g, e){
  from_node <- V(g)[ends(g, e)[1]]
  E(g)[to(from_node)]
}
propagate_orientation <- function(g, e){
  upstream_edge <- get_upstream_edges(g, e)
  if(length(upstream_edge) == 1){ 
    # If 0 upstream edges do nothing.  
    # There shouldn't be more than 1, because that would be a v-structure.
    if(upstream_edge$fixed) E(g)[as.numeric(e)]$fixed <- TRUE
  }
  g
}
ctsdag <- function(dag, interventions){
  interventions # eval before dplyr for clear error
  arcs_in_vstructs <- dag %>%
    vstructs %>%
    {rbind(.[, c(1, 2)], .[, c(3, 2)])} %>%
    arcs2names
  dag_igraph <- arcs(dag) %>%
    data.frame(stringsAsFactors = FALSE) %>%
    mutate(name = arcs2names(.),
           in_vstruct = ifelse(name %in% arcs_in_vstructs, TRUE, FALSE),
           from_targeted = ifelse(from %in% interventions, TRUE, FALSE),
           to_targeted = ifelse(to %in% interventions, TRUE, FALSE),
           fixed = ifelse(in_vstruct + from_targeted + to_targeted > 0, TRUE, FALSE),
           updated = ifelse(fixed, TRUE, FALSE)) %>%
    graph_from_data_frame(directed = TRUE) 
  dag_igraph <- update_edges(dag_igraph, get_upstream_edges, propagate_orientation)
  dag_arcs  <- as_edgelist(dag_igraph)
  arcs_fixed <- E(dag_igraph)$fixed
  tsdag <- dag
  for(i in 1:nrow(dag_arcs)){
    if(!arcs_fixed[i]){
      dag_arc <- as.character(dag_arcs[i, c(1, 2)])
      tsdag <- set.edge(tsdag, dag_arc[1], dag_arc[2])
    }
  }
  tsdag
}

Trying once again, applying the function to the test cases:

all.equal(ctsdag(net1, "C"), net1 )
all.equal(ctsdag(net4, "C"), net4)
all.equal(ctsdag(net5, "C"), net5)
all.equal(ctsdag(net2, "C"), tsdag_23_C)
all.equal(ctsdag(net3, "C"), tsdag_23_C)
all.equal(ctsdag(net1, "A"), tsdag_1234_A)
all.equal(ctsdag(net2, "A"), tsdag_1234_A)
all.equal(ctsdag(net3, "A"), tsdag_1234_A)
all.equal(ctsdag(net4, "A"), tsdag_1234_A)
all.equal(ctsdag(net5, "A"), net5) 

Success!

Counting with custom.strength

Next, I want to make sure bnlearn's custom.strength method can count undirected edges.

Looking at these two test networks:

graphviz.plot(tsdag_23_C)
graphviz.plot(tsdag_1234_A)

The strengths should be 1 for each edge. Both have D-E undirected. tsdag-1234-A has C-D and B-C undirected, while tsdag-23-C has D->C and C->B. So D-E direction should be .5, D->C, and C->B should be .75.

custom.strength(list(tsdag_23_C, tsdag_1234_A), nodes = nodes(tsdag_1234_A), cpdag = F) %>%
  filter(direction >= .5)

Perfect. So we can get a distribution of edge probabilities given the same intervention on all DAGs of a single equivalence class.

equivalence_class <- lapply(list(net1, net2, net3, net4, net5), ctsdag, interventions = "D")
custom.strength(equivalence_class, nodes = nodes(net1), cpdag = F) %>%
  filter(direction >= .5)

robertness/bninfo documentation built on May 27, 2019, 10:32 a.m.