library(bnlearn) data(learning.test) str(learning.test) bn.gs <- gs(learning.test) bn.gs bn2 <- iamb(learning.test) bn3 <- fast.iamb(learning.test) bn4 <- inter.iamb(learning.test) compare(bn.gs, bn2) compare(bn.gs, bn3) compare(bn.gs, bn4) bn.hc <- hc(learning.test, score = "aic") bn.hc compare(bn.hc, bn.gs) par(mfrow = c(1,2)) plot(bn.gs, main = "Constraint-based algorithms", highlight = c("A", "B")) plot(bn.hc, main = "Hill-Climbing", highlight = c("A", "B")) par(mfrow = c(1,2)) highlight.opts <- list(nodes = c("A", "B"), arcs = c("A", "B"), col = "red", fill = "grey") graphviz.plot(bn.hc, highlight = highlight.opts) graphviz.plot(bn.gs, highlight = highlight.opts) bn.AB <- gs(learning.test, blacklist = c("B", "A")) compare(bn.AB, bn.hc) score(bn.AB, learning.test, type = "bde") bn.BA <- gs(learning.test, blacklist = c("A", "B")) score(bn.BA, learning.test, type = "bde") modelstring(bn.hc) undirected.arcs(bn.gs) bn.dag <- set.arc(bn.gs, "A", "B") modelstring(bn.dag) compare(bn.dag, bn.hc) try( set.arc(bn.hc, "E", "A") ) other <- empty.graph(nodes = nodes(bn.hc)) arcs(other) <- data.frame( from = c("A", "A", "B", "D"), to = c("E", "F", "C", "E")) other score(other, data = learning.test, type = "aic") score(bn.hc, data = learning.test, type = "aic") gs(learning.test, debug = TRUE) # ALARM data set example data(alarm) alarm.gs <- gs(alarm) alarm.iamb <- iamb(alarm) alarm.fast.iamb <- fast.iamb(alarm) alarm.inter.iamb <- inter.iamb(alarm) alarm.mmpc <- mmpc(alarm) alarm.hc <- hc(alarm, score = "bic") dag <- empty.graph(names(alarm)) modelstring(dag) <- 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.gs <- gs(alarm, test = "x2") alarm.mc <- gs(alarm, test = "mc-x2", B = 10000) par(mfrow = c(1,2), omi = rep(0, 4), mar = c(1, 0, 1, 0)) graphviz.plot(dag, highlight = list(arcs = arcs(alarm.gs))) graphviz.plot(dag, highlight = list(arcs = arcs(alarm.mc))) # MARKS data set example data(marks) str(marks) ci.test("MECH", "ANL", "ALG", data = marks) ci.test("STAT", "VECT", "ALG", data = marks) ci.test("STAT", "VECT", "ALG", data = marks, test = "zf")$p.value ci.test("STAT", "VECT", "ALG", data = marks, test = "mc-cor")$p.value ci.test("STAT", "VECT", "ALG", data = marks, test = "mi-g")$p.value ci.test("STAT", "VECT", "ALG", data = marks, test = "mc-mi-g")$p.value # Figure 1 rm(list = ls()) data(learning.test) bn.gs = gs(learning.test) bn.hc = hc(learning.test) pdf(file = "compare.pdf", width = 12, height = 6) par(mfrow = c(1,2), omi = rep(0, 4), mar = rep(1,4)) plot(bn.gs, main = "Constraint-based algorithms", highlight = c("A", "B")) plot(bn.hc, main = "Hill-Climbing", highlight = c("A", "B")) dev.off() # Figure 2 rm(list = ls()) data(alarm) dag = empty.graph(names(alarm)) modelstring(dag) = 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 = "") dag.gs = gs(alarm) dag.iiamb = inter.iamb(alarm) dag.hc = hc(alarm, score = "bic") pdf(file = "alarm.pdf", width = 12, height = 12) par(mfrow = c(2,2), omi = rep(0, 4), mar = c(1, 0, 1, 0)) graphviz.plot(dag) graphviz.plot(dag.gs) graphviz.plot(dag.iiamb) graphviz.plot(dag.hc) dev.off() # FIgure 3 rm(list = ls()) data(alarm) dag = empty.graph(names(alarm)) modelstring(dag) = 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.gs = gs(alarm, test = "x2") alarm.mc = gs(alarm, test = "mc-x2", B = 10000) pdf(file = "alarm2.pdf", width = 12, height = 4) par(mfrow = c(1,2), omi = rep(0, 4), mar = c(1, 0, 1, 0)) graphviz.plot(dag, highlight = list(arcs = arcs(alarm.gs))) graphviz.plot(dag, highlight = list(arcs = arcs(alarm.mc))) dev.off() # FIgure 4 rm(list = ls()) data(marks) marks.ug = empty.graph(names(marks)) arcs(marks.ug) = matrix( c("MECH", "VECT", "MECH", "ALG", "VECT", "MECH", "VECT", "ALG", "ALG", "MECH", "ALG", "VECT", "ALG", "ANL", "ALG", "STAT", "ANL", "ALG", "ANL", "STAT", "STAT", "ALG", "STAT", "ANL"), ncol = 2, byrow = TRUE, dimnames = list(c(), c("from", "to"))) marks.gs = gs(marks) marks.hc = hc(marks, score = "bic") pdf(file = "marks.pdf", width = 12, height = 4) par(mfrow = c(1,3), omi = rep(0, 4), mar = c(1, 0, 1, 0)) graphviz.plot(marks.ug) graphviz.plot(marks.gs) graphviz.plot(marks.hc) dev.off() # Table 1 rm(list = ls()) data(alarm) print(gs(alarm)) print(iamb(alarm)) print(fast.iamb(alarm)) print(inter.iamb(alarm)) print(hc(alarm, score = "bic")) system.time(for (i in 1:20) gs(alarm))/20 system.time(for (i in 1:20) iamb(alarm))/20 system.time(for (i in 1:20) fast.iamb(alarm))/20 system.time(for (i in 1:20) inter.iamb(alarm))/20 system.time(for (i in 1:20) hc(alarm, score = "bic"))/20