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

Evaluating the effect of interventions

Interventions affect both the strength and direction of edges in the bootstrap inference of Bayesian networks.

Does the approach to boostrap inference affect results?

The following repeats the analysis of the Sachs 2005 T cell data used in Scutari.

isachs <- read.table("/Users/robertness/Downloads/sachs.interventional.txt", header = T)
for(i in 1:11){
  isachs[ , i] <- factor(isachs[ , i], levels = c("1", "2", "3"))
}
INT <- sapply(1:11, function(x) {which(isachs$INT == x) })
nodes <- names(isachs)[1:11]
names(INT) <- nodes
start <- random.graph(nodes = nodes, method = "melancon", 
                      num = 500, burn.in = 10^5, every = 100)
netlist <- lapply(start, function(net) {
  tabu(isachs[, 1:11], score = "mbde", exp = INT,
       iss = 1, start = net, tabu = 50)
})
intscore <- custom.strength(netlist, nodes = nodes, cpdag = FALSE)
reduced <- reduce_averaging(intscore)

This analysis used randomly generated networks as starting networks and using the custom.strength method. Notably, the way Scutari implements this, there is no resampling. Does one get similar results if using just the boot.strength method? I expect yes, and if not, it would imply a worrisome case where the nature of the greedy search affects how interventions are used. It is worrisome since there is no default method for incorporating whitelists with random start graphs. I compare using by plotting strengths in the random start net case with strengths in the case of no random start nets and resample (boot.strap). I also compare the direction values in edges with high strength.

netlist_standard <- boot.strength(isachs[, 1:11], R = 500, algorithm = "tabu",
                                  algorithm.args = list(score = "mbde", exp = INT,
                                                        iss = 1, tabu = 50)) %>%
reduced_standard <-  reduce_averaging(netlist_standard)
plot(reduced$strength, reduced_standard$strength, 
     xlab = "Strength from custom.strength", 
     ylab = "Strength from standard")

plot(reduced$direction[reduced$strength > .95], 
     reduced_standard$direction[reduced$strength > .95],
     xlab = "custom inference", 
     ylab = "standard inference",
     main = "Direction vals of edges with strength of more than .95")

These are similar enough I believe.

Can one achieve results at least as good with simulated data?

I created a function do_int that creates a targeted intervention using the mutilated network approach described in Koller.

Suppose I were to simulate data with the exact same intervention scheme as the above T cell dataset.

experiment_details <- table(isachs$INT)
names(experiment_details)[2:6] <- names(isachs)[as.numeric(names(experiment_details)[2:6])]
experiment_details

There are 1800 observational datapoints. Interventions on Mek, PIP2, AKT, PKA, and PKC have 1800, 600, 600, 600, and 1200 datapoints.

Mek, PIP2, and AKT were all inhibitors, so their levels were set low.

isachs %>%
  filter(INT == which(names(isachs) == "Mek")) %>%
  select(Mek) %>%
  table
isachs %>%
  filter(INT == which(names(isachs) == "PIP2")) %>%
  select(PIP2) %>%
  table
isachs %>%
  filter(INT == which(names(isachs) == "Akt")) %>%
  select(Akt) %>%
  table

PKA has activations, and they are not uniformly all at the highest level of 3.

isachs %>%
  filter(INT == which(names(isachs) == "PKA")) %>%
  select(PKA) %>%
  table

In PKC's 1200 observations, half are inhibition, the other half are activations, and again the activations split between level 2 and level 3.

isachs %>%
  filter(INT == which(names(isachs) == "PKC")) %>%
  select(PKC) %>%
  table

Scutari provides a fitted version of the network.

load("/Users/robertness/Downloads/sachs.rda")
graphviz.plot(bn)

Using this network, I will simulate data exactly this same intervention scheme.

do_int_named <- function(fitted_net, node, level, n){
  do_int(fitted_net, node, level, n) %>%
    mutate(INT = node)
}
isachs_new <- rbind(
      do_int_named(bn, "PKA", "1", 1), 
      do_int_named(bn, "PKA", "2", 36), 
      do_int_named(bn, "PKA", "3", 563),
      mutate(rbn(bn, 600), INT = 0),
      do_int_named(bn, "Akt", "1", 600),
      do_int_named(bn, "PKC", "1", 600),
      do_int_named(bn, "PIP2", "1", 600),
      do_int_named(bn, "Mek", "1", 600),
      mutate(rbn(bn, 600), INT = 0),
      do_int_named(bn, "PKC", "2", 86),
      do_int_named(bn, "PKC", "3", 514),
      mutate(rbn(bn, 600), INT = 0)
  ) %>%
      {.[, names(isachs)]} %>%
      mutate(INT = sapply(INT, function(item){
        if(item == 0) return(0)
        which(names(.) == item)
        }, USE.NAMES = F)) %>%
      mutate(INT = as.integer(INT))

# The same intervention
identical(isachs$INT, isachs_new$INT)

Now I repeat the custom.strength inference using this simulated data.

netlist_new <- lapply(start, function(net) {
  tabu(isachs_new[, 1:11], score = "mbde", exp = INT,
       iss = 1, start = net, tabu = 50)
})
intscore_new <- custom.strength(netlist_new, nodes = nodes, cpdag = FALSE)
reduced_new <- reduce_averaging(intscore_new)

library(ggplot)
library(tidyr)
plot_data <- data.frame(edges = arcs2names(reduced, directed_edges = F),
                        original = reduced$strength,
                        mirror = reduced_new$strength) %>%
  gather(set, strength, -edges)
centropy_data <- data.frame(edges = arcs2names(reduced, directed_edges = F),
                            original = orientation_entropy(reduced),
                            mirror = orientation_entropy(reduced_new)) %>%
  gather(set, c_entropy, -edges)
plot_data$c_entropy <- centropy_data$c_entropy 
plot_data <- mutate(plot_data, 
                    validation = ifelse(edges %in% arcs2names(arcs(bn), 
                                        directed_edges = F), "TP", "FP"),
                    index = rep(order(order(strength[set == "original"])), 2))

First evaluating strength results, in this case I don't expect the strengths to be the same. One expects the greater faithfulness of the simulated intervention data to improve in detection of true edges and elimination of poor edges. The following plots order the strength values, forming the S curve. This means a stronger S curve in the simulated data.

plot_data %>%
  filter(set == "original") %>%
  ggplot(aes(x = index, y = strength, colour = validation)) +
  geom_point() +
  ggtitle("Original Scutari Dataset")

Below is the same plot with the simlated data, preserving the same edge index ordering of the previous plot.

plot_data %>%
  filter(set == "mirror") %>%
  ggplot(aes(x = index, y = strength, colour = validation)) +
  geom_point() +
  ggtitle("Simulated Mirror Dataset")

With orientation entropy, we expect low probability edges to have low entropy. It would be nice if true positives had less entropy than false positives. We expect the simulated data to have less entropy than the real data.

plot_data %>%
  filter(set == "original") %>%
  ggplot(aes(x = index, y = c_entropy, colour = validation)) +
  geom_point() +
  ylim(c(0, 1)) +
  ggtitle("Original Scutari Dataset")

This is a nice plot that shows that entropy increases in strength (since it is conditioned on strength). The interesting comparison comes from comparing to the simulated dataset.

plot_data %>%
  filter(set == "mirror") %>%
  ggplot(aes(x = index, y = c_entropy, colour = validation)) +
  geom_point() +
  ylim(c(0, 1)) +
  ggtitle("Simulated Mirror Dataset")

The high value comes from PIP3--Plcg, which is an edge that is not affected by the interventions. That is a bonus result that I did not expect to see.

I expect the simulated dataset to be a better fit than the original.

bn_net <- construct_bn(arcs(bn))
score(bn_net, isachs[, 1:11], type = "mbde", exp = INT)

score(bn_net, isachs_new[, 1:11], type = "mbde", exp = INT)

That is not true. That is troubling.

Can I get more resolution in strength values with "perfect" data?

I can simulate a "perfect"" dataset by making sure that each node has a representative number of samples in each one of its states.

To understand why this might be neccessary, we examine the observational values of the sachs data.

isachs %>%
  filter(INT == 0) %>%
  {lapply(.[, 1:11], table)} %>%
  {do.call("rbind", .)}

Plcg, PIP2, PKC, P38, and Jnk, are missing level 3 values in this case. This suprises me, I would have expected level 1 values to be missing. Even with the nodes where level 3 is present, it is usually a small number. Normalization biase?

Continuing, using the rbn function I can simulate 1800 observations and compare to the above table. Looking at the fitted model, are their any 0 probabilities?

#lapply(bn, function(item) item$prob)

Not running that code because the output is ugly, so the answer is yes. So I expect simulated data will have a similar emperical distribution as the results above, including many zeros.

rbn(bn, 1800) %>%
  lapply(table) %>%
  {do.call("rbind", .)} %>%
  {.[names(isachs[1:11]), ]}

It is completely different. So I expect there is some serious issues with faithfulness of the original observational data, with respect to the fitted network. It implies that you needed the interventions to detect edges as well as directions. To prove this I compare the rbn inference to the observational inference. I expect more false positives, and more entropy in strength.

library(parallel)
cl <-  makeCluster(4)
original_obs <- filter(isachs, INT == 0)[, 1:11] %>%
  boot.strength( R = 500, algorithm = "tabu", cluster = cl, cpdag = TRUE,
                algorithm.args = list(score = "bde", iss = 1, tabu = 50))
rbn_obs <- boot.strength(rbn(bn, 1800), R = 500, algorithm = "tabu", cluster = cl, cpdag = TRUE, 
                algorithm.args = list(score = "bde", iss = 1, tabu = 50))

plot_data <- data.frame(edges = arcs2names(original_obs, directed_edges = F),
                        original = original_obs$strength,
                        mirror = rbn_obs$strength) %>%
  gather(set, strength, -edges) %>%
  mutate(validation = ifelse(edges %in% arcs2names(arcs(bn), 
                                        directed_edges = F), "T", "F"),
         index = rep(order(order(strength[set == "original"])), 2))

plot_data %>%
  filter(set == "original") %>%
  ggplot(aes(x = index, y = strength, colour = validation)) +
  geom_point() +
  ggtitle("Original Observational Data")

It would seem to me that entropy in strength is not much of a problem. Indeed I know that by increasing the imaginary sample size, I can increase the amout of entropy. Perhaps, Karen is wrong. Perhaps it is follow to use interventions to try to resolve the presence of an edge since you can artificially increase the probability of a false edge merely by making your network inference method more sensitive, in this case by increasing the weight of a uniform prior. Contrary to what I expected, I only see one false positive. Rather I see many false negatives. Again, this is a problem of sensitivity, not entropy. In the next graph, I still expect strength to decrease for false positives (entropy reduction), but I think the real impact of the more faithful data will be to decrease the amount of false negatives.

plot_data %>%
  filter(set == "mirror") %>%
  ggplot(aes(x = index, y = strength, colour = validation)) +
  geom_point() +
  ggtitle("Simulated Mirror Dataset")

This is a tremendous shock. Strength entropy certainly went down, but false positives went up, as did false negatives. Perhaps the observational data was too ideal. I believe I recall using other observational data from this overall dataset and it didn't do as well as this particular set of observational data.

Why am I getting so much seperation in this simulated data? Trying a greater iss, and smaller sample size.

library(parallel)
cl <-  makeCluster(4)
small_obs <- boot.strength(rbn(bn, 100), R = 500, 
                           algorithm = "tabu", cluster = cl, cpdag = TRUE,
                algorithm.args = list(score = "bde", iss = 1, tabu = 50))
large_obs <- boot.strength(rbn(bn, 500), R = 500, 
                           algorithm = "tabu", cluster = cl, cpdag = TRUE, 
                algorithm.args = list(score = "bde", iss = 1, tabu = 50))
plot_data <- data.frame(edges = arcs2names(small_obs, directed_edges = F),
                        small = small_obs$strength,
                        large = large_obs$strength) %>%
  gather(set, strength, -edges) %>%
  mutate(validation = ifelse(edges %in% arcs2names(arcs(bn), 
                                        directed_edges = F), "T", "F"),
         index = rep(order(order(strength[set == "small"])), 2))

Hopefully, more strength entropy in small sample size.

plot_data %>%
  filter(set == "small") %>%
  ggplot(aes(x = index, y = strength, colour = validation)) +
  geom_point() +
  ggtitle("Small Simulated Dataset")

Less in larger sample size.

plot_data %>%
  filter(set == "large") %>%
  ggplot(aes(x = index, y = strength, colour = validation)) +
  geom_point() +
  ggtitle("Large Simulated Dataset")

We've seen this in my previous paper. More samples decreases entropy in strength, but only if the sampling of the space is full -- there are nonzero counts for each state of each node.

If there were non-zero counts for certain states, how could a targeted intervention play a role?

score(construct_bn(arcs(bn)), filter(isachs, , type = "bic") ideal_observational_data <- function(fitted_net, n){ lapply(fitted_net, function(item){ dimnames(item$prob)[[1]] %>% # item_levels lapply(function(item_level){ ev_list <- structure(list(item_level), names = item$node) sampling <- cpdist(fitted_net, names(fitted_net), evidence = ev_list, method= "lw") %>% sample_n(n) }) %>% {do.call("rbind", .)} %>% set_rownames(NULL) }) %>% {do.call("rbind", .)} %>% set_rownames(NULL) }

ideal_intervention_data <- function(fitted_net, n){ lapply(fitted_net, function(item){ dimnames(item$prob)[[1]] %>% # item_levels lapply(function(level) mutate(do_int(fitted_net, item$node, level, n), INT = node) %>% {do.call("rbind", .)} %>% set_rownames(NULL) }) %>% {do.call("rbind", .)} %>% set_rownames(NULL) } ```

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