Probability of PE before and after V/Q scan
post.test.pe <- function(p) {
lr.1 <- (102/251)/(14/480)
lr.2 <- (105/251)/(217/480)
lr.3 <- (39/251)/(199/480)
lr.4 <- (5/251)/(50/480)
pre <- p
pre.odds <- p/(1 - p)
o.t1 <- pre.odds * lr.1
o.t2 <- pre.odds * lr.2
o.t3 <- pre.odds * lr.3
o.t4 <- pre.odds * lr.4
data.frame(p = p, high = o.t1/(1 + o.t1), intermediate = o.t2/(1 + o.t2),
low = o.t3/(1 + o.t3), normal = o.t4/(1 + o.t4))
}
dat.post <- post.test.pe(seq(0, 1, 0.001))
library(reshape)
dat.post.melt <- melt(dat.post, "p")
library(ggplot2)
ggplot(dat.post.melt, aes(p, value, color = variable)) + geom_line(lwd = 2) +
scale_color_discrete(name = "result") + xlab("Pretest probability") + ylab("Posttest probability") +
geom_hline(y = 0.05)
Thresholds each test result. If pretest probability is lowe than this, the result does not change the posttest probability high enough to justify treatment.
lr.list <- list(high = (102/251)/(14/480), intermediate = (105/251)/(217/480),
low = (39/251)/(199/480), normal = (5/251)/(50/480))
post.test.prob <- function(pre.p, lr) {
pre.odds <- pre.p/(1 - pre.p)
post.odds <- pre.odds * lr # Bayes' formula
post.p <- post.odds/(1 + post.odds)
post.p
}
sapply(simplify = FALSE, lr.list, function(lr) {
f.cutoff <- function(pre.p) post.test.prob(pre.p, lr) - 0.05
uniroot(f.cutoff, c(0, 0.99999999))$root
})
## $high
## [1] 0.003763
##
## $intermediate
## [1] 0.05382
##
## $low
## [1] 0.1231
##
## $normal
## [1] 0.2158
##