inst/extdata/scripts/entropy_trajectory_simulation_no_bayes.R

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

# In this analysis, I generate model averaging entropies for the isachs data.
# I generate entropy in the observation data (no inhibitions), then
# iteratively add data in and getting new bootstrap results, and using those
# results to calculate the prior for the next iteration.
# Then I apply the proposed method to predict.
library(dplyr)
library(bninfo)
library(parallel)
library(combinat)
library(stringr)
cl <- makeCluster(4)
cl <- makeCluster(16)
cl <- makeCluster(40)
#zz <- file("entropy_sim_log_2-8-16.txt", open = "wt")
#sink(zz, type = "message")
d_file <- system.file("extdata", "sachs.discretised.txt", package = "bninfo")
obs_data <- read.delim(d_file, header = TRUE,sep = " ")
data(isachs)
node_names <- names(isachs)[1:11]
levels(obs_data$PIP2) <- c("1", "3", "2")
levels(obs_data$PKC) <- c("1", "3", "2")
for(node in setdiff(node_names, c("PIP2", "PKC"))){
  levels(obs_data[, node]) <- c("1", "2", "3")
}
for(node in node_names){
  obs_data[, node] <- factor(obs_data[, node], levels = sort(levels(obs_data[, node])))
}
int_targets <- c("PKC", "Akt", "PKA", "PIP2", "Mek")
#################################################################
## All permutations of the interventions
#################################################################
permutations <- permn(int_targets)
#################################################################
## This is to create all the combinations of interventions
#################################################################
cbn <- NULL
for(i in 1:(length(int_targets) - 1)){
  cbn_list_names <- combn(int_targets, i) %>%
    apply(., 2, list) %>%
    lapply(unlist) %>%
    sapply(function(item) paste0(sort(item), collapse="-"))
  cbn_list <- as.list(rep(NA, length(cbn_list_names)))
  names(cbn_list) <- cbn_list_names
  cbn <- c(cbn, cbn_list)
}
cbn <- c(cbn, list(NA))
names(cbn)[length(cbn)] <- paste0(sort(int_targets), collapse = "-")

################################################################################
## Algo for calculating the ordering
################################################################################
.new_ordering <- function(ordering, starting_boot, wl = NULL, bl = NULL){
  entropies <- mean(orientation_entropy(starting_boot))
  for(i in 1:length(ordering)){
    sub_order_name <- paste0(sort(ordering[1:i]), collapse = "-")
    print(sub_order_name)
    # To avoid repeats, only compute new entropies. This should take time
    prior_entropy <- cbn[[sub_order_name]]
    if(is.na(prior_entropy)){
      int_dex <- c(which(node_names %in% ordering[1:i]))
      int_data <- filter(isachs, INT %in% int_dex) %>%
        rbind(mutate(sample_n(obs_data, nrow(obs_data), replace = TRUE), INT = NA)) # Add some of original data, resampled.
      INT <- lapply(1:11, function(x) {which(int_data$INT == x)})
      names(INT) <- node_names
      boot_i <- boot.strength(int_data[, 1:11], cl, algorithm = "tabu",
                              algorithm.args = list(score = "mbde",
                                                    iss = 2500,
                                                    exp = INT,
                                                    whitelist = wl,
                                                    blacklist = bl,
                                                    tabu = 50), cpdag = FALSE)
      causal_entropy_i <- mean(orientation_entropy(boot_i))
      cbn[[sub_order_name]] <<- causal_entropy_i
    } else {
      causal_entropy_i <- prior_entropy
    }
    entropies <- c(entropies, causal_entropy_i)
  }
  names(entropies) <- c("obs", ordering)
  entropies
}

new_ordering <- function(ordering, starting_boot, wl = NULL, bl = NULL){
  # This little hack should send errors out as messages,
  # and just retry the algorithm until successful rather than
  # stopping.
  entropies <- "error"
  while(entropies == "error"){
    entropies <- tryCatch(.new_ordering(ordering, starting_boot,
                                        wl = wl, bl = bl),
                          error = function(e){
                            message(e)
                            "error"
                          })
  }
  entropies
}
################################################################################
## Simulation with just observational data
################################################################################
cl <- makeCluster(4)
cl <- makeCluster(16)
cl <- makeCluster(40)
algo_args <- list(score = "bde", prior = "uniform", tabu = 50, whitelist = NULL)
boot_vanilla <- boot.strength(obs_data, cluster = cl, R = 500,
                                algorithm = "tabu",
                              algorithm.args = c(algo_args, list(iss = 1500)), cpdag = TRUE)
averaging_entropies_vanilla <- lapply(permutations, .new_ordering,
                                      starting_boot = boot_vanilla)
save(averaging_entropies_vanilla, file = "averaging_entropies_vanilla_3-12-16_no_bayes_2500iss.Rdata")
permutations1 <- permutations
################################################################################
## Simulation with known directed edges (PKC -> Raf, PKA -> Raf -> Mek -> Erk)
################################################################################
cl <- makeCluster(4)
cl <- makeCluster(16)
cl <- makeCluster(40)
for(i in 1:length(cbn)) cbn[[i]] <- NA
## All permutations of the interventions
permutations <- permn(int_targets)
wl <- matrix(c("PKC", "Raf",
               "PKA", "Raf",
               "Raf", "Mek",
               "Mek", "Erk"),
             ncol = 2, byrow = T)
algo_args <- list(score = "bde", prior = "uniform", tabu = 50, whitelist = wl)
boot_directed <- boot.strength(obs_data, cluster = cl, algorithm = "tabu", R = 1000,
                               algorithm.args = c(algo_args, list(iss = 2500)), cpdag = FALSE)
averaging_entropies_directed <- lapply(permutations, .new_ordering,
                                       starting_boot = boot_directed, wl = wl)
save(averaging_entropies_directed,
     file = "averaging_entropies_directed_3-12-16_no_bayes_2500iss.Rdata")
################################################################################
## Simulation with all undirected edges.
################################################################################
data(tcell_examples)
true_arcs <- arcs(tcell_examples$net)
wl <- rbind(true_arcs, true_arcs[, c(2, 1)]
true_arc_names <- arcs2names(true_arcs)
all_arcs <- boot.strength(obs_data, algorithm = "hc", R = 1) %>%
  select(from, to)
false_arc_names <- setdiff(arcs2names(boot), true_arc_names)
false_arcs <- all_arcs %>%
  mutate(arc_name = arcs2names(boot)) %>%
  filter(arc_name %in% false_arc_names) %>%
  select(from, to)
cl <- makeCluster(4)
cl <- makeCluster(16)
cl <- makeCluster(40)
for(i in 1:length(cbn)) cbn[[i]] <- NA
algo_args <- list(score = "bde", prior = "uniform", tabu = 50, whitelist = wl, blacklist = bl)
boot_undirected <- boot.strength(obs_data, cluster = cl, algorithm = "tabu", R = 1000,
                               algorithm.args = algo_args, cpdag = FALSE)
averaging_entropies_undirected <- lapply(permutations, .new_ordering,
                                         starting_boot = boot_undirected, wl = wl, bl = bl)
save(averaging_entropies_undirected,
     file = "averaging_entropies_undirected_3-8-16_no_bayes_1000iss.Rdata")
################################################################################
## Simulation with all blacklisted missing edges
################################################################################
# bl0 <- t(combn(names(obs_data), 2)) %>%
#   set_colnames(c("from", "to")) %>%
#   arcs2names(F, T) %>%
#   setdiff(arcs2names(wl, F, T)) %>%
#   str_split("--") %>%
#   {do.call("rbind", .)}
# set.seed(10)
# bl <- bl0[sample(1:nrow(bl0), 19), ]
# bl <- rbind(bl, bl[,c(2, 1)])
# algo_args <- list(score = "bde", prior = "uniform", tabu = 50,
#                   blacklist = bl)
# boot_bl <- boot.strength(obs_data, cluster = cl, algorithm = "tabu",
#                                algorithm.args = algo_args, cpdag = FALSE)
# averaging_entropies_with_bl <- lapply(permutations, new_ordering,
#                                       starting_boot = boot_bl, bl = bl)
# save(averaging_entropies_with_bl,
#      file = "averaging_entropies_with_bl_2-8-16.Rdata")
################################################################################
## Simulation with all undirected edges and blacklisted missing edges
################################################################################
# wl <- arcs(moral(tcell_examples$net))
# algo_args <- list(score = "bde", prior = "uniform", tabu = 50,
#                   whitelist = wl, blacklist = bl)
# boot_undirected_bl <- boot.strength(obs_data, cluster = cl, algorithm = "tabu",
#                                algorithm.args = algo_args, cpdag = FALSE)
# averaging_entropies_undirected_bl <- lapply(permutations, new_ordering,
#                                          starting_boot = boot_undirected_bl,
#                                          wl = wl, bl = bl)
# save(averaging_entropies_undirected_bl,
#      file = "averaging_entropies_undirected_bl_2-8-16.Rdata")
################################################################################
## Save the data
################################################################################
results <- list(vanilla = averaging_entropies_vanilla,
                directed = averaging_entropies_directed,
                undirected = averaging_entropies_undirected)
save(results, file = "entropy_trajectory_simulations_3-8-16_no_bayes_1000iss.Rdata")