inst/extdata/scripts/intervention_sim_updated.R

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

library(dplyr)
library(bninfo)
library(parallel)
library(combinat)
library(stringr)
cl <- makeCluster(4)
cl <- makeCluster(16)
cl <- makeCluster(16)
#zz <- file("simulation_log.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 = " ")
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 = c("1", "2", "3"))
}
obs_data <- read.delim(d_file, header = TRUE, sep = " ")
data(isachs)
int_targets <- c("PKC", "Akt", "PKA", "PIP2", "Mek")


################################################################################
## Algo for the sim.
################################################################################
do_sim <- function(model_averaging, algo_args = NULL, base_wl = NULL){
  candidates <- int_targets
  selected <- NULL
  entropies <- NULL
  info_gains <- NULL
  model_averaging$entropy <- orientation_entropy(model_averaging)
  all_results <- list()
  #p_vals <- NULL
  while(length(candidates) > 0){
    #null_dist <- null_IG_distribution(obs_data, model_averaging, selected,
    #                                  k = 30, debug = TRUE)
    sim_results <- select_next_inhibition(obs_data, model_averaging, selected,
                                          candidates, algo_args = algo_args,
                                          base_wl = base_wl, k = 30, debug = TRUE)
    next_inh <- sim_results$top_candidate
    next_entropy <- sim_results$predicted_entropy
    next_gain <- sim_results$predicted_gain
    entropies <- c(entropies, next_entropy)
    info_gains <- c(info_gains, next_gain)
    #p_val <- right_side_pval(null_dist, next_gain)
    #p_vals <- c(p_vals, p_val)
    candidates <- setdiff(candidates, next_inh)
    selected <- c(selected, next_inh)
    all_results <- c(all_results, sim_results)
    names(all_results)[length(all_results)] <- paste0(selected, collapse = "-")
    #message("null dist when ", selected, " is selected looks like ",
    #        paste(round(head(null_dist), 4), collapse = ", "))
    #print(p_val)
    print(selected)
    if(length(candidates) == 0) break
  }
  list(selected = selected, entropies = entropies,
       info_gains = info_gains, all_results = all_results)
}
################################################################################
## Simulation with just observational data
################################################################################
cl <- makeCluster(4)
cl <- makeCluster(16)
cl <- makeCluster(16)
algo_args <- list(score = "bde", prior = "uniform", tabu = 50, whitelist = NULL, iss = 2500)
boot_vanilla <- boot.strength(obs_data, cluster = cl, algorithm = "tabu", R = 1000,
                              algorithm.args = algo_args, cpdag = FALSE)
npr_time <- system.time(no_prior_result <- do_sim(boot_vanilla, algo_args = algo_args))
save(no_prior_result, file = "no_prior_result_3-13-16_iss2500.Rdata")

################################################################################
## Simulation with known directed edges (PKC -> Raf, PKA -> Raf -> Mek -> Erk)
################################################################################
cl <- makeCluster(4)
cl <- makeCluster(16)
cl <- makeCluster(16)
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, iss = 2500)
boot_directed <- boot.strength(obs_data, cluster = cl, algorithm = "tabu", R = 1000,
                               algorithm.args = algo_args, cpdag = FALSE)
dpr_time <- system.time(directed_prior_result <- do_sim(boot_directed, algo_args, base_wl = wl))
save(directed_prior_result, file = "directed_prior_result_3-13-16_iss2500.Rdata")
################################################################################
## Simulation with all undirected edges.
################################################################################
data(tcell_examples)
wl <- arcs(moral(tcell_examples$net))
algo_args <- list(score = "bde", prior = "uniform", tabu = 50, whitelist = wl, iss = 1000)
boot_undirected <- boot.strength(obs_data, cluster = cl, algorithm = "tabu", R = 1000,
                               algorithm.args = algo_args, cpdag = FALSE)
upr_time <- system.time(undirected_prior_result  <- do_sim(boot_undirected))
save(undirected_prior_result, file = "undirected_prior_result-3-8-16_iss1000.Rdata")
################################################################################
## Simulation with all undirected edges and blacklisted missing edges
################################################################################
wl <- arcs(moral(tcell_examples$net))
bl <- t(combn(names(obs_data), 2)) %>%
  set_colnames(c("from", "to")) %>%
  arcs2names(F, T) %>%
  setdiff(arcs2names(wl, F, T)) %>%
  str_split("--") %>%
  {do.call("rbind", .)} %>%
  rbind( .[,c(2, 1)])
algo_args <- list(score = "bde", prior = "uniform", tabu = 50,
                  whitelist = wl, blacklist = bl)
boot_bl_wl <- boot.strength(obs_data, cluster = cl, algorithm = "tabu", R = 1000,
                                 algorithm.args = algo_args, cpdag = FALSE)
bl_wl_prior_result  <- do_sim(averaging_entropy(boot_bl_wl))
save(bl_wl_prior_result, file = "bl_wl_prior_result-3-4-16.Rdata")


################################################################################
## Save the data
################################################################################
results <- list(no_prior_result,
                directed_prior_result,
                undirected_prior_result,
                bl_wl_prior_result)
#save(results, file = "inst/extdata/entropy_trajectory_simulations.Rdata")
save(results, file = "inhibition_selection_simulations_3-4-16.Rdata")