R/transition_sequence.R

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

#' Get upstream edges
#'
#' Used in finding transition sequence equivalence class.
get_upstream_edges <- function(g, e){
  require(igraph)
  from_node <- igraph::V(g)[igraph::ends(g, e)[1]]
  igraph::incident(g, from_node, mode = "in")
}
#' Propagate orientation
#'
#' Callback for use in finding transition sequence equivalence class
propagate_orientation <- function(g, e){
  require(igraph)
  upstream_edges <- get_upstream_edges(g, e)
  if(length(upstream_edges) >= 1){
    # If 0 upstream incoming edges do nothing.
    # If any of the upstream incoming edges are fixed reversing the edge would introduce a V structure.
    if(any(upstream_edges$fixed)) igraph::E(g)[as.numeric(e)]$fixed <- TRUE
  }
  g
}

#' Skewed Edge Priors
#'
#' Identifies which edges in an edge-wise priors (beta argument in bnlearn's score())
#' have different prior orientation probabilities.
skewed_beta <- function(beta){
  matching_prob <- apply(beta, 1, function(row1){
    match <- apply(beta, 1, function(row2){
      same_edge <- identical(as.character(row1[c("from", "to")]),
                             as.character(row2[c("to", "from")]))
      same_prob <- row1["prob"] == row2["prob"]
      same_edge & same_prob
    })
    any(match)
  })
  !matching_prob
}

#' Transition sequence equivalent class
#'
#' Convert a dag to its transition sequence equivalence class, given a interventions.
#' If a beta argument is provided, the algorithm assumes you working with a Bayesian score that does
#' not have uniform prior.  In this case, the members of the traditional PDAG equivalence class do
#' not have the same posterior.  So to preserve posterior equivalence, edges with skewed prior
#' direction probability (probability of one direction greater than another) remain directed in the PDAG.
#' @param a dag to be converted to a tsdag
#' @param The interventions used in calculating the tsdag, defaults to NULL, which is the same as cpdag().
#' @param beta. beta is a data frame with columns from, to and prob specifying the prior probability for a set of arcs. A uniform probability distribution is assumed for the remaining arcs.  This is the same argument of score-related functions in bnlearn.
#' @export
ctsdag <- function(dag, interventions = NULL, beta = NULL){
  require(igraph)
  interventions # eval before dplyr for clear error
  # Construct an edge list of arcs belonging to amoral v-stuctures
  v_structs <-  bnlearn::vstructs(dag, moral = FALSE)
  arcs_in_vstructs <- arcs2names(rbind(v_structs[, c(1, 2)],
                                       v_structs[, c(3, 2)]))
  dag_igraph <- data.frame(bnlearn::arcs(dag), stringsAsFactors = FALSE)
  # Label the fixed arcs as those that are in v-structures, or coerced by intervention, or have skewed prior orientation probability.
  dag_igraph <-  dplyr::mutate(dag_igraph, name = arcs2names(dag_igraph),
                               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),
                               skewed_prior = FALSE)
  if(!is.null(beta)){
    skewed_arcs <- NULL
    beta_idx <- skewed_beta(beta)
    if(any(beta_idx)){
      skewed_arcs <- arcs2names(beta)[beta_idx]
      dag_igraph$skewed_prior <- ifelse(dag_igraph$name %in% skewed_arcs, TRUE, FALSE)
    }
  }
  dag_igraph <-  dplyr::mutate(dag_igraph,
                               fixed = ifelse(in_vstruct + from_targeted + to_targeted + skewed_prior > 0, TRUE, FALSE),
                               updated = ifelse(fixed, TRUE, FALSE))
  # Propogate orientation
  dag_igraph <- igraph::graph_from_data_frame(dag_igraph, directed = TRUE)
  dag_igraph <- lucy::update_edges(dag_igraph, get_upstream_edges, propagate_orientation)
  dag_arcs  <- igraph::as_edgelist(dag_igraph)
  arcs_fixed <- igraph::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 <- bnlearn::set.edge(tsdag, dag_arc[1], dag_arc[2])
    }
  }
  tsdag
}