Skip to contents
devtools::load_all()
#>  Loading testerror
library(tidyverse)
#> ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
#>  dplyr     1.1.2      readr     2.1.4
#>  forcats   1.0.0      stringr   1.5.0
#>  ggplot2   3.5.0      tibble    3.2.1
#>  lubridate 1.9.3      tidyr     1.3.0
#>  purrr     1.0.1     
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#>  readr::edition_get()   masks testthat::edition_get()
#>  dplyr::filter()        masks stats::filter()
#>  purrr::is_null()       masks testthat::is_null()
#>  dplyr::lag()           masks stats::lag()
#>  readr::local_edition() masks testthat::local_edition()
#>  dplyr::matches()       masks tidyr::matches(), testthat::matches()
#>  Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(testerror)
options(mc.cores = 2) # parallel::detectCores())
rstan::rstan_options(auto_write = FALSE)
here::i_am("vignettes/testerror.Rmd")
#> here() starts at /home/vp22681/Dropbox/Git/testerror
source(here::here("vignettes/formatting.R"))

BinaxNOW test results


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"
)
#> Warning: There were 11 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems

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()
Test Positivity Estimated prevalence Sensitivity Specificity Method
Binax 26/786 (3.31%) 0.32% [0.00% — 3.25%] 75.13% [53.83% — 89.57%] 97.16% [95.70% — 98.65%] bayes (logit)

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.


# 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


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"
)
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)
Test Positivity Estimated prevalence Sensitivity Specificity Method
1 2/796 (0.25%) 0.02% [0.00% — 0.43%] 97.65% [78.00% — 99.93%] 99.81% [99.50% — 99.97%] bayes (logit)
3 30/796 (3.77%) 5.70% [3.08% — 20.06%] 60.80% [16.48% — 96.17%] 99.80% [99.30% — 99.97%] bayes (logit)
4 2/796 (0.25%) 0.03% [0.00% — 0.45%] 96.04% [59.54% — 99.90%] 99.81% [99.50% — 99.96%] bayes (logit)
5 1/796 (0.13%) 0.01% [0.00% — 0.30%] 93.53% [32.94% — 99.86%] 99.79% [99.48% — 99.94%] bayes (logit)
6A 4/796 (0.50%) 0.14% [0.00% — 1.00%] 93.94% [37.82% — 99.86%] 99.75% [99.36% — 99.96%] bayes (logit)
6B 0/796 (0.00%) 0.01% [0.00% — 0.29%] 89.30% [6.89% — 99.84%] 99.90% [99.69% — 99.99%] bayes (logit)
7F 7/796 (0.88%) 0.37% [0.00% — 1.30%] 96.68% [66.50% — 99.92%] 99.61% [99.09% — 99.92%] bayes (logit)
9V 1/796 (0.13%) 0.02% [0.00% — 0.31%] 95.18% [48.15% — 99.89%] 99.86% [99.60% — 99.97%] bayes (logit)
14 1/796 (0.13%) 0.02% [0.00% — 0.29%] 97.57% [76.52% — 99.93%] 99.86% [99.59% — 99.97%] bayes (logit)
18C 3/796 (0.38%) 0.02% [0.00% — 0.60%] 89.52% [8.35% — 99.86%] 99.62% [99.24% — 99.85%] bayes (logit)
19A 2/796 (0.25%) 0.03% [0.00% — 0.41%] 97.28% [74.03% — 99.91%] 99.81% [99.51% — 99.96%] bayes (logit)
19F 3/796 (0.38%) 0.07% [0.00% — 0.97%] 88.86% [10.17% — 99.81%] 99.78% [99.43% — 99.96%] bayes (logit)
23F 4/796 (0.50%) 0.15% [0.00% — 0.90%] 95.52% [47.99% — 99.91%] 99.75% [99.39% — 99.96%] bayes (logit)
PCV13 59/796 (7.41%) 6.91% [4.24% — 21.10%] 67.01% [23.96% — 94.51%] 96.92% [95.91% — 97.82%] bayes (logit)
Binax 26/786 (3.31%) 0.32% [0.00% — 3.25%] 75.13% [53.83% — 89.57%] 97.16% [95.70% — 98.65%] bayes (logit)

Sensitivity analysis with different UAD test parameters


# 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"
)
#> Warning: There were 12 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems

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)
Test Positivity Estimated prevalence Sensitivity Specificity Method
1 2/796 (0.25%) 0.02% [0.00% — 0.41%] 97.50% [76.57% — 99.92%] 99.80% [99.48% — 99.96%] bayes (logit)
3 30/796 (3.77%) 5.01% [3.13% — 8.30%] 67.25% [42.36% — 88.71%] 99.77% [99.27% — 99.97%] bayes (logit)
4 2/796 (0.25%) 0.02% [0.00% — 0.43%] 96.27% [59.57% — 99.92%] 99.79% [99.49% — 99.96%] bayes (logit)
5 1/796 (0.13%) 0.01% [0.00% — 0.30%] 93.68% [32.70% — 99.88%] 99.78% [99.47% — 99.94%] bayes (logit)
6A 4/796 (0.50%) 0.14% [0.00% — 0.93%] 93.59% [33.82% — 99.88%] 99.73% [99.33% — 99.95%] bayes (logit)
6B 0/796 (0.00%) 0.01% [0.00% — 0.28%] 88.34% [7.84% — 99.82%] 99.90% [99.66% — 99.99%] bayes (logit)
7F 7/796 (0.88%) 0.30% [0.00% — 1.19%] 96.62% [64.59% — 99.91%] 99.55% [99.07% — 99.89%] bayes (logit)
9V 1/796 (0.13%) 0.02% [0.00% — 0.31%] 95.15% [47.46% — 99.89%] 99.85% [99.55% — 99.97%] bayes (logit)
14 1/796 (0.13%) 0.01% [0.00% — 0.29%] 97.42% [76.91% — 99.93%] 99.85% [99.57% — 99.97%] bayes (logit)
18C 3/796 (0.38%) 0.02% [0.00% — 0.61%] 88.23% [8.68% — 99.83%] 99.59% [99.20% — 99.84%] bayes (logit)
19A 2/796 (0.25%) 0.02% [0.00% — 0.41%] 97.29% [75.59% — 99.91%] 99.80% [99.49% — 99.96%] bayes (logit)
19F 3/796 (0.38%) 0.06% [0.00% — 0.88%] 89.31% [8.65% — 99.87%] 99.76% [99.39% — 99.95%] bayes (logit)
23F 4/796 (0.50%) 0.13% [0.00% — 0.88%] 95.35% [48.50% — 99.91%] 99.73% [99.36% — 99.96%] bayes (logit)
PCV13 59/796 (7.41%) 6.01% [4.28% — 8.97%] 72.12% [51.75% — 86.58%] 96.67% [95.69% — 97.50%] bayes (logit)
Binax 26/786 (3.31%) 0.32% [0.00% — 3.25%] 75.13% [53.83% — 89.57%] 97.16% [95.70% — 98.65%] bayes (logit)