library(Rgraphviz)
## Loading required package: graph
## Loading required package: BiocGenerics
## Loading required package: parallel
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which, which.max, which.min
## Loading required package: grid
library(gRain)
## Loading required package: gRbase
yn <- c("yes", "no")
B <- cptable(~BatCave, values=c(1, 99), levels = yn)
L <- cptable(~LabModified, values = c(1,1), levels = yn)
J.BL <- cptable(~JumpToHumans|BatCave:LabModified, values = c(1, 1, 98, 10, 2, 88, 1, 1, 1, 2, 2, 1), levels =c("Other", "WetMarket", "HunanLab"))
V.J <- cptable(~VirusFoundAtMarket|JumpToHumans, values = c(1, 1, 9, 1, 2, 8), levels = yn)
P.B <- cptable(~BatCavePreprint|BatCave:JumpToHumans, values = c(9, 1, 1, 9, 9, 1, 1, 9, 99, 1, 1, 99), levels = yn)
W.JL <- cptable(~WhoClearsLab|JumpToHumans:LabModified, values = c(1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 3), levels = yn)
R.L <- cptable(~RnaModified|LabModified, values = c(2, 1), levels = yn)
corona_cpt <- compileCPT(B, L, J.BL, V.J, P.B, W.JL, R.L)
corona_cpt
## cpt_spec with probabilities:
## P( BatCave )
## P( LabModified )
## P( JumpToHumans | BatCave LabModified )
## P( VirusFoundAtMarket | JumpToHumans )
## P( BatCavePreprint | BatCave JumpToHumans )
## P( WhoClearsLab | JumpToHumans LabModified )
## P( RnaModified | LabModified )
corona_cpt$BatCave
## BatCave
## yes no
## 0.01 0.99
corona_cpt$LabModified
## LabModified
## yes no
## 0.5 0.5
corona_cpt$JumpToHumans
## , , LabModified = yes
##
## BatCave
## JumpToHumans yes no
## Other 0.01 0.10
## WetMarket 0.01 0.02
## HunanLab 0.98 0.88
##
## , , LabModified = no
##
## BatCave
## JumpToHumans yes no
## Other 0.3333333 0.4
## WetMarket 0.3333333 0.4
## HunanLab 0.3333333 0.2
corona_cpt$VirusFoundAtMarket
## JumpToHumans
## VirusFoundAtMarket Other WetMarket HunanLab
## yes 0.5 0.9 0.2
## no 0.5 0.1 0.8
corona_cpt$BatCavePreprint
## , , JumpToHumans = Other
##
## BatCave
## BatCavePreprint yes no
## yes 0.9 0.1
## no 0.1 0.9
##
## , , JumpToHumans = WetMarket
##
## BatCave
## BatCavePreprint yes no
## yes 0.9 0.1
## no 0.1 0.9
##
## , , JumpToHumans = HunanLab
##
## BatCave
## BatCavePreprint yes no
## yes 0.99 0.01
## no 0.01 0.99
corona_cpt$WhoClearsLab
## , , LabModified = yes
##
## JumpToHumans
## WhoClearsLab Other WetMarket HunanLab
## yes 0.5 0.5 0.3333333
## no 0.5 0.5 0.6666667
##
## , , LabModified = no
##
## JumpToHumans
## WhoClearsLab Other WetMarket HunanLab
## yes 0.5 0.5 0.25
## no 0.5 0.5 0.75
corona_cpt$RnaModified
## LabModified
## RnaModified yes no
## yes 0.6666667 0.6666667
## no 0.3333333 0.3333333
corona_bn <- grain(corona_cpt)
corona_bn <- propagate(corona_bn)
corona_bn
## Independence network: Compiled: TRUE Propagated: TRUE
## Nodes: chr [1:7] "BatCave" "LabModified" "JumpToHumans" "VirusFoundAtMarket" ...
querygrain(corona_bn, nodes = c("BatCave", "JumpToHumans"), type = "joint")
## JumpToHumans
## BatCave Other WetMarket HunanLab
## yes 0.001716667 0.001716667 0.006566667
## no 0.247500000 0.207900000 0.534600000
corona_bn2 <- setEvidence(corona_bn, evidence = list(BatCavePreprint = "yes", VirusFoundAtMarket = "no", WhoClearsLab = "yes", RnaModified = "no"))
corona_bn2 <- propagate(corona_bn2)
corona_bn2
## Independence network: Compiled: TRUE Propagated: TRUE
## Nodes: chr [1:7] "BatCave" "LabModified" "JumpToHumans" "VirusFoundAtMarket" ...
## Evidence:
## nodes is.hard.evidence hard.state
## 1 BatCavePreprint TRUE yes
## 2 VirusFoundAtMarket TRUE no
## 3 WhoClearsLab TRUE yes
## 4 RnaModified TRUE no
## pEvidence: 0.003558
querygrain(corona_bn2, nodes = c("BatCave", "JumpToHumans"), type = "joint")
## JumpToHumans
## BatCave Other WetMarket HunanLab
## yes 0.03618708 0.007237415 0.1521122
## no 0.57969589 0.097388909 0.1273785