devtools::load_all()
#> ℹ Loading testerror
#> ℹ Re-compiling testerror (debug build)
#> ── R CMD INSTALL ───────────────────────────────────────────────────────────────
#> * installing *source* package ‘testerror’ ...
#> ** this is package ‘testerror’ version ‘0.1.0’
#> ** using staged installation
#> * DONE (testerror)
library(tidyverse)
#> ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
#> ✔ dplyr 1.2.1 ✔ readr 2.2.0
#> ✔ forcats 1.0.1 ✔ stringr 1.6.0
#> ✔ ggplot2 4.0.3 ✔ tibble 3.3.1
#> ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
#> ✔ purrr 1.2.2
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ readr::edition_get() masks testthat::edition_get()
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag() masks stats::lag()
#> ✖ readr::local_edition() masks testthat::local_edition()
#> ℹ 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 /tmp/Rtmpk0EDlu/Rbuild19607ab3e45/testerror
source(here::here("vignettes/formatting.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"
)
#> recompiling to avoid crashing R session#> Warning: There were 59 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
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
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.38% [0.00% — 3.43%] | 75.94% [55.55% — 88.54%] | 97.25% [95.72% — 99.32%] | bayes (logit) |
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"
)
#> recompiling to avoid crashing R session
#> : 4.865 seconds (Total)
#> Chain 2: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.03% [0.00% — 0.43%] | 97.48% [77.67% — 99.92%] | 99.81% [99.48% — 99.96%] | bayes (logit) |
| 3 | 30/796 (3.77%) | 5.62% [3.10% — 21.46%] | 61.73% [15.85% — 96.45%] | 99.80% [99.28% — 99.97%] | bayes (logit) |
| 4 | 2/796 (0.25%) | 0.02% [0.00% — 0.44%] | 96.09% [61.41% — 99.89%] | 99.81% [99.51% — 99.97%] | bayes (logit) |
| 5 | 1/796 (0.13%) | 0.01% [0.00% — 0.32%] | 93.69% [33.37% — 99.88%] | 99.79% [99.49% — 99.94%] | bayes (logit) |
| 6A | 4/796 (0.50%) | 0.15% [0.00% — 0.96%] | 93.94% [35.36% — 99.87%] | 99.75% [99.39% — 99.96%] | bayes (logit) |
| 6B | 0/796 (0.00%) | 0.01% [0.00% — 0.31%] | 87.68% [6.73% — 99.84%] | 99.91% [99.67% — 99.99%] | bayes (logit) |
| 7F | 7/796 (0.88%) | 0.35% [0.00% — 1.27%] | 96.59% [65.45% — 99.92%] | 99.60% [99.12% — 99.91%] | bayes (logit) |
| 9V | 1/796 (0.13%) | 0.02% [0.00% — 0.32%] | 95.30% [49.03% — 99.88%] | 99.86% [99.58% — 99.98%] | bayes (logit) |
| 14 | 1/796 (0.13%) | 0.02% [0.00% — 0.29%] | 97.56% [76.99% — 99.93%] | 99.85% [99.59% — 99.97%] | bayes (logit) |
| 18C | 3/796 (0.38%) | 0.02% [0.00% — 0.67%] | 89.24% [8.77% — 99.84%] | 99.62% [99.23% — 99.85%] | bayes (logit) |
| 19A | 2/796 (0.25%) | 0.02% [0.00% — 0.43%] | 97.27% [74.20% — 99.92%] | 99.81% [99.50% — 99.96%] | bayes (logit) |
| 19F | 3/796 (0.38%) | 0.07% [0.00% — 1.09%] | 89.86% [9.81% — 99.84%] | 99.78% [99.41% — 99.96%] | bayes (logit) |
| 23F | 4/796 (0.50%) | 0.16% [0.00% — 0.88%] | 95.29% [50.11% — 99.88%] | 99.75% [99.39% — 99.96%] | bayes (logit) |
| PCV13 | 59/796 (7.41%) | 6.77% [4.26% — 22.10%] | 67.35% [22.47% — 94.77%] | 96.91% [95.91% — 97.79%] | bayes (logit) |
| Binax | 26/786 (3.31%) | 0.38% [0.00% — 3.43%] | 75.94% [55.55% — 88.54%] | 97.25% [95.72% — 99.32%] | bayes (logit) |
# 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"
)
#> recompiling to avoid crashing R session
#> recompiling to avoid crashing R session
#> 1: Elapsed Time: 1.25 seconds (Warm-up)
#> Chain 1: 1.004 seconds (Sampling)
#> Chain 1: 2.254 seconds (Total)
#> Chain 1:
#> Warning: There were 2 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.42%] | 97.53% [76.82% — 99.92%] | 99.80% [99.47% — 99.96%] | bayes (logit) |
| 3 | 30/796 (3.77%) | 5.02% [3.14% — 8.30%] | 66.88% [42.33% — 87.26%] | 99.77% [99.25% — 99.97%] | bayes (logit) |
| 4 | 2/796 (0.25%) | 0.02% [0.00% — 0.48%] | 96.16% [59.08% — 99.91%] | 99.80% [99.48% — 99.96%] | bayes (logit) |
| 5 | 1/796 (0.13%) | 0.01% [0.00% — 0.29%] | 93.92% [33.01% — 99.87%] | 99.78% [99.45% — 99.94%] | bayes (logit) |
| 6A | 4/796 (0.50%) | 0.13% [0.00% — 0.95%] | 93.59% [34.60% — 99.89%] | 99.73% [99.35% — 99.95%] | bayes (logit) |
| 6B | 0/796 (0.00%) | 0.01% [0.00% — 0.33%] | 87.53% [6.54% — 99.85%] | 99.90% [99.65% — 99.99%] | bayes (logit) |
| 7F | 7/796 (0.88%) | 0.28% [0.00% — 1.24%] | 96.68% [66.28% — 99.92%] | 99.55% [99.05% — 99.90%] | bayes (logit) |
| 9V | 1/796 (0.13%) | 0.02% [0.00% — 0.31%] | 95.37% [51.51% — 99.87%] | 99.85% [99.56% — 99.97%] | bayes (logit) |
| 14 | 1/796 (0.13%) | 0.01% [0.00% — 0.29%] | 97.48% [74.43% — 99.93%] | 99.85% [99.56% — 99.97%] | bayes (logit) |
| 18C | 3/796 (0.38%) | 0.02% [0.00% — 0.63%] | 88.98% [9.14% — 99.83%] | 99.59% [99.20% — 99.82%] | bayes (logit) |
| 19A | 2/796 (0.25%) | 0.02% [0.00% — 0.43%] | 97.18% [75.18% — 99.91%] | 99.80% [99.49% — 99.96%] | bayes (logit) |
| 19F | 3/796 (0.38%) | 0.06% [0.00% — 0.89%] | 88.65% [8.91% — 99.86%] | 99.76% [99.40% — 99.95%] | bayes (logit) |
| 23F | 4/796 (0.50%) | 0.14% [0.00% — 0.82%] | 95.32% [51.51% — 99.87%] | 99.73% [99.33% — 99.96%] | bayes (logit) |
| PCV13 | 59/796 (7.41%) | 6.03% [4.28% — 9.03%] | 72.09% [50.95% — 86.26%] | 96.67% [95.73% — 97.52%] | bayes (logit) |
| Binax | 26/786 (3.31%) | 0.38% [0.00% — 3.43%] | 75.94% [55.55% — 88.54%] | 97.25% [95.72% — 99.32%] | bayes (logit) |