1 Example

1.1 bnlearn package

learning sturcture
inference

library(bnlearn)

## 
## Attaching package: 'bnlearn'

## The following object is masked from 'package:stats':
## 
##     sigma

# coronary all columns are factors
#   P. Work (strenuous physical work): a two-level factor with levels no and yes.
#   Pressure (systolic blood pressure): a two-level factor with levels <140 and >140.
coronary[1:5,]

##   Smoking M. Work P. Work Pressure Proteins Family
## 1      no      no      no     <140       <3    neg
## 2      no      no      no     <140       <3    neg
## 3      no      no      no     <140       <3    neg
## 4      no      no      no     <140       <3    neg
## 5      no      no      no     <140       <3    neg

# ---- Learn structure --------

bn_df <- data.frame(coronary[1:700,])

# constraint-based (gs, iamb, fast.iamb and inter.iamb) and score-based (hc) algorithms
par(mfrow = c(1,2), omi = rep(0, 4), mar = c(1, 1, 1, 1))
res <- gs(bn_df) # onstraint-based algorithms
plot(res)
res <- hc(bn_df) # Score-based algorithms:
plot(res)

# it does not make sense to have Family as a variable condition on M.Work. 
# Therefore, we need to modify the derived structure. 
# remove the link between M.Work and Family.
res$arcs <- res$arcs[-which((res$arcs[,'from'] == "M..Work" & res$arcs[,'to'] == "Family")),]
plot(res)


# After learning the structure, we need to find out the conditional probability tables (CPTs) at each node. The bn.fit function runs the EM algorithm to learn CPT for different nodes in the above graph.
fittedbn <- bn.fit(res, data = bn_df)


print(fittedbn$Proteins)

## 
##   Parameters of node Proteins (multinomial distribution)
## 
## Conditional probability table:
##  <3 >3 
##  1  0

# ------ Inference ----------
# Perform conditional probability queries (CPQs)
cpquery(fittedbn, event=(Proteins=="<3"), evidence=(Smoking=="no"))

## [1] 1

cpquery(fittedbn, event=(Proteins=="<3"), evidence=(Smoking=="no"&Pressure==">140"))

## [1] 1

cpquery(fittedbn, event=(Proteins=="<3"), evidence=(Pressure==">140"))

## [1] 1

cpquery(fittedbn, event=(Smoking=="no"&Pressure==">140"), evidence=(Proteins=="<3"))

## [1] 0.164

str <- "Proteins=='<3'"
cpquery(fittedbn, eval(parse(text = str)), evidence=(Smoking=="no"))

## [1] 1

1.2 gRain package

library(gRain)

## Loading required package: gRbase

## 
## Attaching package: 'gRbase'

## The following objects are masked from 'package:bnlearn':
## 
##     children, parents

yn   <- c("yes","no")
ynm  <- c("yes","no","maybe")
# To specify P(a|b,c) one may write ~a|b:c, ~a:b:c, ~a|b+c, ~a+b+c or c("a","b","c").
# for a, yes:1, No:99
a    <- cptable( ~ asia, values=c(1,99), levels=yn)
# for tub, yes:5, no:95, maybe 1. 
# for asia, yes:1, no:99, maybe 999
t.a  <- cptable( ~ tub : asia, values=c(5, 95, 1, 99, 1, 999),  levels=ynm)
d.a  <- cptable( ~ dia : asia, values=c(5, 5, 1, 99, 100, 999), levels=ynm)
cptlist <- compileCPT(list(a, t.a, d.a))
grn <- grain(cptlist)
## Bayesian Network
summary (grn)

## Independence network: Compiled: FALSE Propagated: FALSE 
##  Nodes : chr [1:3] "asia" "tub" "dia"

plot (grn)

q1 <- setEvidence(grn, c("asia"), c("yes"))
querygrain(q1)

## $asia
## asia
## yes  no 
##   1   0 
## 
## $tub
## tub
##        yes         no      maybe 
## 0.04950495 0.94059406 0.00990099 
## 
## $dia
## dia
##        yes         no      maybe 
## 0.45454545 0.45454545 0.09090909

q2 <- setEvidence(grn, c("asia"), c("no"))
querygrain(q2, nodes="tub")$tub

## tub
##          yes           no        maybe 
## 0.0900818926 0.0009099181 0.9090081893

library(gRain)
## Example: Specifying conditional probabilities as a matrix
bayes.levels  <- c('Enzyme', 'Keratine', 'unknown')
root.node     <- cptable( ~R, values=c( 1, 1, 1 ), levels=bayes.levels)
cond.prob.tbl <- t(matrix( c( 1, 0, 0, 0, 1, 0, 0.5, 0.5, 0 ),
   nrow=3, ncol=3, byrow=TRUE, dimnames=list(bayes.levels, bayes.levels)))
cond.prob.tbl

##          Enzyme Keratine unknown
## Enzyme        1        0     0.5
## Keratine      0        1     0.5
## unknown       0        0     0.0

## Notice above: Columns represent parent states; rows represent child states
query.node    <- cptable( ~ Q | R, values=cond.prob.tbl, levels=bayes.levels )
sister.node   <- cptable( ~ S | R, values=cond.prob.tbl, levels=bayes.levels )
## Testing 
compile(grain(compileCPT(list( root.node, query.node, sister.node ))), propagate=TRUE)

## Independence network: Compiled: TRUE Propagated: TRUE 
##   Nodes: chr [1:3] "R" "Q" "S"

dat <- read.table(text="A B 
1 2
1 3
1 4
2 5
3 7", header=TRUE)

x <- table(dat)
x %*% t(x)

##    A
## A   1 2 3
##   1 3 0 0
##   2 0 1 0
##   3 0 0 1

library(gRain)
data(chestSim1000, package="gRbase")
head(chestSim1000)

##   asia tub smoke lung bronc either xray dysp
## 1   no  no    no   no   yes     no   no  yes
## 2   no  no   yes   no   yes     no   no  yes
## 3   no  no   yes   no    no     no   no   no
## 4   no  no    no   no    no     no   no   no
## 5   no  no   yes   no   yes     no   no  yes
## 6   no  no   yes  yes   yes    yes  yes  yes

s <- xtabs(~smoke, chestSim1000)
d.b <- xtabs(~dysp+bronc, chestSim1000)

chestSim1000 <- chestSim1000[,c("smoke", "bronc", "dysp")]

s <- xtabs(~smoke, chestSim1000)
s <- as.parray(s, normalize="first")
b <- xtabs(~bronc+., chestSim1000)
b <- as.parray(b, normalize="first")
d <- xtabs(~dysp+., chestSim1000)
d <- as.parray(d, normalize="first")
# cpt.list <- compileCPT(list(s, b, d))
# cpt.list

Bayesian Networks

Updated: 18 October, 2016

1 Example

1.1 bnlearn package

1.2 gRain package