This is a little Bayesian analysis of the Daniela Poggiali case.

According to a statistical report by an Italian forensic pathologist and an Italian epidemiologist, called “TM” in the sequel, Daniela killed about 80 patients in a two year period. They have two sources of evidence. One is data on deaths during Daniela’s shifts and not during her shifts. It seems that the death rate is elevated by a factor of about 2 when she is on duty. The p-value for an elevated rate is 1 in a thousand. The other source of evidence is the post-mortem concentration of Potassium in the ocular fluid of an elderly and terminally ill patient, Mrs Calderoni, a number of hours after her otherwise not unexpected decease.

In this R markdown document I only represent a subjective Bayesian account of part of this case.

If TM are right then Daniela Poggiali is actually the cause of about 80 excess deaths. Of course, maybe she killed out of pity for the suffering of terminally ill patients, but anyway, since this was not authorised by medical specialists or family members or indeed the patient, it is a terrible crime. It is a crime which would make her the almost the most deadly health-care serial killer (HCSK) ever. That would not only be horrfic. It would also be amazing. I estimate the chance that a particular nurse in Italy is the most deadly HCSK ever at about 1 in a million.

Of course, nurses also do occasionally kill patients, whether our of pity, personal animosity, or for personal gain (the opportunity to steal some gold trinkets), or whatever. I estimate the rate of killings of this kind at perhaps 1 per 1000 deaths. Note that it is certainly the case that most deaths in hospitals are due to medical errors, though mainly in situations where the patient was expected to die soon, anyway. For various reasons, mainly economic, medical errors are hardly ever admitted and even less often registered as such. The legal and economic consequences are too great – and some rate of error is of course unavoidable.

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")

DS <-cptable(~DanielaSerialKiller, levels = yn, values = c(1, 10^6 - 1))
DC.DS <- cptable(~DanielaKilledCalderoni|DanielaSerialKiller, levels = yn, values = c(10^2 - 1, 1, 1, 10^3 - 1))
TM.DS <- cptable(~TMconcludeDSK|DanielaSerialKiller, levels = yn, values = c(10^3 - 1, 1, 1, 10^3 - 1))
daniela_cpt <- compileCPT(DS, DC.DS, TM.DS)
daniela_cpt
## cpt_spec with probabilities:
##  P( DanielaSerialKiller )
##  P( DanielaKilledCalderoni | DanielaSerialKiller )
##  P( TMconcludeDSK | DanielaSerialKiller )
daniela_cpt$DanielaSerialKiller
## DanielaSerialKiller
##      yes       no 
## 0.000001 0.999999
daniela_cpt$DanielaKilledCalderoni
##                       DanielaSerialKiller
## DanielaKilledCalderoni  yes    no
##                    yes 0.99 0.001
##                    no  0.01 0.999
daniela_cpt$TMconcludeDSK
##              DanielaSerialKiller
## TMconcludeDSK   yes    no
##           yes 0.999 0.001
##           no  0.001 0.999
daniela_bn <- grain(daniela_cpt)
daniela_bn <- propagate(daniela_bn)
querygrain(daniela_bn, nodes = c("DanielaSerialKiller", "DanielaKilledCalderoni"), type = "joint")
##                    DanielaKilledCalderoni
## DanielaSerialKiller         yes         no
##                 yes 0.000000990 0.00000001
##                 no  0.000999999 0.99899900
daniela_bn2 <- setEvidence(daniela_bn, evidence = list(TMconcludeDSK = "yes"))
querygrain(daniela_bn2, nodes = c("DanielaSerialKiller", "DanielaKilledCalderoni"), type = "joint")
##                    DanielaKilledCalderoni
## DanielaSerialKiller         yes          no
##                 yes 0.000988024 9.98004e-06
##                 no  0.000999002 9.98003e-01
saveHuginNet(daniela_bn, "daniela.net")