# 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.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.