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