---
title: "Bayesian adjustment of panel test error" 
output: html_document
vignette: > 
  %\VignetteIndexEntry{Bayesian adjustment of panel test error} 
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
devtools::load_all()
library(tidyverse)
library(testerror)
options(mc.cores = 2) # parallel::detectCores())
rstan::rstan_options(auto_write = FALSE)
here::i_am("vignettes/testerror.Rmd")
source(here::here("vignettes/formatting.R"))
```

# BinaxNOW test results

```{r}

binax_pos = 26
binax_n = 786

# Fit beta distributions to sinclair 2013 data:
# We apply a widening to increase uncertainty in sensitivity
binax_sens = beta_params(median = 0.740, lower = 0.666, upper = 0.823, widen = 5)
binax_spec = beta_params(median = 0.972, lower=  0.927, upper = 0.998)

tmp1 = testerror::bayesian_true_prevalence_model(
  pos_obs = binax_pos,
  n_obs = binax_n,
  sens = binax_sens,
  spec = binax_spec,
  model_type = "logit"
)

binax_summ = tmp1$summary %>%
  dplyr::transmute(
  Test = "Binax",
  Positivity = sprintf("%d/%d (%1.2f%%)", binax_pos, binax_n, binax_pos/binax_n*100),
  `Estimated prevalence` = prevalence.label,
  `Sensitivity` = sens.label,
  `Specificity` = spec.label,
  `Method` = prevalence.method
  )

binax_summ %>% default_table()
```

# UAD1 results

For this description we are going to assume we have only count data for panel
and components, and we have limited information about component tests.

```{r}

# data from Forstner et al (2019)
components = tibble::tibble(
  serotype = factor(c("1", "3", "4", "5", "6A", "6B", "7F", "9V", "14", "18C", "19A", "19F", "23F")),
  pos = c(2, 30, 2, 1, 4, 0, 7, 1, 1, 3, 2, 3, 4),
  n = 796
)

uad1_pos = 59
uad1_n = 796

# control group data
# negatives at 
uad1_spec = spec_prior() %>% 
  # data from Pride et al
  update_posterior(0,17)

# Datd from Pride et al
uad1_controls = tibble::tibble(
  serotype = factor(c("1", "3", "4", "5", "6A", "6B", "7F", "9V", "14", "18C", "19A", "19F", "23F")),
  false_neg_diseased = c(0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L),
  n_diseased = c(7, 1, 3, 1, 1, 0, 4, 2, 7, 0, 6, 0, 2),
  # In the design of UAD1 the cut off is calibrated on 400 negative samples and 
  # set to make 2 positive
  
  # Forstner et al (2019) control group
  # In the control group of 397 non-CAP patients, 1 patient had to
  # be excluded because of indeterminable results (0.3%, see Fig. 1)
  # and SSUAD was positive in 3 patients (0.8%). Of those 3 non-CAP
  # patients, 2 were positive for serotype 18C and 1 person for two
  # pneumococcal serotypes (5, 7F).

  false_pos_controls = 2+c(0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 2L, 0L, 0L, 0L),
  n_controls = 400+396
)


uad1_sens = sens_prior() %>% 
  # data from Pride et al
  update_posterior(17,17)

```

With this we can use a bayesian model to construct adjusted estimates

```{r}

corr3 = testerror::bayesian_panel_true_prevalence_model(
  panel_pos_obs = uad1_pos,
  panel_n_obs = uad1_n,
  panel_name = "PCV13",
  
  pos_obs = components$pos,
  n_obs = components$n,
  test_names = components$serotype,
  
  false_pos_controls = uad1_controls$false_pos_controls,
  n_controls = uad1_controls$n_controls,
  false_neg_diseased = uad1_controls$false_neg_diseased,
  n_diseased = uad1_controls$n_diseased,
  
  panel_sens = uad1_sens,
  model_type = "logit"
)
```

```{r}
corr3$summary %>% 
  mutate(
    pos_obs = c(components$pos, uad1_pos),
    n_obs = c(components$n, uad1_n)
  ) %>%
  dplyr::transmute(
  Test = test,
  Positivity = sprintf("%d/%d (%1.2f%%)", pos_obs, n_obs, pos_obs/n_obs*100),
  `Estimated prevalence` = prevalence.label,
  `Sensitivity` = sens.label,
  `Specificity` = spec.label,
  `Method` = prevalence.method
) %>% 
  bind_rows(binax_summ) %>%
  default_table() %>% huxtable::set_top_border(row = huxtable::final(2), col=huxtable::everywhere, value = 1)

```

# Sensitivity analysis with different UAD test parameters


```{r}

# senstivity / specificity from Kakiuchi et al 

kak_uad1_sens = beta_params(median = 0.741, lower = 0.537, upper = 0.889)
kak_uad1_spec = beta_params(median = 0.954, lower = 0.917, upper = 0.978) %>% 
  # data from Pride et al
  update_posterior(17,17)

corr4 = testerror::bayesian_panel_true_prevalence_model(
  panel_pos_obs = uad1_pos,
  panel_n_obs = uad1_n,
  panel_name = "PCV13",
  
  pos_obs = components$pos,
  n_obs = components$n,
  test_names = components$serotype,
  
  false_pos_controls = uad1_controls$false_pos_controls,
  n_controls = uad1_controls$n_controls,
  false_neg_diseased = uad1_controls$false_neg_diseased,
  n_diseased = uad1_controls$n_diseased,
  
  panel_sens = kak_uad1_sens,
  panel_spec = kak_uad1_spec,
  model_type="logit"
)

corr4$summary %>% 
  mutate(
    pos_obs = c(components$pos, uad1_pos),
    n_obs = c(components$n, uad1_n)
  ) %>%
  dplyr::transmute(
  Test = test,
  Positivity = sprintf("%d/%d (%1.2f%%)", pos_obs, n_obs, pos_obs/n_obs*100),
  `Estimated prevalence` = prevalence.label,
  `Sensitivity` = sens.label,
  `Specificity` = spec.label,
  `Method` = prevalence.method
) %>% bind_rows(
  binax_summ
) %>% 
  default_table() %>% huxtable::set_top_border(row = huxtable::final(2), col=huxtable::everywhere, value = 1)

```