---

In this analysis I am going to build a Bayesian Network in R based on medical data to identify the likelihood of six possible diseases based on clinical presentations in infants.

The idea and data for this analysis came from the CHILD network created by David J. Spiegelhalter, A. Philip Dawid, Steffen L. Lauritzen, Robert G. Cowell. Here is a citation to an article on their work and the structure of a Bayesian Network they created:

David J. Spiegelhalter, A. Philip Dawid, Steffen L. Lauritzen, Robert G. Cowell “Bayesian Analysis in Expert Systems,” Statistical Science, Statist. Sci. 8(3), 219–247, (August, 1993)

I am going to use the Spiegelhalter et al. model as guidance to build the structure of a network. I will then parameterize with medical data, visualize the model, and perform queries. I will also visualize evidence based queries.

Load necessary packages.

library(bnlearn)

if (!require(“BiocManager”, quietly = TRUE))
 install.packages(“BiocManager”)

BiocManager::install(“Rgraphviz”)

if (!require(“BiocManager”, quietly = TRUE))
 install.packages(“BiocManager”)

BiocManager::install(“graph”)

library(Rgraphviz)
library(gRbase)
library(graph)

Construct the BN structure in R.

# create an empty graph
dag <- empty.graph(nodes = c(“BirthAsphyxia”, “Disease”, “Age”, “LVH”, “DuctFlow”, “CardiacMixing”, “LungParench”, “LungFlow”, “Sick”, “HypDistrib”, “HypoxiaInO2”, “CO2”,
 “ChestXray”, “Grunting”, “LVHreport”, “LowerBodyO2”, “RUQO2”, “CO2Report”, “XrayReport”, “GruntingReport”))

quartz()
plot(dag)

dag <- set.arc(dag, from = “BirthAsphyxia”, to = “Disease”) #birth asphyxia -> disease
dag <- set.arc(dag, from = “Disease”, to = “Age”) #disease -> age at presentation 
dag <- set.arc(dag, from = “Disease”, to = “LVH”) #disease -> LVH
dag <- set.arc(dag, from = “Disease”, to = “DuctFlow”) #disease -> duct flow
dag <- set.arc(dag, from = “Disease”, to = “CardiacMixing”) #disease -> cardiac mixing
dag <- set.arc(dag, from = “Disease”, to = “LungParench”) #disease -> lung parenchema 
dag <- set.arc(dag, from = “Disease”, to = “LungFlow”) #disease -> lung flow
dag <- set.arc(dag, from = “Disease”, to = “Sick”) #disease -> sick
dag <- set.arc(dag, from = “LVH”, to = “LVHreport”) #LVH-> LVH report
dag <- set.arc(dag, from = “DuctFlow”, to = “HypDistrib”) #duct flow -> hypoxia distribution
dag <- set.arc(dag, from = “CardiacMixing”, to = “HypDistrib”) #cardiac mixing -> hypoxia distribution
dag <- set.arc(dag, from = “CardiacMixing”, to = “HypoxiaInO2”) #cardiac mixing -> hypoxia in O2
dag <- set.arc(dag, from = “LungParench”, to = “HypoxiaInO2”) #lung parenchema -> hypoxia in O2
dag <- set.arc(dag, from = “LungParench”, to = “CO2”) #lung parenchema -> CO2
dag <- set.arc(dag, from = “LungParench”, to = “ChestXray”) #lung parenchema -> chest x-ray
dag <- set.arc(dag, from = “LungParench”, to = “Grunting”) #lung parenchema -> grunting 
dag <- set.arc(dag, from = “LungFlow”, to = “ChestXray”) #lung flow -> chest x-ray
dag <- set.arc(dag, from = “Sick”, to = “Grunting”) #sick -> grunting
dag <- set.arc(dag, from = “HypDistrib”, to = “LowerBodyO2”) #hypoxia distribution -> lower body O2
dag <- set.arc(dag, from = “HypoxiaInO2”, to = “LowerBodyO2”) #hypoxia in O2 -> lower body O2
dag <- set.arc(dag, from = “HypoxiaInO2”, to = “RUQO2”) #hypoxia in O2 -> right up quad O2
dag <- set.arc(dag, from = “CO2”, to = “CO2Report”) #CO2 -> CO2 report
dag <- set.arc(dag, from = “ChestXray”, to = “XrayReport”) #chest x-ray -> chest x-ray report
dag <- set.arc(dag, from = “Grunting”, to = “GruntingReport”) #grunting -> grunting report

quartz()
plot(dag)

The joint distribution:

modelstring(dag)

“[n1][n2|n1][n3|n2][n4|n2][n5|n2][n6|n2][n7|n2][n8|n2][n9|n2][n10|n5:n6][n11|n6:n7][n12|n7][n13|n7:n8][n14|n7:n9][n15|n4][n16|n10:n11][n17|n11][n18|n12][n19|n13][n20|n14]”

“[BirthAsphyxia][Disease|BirthAsphyxia][Age|Disease][LVH|Disease][DuctFlow|Disease][CardiacMixing|Disease][LungParench|Disease][LungFlow|Disease][Sick|Disease][HypDistrib|DuctFlow:CardiacMixing][HypoxiaInO2|CardiacMixing:LungParench][CO2|LungParench][ChestXray|LungParench:LungFlow][Grunting|LungParench:Sick][LVHreport|LVH][LowerBodyO2|HypDistrib:HypoxiaInO2][RUQO2|HypoxiaInO2][CO2Report|CO2][XrayReport|ChestXray][GruntingReport|Grunting]”

nodes(dag)
arcs(dag)

Parametrize the network.

child<- read.csv(“child_network.csv”, header=TRUE, colClasses=”factor”)
dim(child)
head(child)

bn.mle <- bn.fit(dag, data = child, method = “mle”) #maximum likelihood estimation
bn.mle

quartz()
graphviz.chart(bn.mle, grid = TRUE, main = “BN”)

Explore the conditional probability tables;

for example, what is the conditional probability table for n1?
bn.mle$BirthAsphyxia

What is the conditional probability table for n16?
bn.mle$LowerBodyO2

Perform queries on the data.

Suppose “lower body O2 < 5” and “X-ray report = Plethoric”, what can we deduce about the disease?

cpquery(bn.mle, event = (Disease == “Fallot”), evidence = ((LowerBodyO2 == “<5” & XrayReport == “Plethoric”))) #0.09681228

cpquery(bn.mle, event = (Disease == “Lung”), evidence = ((LowerBodyO2 == “<5” & XrayReport == “Plethoric”))) #0.02466368

cpquery(bn.mle, event = (Disease == “PAIVS”), evidence = ((LowerBodyO2 == “<5” & XrayReport == “Plethoric”))) #0.08390805

cpquery(bn.mle, event = (Disease == “PFC”), evidence = ((LowerBodyO2 == “<5” & XrayReport == “Plethoric”))) #0.01668521

cpquery(bn.mle, event = (Disease == “TAPVD”), evidence = ((LowerBodyO2 == “<5” & XrayReport == “Plethoric”))) #0.05727377

cpquery(bn.mle, event = (Disease == “TGA”), evidence = ((LowerBodyO2 == “<5” & XrayReport == “Plethoric”))) #0.7070707

library(gRain)
junction <- compile(as.grain(bn.mle))

jedu <- setEvidence(junction, nodes = c(“LowerBodyO2”, “XrayReport”), states = c(“<5”, “Plethoric”)) #set the evidence

quartz()
graphviz.chart(bn.mle, grid = TRUE, main = “Original BN”)
graphviz.chart(as.bn.fit(jedu, including.evidence = TRUE), grid = TRUE,
 bar.col = c(BirthAsphyxia = “black”, Disease = “red”, Age = “black”, LVH = “black”,
 DuctFlow = “black”, CardiacMixing = “black”, LungParench=”black”,
 LungFlow=”black”, Sick=”black”, HypDistrib=”black”, HypoxiaInO2=”black”,
 CO2=”black”, ChestXray=”black”, Grunting=”black”, LVHreport=”black”,
 LowerBodyO2=”red”, RUQO2=”black”, CO2Report=”black”, XrayReport=”red”,
 GruntingReport=”black”),
 strip.bg = c(BirthAsphyxia = “transparent”, Disease = “red”, Age = “transparent”, LVH = “transparent”,
 DuctFlow = “transparent”, CardiacMixing = “transparent”, LungParench=”transparent”,
 LungFlow=”transparent”, Sick=”transparent”, HypDistrib=”transparent”, HypoxiaInO2=”transparent”,
 CO2=”transparent”, ChestXray=”transparent”, Grunting=”transparent”, LVHreport=”transparent”,
 LowerBodyO2=”red”, RUQO2=”transparent”, CO2Report=”transparent”, XrayReport=”red”,
 GruntingReport=”transparent”), main = “BN with Evidence”)

We can deduce the disease is likely TGA; there is a 70% probability.

Suppose “lower body O2 < 5” and “X-ray report = Oligaemic” but the child is “not grunting”, what can we deduce about the disease? 

cpquery(bn.mle, event = (Disease == “Fallot”), evidence = ((LowerBodyO2 == “<5”) & (XrayReport == “Oligaemic”) & (Grunting == “no”))) #0.4626866

cpquery(bn.mle, event = (Disease == “Lung”), evidence = ((LowerBodyO2 == “<5”) & (XrayReport == “Oligaemic”) & (Grunting == “no”))) #00.01447527

cpquery(bn.mle, event = (Disease == “PAIVS”), evidence = ((LowerBodyO2 == “<5”) & (XrayReport == “Oligaemic”) & (Grunting == “no”))) #0.3218263

cpquery(bn.mle, event = (Disease == “PFC”), evidence = ((LowerBodyO2 == “<5”) & (XrayReport == “Oligaemic”) & (Grunting == “no”))) # 0.04091456

cpquery(bn.mle, event = (Disease == “TAPVD”), evidence = ((LowerBodyO2 == “<5”) & (XrayReport == “Oligaemic”) & (Grunting == “no”))) #0.009888752

cpquery(bn.mle, event = (Disease == “TGA”), evidence = ((LowerBodyO2 == “<5”) & (XrayReport == “Oligaemic”) & (Grunting == “no”))) #0.08656036

jedu <- setEvidence(junction, nodes = c(“LowerBodyO2”, “XrayReport”, “Grunting”), states = c(“<5”, “Oligaemic”, “no”)) #set the evidence

quartz()
graphviz.chart(bn.mle, grid = TRUE, main = “Original BN”)
graphviz.chart(as.bn.fit(jedu, including.evidence = TRUE), grid = TRUE,
 bar.col = c(BirthAsphyxia = “black”, Disease = “red”, Age = “black”, LVH = “black”,
 DuctFlow = “black”, CardiacMixing = “black”, LungParench=”black”,
 LungFlow=”black”, Sick=”black”, HypDistrib=”black”, HypoxiaInO2=”black”,
 CO2=”black”, ChestXray=”black”, Grunting=”red”, LVHreport=”black”,
 LowerBodyO2=”red”, RUQO2=”black”, CO2Report=”black”, XrayReport=”red”,
 GruntingReport=”black”),
 strip.bg = c(BirthAsphyxia = “transparent”, Disease = “red”, Age = “transparent”, LVH = “transparent”,
 DuctFlow = “transparent”, CardiacMixing = “transparent”, LungParench=”transparent”,
 LungFlow=”transparent”, Sick=”transparent”, HypDistrib=”transparent”, HypoxiaInO2=”transparent”,
 CO2=”transparent”, ChestXray=”transparent”, Grunting=”red”, LVHreport=”transparent”,
 LowerBodyO2=”red”, RUQO2=”transparent”, CO2Report=”transparent”, XrayReport=”red”,
 GruntingReport=”transparent”), main = “BN with Evidence”)

We can deduce that the disease is most likely Fallot with a 46% probability; alternatively, there is a 32% probability it is PAIVS.

Baby Bea is “grunting” with “mild cardiac mixing”. Baby Hank is “not grunting” with “complete cardiac mixing”. What can we learn about these babies?

#Bea
cpquery(bn.mle, event = (Disease == “Fallot”), evidence = ((Grunting == “yes”) & (CardiacMixing == “Mild”))) #0.1809291

cpquery(bn.mle, event = (Disease == “Lung”), evidence = ((Grunting == “yes”) & (CardiacMixing == “Mild”))) #0.4311224

cpquery(bn.mle, event = (Disease == “PAIVS”), evidence = ((Grunting == “yes”) & (CardiacMixing == “Mild”))) #0.01578947

cpquery(bn.mle, event = (Disease == “PFC”), evidence = ((Grunting == “yes”) & (CardiacMixing == “Mild”))) #0.1182796

cpquery(bn.mle, event = (Disease == “TAPVD”), evidence = ((Grunting == “yes”) & (CardiacMixing == “Mild”))) #0.0234375

cpquery(bn.mle, event = (Disease == “TGA”), evidence = ((Grunting == “yes”) & (CardiacMixing == “Mild”))) #0.1741573

jedu <- setEvidence(junction, nodes = c(“Grunting”, “CariacMixing”), states = c(“yes”, “Mild”)) #set the evidence

quartz()
graphviz.chart(bn.mle, grid = TRUE, main = “Original BN”)
graphviz.chart(as.bn.fit(jedu, including.evidence = TRUE), grid = TRUE,
 bar.col = c(BirthAsphyxia = “black”, Disease = “red”, Age = “black”, LVH = “black”,
 DuctFlow = “black”, CardiacMixing = “red”, LungParench=”black”,
 LungFlow=”black”, Sick=”black”, HypDistrib=”black”, HypoxiaInO2=”black”,
 CO2=”black”, ChestXray=”black”, Grunting=”red”, LVHreport=”black”,
 LowerBodyO2=”black”, RUQO2=”black”, CO2Report=”black”, XrayReport=”black”,
 GruntingReport=”black”),
 strip.bg = c(BirthAsphyxia = “transparent”, Disease = “red”, Age = “transparent”, LVH = “transparent”,
 DuctFlow = “transparent”, CardiacMixing = “red”, LungParench=”transparent”,
 LungFlow=”transparent”, Sick=”transparent”, HypDistrib=”transparent”, HypoxiaInO2=”transparent”,
 CO2=”transparent”, ChestXray=”transparent”, Grunting=”red”, LVHreport=”transparent”,
 LowerBodyO2=”transparent”, RUQO2=”transparent”, CO2Report=”transparent”, XrayReport=”transparent”,
 GruntingReport=”transparent”), main = “BN with Evidence”)

The most likely disease that Bea has is lung disease with a 43% probability.

#Hank 
cpquery(bn.mle, event = (Disease == “Fallot”), evidence = ((Grunting == “no”) & (CardiacMixing == “Complete”))) #0.4435653

cpquery(bn.mle, event = (Disease == “Lung”), evidence = ((Grunting == “no”) & (CardiacMixing == “Complete”))) #0.001911132

cpquery(bn.mle, event = (Disease == “PAIVS”), evidence = ((Grunting == “no”) & (CardiacMixing == “Complete”))) #0.4176822

cpquery(bn.mle, event = (Disease == “PFC”), evidence = ((Grunting == “no”) & (CardiacMixing == “Complete”))) #0.01378327

cpquery(bn.mle, event = (Disease == “TAPVD”), evidence = ((Grunting == “no”) & (CardiacMixing == “Complete”))) #0.06555819

cpquery(bn.mle, event = (Disease == “TGA”), evidence = ((Grunting == “no”) & (CardiacMixing == “Complete”))) #0.05988727

jedu <- setEvidence(junction, nodes = c(“Grunting”, “CariacMixing”), states = c(“no”, “Complete”)) #set the evidence

quartz()
graphviz.chart(bn.mle, grid = TRUE, main = “Original BN”)
graphviz.chart(as.bn.fit(jedu, including.evidence = TRUE), grid = TRUE,
 bar.col = c(BirthAsphyxia = “black”, Disease = “red”, Age = “black”, LVH = “black”,
 DuctFlow = “black”, CardiacMixing = “red”, LungParench=”black”,
 LungFlow=”black”, Sick=”black”, HypDistrib=”black”, HypoxiaInO2=”black”,
 CO2=”black”, ChestXray=”black”, Grunting=”red”, LVHreport=”black”,
 LowerBodyO2=”black”, RUQO2=”black”, CO2Report=”black”, XrayReport=”black”,
 GruntingReport=”black”),
 strip.bg = c(BirthAsphyxia = “transparent”, Disease = “red”, Age = “transparent”, LVH = “transparent”,
 DuctFlow = “transparent”, CardiacMixing = “red”, LungParench=”transparent”,
 LungFlow=”transparent”, Sick=”transparent”, HypDistrib=”transparent”, HypoxiaInO2=”transparent”,
 CO2=”transparent”, ChestXray=”transparent”, Grunting=”red”, LVHreport=”transparent”,
 LowerBodyO2=”transparent”, RUQO2=”transparent”, CO2Report=”transparent”, XrayReport=”transparent”,
 GruntingReport=”transparent”), main = “BN with Evidence”)

The most likely disease that Hank has is Fallot with a 44% probability.

In conclusion, this was an example of building a Bayesian Network in R. I started with building the structure of the network, and then showed how to parameterize using data. Visualizations and queries were then performed.

By Aspen Gulley on October 18, 2022.