tests/testthat/test-performance.R

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

context("performance")

data(alarm)
alarm_truth <- empty.graph(names(alarm))
modelstring(alarm_truth) <- 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 = "")
alarm_infer <- hc(alarm)
small_alarm_truth <- arcs(alarm_truth)[1:6, ] %>% construct_bn
inferred_small_alarm_truth <- alarm[, nodes(small_alarm_truth)] %>% hc

gauss_truth <- empty.graph(names(gaussian.test))
modelstring(gauss_truth) <- "[A][B][E][G][C|A:B][D|B][F|A:D:E:G]"
gauss_infer <- hc(gaussian.test)


expect_that("performance_outcomes returns correct edge name sets.", {
            # Use a small version of the alarm network as a test case
              small_alarm_truth <- arcs(alarm_truth)[1:6, ] %>% construct_bn
              small_alarm_inferred <- alarm[, nodes(small_alarm_truth)] %>% hc
              inf <- c("LVV->PCWP", "LVV->CVP", "LVV->HYP", "LVF->HIST", "LVF->LVV", "LVV->STKV", "LVF->STKV", "STKV->HYP")
              truth <- c("LVF->HIST", "LVF->LVV", "HYP->LVV", "HYP->STKV", "LVV->CVP", "LVV->PCWP")
              list(tp = intersect(inf, truth),
                   fp = setdiff(inf, truth),
                   fn = setdiff(truth, inf)) %>%
                expect_identical(performance_outcomes(small_alarm_inferred, small_alarm_truth))
              # test undirected case
              inf <- unique(arcs2names(arcs(moral(small_alarm_inferred)), directed = FALSE))
              truth <- unique(arcs2names(arcs(moral(small_alarm_truth)), directed = FALSE))
              list(tp = intersect(inf, truth),
                   fp = setdiff(inf, truth),
                   fn = setdiff(truth, inf)) %>%
                expect_identical(performance_outcomes(moral(small_alarm_inferred),
                                                      moral(small_alarm_truth)))
              # Koller style cpdag
              name_edge_df(cpdag(small_alarm_inferred, moral = FALSE))$edge_name
              inf <- name_edge_df(cpdag(small_alarm_inferred, moral = FALSE))$edge_name
              truth <- name_edge_df(cpdag(small_alarm_truth, moral = FALSE))$edge_name
              list(tp = intersect(inf, truth),
                   fp = setdiff(inf, truth),
                   fn = setdiff(truth, inf)) %>%
                expect_identical(performance_outcomes(cpdag(small_alarm_inferred, moral = FALSE),
                                                      cpdag(small_alarm_truth, moral = FALSE)))

})



library(ROCR)
set.seed(5)
net <- simGaussianNet(8)
net.structure <- bn.net(net)
sim.data <- rbn(net, 1000)
ma.results <- boot.strength(sim.data, R = 100, m = 1000, algorithm = "tabu",
                            algorithm.args = list(score = "bic-g"))
output <- get_predictions(net_structure, ma_results)

## Want to superimpose two plots

inferred_cpdag <- empty.graph(nodes(ground_truth))
modelstring(inferred_cpdag) <- "[PAP][FIO2][APL][ANES][INT][KINK][DISC][LVV][STKV][HR][ERCA][ACO2][VMCH][CVP|LVV][PCWP|LVV][CO|STKV:HR][HREK|HR:ERCA][HRSA|HR:ERCA][PRSS|KINK][MVS|VMCH][HYP|LVV:STKV][LVF|LVV][ERLO|HR][SHNT|INT][VALV|INT][VTUB|DISC:VMCH][HIST|LVF][TPR|CO][HRBP|ERLO:HR][PMB|PAP:SHNT][PVS|FIO2:VALV][VLNG|INT:ACO2:VALV][BP|TPR:CO][SAO2|SHNT:PVS][ECO2|ACO2:VLNG][MINV|INT:VLNG][CCHL|TPR:SAO2:HR]"
inferred_cpdag <- cpdag(inferred_cpdag, moral = FALSE)
# Needs testing
performance_arc_list(inferred, cpdag(ground_truth, moral=FALSE))



    highlight_ground_truth <- list(
      detected = list(detected_edges, "green", "solid"),
      fn = list(false_negatives, "black", "dashed")
    ) %>%
    lapply(function(item){
      if(length(item[[1]]) > 0){
        data.frame(edge_name = item[[1]], col = item[[2]], lty = item[[3]])
      } else {
        NULL
      }
    }) %>%
    lapply(function(item){
      if(!is.null(item)){
        merge(item, ground_truth_df, by = "edge_name") %>%
          select(from, to, col, lty)
      } else {
        NULL
      }
    }) %>%
{do.call("rbind", .)} %>%
{list(arcs = as.matrix(dplyr::select(., from, to)),
      col = as.character(.$col),
      lty = as.character(.$lty))}
##
highlight_inferred <- list(
  detected = list(detected_edges, "green", "solid"),
  fp = list(false_positives, "darkred", "solid")
) %>%
  lapply(function(item){
    if(length(item[[1]]) > 0){
      data.frame(edge_name = item[[1]], col = item[[2]], lty = item[[3]])
    } else {
      NULL
    }
  }) %>%
  lapply(function(item){
    if(!is.null(item)){
      merge(item, inferred_df, by = "edge_name") %>%
        select(from, to, col, lty)
    } else {
      NULL
    }
  }) %>%
{do.call("rbind", .)} %>%
{list(arcs = as.matrix(dplyr::select(., from, to)),
      col = as.character(.$col),
      lty = as.character(.$lty))}
list(ground_truth = highlight_ground_truth,
     inferred = highlight_inferred)

}


x <- performance_plot_list(inferred, groundtruth)
test_that("performance plot, used to compare an inferred structure to a ground truth reference,
          works on partially directed networks", {
            gt <- arcs(groundtruth)[1:10, ] %>% construct_bn %>% cpdag
            inf <- alarm[, nodes(gt)] %>% hc %>% cpdag
            # Have to eyeball this test.

            performance_plot(gt, inf, plot_truth = TRUE)
            performance_plot(gt, inf, plot_truth = FALSE)

          })

##

graphviz.plot(groundtruth, highlight = highlight_gt)
graphviz.plot(inferred, highlight = highlight_inferred)

test_that("when working with CPDAGs we use the Koller (as opposed to Pearl) approach,
          which coerces edges in v-structures", {
            expect_true(FALSE)
          })





duplicated_row <- function(mat){
  mat <- as.matrix(mat)
  apply(mat, 1, function(row){
    paste0(sort(as.character(row)), collapse="_")
  }) %>%
    duplicated
}