inst/extdata/scripts/alarm_analysis_1.R

# Repeating the experiments on alarm data
library(bninfo)
library(dplyr)
library(magrittr)
data(alarm)
####################################################
# Use alarm because it can have a causal interpretation.
# But take it down to 9 interconnected nodes so analysis is not
# plagued by dimensionality.  Also, limit the amount of V-structures.
####################################################
full_net <- empty.graph(names(alarm))
modelstring(full_net) <- paste("[HIST|LVF][CVP|LVV][PCWP|LVV][HYP][LVV|HYP:LVF]",
                         "[LVF][STKV|HYP:LVF][ERLO][HRBP|ERLO:HR][HREK|ERCA:HR][ERCA]",
                         "[HRSA|ERCA:HR][ANES][APL][TPR|APL][ECO2|ACO2:VLNG][KINK]",
                         "[MINV|INT:VLNG][FIO2][PVS|FIO2:VALV][SAO2|PVS:SHNT][PAP|PMB][PMB]",
                         "[SHNT|INT:PMB][INT][PRSS|INT:KINK:VTUB][DISC][MVS][VMCH|MVS]",
                         "[VTUB|DISC:VMCH][VLNG|INT:KINK:VTUB][VALV|INT:VLNG][ACO2|VALV]",
                         "[CCHL|ACO2:ANES:SAO2:TPR][HR|CCHL][CO|HR:STKV][BP|CO:TPR]", sep = "")
.subset <- c("PMB", "PAP", "SHNT", "INT", "MINV", "VLNG", "KINK", "PRSS", "VTUB")
ground_truth <- subgraph(full_net, .subset)
####################################################
# Use simulated from the network structure so data is
# definately faithful to the graph structure.
####################################################
.data <- bn.fit(ground_truth, alarm[, .subset], "bayes") %>%
  rbn(100000)


###############################################################
# Make an isachs-like dataset with faithful interventions.
###############################################################
obs <- cbind(rbn(ground_truth, 500, .data, fit = "bayes"), intervention = 0)
ialarm <- lapply(names(.data), function(node){
  level <- levels(.data$node)[1]
  ground_truth %>%
    mutilated(evidence = structure(list(level), names = node)) %>%
    rbn(500, data = .data,fit = "bayes") %>%
    mutate(intervention = which(names(.) == node))
  }) %>%
  {do.call("rbind", .)} %>%
  {rbind(., obs)}


###############################################################
# Check the strength of edges in the bootstrap results
# There should be less entropy in the model averaging results.
###############################################################
# Set cpdag to FALSE in case V-structures provide causal direction
# information

boot <- boot.strength(.data, algorithm = "hc",
                      algorithm.args = list(score = "bde"),
                      cpdag = FALSE)

boot %>%
  select(from, to, strength) %>%
  reduce_averaging

##################
# Ignoring direction, true edges all have high strength in
# model averaging results.
# Few or no false positive edges given bootstrap threshold of .53.

# 1   PMB  PAP    1.000
# 2   PMB SHNT    1.000
# 16 SHNT  INT    1.000
# 22  INT MINV    1.000
# 25  INT PRSS    0.885
# 23  INT VLNG    0.885
# 34 KINK PRSS    0.945
# 33 VLNG VTUB    0.985
# 36 PRSS VTUB    1.000
# 22  INT MINV    1.000
# 27 MINV VLNG    1.000

##################
# The only real false negative, given the threshold being
# VLNG to KINK.  It's close anyway.

# 31 VLNG KINK    0.480

##################
# Not too many false negatives.

# from   to strength
# 3   PMB  INT    0.985
# 9   PAP SHNT    0.530
# 13  PAP KINK    0.640
# 26  INT VTUB    0.660
# 32 VLNG PRSS    0.125
# 35 KINK VTUB    0.785


###############################################################
# Sanity check to make sure mbde works with CS prior.
# There should be less entropy in the model averaging results.
###############################################################

mean(averaging_entropy(boot)$entropy)
# initial entropy is .6651

beta <- boot %>%
  mutate(prob = strength * direction) %>%
  select(from, to, prob)

INT <- lapply(1:9, function(x) {which(ialarm$intervention == x)})
names(INT) <- names(ialarm)[1:9]
library(parallel)
cl <- makeCluster(4)
boot_mbde <- boot.strength(ialarm[, 1:9], cl, algorithm = "tabu",
                        algorithm.args = list(score = "mbde",
                                              exp = INT, prior = "cs",
                                              beta = beta,
                                              iss = nrow(.data),
                                              tabu = 50), cpdag = FALSE)
# The following call to boot.strength is an intentional error.
# It is called with all the same arguments as would be called with
# mbde.  If the results are the same, then something is wrong.
boot_bde <- boot.strength(ialarm[, 1:9], cl, algorithm = "tabu",
                          algorithm.args = list(score = "bde",
                                                exp = INT, prior = "cs",
                                                beta = beta,
                                                iss = nrow(.data),
                                                tabu = 50), cpdag = FALSE)
mean(averaging_entropy(boot)$entropy) > mean(averaging_entropy(boot_mbde)$entropy)
# The prior was constructed on a dataset of 100000 observations.
# So additional data with no interventional information should
# not be more informative than interventional data.
obs <- rbn(ground_truth, nrow(ialarm), .data, fit = "bayes")
boot_bde <- boot.strength(obs, cl, algorithm = "tabu",
                          algorithm.args = list(score = "bde",
                                                prior = "cs",
                                                beta = beta,
                                                iss = nrow(.data),
                                                tabu = 50), cpdag = FALSE)
abs(mean(averaging_entropy(boot_bde)$entropy) -
      mean(averaging_entropy(boot_mbde)$entropy)) > .01

boot_mbde %>%
  reduce_averaging


#    from   to  boot_s boot_d bde_s bde_d mbde_s mbde_d
# 1   PMB  PAP  1.000   0.905     1 1.000  0.905  0.960
# 2   PMB SHNT  1.000   1.000     1 1.000  0.915  0.906
# 16 SHNT  INT  1.000   0.000     1 0.000  0.920  0.071
# 22  INT MINV  1.000   1.000     1 1.000  0.930  0.929
# 25  INT PRSS  1.000   1.000     1 1.000  0.900  0.927
# 23  INT VLNG  1.000   1.000     1 1.000  0.990  0.989
# 34 KINK PRSS  1.000   1.000     1 1.000  0.870  0.959
# 33 VLNG VTUB  1.000   1.000     1 1.000  0.995  0.914
# 36 PRSS VTUB  1.000   0.550     1 1.000  0.910  0.900
# 22  INT MINV  1.000   1.000     1 1.000  0.930  0.929
# 27 MINV VLNG  1.000   1.000     1 0.000  0.995  0.281

# Not much change from the boot results to the mbde results.
# The initial boot results have much certainty in terms of
# direction.  I think this is because of the v-structures.
# No, it is because I used hc.  Use tabu.
robertness/bninfo documentation built on May 27, 2019, 10:32 a.m.