# A simple R function to compute pathogenicity probabilities, as in
# the recent paper by M. Gollob, J.S. Rosenthal, and K.E. Thorpe
# available at: http://probability.ca/jeff/ftpdir/pathprob.pdf
# To download R for free, see: www.r-project.org
# To calculate the probabilities online, see: www.probability.ca/pathprob
pathprob = function(p,q,r, verbosity=2) {
y = (r-q)/(1-q)
z = p*(1-y)/(1-p*y)
w = p * y / (p * r + (1-p) * q)
u = (1-z) * w
if (verbosity >= 2) {
cat("\nIf the disease prevalence (p) is equal to", p)
cat("\nand the healthy variant rate (q) is equal to", q)
cat("\nand the disease variant rate (r) is equal to", r)
cat("\nthen:\n")
cat("\tThe unconditional probability that the variant is disease-causing is", w, "\n")
cat("\tThe unconditional probability that a patient with the variant will get the disease is", z+u, "\n")
cat("\tThe probability, conditional on disease, that the variant was the *SOLE* cause of the disease is", u/(z+u), "\n")
cat("\tAnd, conditional on disease, the pathogenicity probability that the variant was *A* cause of the disease is", w/(z+u), "\n\n")
} else if (verbosity >= 1) {
cat("pathogenicity probability =", w / (z+u), "\n")
}
return( w / (z+u) )
}
# Examples:
pathprob(1/4000, 0.02, 0.4)
pathprob(1/400, 0.01, 0.4, 1)
pathprob(1/50, 0.1, 0.2, 1)
pathprob(1/400, 0.1, 0.15, 1)
pathprob(1/400, 0.1, 1, 1)
pathprob(1/400, 0, 0.5, 1)
pathprob(1/400, 0.1, 0.1, 1)