false-positive.R

Sys.setlocale(locale = "cht") # for windows os
library(magrittr)
library(ggpubr)
library(data.table)

p.sensitive <- c(0.8, 0.9, 0.999999)
p.specific <- c(0.9, 0.99, 0.999, 0.9999)
p.prevalence.rate <-
  c(
    seq(0.000001, 0.000009, 0.0000001),
    seq(0.00001, 0.00009, 0.000001),
    seq(0.0001, 0.0009, 0.00001),
    seq(0.001, 0.009, 0.0001),
    seq(0.01, 0.09, 0.001),
    seq(0.1, 0.9, 0.01),
    1
  )

myF <-
  function(p.sensitive,
           p.specific,
           p.prevalence.rate) {
    p.sensitive * p.prevalence.rate /
      (p.sensitive * p.prevalence.rate + (1 - p.specific) * (1 - p.prevalence.rate))
  }

dt <-
  expand.grid(p.sensitive, p.specific, p.prevalence.rate) %>%
  as.data.table %>%
  set_names(c("p.sensitive", "p.specific", "p.prevalence.rate")) %>%
  .[, true.positive := myF(p.sensitive, p.specific, p.prevalence.rate)] %>%
  .[, false.positive := 1 - true.positive] %>%
  .[, `靈敏度` := p.sensitive %>% as.factor] %>%
  .[, `特異度` := p.specific %>% as.factor] %>%
  setkey(p.prevalence.rate, p.sensitive, p.specific)

ggplot(dt,
       aes(p.prevalence.rate,
           false.positive)) +
  geom_line(aes(
    size = `靈敏度`,
    color = `特異度`,
    group = interaction(`靈敏度`, `特異度`)
  )) +
  scale_size_manual(values = c(0.8, 1.25, 2)) +
  scale_color_manual(values = c('#e66101', '#fdb863', '#b2abd2', '#5e3c99')) +
  scale_x_continuous(
    "疾病盛行率（%）",
    trans = "log10",
    breaks = c(0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1, 1),
    labels = (c(
      0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1, 1
    ) * 100) %>% formatC(
      format = "f",
      digits = 4,
      drop0trailing = T
    )
  ) +
  scale_y_continuous(
    "偽陽性機率\n（某人測得陽性反應但實際並無感染的機率；越低越好）",
    breaks = seq(0, 1, 0.2),
    labels = seq(0, 1, 0.2)
  ) +
  theme_pubr(16,
             border = T,
             legend = "right") +
  theme(plot.title = element_text(size = 24, face = "bold"), 
        plot.subtitle = element_text(size = 16),
        plot.caption = element_text(size = 16), plot.caption.position = "plot"
        ) +
  labs(title  = "盛行率過低時的偽陽性膨脹現象",
       subtitle = "即使有超級準確的神之測試方法（深紫色粗線），\n在盛行過低時偽陽性結果的機率仍過高。",
       caption = "靈敏度：測試方法對實際是陽性樣本正確地判斷為陽性的機率\n特異度：測試方法對實際是陰性樣本正確地判斷為陰性的機率") +
  windows(12, 8, antialias = "cleartype")
ggsave("false-positive.png")