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.max, which.min
## Loading required package: grid
library(gRain)
## Loading required package: gRbase
yn <- c("yes", "no")
S <- cptable(~SerialKiller, values=c(1, 999999), levels = yn)
U <- cptable(~UnknownCause, values=c(1, 99), levels = yn)
E.S <- cptable(~EvidenceOneMurder:SerialKiller, values = c(99, 1, 1, 999), levels = yn)
C.SU <- cptable(~ClusterofEvents:SerialKiller:UnknownCause, values = c(9999, 1, 1, 99, 9999, 1, 1, 99999), levels = yn)
killerNurse <- compileCPT(S, U, C.SU, E.S)
killerNurse
## cpt_spec with probabilities:
## P( SerialKiller )
## P( UnknownCause )
## P( ClusterofEvents | SerialKiller UnknownCause )
## P( EvidenceOneMurder | SerialKiller )
killerNurse$SerialKiller
## SerialKiller
## yes no
## 0.000001 0.999999
killerNurse$UnknownCause
## UnknownCause
## yes no
## 0.01 0.99
killerNurse$EvidenceOneMurder
## SerialKiller
## EvidenceOneMurder yes no
## yes 0.99 0.001
## no 0.01 0.999
killerNurse$ClusterofEvents
## , , UnknownCause = yes
##
## SerialKiller
## ClusterofEvents yes no
## yes 0.9999 0.01
## no 0.0001 0.99
##
## , , UnknownCause = no
##
## SerialKiller
## ClusterofEvents yes no
## yes 0.9999 0.00001
## no 0.0001 0.99999
killerNurse_bn <- grain(killerNurse)
killerNurse_bn <- propagate(killerNurse_bn)
killerNurse_bn
## Independence network: Compiled: TRUE Propagated: TRUE
## Nodes: chr [1:4] "SerialKiller" "UnknownCause" "ClusterofEvents" ...
##################################################################################
# No instantiation of evidence
querygrain(killerNurse_bn, nodes = c("SerialKiller", "UnknownCause"), type = "joint")
## UnknownCause
## SerialKiller yes no
## yes 0.00000001 0.00000099
## no 0.00999999 0.98999901
round(100*querygrain(killerNurse_bn, nodes = c("SerialKiller", "UnknownCause"), type = "joint"))
## UnknownCause
## SerialKiller yes no
## yes 0 0
## no 1 99
##################################################################################
# We observe no mysteriouus cluster and no evidence of a murder
killerNurse_bn0 <- setEvidence(killerNurse_bn, evidence = list(ClusterofEvents="no", EvidenceOneMurder="no"))
querygrain(killerNurse_bn0, nodes = c("SerialKiller", "UnknownCause"), type = "joint")
## UnknownCause
## SerialKiller yes no
## yes 1.001112e-14 9.911009e-13
## no 9.901088e-03 9.900989e-01
round(100*querygrain(killerNurse_bn0, nodes = c("SerialKiller", "UnknownCause"), type = "joint"))
## UnknownCause
## SerialKiller yes no
## yes 0 0
## no 1 99
# 99% chance nothing going on at all. 1% chance that things have changed but we did not see evidence of it
##################################################################################
# We just observe a suspicious cluster of events involving our nurse
killerNurse_bn1 <- setEvidence(killerNurse_bn, evidence = list(ClusterofEvents="yes"))
querygrain(killerNurse_bn1, nodes = c("SerialKiller", "UnknownCause"), type = "joint")
## UnknownCause
## SerialKiller yes no
## yes 9.016248e-05 0.008926085
## no 9.017141e-01 0.089269692
round(100*querygrain(killerNurse_bn1, nodes = c("SerialKiller", "UnknownCause"), type = "joint"))
## UnknownCause
## SerialKiller yes no
## yes 0 1
## no 90 9
# 90% chance something changed but our nurse is not a serial killer.
# 9% chance that actually nothing changed at all
# 1% chance we have a serial killer on our hands
##################################################################################
# We observe a suspicious cluster and strong evidence of wrong doing by our nurse on one particular patient
killerNurse_bn11 <- setEvidence(killerNurse_bn, evidence = list(ClusterofEvents="yes", EvidenceOneMurder="yes"))
querygrain(killerNurse_bn11, nodes = c("SerialKiller", "UnknownCause"), type = "joint")
## UnknownCause
## SerialKiller yes no
## yes 0.009000729 0.89107219
## no 0.090925458 0.00900162
round(100*querygrain(killerNurse_bn11, nodes = c("SerialKiller", "UnknownCause"), type = "joint"))
## UnknownCause
## SerialKiller yes no
## yes 1 89
## no 9 1
# 90% chance our nurse is a serial killer
# 9% chance she is not a serial killer but something else caused the cluster of events
# 1% chance nothing is going on at all; 1% chance she's a serial killer and things changed
##################################################################################
# We observe a suspicious cluster but no strong evidence of wrong doing by our nurse on any one patient
killerNurse_bn12 <- setEvidence(killerNurse_bn, evidence = list(ClusterofEvents="yes", EvidenceOneMurder="no"))
querygrain(killerNurse_bn12, nodes = c("SerialKiller", "UnknownCause"), type = "joint")
## UnknownCause
## SerialKiller yes no
## yes 9.106558e-07 9.015493e-05
## no 9.098352e-01 9.007369e-02
round(100*querygrain(killerNurse_bn12, nodes = c("SerialKiller", "UnknownCause"), type = "joint"))
## UnknownCause
## SerialKiller yes no
## yes 0 0
## no 91 9
# 91% chance something changed, the nurse is not a serial killer
# 9% chance nothing going on at all
###################################################################################
## Experimental graphical representation ##
###################################################################################
library(png)
library(igraph)
##
## Attaching package: 'igraph'
## The following objects are masked from 'package:gRbase':
##
## is_dag, topo_sort
## The following objects are masked from 'package:graph':
##
## degree, edges, intersection, union
## The following objects are masked from 'package:BiocGenerics':
##
## normalize, path, union
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
png("barsS.png")
barplot(height = c(30, 70), horiz = "TRUE", col = c("orange", "blue"), xlim = c(0, 100), axes = FALSE)
graphics.off()
PictureS <- readPNG("barsS.png")
png("barsC.png")
barplot(height = c(10, 90), horiz = "TRUE", col = c("orange", "blue"), xlim = c(0, 100), axes = FALSE)
graphics.off()
PictureC <- readPNG("barsC.png")
png("barsE.png")
barplot(height = c(70, 30), horiz = "TRUE", col = c("orange", "blue"), xlim = c(0, 100), axes = FALSE)
graphics.off()
PictureE <- readPNG("barsE.png")
png("barsU.png")
barplot(height = c(90, 10), horiz = "TRUE", col = c("orange", "blue"), xlim = c(0, 100), axes = FALSE)
graphics.off()
PictureU <- readPNG("barsU.png")
g <- graph(c("S", "C", "U", "C", "S", "E"))
l <- matrix(data = c(20, 100, 50, 30, 110, 90, 150, 20), 4, 2, byrow = TRUE, dimnames = list(nodes = c("S", "C", "U", "E"), coord = c("x", "y")))
plot(g, vertex.shape = "rectangle", vertex.size = 62, vertex.size2 = 31, layout = l, vertex.label = c("S", "C", "U", "E"), vertex.color = "white")
rasterImage(PictureC, -0.845, -0.9, -0.245, -0.6)
rasterImage(PictureE, 0.7, -1.15, 1.3, -0.85)
rasterImage(PictureS, -1.3, 0.85, -0.7, 1.15)
rasterImage(PictureU, 0.08, 0.6, 0.68, 0.9)
text(-1, 1.25, "S: Serial Killer?")
text(-1.55, 1.05, "Yes: 70%")
text(-1.55, 0.95, " No: 30%")
text(+1, -1.25, "E: Evidence of one Murder")
text(0.45, -0.95, "Yes: 30%")
text(0.45, -1.05, " No: 70%")
text(0.4, 1, "U: Unknown unknowns")
text(-0.15, 0.8, "Yes: 10%")
text(-0.15, 0.7, " No: 90%")
text(-0.5, -1, "C: Cluster of events")
text(-1.075, -0.7, "Yes: 90%")
text(-1.075, -0.8, " No: 10%")