# The problem of false positives: Antibody tests in times of Corona
# R code from https://philippmasur.de/blog/2020-04-14-false-positives-antibody-tests/

library(tidyverse)   # for data wrangling
library(MonteCarlo)  # for the simulation

test_sim <- function(rate = 0.001, 
                     pop = 100000, 
                     sensitivity = .95, 
                     specificity = .95) {

  # Create population with an infected subpopulation
  data <- data.frame(id = 1:pop, corona = c(rep(1, pop*rate), rep(0, pop-rate*pop)))

  # Sample individual person randomly
  person <- sample_n(data, 1)

  # "test" the person
  if(person$corona == 1) {

    prob_sens <- sample(c(rep("correct", (sensitivity*100)), 
                     rep("false", (100-sensitivity*100))), 1)

    if(prob_sens == "correct") {
      person$test_result <- 1
    } else {
      person$test_result <- 0
    }

  } else {

    prob_spec <- sample(c(rep("correct", (specificity*100)), 
                     rep("false", (100-specificity*100))), 1)

    if(prob_spec == "correct") {
      person$test_result <- 0
    } else {
      person$test_result <- 1
    }
  }

  # return "true" infection and test result
  return(list("infected" = person$corona,
              "test" = person$test_result))
}

set.seed(42)
test_sim()

# Define simulation parameters
rate_grid <- c(0.001, 0.010, 0.050, 0.100, 0.200) # Rate of infection 
pop_grid <- c(100000)                             # Population size    
sens_grid <- c(.85, .90, .95, .99)                # Sensitivity of the test
spec_grid <- c(.85, .90, .95, .99)                # Specificity of the test

# Collect simulation parameters in list
(param_list <- list("rate" = rate_grid, 
                   "pop" = pop_grid,
                   "sensitivity" = sens_grid,
                   "specificity" = spec_grid))

# Run simulation
result <- MonteCarlo(func = test_sim,            # pass test function
                     nrep = 1000,                # number of tests
                     param_list = param_list)    # provide parameters

# Create result data frame
df <- MakeFrame(result) %>% as_tibble %>%
  mutate(test_result = ifelse(infected == 0 & test == 0, "correct negative",
                              ifelse(infected == 1 & test == 1, "correct positive",
                                     ifelse(infected == 0 & test == 1, "false positive", 
                                            ifelse(infected == 1 & test == 0, "false negative", NA)))))

head(df)

# plot
df %>%
  filter(test_result != "correct negative") %>%
  group_by(rate, sensitivity, specificity) %>%
  summarise("correct positive" = length(which(test_result == "correct positive")),
            "false negative" = length(which(test_result == "false negative")),
            "false positive" = length(which(test_result == "false positive"))) %>%
  gather(key, value, -rate, -sensitivity, -specificity) %>%
  mutate(value = value/1000) %>%     # share of the 1,000 simulated tests per condition
  ggplot(aes(x = specificity, y = value, color = key)) +
  geom_point() +
  geom_line() +
  scale_color_brewer(palette = "Dark2") +
  facet_grid(rate ~ sensitivity, labeller = label_both) +
  scale_y_continuous(labels = scales::percent) +
  theme_bw() +
  labs(y = "share of all tests", color = "test result")