diff --git a/_pkgdown.yml b/_pkgdown.yml
index 082351f..ecae09f 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -36,6 +36,7 @@ articles:
- usability_limits
- check_poisson
- visits_or_days
+ - gsm_perf
reference:
- title: S3 Interface
diff --git a/docs/404.html b/docs/404.html
index 72bd360..f31c8d5 100644
--- a/docs/404.html
+++ b/docs/404.html
@@ -45,7 +45,7 @@
Test Data and Standard Processingdf_sim_sites <- sim_sites(df_site, df_visit)
df_visit
-## # A tibble: 59,221 × 9
-## patnum site_number is_ur max_visit_mean max_vi…¹ ae_pe…² visit n_ae study…³
-## <chr> <chr> <lgl> <dbl> <dbl> <dbl> <int> <int> <chr>
-## 1 P000001 S0001 TRUE 20 4 0.25 1 0 ae_per…
-## 2 P000001 S0001 TRUE 20 4 0.25 2 1 ae_per…
-## 3 P000001 S0001 TRUE 20 4 0.25 3 2 ae_per…
-## 4 P000001 S0001 TRUE 20 4 0.25 4 2 ae_per…
-## 5 P000001 S0001 TRUE 20 4 0.25 5 2 ae_per…
-## 6 P000001 S0001 TRUE 20 4 0.25 6 2 ae_per…
-## 7 P000001 S0001 TRUE 20 4 0.25 7 2 ae_per…
-## 8 P000001 S0001 TRUE 20 4 0.25 8 4 ae_per…
-## 9 P000001 S0001 TRUE 20 4 0.25 9 4 ae_per…
-## 10 P000001 S0001 TRUE 20 4 0.25 10 4 ae_per…
-## # … with 59,211 more rows, and abbreviated variable names ¹max_visit_sd,
-## # ²ae_per_visit_mean, ³study_id
+## # A tibble: 58,367 × 9
+## patnum site_number is_ur max_visit_mean max_visit_sd ae_per_visit_mean visit
+## <chr> <chr> <lgl> <dbl> <dbl> <dbl> <int>
+## 1 P000001 S0001 TRUE 20 4 0.25 1
+## 2 P000001 S0001 TRUE 20 4 0.25 2
+## 3 P000001 S0001 TRUE 20 4 0.25 3
+## 4 P000001 S0001 TRUE 20 4 0.25 4
+## 5 P000001 S0001 TRUE 20 4 0.25 5
+## 6 P000001 S0001 TRUE 20 4 0.25 6
+## 7 P000001 S0001 TRUE 20 4 0.25 7
+## 8 P000001 S0001 TRUE 20 4 0.25 8
+## 9 P000001 S0001 TRUE 20 4 0.25 9
+## 10 P000001 S0001 TRUE 20 4 0.25 10
+## # ℹ 58,357 more rows
+## # ℹ 2 more variables: n_ae <int>, study_id <chr>
df_site
## # A tibble: 300 × 6
-## study_id site_number n_pat n_pat_with_med75 visit_med75 mean_ae_s…¹
-## <chr> <chr> <int> <int> <dbl> <dbl>
-## 1 ae_per_visit: 0.05 S0001 10 8 16 0
-## 2 ae_per_visit: 0.05 S0002 10 9 16 0.333
-## 3 ae_per_visit: 0.05 S0003 10 10 15 0.3
-## 4 ae_per_visit: 0.05 S0004 10 9 16 0.111
-## 5 ae_per_visit: 0.05 S0005 10 9 18 0.556
-## 6 ae_per_visit: 0.05 S0006 10 8 17 0
-## 7 ae_per_visit: 0.05 S0007 10 10 16 0.6
-## 8 ae_per_visit: 0.05 S0008 10 10 15 0.5
-## 9 ae_per_visit: 0.05 S0009 10 9 17 0.111
-## 10 ae_per_visit: 0.05 S0010 10 10 16 0.5
-## # … with 290 more rows, and abbreviated variable name ¹mean_ae_site_med75
+## study_id site_number n_pat n_pat_with_med75 visit_med75 mean_ae_site_med75
+## <chr> <chr> <int> <int> <dbl> <dbl>
+## 1 ae_per_vis… S0001 10 9 16 0.556
+## 2 ae_per_vis… S0002 10 10 14 0.3
+## 3 ae_per_vis… S0003 10 8 15 0.375
+## 4 ae_per_vis… S0004 10 8 14 0.375
+## 5 ae_per_vis… S0005 10 9 17 0.333
+## 6 ae_per_vis… S0006 10 7 15 0.143
+## 7 ae_per_vis… S0007 10 8 16 0.375
+## 8 ae_per_vis… S0008 10 7 19 0.714
+## 9 ae_per_vis… S0009 10 9 15 0.333
+## 10 ae_per_vis… S0010 10 7 15 0.714
+## # ℹ 290 more rows
df_sim_sites
## # A tibble: 300 × 10
-## study…¹ site_…² n_pat n_pat…³ visit…⁴ mean_…⁵ mean_…⁶ n_pat…⁷ pval prob_…⁸
-## <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
-## 1 ae_per… S0001 10 8 16 0 0.699 860 0.00834 0.004
-## 2 ae_per… S0002 10 9 16 0.333 0.696 859 0.307 0.132
-## 3 ae_per… S0003 10 10 15 0.3 0.643 898 0.230 0.106
-## 4 ae_per… S0004 10 9 16 0.111 0.698 859 0.0252 0.016
-## 5 ae_per… S0005 10 9 18 0.556 0.791 704 0.570 0.28
-## 6 ae_per… S0006 10 8 17 0 0.742 787 0.00559 0.007
-## 7 ae_per… S0007 10 10 16 0.6 0.693 858 1 0.44
-## 8 ae_per… S0008 10 10 15 0.5 0.640 898 0.840 0.358
-## 9 ae_per… S0009 10 9 17 0.111 0.742 786 0.0178 0.013
-## 10 ae_per… S0010 10 10 16 0.5 0.695 858 0.699 0.304
-## # … with 290 more rows, and abbreviated variable names ¹study_id, ²site_number,
-## # ³n_pat_with_med75, ⁴visit_med75, ⁵mean_ae_site_med75, ⁶mean_ae_study_med75,
-## # ⁷n_pat_with_med75_study, ⁸prob_low
+## study_id site_number n_pat n_pat_with_med75 visit_med75 mean_ae_site_med75
+## <chr> <chr> <int> <int> <dbl> <dbl>
+## 1 ae_per_vis… S0001 10 9 16 0.556
+## 2 ae_per_vis… S0002 10 10 14 0.3
+## 3 ae_per_vis… S0003 10 8 15 0.375
+## 4 ae_per_vis… S0004 10 8 14 0.375
+## 5 ae_per_vis… S0005 10 9 17 0.333
+## 6 ae_per_vis… S0006 10 7 15 0.143
+## 7 ae_per_vis… S0007 10 8 16 0.375
+## 8 ae_per_vis… S0008 10 7 19 0.714
+## 9 ae_per_vis… S0009 10 9 15 0.333
+## 10 ae_per_vis… S0010 10 7 15 0.714
+## # ℹ 290 more rows
+## # ℹ 4 more variables: mean_ae_study_med75 <dbl>, n_pat_with_med75_study <int>,
+## # pval <dbl>, prob_low <dbl>
Simulate Studies
@@ -257,17 +260,9 @@ Simulate Studies
set.seed(1)
-future::plan(multiprocess)
-## Warning: Strategy 'multiprocess' is deprecated in future (>= 1.20.0)
-## [2020-10-30]. Instead, explicitly specify either 'multisession' (recommended) or
-## 'multicore'. In the current R session, 'multiprocess' equals 'multisession'.
-## Warning in supportsMulticoreAndRStudio(...): [ONE-TIME WARNING] Forked
-## processing ('multicore') is not supported when running R from RStudio
-## because it is considered unstable. For more details, how to control forked
-## processing or not, and how to silence this warning in future R sessions,
-## see ?parallelly::supportsMulticore
-
-df_sim_studies <- sim_studies(df_visit, df_site,
+plan(multisession, workers = 6)
+
+df_sim_studies <- sim_studies(df_visit, df_site,
r = 100,
parallel = TRUE,
poisson_test = TRUE,
@@ -276,20 +271,20 @@ Simulate Studies
df_sim_studies
## # A tibble: 30,000 × 8
-## r study_id site_number visit_me…¹ n_pat…² n_pat…³ pval prob_…⁴
-## <dbl> <chr> <chr> <dbl> <int> <dbl> <dbl> <dbl>
-## 1 1 ae_per_visit: 0.05 S0001 16 8 860 0.389 0.208
-## 2 1 ae_per_visit: 0.05 S0002 16 9 859 0.841 0.362
-## 3 1 ae_per_visit: 0.05 S0003 15 10 898 1 1
-## 4 1 ae_per_visit: 0.05 S0004 16 9 859 0.672 0.363
-## 5 1 ae_per_visit: 0.05 S0005 18 9 704 1 1
-## 6 1 ae_per_visit: 0.05 S0006 17 8 787 0.543 0.26
-## 7 1 ae_per_visit: 0.05 S0007 16 10 858 1 0.48
-## 8 1 ae_per_visit: 0.05 S0008 15 10 898 0.105 0.047
-## 9 1 ae_per_visit: 0.05 S0009 17 9 786 1 0.562
-## 10 1 ae_per_visit: 0.05 S0010 16 10 858 1 0.576
-## # … with 29,990 more rows, and abbreviated variable names ¹visit_med75,
-## # ²n_pat_with_med75, ³n_pat_study, ⁴prob_low
+## r study_id site_number visit_med75 n_pat_with_med75 n_pat_study pval
+## <dbl> <chr> <chr> <dbl> <int> <dbl> <dbl>
+## 1 1 ae_per_vis… S0001 16 9 819 1
+## 2 1 ae_per_vis… S0002 14 10 911 1
+## 3 1 ae_per_vis… S0003 15 8 872 0.0511
+## 4 1 ae_per_vis… S0004 14 8 913 0.827
+## 5 1 ae_per_vis… S0005 17 9 750 1
+## 6 1 ae_per_vis… S0006 15 7 873 1
+## 7 1 ae_per_vis… S0007 16 8 820 0.391
+## 8 1 ae_per_vis… S0008 19 7 564 1
+## 9 1 ae_per_vis… S0009 15 9 871 1
+## 10 1 ae_per_vis… S0010 15 7 873 1
+## # ℹ 29,990 more rows
+## # ℹ 1 more variable: prob_low <dbl>
Check p-value Probabilities
@@ -302,28 +297,27 @@ Check p-value Probabilities
+
df_check_pval <- get_ecd_values(df_sim_studies, df_sim_sites, val_str = "pval")
df_check_pval
## # A tibble: 300 × 11
-## study…¹ site_…² n_pat n_pat…³ visit…⁴ mean_…⁵ mean_…⁶ n_pat…⁷ pval prob_…⁸
-## <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
-## 1 ae_per… S0001 10 8 16 0 0.699 860 0.00834 0.004
-## 2 ae_per… S0002 10 9 16 0.333 0.696 859 0.307 0.132
-## 3 ae_per… S0003 10 10 15 0.3 0.643 898 0.230 0.106
-## 4 ae_per… S0004 10 9 16 0.111 0.698 859 0.0252 0.016
-## 5 ae_per… S0005 10 9 18 0.556 0.791 704 0.570 0.28
-## 6 ae_per… S0006 10 8 17 0 0.742 787 0.00559 0.007
-## 7 ae_per… S0007 10 10 16 0.6 0.693 858 1 0.44
-## 8 ae_per… S0008 10 10 15 0.5 0.640 898 0.840 0.358
-## 9 ae_per… S0009 10 9 17 0.111 0.742 786 0.0178 0.013
-## 10 ae_per… S0010 10 10 16 0.5 0.695 858 0.699 0.304
-## # … with 290 more rows, 1 more variable: pval_ecd <dbl>, and abbreviated
-## # variable names ¹study_id, ²site_number, ³n_pat_with_med75, ⁴visit_med75,
-## # ⁵mean_ae_site_med75, ⁶mean_ae_study_med75, ⁷n_pat_with_med75_study,
-## # ⁸prob_low
-
+## study_id site_number n_pat n_pat_with_med75 visit_med75 mean_ae_site_med75
+## <chr> <chr> <int> <int> <dbl> <dbl>
+## 1 ae_per_vis… S0001 10 9 16 0.556
+## 2 ae_per_vis… S0002 10 10 14 0.3
+## 3 ae_per_vis… S0003 10 8 15 0.375
+## 4 ae_per_vis… S0004 10 8 14 0.375
+## 5 ae_per_vis… S0005 10 9 17 0.333
+## 6 ae_per_vis… S0006 10 7 15 0.143
+## 7 ae_per_vis… S0007 10 8 16 0.375
+## 8 ae_per_vis… S0008 10 7 19 0.714
+## 9 ae_per_vis… S0009 10 9 15 0.333
+## 10 ae_per_vis… S0010 10 7 15 0.714
+## # ℹ 290 more rows
+## # ℹ 5 more variables: mean_ae_study_med75 <dbl>, n_pat_with_med75_study <int>,
+## # pval <dbl>, prob_low <dbl>, pval_ecd <dbl>
+
df_check_pval %>%
ggplot(aes(log(pval, base = 10), log(pval_ecd, base = 10))) +
geom_point(alpha = 0.5, size = 2) +
@@ -342,28 +336,27 @@ Perform Same Check on
Note: Dots on the edge of the graph that are cut off have zero value
for the corresponding axis. For the tie simulations using 1000 repeats
the smallest value greater zero we get for prob_low is 0.001 (1e-3).
-
+
df_check_prob <- get_ecd_values(df_sim_studies, df_sim_sites, val_str = "prob_low")
df_check_prob
## # A tibble: 300 × 11
-## study…¹ site_…² n_pat n_pat…³ visit…⁴ mean_…⁵ mean_…⁶ n_pat…⁷ pval prob_…⁸
-## <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
-## 1 ae_per… S0001 10 8 16 0 0.699 860 0.00834 0.004
-## 2 ae_per… S0002 10 9 16 0.333 0.696 859 0.307 0.132
-## 3 ae_per… S0003 10 10 15 0.3 0.643 898 0.230 0.106
-## 4 ae_per… S0004 10 9 16 0.111 0.698 859 0.0252 0.016
-## 5 ae_per… S0005 10 9 18 0.556 0.791 704 0.570 0.28
-## 6 ae_per… S0006 10 8 17 0 0.742 787 0.00559 0.007
-## 7 ae_per… S0007 10 10 16 0.6 0.693 858 1 0.44
-## 8 ae_per… S0008 10 10 15 0.5 0.640 898 0.840 0.358
-## 9 ae_per… S0009 10 9 17 0.111 0.742 786 0.0178 0.013
-## 10 ae_per… S0010 10 10 16 0.5 0.695 858 0.699 0.304
-## # … with 290 more rows, 1 more variable: prob_low_ecd <dbl>, and abbreviated
-## # variable names ¹study_id, ²site_number, ³n_pat_with_med75, ⁴visit_med75,
-## # ⁵mean_ae_site_med75, ⁶mean_ae_study_med75, ⁷n_pat_with_med75_study,
-## # ⁸prob_low
-
+## study_id site_number n_pat n_pat_with_med75 visit_med75 mean_ae_site_med75
+## <chr> <chr> <int> <int> <dbl> <dbl>
+## 1 ae_per_vis… S0001 10 9 16 0.556
+## 2 ae_per_vis… S0002 10 10 14 0.3
+## 3 ae_per_vis… S0003 10 8 15 0.375
+## 4 ae_per_vis… S0004 10 8 14 0.375
+## 5 ae_per_vis… S0005 10 9 17 0.333
+## 6 ae_per_vis… S0006 10 7 15 0.143
+## 7 ae_per_vis… S0007 10 8 16 0.375
+## 8 ae_per_vis… S0008 10 7 19 0.714
+## 9 ae_per_vis… S0009 10 9 15 0.333
+## 10 ae_per_vis… S0010 10 7 15 0.714
+## # ℹ 290 more rows
+## # ℹ 5 more variables: mean_ae_study_med75 <dbl>, n_pat_with_med75_study <int>,
+## # pval <dbl>, prob_low <dbl>, prob_low_ecd <dbl>
+
df_check_prob %>%
ggplot(aes(log(prob_low, base = 10), log(prob_low_ecd, base = 10))) +
geom_point(alpha = 0.5, size = 2) +
@@ -387,6 +380,8 @@ Conclusion
+plan(sequential)
@@ -408,7 +403,7 @@ Conclusion
-Site built with pkgdown 2.0.5.
+Site built with pkgdown 2.0.7.
diff --git a/docs/articles/check_poisson_files/figure-html/unnamed-chunk-4-1.png b/docs/articles/check_poisson_files/figure-html/unnamed-chunk-4-1.png
index 4eb24a0..e25f87d 100644
Binary files a/docs/articles/check_poisson_files/figure-html/unnamed-chunk-4-1.png and b/docs/articles/check_poisson_files/figure-html/unnamed-chunk-4-1.png differ
diff --git a/docs/articles/check_poisson_files/figure-html/unnamed-chunk-5-1.png b/docs/articles/check_poisson_files/figure-html/unnamed-chunk-5-1.png
index deeab82..8a3a9bb 100644
Binary files a/docs/articles/check_poisson_files/figure-html/unnamed-chunk-5-1.png and b/docs/articles/check_poisson_files/figure-html/unnamed-chunk-5-1.png differ
diff --git a/docs/articles/gsm_perf.html b/docs/articles/gsm_perf.html
new file mode 100644
index 0000000..f8ae437
--- /dev/null
+++ b/docs/articles/gsm_perf.html
@@ -0,0 +1,1643 @@
+
+
+
+
+
+
+
+Comparing {simaerep} and {gsm} Performance • simaerep
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Comparing {simaerep} and {gsm} Performance
+
+
+ Source: vignettes/gsm_perf.Rmd
+ gsm_perf.Rmd
+
+
+
+
+
+
+Load
+
+
+suppressPackageStartupMessages(library(tidyverse))
+suppressPackageStartupMessages(library(knitr))
+suppressPackageStartupMessages(library(furrr))
+suppressPackageStartupMessages(library(future))
+suppressPackageStartupMessages(library(simaerep))
+
+plan(multisession, workers = 6)
+
+
+Install {gsm}
+
+
+devtools::install_github("Gilead-BioStats/gsm", ref = "main")
+
+
+Introduction
+
+The {gsm} R
+package provides a standardized Risk Based Quality Monitoring (RBQM)
+framework for clinical trials that pairs a flexible data pipeline with
+robust reports. It also uses Funnel Plots to flag outliers which provide
+broader tolerance limits for sites with low exposure and narrower limits
+for sites with higher exposure. This method is different to the event
+rate based limits we have used in previous heuristics to measure
+{simaerep} performance. Funnel plots are discussed in greater detail by
+Zink et
+al. 2018
+One of the draw backs of using funnel plots for flagging is that they
+assume that the AE rate remains constant over the course of the
+study.
+
+
+Prepare Data
+
+
+Load Portfolio Configurations
+
+We have prepared a snapshot of the AE reporting configuration of our
+current portfolio. For each study we have also measured a visit-specific
+AE rate which allows us to generate a synthetic portfolio with flexible
+AE rates across a study.
+
+df_config <- readr::read_csv("ae_conf_20240220.csv")
+df_ae_rates <- readr::read_csv("ae_rates_20240220.csv")
+
+df_config %>%
+ head(25) %>%
+ knitr::kable()
+
+
+
+
+study_id
+ae_per_visit_mean
+site_number
+max_visit_sd
+max_visit_mean
+n_pat
+
+
+
+0001
+0.2973806
+4746
+3.5355339
+16.50000
+2
+
+
+0001
+0.2973806
+4747
+0.7071068
+27.50000
+2
+
+
+0001
+0.2973806
+4750
+0.5000000
+19.25000
+4
+
+
+0001
+0.2973806
+4815
+8.3765546
+22.16667
+6
+
+
+0001
+0.2973806
+4816
+0.0000000
+27.00000
+1
+
+
+0001
+0.2973806
+4817
+0.0000000
+31.00000
+1
+
+
+0001
+0.2973806
+4818
+0.0000000
+30.00000
+1
+
+
+0001
+0.2973806
+4891
+2.3804761
+22.00000
+7
+
+
+0001
+0.2973806
+4893
+0.0000000
+20.00000
+2
+
+
+0001
+0.2973806
+4932
+1.0954451
+26.80000
+5
+
+
+0001
+0.2973806
+4941
+1.7320508
+27.00000
+3
+
+
+0001
+0.2973806
+4942
+1.4142136
+25.00000
+2
+
+
+0001
+0.2973806
+4943
+0.0000000
+25.00000
+1
+
+
+0001
+0.2973806
+4966
+0.0000000
+18.00000
+1
+
+
+0001
+0.2973806
+4967
+0.8164966
+18.00000
+4
+
+
+0001
+0.2973806
+4968
+0.0000000
+24.00000
+1
+
+
+0001
+0.2973806
+4969
+1.5275252
+21.33333
+3
+
+
+0001
+0.2973806
+4984
+1.7078251
+30.75000
+4
+
+
+0001
+0.2973806
+4985
+8.9087972
+21.16667
+6
+
+
+0001
+0.2973806
+4986
+0.5773503
+20.33333
+3
+
+
+0001
+0.2973806
+4988
+1.4142136
+19.00000
+2
+
+
+0001
+0.2973806
+5079
+2.1602469
+23.00000
+4
+
+
+0001
+0.2973806
+5080
+0.0000000
+23.00000
+1
+
+
+0001
+0.2973806
+5081
+0.0000000
+23.00000
+1
+
+
+0001
+0.2973806
+5082
+1.4142136
+22.00000
+2
+
+
+
+
+
+
+study_id
+cum_visit
+ae_rate
+n_pat
+
+
+
+0001
+1
+0.017699
+113
+
+
+0001
+2
+0.088496
+113
+
+
+0001
+3
+0.230088
+113
+
+
+0001
+4
+0.223214
+112
+
+
+0001
+5
+0.187500
+112
+
+
+0001
+6
+0.196429
+112
+
+
+0001
+7
+0.396396
+111
+
+
+0001
+8
+0.207207
+111
+
+
+0001
+9
+0.216216
+111
+
+
+0001
+10
+0.315315
+111
+
+
+0001
+11
+0.263636
+110
+
+
+0001
+12
+0.318182
+110
+
+
+0001
+13
+0.345455
+110
+
+
+0001
+14
+0.379630
+108
+
+
+0001
+15
+0.401869
+107
+
+
+0001
+16
+0.485981
+107
+
+
+0001
+17
+0.383178
+107
+
+
+0001
+18
+0.396226
+106
+
+
+0001
+19
+0.370000
+100
+
+
+0001
+20
+0.439560
+91
+
+
+0001
+21
+0.369863
+73
+
+
+0001
+22
+0.370968
+62
+
+
+0001
+23
+0.345455
+55
+
+
+0001
+24
+0.446809
+47
+
+
+0001
+25
+0.292683
+41
+
+
+
+
+
+Simulate Portfolio
+
+We generate two synthetic portfolios with no AE under-reporting
+sites. One portfolio with a fixed AE rate for all visits and another one
+with a flexible visit-specific AE rate.
+
+df_portf_fix <- sim_test_data_portfolio(df_config, parallel = TRUE, progress = TRUE)
+
+df_portf_fix %>%
+ head(25) %>%
+ knitr::kable()
+
+
+
+
+study_id
+ae_per_visit_mean
+site_number
+max_visit_sd
+max_visit_mean
+patnum
+visit
+n_ae
+
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+1
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+2
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+3
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+4
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+5
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+6
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+7
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+8
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+9
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+10
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+11
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+1
+1
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+2
+1
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+3
+1
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+4
+2
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+5
+2
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+6
+3
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+7
+4
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+8
+4
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+9
+5
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+10
+5
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+11
+6
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+12
+6
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+13
+6
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+14
+6
+
+
+
+
+df_portf_flex <- sim_test_data_portfolio(df_config, df_ae_rates = df_ae_rates, parallel = TRUE, progress = TRUE)
+
+
+Compare AE rates
+
+Next we confirm the different AE rates in our two synthetic
+portfolios.
+
+df_rate_fix <- df_portf_fix %>%
+ mutate(ae_rate = coalesce(n_ae - lag(n_ae), n_ae), .by = c("study_id", "patnum")) %>%
+ summarise(ae_rate = mean(ae_rate), .by = c("study_id", "visit")) %>%
+ mutate(rate = "fix")
+
+
+df_rate_flex <- df_portf_flex %>%
+ mutate(ae_rate = coalesce(n_ae - lag(n_ae), n_ae), .by = c("study_id", "patnum")) %>%
+ summarise(ae_rate = mean(ae_rate), .by = c("study_id", "visit")) %>%
+ mutate(rate = "flex")
+
+bind_rows(df_rate_flex, df_rate_fix) %>%
+ ggplot(aes(visit, ae_rate)) +
+ geom_line(aes(group = study_id), alpha = 0.2) +
+ geom_smooth() +
+ facet_wrap(~ rate) +
+ labs(title = "Average AE rates per Study")
+
+
+bind_rows(df_rate_flex, df_rate_fix) %>%
+ filter(dense_rank(study_id) <= 16) %>%
+ ggplot(aes(visit, ae_rate)) +
+ geom_line(aes(group = rate, color = rate)) +
+ facet_wrap(~ study_id, scales = "free") +
+ labs(title = "Average AE rates for Selected Studies")
+
+We can confirm that the AE rates in the “flexible” portfolio are not
+constant. Moreover we see that the AE rate profile is very unique for
+each study.
+
+
+
+Apply {gsm}
+
+
+Example
+
+Here we demonstrate how to use the {gsm} package on our simulated
+portfolios.
+
+get_SUBJ <- function(df_portf) {
+ df_portf %>%
+ select(study_id, siteid = site_number, subjid = patnum, timeonstudy = visit) %>%
+ summarise(timeonstudy = max(timeonstudy), .by = c(study_id, siteid, subjid)) %>%
+ group_by(study_id) %>%
+ nest()
+}
+
+
+get_AE <- function(df_portf) {
+ df_portf_fix %>%
+ select(study_id, subjid = patnum, n_ae) %>%
+ summarise(n_ae = max(n_ae), .by = c(study_id, subjid)) %>%
+ filter(n_ae > 0) %>%
+ mutate(n_ae = map(n_ae, ~ tibble(n = seq(1, .)), .progress = TRUE)) %>%
+ unnest(n_ae) %>%
+ select(- n) %>%
+ group_by(study_id) %>%
+ nest()
+}
+
+dfSUBJ_fix <- get_SUBJ(df_portf_fix)
+dfAE_fix <- get_AE(df_portf_fix)
+
+dfInput <- gsm::AE_Map_Raw(list(dfSUBJ = dfSUBJ_fix$data[[1]], dfAE = dfAE_fix$data[[1]]))
+dfInput
+## # A tibble: 113 × 5
+## SubjectID SiteID Exposure Count Rate
+## <chr> <chr> <int> <int> <dbl>
+## 1 0001 4746 11 0 0
+## 2 0002 4746 17 6 0.353
+## 3 0003 4747 27 6 0.222
+## 4 0004 4747 27 5 0.185
+## 5 0005 4750 18 10 0.556
+## 6 0006 4750 19 8 0.421
+## 7 0007 4750 19 8 0.421
+## 8 0008 4750 18 7 0.389
+## 9 0009 4815 15 3 0.2
+## 10 0010 4815 31 2 0.0645
+## # ℹ 103 more rows
+
+dfTransformed <- gsm::Transform_Rate(
+ dfInput,
+ strNumeratorCol = "Count",
+ strDenominatorCol = "Exposure"
+ )
+dfTransformed
+## # A tibble: 44 × 4
+## GroupID Numerator Denominator Metric
+## <chr> <int> <int> <dbl>
+## 1 4746 6 28 0.214
+## 2 4747 11 54 0.204
+## 3 4750 33 74 0.446
+## 4 4815 33 145 0.228
+## 5 4816 8 27 0.296
+## 6 4817 5 31 0.161
+## 7 4818 5 30 0.167
+## 8 4891 45 144 0.312
+## 9 4893 13 40 0.325
+## 10 4932 36 132 0.273
+## # ℹ 34 more rows
+
+dfAnalyzed <- gsm::Analyze_NormalApprox(dfTransformed)
+dfAnalyzed
+## # A tibble: 44 × 7
+## GroupID Numerator Denominator Metric OverallMetric Factor Score
+## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
+## 1 4815 33 145 0.228 0.300 1.31 -1.65
+## 2 4941 17 80 0.212 0.300 1.31 -1.48
+## 3 4817 5 31 0.161 0.300 1.31 -1.47
+## 4 4818 5 30 0.167 0.300 1.31 -1.39
+## 5 5079 20 89 0.225 0.300 1.31 -1.35
+## 6 4747 11 54 0.204 0.300 1.31 -1.34
+## 7 5311 3 20 0.15 0.300 1.31 -1.27
+## 8 5084 4 21 0.190 0.300 1.31 -0.953
+## 9 4943 5 25 0.2 0.300 1.31 -0.949
+## 10 4746 6 28 0.214 0.300 1.31 -0.860
+## # ℹ 34 more rows
+
+dfFlagged <- gsm::Flag_NormalApprox(dfAnalyzed, vThreshold = c(-3, -2, 2, 3))
+dfFlagged
+## # A tibble: 44 × 8
+## GroupID Numerator Denominator Metric OverallMetric Factor Score Flag
+## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
+## 1 5168 9 15 0.6 0.300 1.31 2.22 1
+## 2 5194 76 197 0.386 0.300 1.31 2.31 1
+## 3 4750 33 74 0.446 0.300 1.31 2.40 1
+## 4 4815 33 145 0.228 0.300 1.31 -1.65 0
+## 5 4941 17 80 0.212 0.300 1.31 -1.48 0
+## 6 4817 5 31 0.161 0.300 1.31 -1.47 0
+## 7 4818 5 30 0.167 0.300 1.31 -1.39 0
+## 8 5079 20 89 0.225 0.300 1.31 -1.35 0
+## 9 4747 11 54 0.204 0.300 1.31 -1.34 0
+## 10 5311 3 20 0.15 0.300 1.31 -1.27 0
+## # ℹ 34 more rows
+
+dfSummary <- gsm::Summarize(dfFlagged)
+dfSummary
+## # A tibble: 44 × 6
+## GroupID Numerator Denominator Metric Score Flag
+## <chr> <int> <int> <dbl> <dbl> <dbl>
+## 1 5168 9 15 0.6 2.22 1
+## 2 4750 33 74 0.446 2.40 1
+## 3 5194 76 197 0.386 2.31 1
+## 4 5083 18 44 0.409 1.39 0
+## 5 5265 16 42 0.381 1.01 0
+## 6 5130 52 140 0.371 1.62 0
+## 7 4986 21 58 0.362 0.908 0
+## 8 5131 8 23 0.348 0.442 0
+## 9 4985 38 111 0.342 0.860 0
+## 10 5224 8 24 0.333 0.316 0
+## # ℹ 34 more rows
+
+dfBounds <- gsm::Analyze_NormalApprox_PredictBounds(dfTransformed, vThreshold = c(-3, -2, 2, 3))
+dfBounds
+## # A tibble: 1,246 × 5
+## Threshold Denominator LogDenominator Numerator Metric
+## <dbl> <dbl> <dbl> <dbl> <dbl>
+## 1 -3 28.1 3.34 0.0765 0.00272
+## 2 -3 28.8 3.36 0.187 0.00649
+## 3 -3 29.6 3.39 0.299 0.0101
+## 4 -3 30.3 3.41 0.413 0.0136
+## 5 -3 31.0 3.43 0.527 0.0170
+## 6 -3 31.7 3.46 0.643 0.0203
+## 7 -3 32.5 3.48 0.760 0.0234
+## 8 -3 33.2 3.50 0.878 0.0264
+## 9 -3 33.9 3.52 0.997 0.0294
+## 10 -3 34.7 3.55 1.12 0.0322
+## # ℹ 1,236 more rows
+
+chart <- gsm::Visualize_Scatter(dfFlagged, dfBounds)
+chart
+
+
+
+Simulate UR
+
+We write a function that removes a given ratio of AEs from one site
+in the data set and returns its z-score.
+
+sim_site_ur_gsm <- function(site, ur_rate, dfTransformed) {
+ dfTransformed <- dfTransformed %>%
+ mutate(
+ Numerator = ifelse(GroupID == site, Numerator * (1 - ur_rate), Numerator),
+ Metric = Numerator / Denominator
+ )
+
+ gsm::Analyze_NormalApprox(dfTransformed) %>%
+ filter(GroupID == site) %>%
+ pull(Score)
+}
+
+sim_site_ur_gsm("4747", ur_rate = 0.75, dfTransformed)
+## [1] -3.095796
+We write another function that systematically applies this
+sim_site_ur_gsm across all sites in all studies across a
+range of under-reporting ratios.
+
+sim_ur_gsm <- function(dfSUBJ, dfAE) {
+ dfSUBJ %>%
+ inner_join(dfAE, by = "study_id") %>%
+ ungroup() %>%
+ mutate(
+ trans = map2(data.x, data.y, ~ gsm::AE_Map_Raw(list(dfSUBJ = .x, dfAE = .y))),
+ trans = map(trans, ~ gsm::Transform_Rate(., strNumeratorCol = "Count", strDenominatorCol = "Exposure")),
+ sites = map(data.x, ~ distinct(., siteid))
+ ) %>%
+ select(- starts_with("data.")) %>%
+ unnest(sites) %>%
+ mutate(ur = list(tibble(ur_rate = c(0, 0.1, 0.25, 0.5, 0.75, 1)))) %>%
+ unnest(ur) %>%
+ mutate(
+ score = pmap_dbl(list(siteid, ur_rate, trans), sim_site_ur_gsm, .progress = TRUE)
+ )
+}
+
+df_sim_gsm_fix <- sim_ur_gsm(dfSUBJ_fix, dfAE_fix)
+## Warning: There were 3116 warnings in `mutate()`.
+## The first warning was:
+## ℹ In argument: `z_i = ifelse(...)`.
+## Caused by warning:
+## ! There was 1 warning in `mutate()`.
+## ℹ In argument: `z_0 = ifelse(...)`.
+## Caused by warning in `sqrt()`:
+## ! NaNs produced
+## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3115 remaining warnings.
+
+df_sim_gsm_fix
+## # A tibble: 129,966 × 5
+## study_id trans siteid ur_rate score
+## <chr> <list> <chr> <dbl> <dbl>
+## 1 0001 <tibble [44 × 4]> 4746 0 -0.860
+## 2 0001 <tibble [44 × 4]> 4746 0.1 -1.07
+## 3 0001 <tibble [44 × 4]> 4746 0.25 -1.38
+## 4 0001 <tibble [44 × 4]> 4746 0.5 -1.87
+## 5 0001 <tibble [44 × 4]> 4746 0.75 -2.33
+## 6 0001 <tibble [44 × 4]> 4746 1 -2.75
+## 7 0001 <tibble [44 × 4]> 4747 0 -1.34
+## 8 0001 <tibble [44 × 4]> 4747 0.1 -1.61
+## 9 0001 <tibble [44 × 4]> 4747 0.25 -1.99
+## 10 0001 <tibble [44 × 4]> 4747 0.5 -2.57
+## # ℹ 129,956 more rows
+We repeat the same steps for the portfolio with the flexible AE
+rates.
+
+dfSUBJ_flex <- get_SUBJ(df_portf_flex)
+dfAE_flex<- get_AE(df_portf_flex)
+
+
+df_sim_gsm_flex <- sim_ur_gsm(dfSUBJ_flex, dfAE_flex)
+## Warning: There were 3085 warnings in `mutate()`.
+## The first warning was:
+## ℹ In argument: `z_i = ifelse(...)`.
+## Caused by warning:
+## ! There was 1 warning in `mutate()`.
+## ℹ In argument: `z_0 = ifelse(...)`.
+## Caused by warning in `sqrt()`:
+## ! NaNs produced
+## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3084 remaining warnings.
+
+
+
+UR {simaerep}
+
+We simulate under-reporting for both portfolios using {simaerep}
+
+df_sim_simaerep_fix <- sim_ur_scenarios(
+ df_portf_fix,
+ extra_ur_sites = 0,
+ ur_rate = c(0, 0.1, 0.25, 0.5, 0.75, 1),
+ parallel = TRUE,
+ poisson = TRUE,
+ prob_lower = TRUE,
+ progress = TRUE
+)
+
+df_sim_simaerep_flex <- sim_ur_scenarios(
+ df_portf_flex,
+ extra_ur_sites = 0,
+ ur_rate = c(0, 0.1, 0.25, 0.5, 0.75, 1),
+ parallel = TRUE,
+ poisson = TRUE,
+ prob_lower = TRUE,
+ progress = TRUE
+)
+
+
+Evaluate
+
+
+Combine Results
+
+
+df_sim_gsm_fix$ae_rate <- "AE rate: fix"
+df_sim_gsm_flex$ae_rate <- "AE rate: flexible"
+df_sim_simaerep_fix$ae_rate <- "AE rate: fix"
+df_sim_simaerep_flex$ae_rate <- "AE rate: flexible"
+
+df_sim_thresh2 <- bind_rows(df_sim_gsm_fix, df_sim_gsm_flex) %>%
+ mutate(
+ is_ur = score <= -2,
+ type = "{gsm} - thresh: -2",
+ site_number = siteid
+ ) %>%
+ select(type, ae_rate, study_id, site_number, ur_rate, is_ur, score)
+
+
+df_sim_simaerep_threshp95 <- bind_rows(df_sim_simaerep_fix, df_sim_simaerep_flex) %>%
+ mutate(
+ is_ur = prob_low_prob_ur >= 0.95,
+ type = "{simaerep} - thresh: 0.95"
+ ) %>%
+ select(type, ae_rate, study_id, site_number, ur_rate, is_ur, score = prob_low_prob_ur)
+
+
+df_eval <- bind_rows(
+ df_sim_thresh2,
+ df_sim_simaerep_threshp95,
+)
+
+
+Aggregate
+
+
+get_prop_test_ci95 <- function(..., ix) {
+
+ stopifnot(ix %in% c(1, 2))
+
+ tryCatch(
+ prop.test(...)$conf.int[ix],
+ error = function(cnd) c(NA, NA)[ix]
+ )
+}
+
+aggr_results <- function(df_eval) {
+
+df_perf <- df_eval %>%
+ summarise(
+ n = n(),
+ .by = c(type, ae_rate, ur_rate, is_ur)
+ ) %>%
+ pivot_wider(
+ names_from = is_ur,
+ values_from = n,
+ names_prefix = "is_ur_",
+ values_fill = 0
+ ) %>%
+ mutate(
+ n_sites = is_ur_TRUE + is_ur_FALSE + is_ur_NA,
+ ratio = is_ur_TRUE / n_sites,
+ ratio_type = ifelse(ur_rate == 0, "fpr", "tpr"),
+ ci95_low = map2_dbl(is_ur_TRUE, n_sites, ~ get_prop_test_ci95(.x, .y, ix = 1)),
+ ci95_high = map2_dbl(is_ur_TRUE, n_sites, ~ get_prop_test_ci95(.x, .y, ix = 2))
+ )
+}
+
+df_perf <- aggr_results(df_eval)
+
+
+
+Results
+
+
+Table
+
+
+
+
+
+
+type
+ae_rate
+ur_rate
+is_ur_FALSE
+is_ur_TRUE
+is_ur_NA
+n_sites
+ratio
+ratio_type
+ci95_low
+ci95_high
+
+
+
+{gsm} - thresh: -2
+AE rate: fix
+0.00
+20830
+299
+532
+21661
+0.014
+fpr
+0.012
+0.015
+
+
+{gsm} - thresh: -2
+AE rate: fix
+0.10
+19851
+1278
+532
+21661
+0.059
+tpr
+0.056
+0.062
+
+
+{gsm} - thresh: -2
+AE rate: fix
+0.25
+15848
+5282
+531
+21661
+0.244
+tpr
+0.238
+0.250
+
+
+{gsm} - thresh: -2
+AE rate: fix
+0.50
+9154
+11987
+520
+21661
+0.553
+tpr
+0.547
+0.560
+
+
+{gsm} - thresh: -2
+AE rate: fix
+0.75
+5214
+15941
+506
+21661
+0.736
+tpr
+0.730
+0.742
+
+
+{gsm} - thresh: -2
+AE rate: fix
+1.00
+2915
+18251
+495
+21661
+0.843
+tpr
+0.838
+0.847
+
+
+{gsm} - thresh: -2
+AE rate: flexible
+0.00
+20899
+241
+521
+21661
+0.011
+fpr
+0.010
+0.013
+
+
+{gsm} - thresh: -2
+AE rate: flexible
+0.10
+20669
+471
+521
+21661
+0.022
+tpr
+0.020
+0.024
+
+
+{gsm} - thresh: -2
+AE rate: flexible
+0.25
+19742
+1399
+520
+21661
+0.065
+tpr
+0.061
+0.068
+
+
+{gsm} - thresh: -2
+AE rate: flexible
+0.50
+16448
+4698
+515
+21661
+0.217
+tpr
+0.211
+0.222
+
+
+{gsm} - thresh: -2
+AE rate: flexible
+0.75
+12320
+8830
+511
+21661
+0.408
+tpr
+0.401
+0.414
+
+
+{gsm} - thresh: -2
+AE rate: flexible
+1.00
+8727
+12437
+497
+21661
+0.574
+tpr
+0.568
+0.581
+
+
+{simaerep} - thresh: 0.95
+AE rate: fix
+0.00
+21367
+294
+0
+21661
+0.014
+fpr
+0.012
+0.015
+
+
+{simaerep} - thresh: 0.95
+AE rate: fix
+0.10
+20366
+1295
+0
+21661
+0.060
+tpr
+0.057
+0.063
+
+
+{simaerep} - thresh: 0.95
+AE rate: fix
+0.25
+17155
+4506
+0
+21661
+0.208
+tpr
+0.203
+0.214
+
+
+{simaerep} - thresh: 0.95
+AE rate: fix
+0.50
+10635
+11026
+0
+21661
+0.509
+tpr
+0.502
+0.516
+
+
+{simaerep} - thresh: 0.95
+AE rate: fix
+0.75
+5845
+15816
+0
+21661
+0.730
+tpr
+0.724
+0.736
+
+
+{simaerep} - thresh: 0.95
+AE rate: fix
+1.00
+4053
+17608
+0
+21661
+0.813
+tpr
+0.808
+0.818
+
+
+{simaerep} - thresh: 0.95
+AE rate: flexible
+0.00
+21340
+321
+0
+21661
+0.015
+fpr
+0.013
+0.017
+
+
+{simaerep} - thresh: 0.95
+AE rate: flexible
+0.10
+20234
+1427
+0
+21661
+0.066
+tpr
+0.063
+0.069
+
+
+{simaerep} - thresh: 0.95
+AE rate: flexible
+0.25
+16661
+5000
+0
+21661
+0.231
+tpr
+0.225
+0.237
+
+
+{simaerep} - thresh: 0.95
+AE rate: flexible
+0.50
+9878
+11783
+0
+21661
+0.544
+tpr
+0.537
+0.551
+
+
+{simaerep} - thresh: 0.95
+AE rate: flexible
+0.75
+5216
+16445
+0
+21661
+0.759
+tpr
+0.753
+0.765
+
+
+{simaerep} - thresh: 0.95
+AE rate: flexible
+1.00
+3595
+18066
+0
+21661
+0.834
+tpr
+0.829
+0.839
+
+
+
+
+
+Plot Performance Metrics
+
+
+- {gsm} has better performance than {simaerep} when the AE rate is
+fixed, while {simaerep} greatly outperforms {gsm} when the AE rate is
+flexible and mimics the AE rates encountered in real study data
+sets.
+
+
+plot_perf <- function(df_perf) {
+
+ df_perf %>%
+ mutate(ur_rate = paste0("under-reporting rate: ", ur_rate, " - ", ratio_type),
+ ur_rate = ifelse(str_detect(ur_rate, "fpr"), "fpr", ur_rate)) %>%
+ group_by(ur_rate) %>%
+ ggplot(aes(type, ratio)) +
+ geom_errorbar(aes(ymin = ci95_low, ymax = ci95_high, color = type), linewidth = 1) +
+ facet_grid(ur_rate ~ ae_rate) +
+ coord_flip() +
+ theme(legend.position = "bottom") +
+ labs(
+ x = "",
+ y = "CI95 Performance Ratio",
+ title = "{simaerep} vs {gsm} Performance"
+ ) +
+ scale_color_manual(values = c("#5491CC", "#F46626"))
+}
+
+plot_perf(df_perf)
+
+
+plan(sequential)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/articles/gsm_perf_files/figure-html/plot-1.png b/docs/articles/gsm_perf_files/figure-html/plot-1.png
new file mode 100644
index 0000000..2feedb9
Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/plot-1.png differ
diff --git a/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-10-1.png b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-10-1.png
new file mode 100644
index 0000000..5098a97
Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-10-1.png differ
diff --git a/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-3-1.png b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-3-1.png
new file mode 100644
index 0000000..91a7a28
Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-3-1.png differ
diff --git a/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-4-1.png b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-4-1.png
new file mode 100644
index 0000000..e27fc8a
Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-4-1.png differ
diff --git a/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-6-1.png b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-6-1.png
new file mode 100644
index 0000000..524a0f5
Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-6-1.png differ
diff --git a/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-9-1.png b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-9-1.png
new file mode 100644
index 0000000..7289510
Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-9-1.png differ
diff --git a/docs/articles/index.html b/docs/articles/index.html
index 8b6151d..d1f3485 100644
--- a/docs/articles/index.html
+++ b/docs/articles/index.html
@@ -23,7 +23,7 @@
@@ -67,6 +67,9 @@
Aggregate AEs by Days or Visit?
+
+ Comparing {simaerep} and {gsm} Performance
+
## # A tibble: 59,221 × 9
-## patnum site_number is_ur max_visit_mean max_vi…¹ ae_pe…² visit n_ae study…³
-## <chr> <chr> <lgl> <dbl> <dbl> <dbl> <int> <int> <chr>
-## 1 P000001 S0001 TRUE 20 4 0.25 1 0 ae_per…
-## 2 P000001 S0001 TRUE 20 4 0.25 2 1 ae_per…
-## 3 P000001 S0001 TRUE 20 4 0.25 3 2 ae_per…
-## 4 P000001 S0001 TRUE 20 4 0.25 4 2 ae_per…
-## 5 P000001 S0001 TRUE 20 4 0.25 5 2 ae_per…
-## 6 P000001 S0001 TRUE 20 4 0.25 6 2 ae_per…
-## 7 P000001 S0001 TRUE 20 4 0.25 7 2 ae_per…
-## 8 P000001 S0001 TRUE 20 4 0.25 8 4 ae_per…
-## 9 P000001 S0001 TRUE 20 4 0.25 9 4 ae_per…
-## 10 P000001 S0001 TRUE 20 4 0.25 10 4 ae_per…
-## # … with 59,211 more rows, and abbreviated variable names ¹max_visit_sd,
-## # ²ae_per_visit_mean, ³study_id## # A tibble: 58,367 × 9
+## patnum site_number is_ur max_visit_mean max_visit_sd ae_per_visit_mean visit
+## <chr> <chr> <lgl> <dbl> <dbl> <dbl> <int>
+## 1 P000001 S0001 TRUE 20 4 0.25 1
+## 2 P000001 S0001 TRUE 20 4 0.25 2
+## 3 P000001 S0001 TRUE 20 4 0.25 3
+## 4 P000001 S0001 TRUE 20 4 0.25 4
+## 5 P000001 S0001 TRUE 20 4 0.25 5
+## 6 P000001 S0001 TRUE 20 4 0.25 6
+## 7 P000001 S0001 TRUE 20 4 0.25 7
+## 8 P000001 S0001 TRUE 20 4 0.25 8
+## 9 P000001 S0001 TRUE 20 4 0.25 9
+## 10 P000001 S0001 TRUE 20 4 0.25 10
+## # ℹ 58,357 more rows
+## # ℹ 2 more variables: n_ae <int>, study_id <chr>
df_site## # A tibble: 300 × 6
-## study_id site_number n_pat n_pat_with_med75 visit_med75 mean_ae_s…¹
-## <chr> <chr> <int> <int> <dbl> <dbl>
-## 1 ae_per_visit: 0.05 S0001 10 8 16 0
-## 2 ae_per_visit: 0.05 S0002 10 9 16 0.333
-## 3 ae_per_visit: 0.05 S0003 10 10 15 0.3
-## 4 ae_per_visit: 0.05 S0004 10 9 16 0.111
-## 5 ae_per_visit: 0.05 S0005 10 9 18 0.556
-## 6 ae_per_visit: 0.05 S0006 10 8 17 0
-## 7 ae_per_visit: 0.05 S0007 10 10 16 0.6
-## 8 ae_per_visit: 0.05 S0008 10 10 15 0.5
-## 9 ae_per_visit: 0.05 S0009 10 9 17 0.111
-## 10 ae_per_visit: 0.05 S0010 10 10 16 0.5
-## # … with 290 more rows, and abbreviated variable name ¹mean_ae_site_med75
df_sim_sites## # A tibble: 300 × 10
-## study…¹ site_…² n_pat n_pat…³ visit…⁴ mean_…⁵ mean_…⁶ n_pat…⁷ pval prob_…⁸
-## <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
-## 1 ae_per… S0001 10 8 16 0 0.699 860 0.00834 0.004
-## 2 ae_per… S0002 10 9 16 0.333 0.696 859 0.307 0.132
-## 3 ae_per… S0003 10 10 15 0.3 0.643 898 0.230 0.106
-## 4 ae_per… S0004 10 9 16 0.111 0.698 859 0.0252 0.016
-## 5 ae_per… S0005 10 9 18 0.556 0.791 704 0.570 0.28
-## 6 ae_per… S0006 10 8 17 0 0.742 787 0.00559 0.007
-## 7 ae_per… S0007 10 10 16 0.6 0.693 858 1 0.44
-## 8 ae_per… S0008 10 10 15 0.5 0.640 898 0.840 0.358
-## 9 ae_per… S0009 10 9 17 0.111 0.742 786 0.0178 0.013
-## 10 ae_per… S0010 10 10 16 0.5 0.695 858 0.699 0.304
-## # … with 290 more rows, and abbreviated variable names ¹study_id, ²site_number,
-## # ³n_pat_with_med75, ⁴visit_med75, ⁵mean_ae_site_med75, ⁶mean_ae_study_med75,
-## # ⁷n_pat_with_med75_study, ⁸prob_lowSimulate Studies
@@ -257,17 +260,9 @@ Simulate Studies
set.seed(1)
-future::plan(multiprocess)
set.seed(1)
-future::plan(multiprocess)## Warning: Strategy 'multiprocess' is deprecated in future (>= 1.20.0)
-## [2020-10-30]. Instead, explicitly specify either 'multisession' (recommended) or
-## 'multicore'. In the current R session, 'multiprocess' equals 'multisession'.## Warning in supportsMulticoreAndRStudio(...): [ONE-TIME WARNING] Forked
-## processing ('multicore') is not supported when running R from RStudio
-## because it is considered unstable. For more details, how to control forked
-## processing or not, and how to silence this warning in future R sessions,
-## see ?parallelly::supportsMulticore
-df_sim_studies <- sim_studies(df_visit, df_site,
+plan(multisession, workers = 6)
+
+df_sim_studies <- sim_studies(df_visit, df_site,
r = 100,
parallel = TRUE,
poisson_test = TRUE,
@@ -276,20 +271,20 @@ Simulate Studies
df_sim_studies
## # A tibble: 30,000 × 8
-## r study_id site_number visit_me…¹ n_pat…² n_pat…³ pval prob_…⁴
-## <dbl> <chr> <chr> <dbl> <int> <dbl> <dbl> <dbl>
-## 1 1 ae_per_visit: 0.05 S0001 16 8 860 0.389 0.208
-## 2 1 ae_per_visit: 0.05 S0002 16 9 859 0.841 0.362
-## 3 1 ae_per_visit: 0.05 S0003 15 10 898 1 1
-## 4 1 ae_per_visit: 0.05 S0004 16 9 859 0.672 0.363
-## 5 1 ae_per_visit: 0.05 S0005 18 9 704 1 1
-## 6 1 ae_per_visit: 0.05 S0006 17 8 787 0.543 0.26
-## 7 1 ae_per_visit: 0.05 S0007 16 10 858 1 0.48
-## 8 1 ae_per_visit: 0.05 S0008 15 10 898 0.105 0.047
-## 9 1 ae_per_visit: 0.05 S0009 17 9 786 1 0.562
-## 10 1 ae_per_visit: 0.05 S0010 16 10 858 1 0.576
-## # … with 29,990 more rows, and abbreviated variable names ¹visit_med75,
-## # ²n_pat_with_med75, ³n_pat_study, ⁴prob_lowCheck p-value Probabilities
@@ -302,28 +297,27 @@ Check p-value Probabilities
+
df_check_pval <- get_ecd_values(df_sim_studies, df_sim_sites, val_str = "pval")
df_check_pval
## # A tibble: 300 × 11
-## study…¹ site_…² n_pat n_pat…³ visit…⁴ mean_…⁵ mean_…⁶ n_pat…⁷ pval prob_…⁸
-## <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
-## 1 ae_per… S0001 10 8 16 0 0.699 860 0.00834 0.004
-## 2 ae_per… S0002 10 9 16 0.333 0.696 859 0.307 0.132
-## 3 ae_per… S0003 10 10 15 0.3 0.643 898 0.230 0.106
-## 4 ae_per… S0004 10 9 16 0.111 0.698 859 0.0252 0.016
-## 5 ae_per… S0005 10 9 18 0.556 0.791 704 0.570 0.28
-## 6 ae_per… S0006 10 8 17 0 0.742 787 0.00559 0.007
-## 7 ae_per… S0007 10 10 16 0.6 0.693 858 1 0.44
-## 8 ae_per… S0008 10 10 15 0.5 0.640 898 0.840 0.358
-## 9 ae_per… S0009 10 9 17 0.111 0.742 786 0.0178 0.013
-## 10 ae_per… S0010 10 10 16 0.5 0.695 858 0.699 0.304
-## # … with 290 more rows, 1 more variable: pval_ecd <dbl>, and abbreviated
-## # variable names ¹study_id, ²site_number, ³n_pat_with_med75, ⁴visit_med75,
-## # ⁵mean_ae_site_med75, ⁶mean_ae_study_med75, ⁷n_pat_with_med75_study,
-## # ⁸prob_low
-
+## study_id site_number n_pat n_pat_with_med75 visit_med75 mean_ae_site_med75
+## <chr> <chr> <int> <int> <dbl> <dbl>
+## 1 ae_per_vis… S0001 10 9 16 0.556
+## 2 ae_per_vis… S0002 10 10 14 0.3
+## 3 ae_per_vis… S0003 10 8 15 0.375
+## 4 ae_per_vis… S0004 10 8 14 0.375
+## 5 ae_per_vis… S0005 10 9 17 0.333
+## 6 ae_per_vis… S0006 10 7 15 0.143
+## 7 ae_per_vis… S0007 10 8 16 0.375
+## 8 ae_per_vis… S0008 10 7 19 0.714
+## 9 ae_per_vis… S0009 10 9 15 0.333
+## 10 ae_per_vis… S0010 10 7 15 0.714
+## # ℹ 290 more rows
+## # ℹ 5 more variables: mean_ae_study_med75 <dbl>, n_pat_with_med75_study <int>,
+## # pval <dbl>, prob_low <dbl>, pval_ecd <dbl>
+
df_check_pval %>%
ggplot(aes(log(pval, base = 10), log(pval_ecd, base = 10))) +
geom_point(alpha = 0.5, size = 2) +
@@ -342,28 +336,27 @@ Perform Same Check on
Note: Dots on the edge of the graph that are cut off have zero value
for the corresponding axis. For the tie simulations using 1000 repeats
the smallest value greater zero we get for prob_low is 0.001 (1e-3).
-
+
df_check_prob <- get_ecd_values(df_sim_studies, df_sim_sites, val_str = "prob_low")
df_check_prob
## # A tibble: 300 × 11
-## study…¹ site_…² n_pat n_pat…³ visit…⁴ mean_…⁵ mean_…⁶ n_pat…⁷ pval prob_…⁸
-## <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
-## 1 ae_per… S0001 10 8 16 0 0.699 860 0.00834 0.004
-## 2 ae_per… S0002 10 9 16 0.333 0.696 859 0.307 0.132
-## 3 ae_per… S0003 10 10 15 0.3 0.643 898 0.230 0.106
-## 4 ae_per… S0004 10 9 16 0.111 0.698 859 0.0252 0.016
-## 5 ae_per… S0005 10 9 18 0.556 0.791 704 0.570 0.28
-## 6 ae_per… S0006 10 8 17 0 0.742 787 0.00559 0.007
-## 7 ae_per… S0007 10 10 16 0.6 0.693 858 1 0.44
-## 8 ae_per… S0008 10 10 15 0.5 0.640 898 0.840 0.358
-## 9 ae_per… S0009 10 9 17 0.111 0.742 786 0.0178 0.013
-## 10 ae_per… S0010 10 10 16 0.5 0.695 858 0.699 0.304
-## # … with 290 more rows, 1 more variable: prob_low_ecd <dbl>, and abbreviated
-## # variable names ¹study_id, ²site_number, ³n_pat_with_med75, ⁴visit_med75,
-## # ⁵mean_ae_site_med75, ⁶mean_ae_study_med75, ⁷n_pat_with_med75_study,
-## # ⁸prob_low
-
+## study_id site_number n_pat n_pat_with_med75 visit_med75 mean_ae_site_med75
+## <chr> <chr> <int> <int> <dbl> <dbl>
+## 1 ae_per_vis… S0001 10 9 16 0.556
+## 2 ae_per_vis… S0002 10 10 14 0.3
+## 3 ae_per_vis… S0003 10 8 15 0.375
+## 4 ae_per_vis… S0004 10 8 14 0.375
+## 5 ae_per_vis… S0005 10 9 17 0.333
+## 6 ae_per_vis… S0006 10 7 15 0.143
+## 7 ae_per_vis… S0007 10 8 16 0.375
+## 8 ae_per_vis… S0008 10 7 19 0.714
+## 9 ae_per_vis… S0009 10 9 15 0.333
+## 10 ae_per_vis… S0010 10 7 15 0.714
+## # ℹ 290 more rows
+## # ℹ 5 more variables: mean_ae_study_med75 <dbl>, n_pat_with_med75_study <int>,
+## # pval <dbl>, prob_low <dbl>, prob_low_ecd <dbl>
+
df_check_prob %>%
ggplot(aes(log(prob_low, base = 10), log(prob_low_ecd, base = 10))) +
geom_point(alpha = 0.5, size = 2) +
@@ -387,6 +380,8 @@ Conclusion
+plan(sequential)
@@ -408,7 +403,7 @@ Conclusion
-Site built with pkgdown 2.0.5.
+Site built with pkgdown 2.0.7.
diff --git a/docs/articles/check_poisson_files/figure-html/unnamed-chunk-4-1.png b/docs/articles/check_poisson_files/figure-html/unnamed-chunk-4-1.png
index 4eb24a0..e25f87d 100644
Binary files a/docs/articles/check_poisson_files/figure-html/unnamed-chunk-4-1.png and b/docs/articles/check_poisson_files/figure-html/unnamed-chunk-4-1.png differ
diff --git a/docs/articles/check_poisson_files/figure-html/unnamed-chunk-5-1.png b/docs/articles/check_poisson_files/figure-html/unnamed-chunk-5-1.png
index deeab82..8a3a9bb 100644
Binary files a/docs/articles/check_poisson_files/figure-html/unnamed-chunk-5-1.png and b/docs/articles/check_poisson_files/figure-html/unnamed-chunk-5-1.png differ
diff --git a/docs/articles/gsm_perf.html b/docs/articles/gsm_perf.html
new file mode 100644
index 0000000..f8ae437
--- /dev/null
+++ b/docs/articles/gsm_perf.html
@@ -0,0 +1,1643 @@
+
+
+
+
+
+
+
+Comparing {simaerep} and {gsm} Performance • simaerep
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Comparing {simaerep} and {gsm} Performance
+
+
+ Source: vignettes/gsm_perf.Rmd
+ gsm_perf.Rmd
+
+
+
+
+
+
+Load
+
+
+suppressPackageStartupMessages(library(tidyverse))
+suppressPackageStartupMessages(library(knitr))
+suppressPackageStartupMessages(library(furrr))
+suppressPackageStartupMessages(library(future))
+suppressPackageStartupMessages(library(simaerep))
+
+plan(multisession, workers = 6)
+
+
+Install {gsm}
+
+
+devtools::install_github("Gilead-BioStats/gsm", ref = "main")
+
+
+Introduction
+
+The {gsm} R
+package provides a standardized Risk Based Quality Monitoring (RBQM)
+framework for clinical trials that pairs a flexible data pipeline with
+robust reports. It also uses Funnel Plots to flag outliers which provide
+broader tolerance limits for sites with low exposure and narrower limits
+for sites with higher exposure. This method is different to the event
+rate based limits we have used in previous heuristics to measure
+{simaerep} performance. Funnel plots are discussed in greater detail by
+Zink et
+al. 2018
+One of the draw backs of using funnel plots for flagging is that they
+assume that the AE rate remains constant over the course of the
+study.
+
+
+Prepare Data
+
+
+Load Portfolio Configurations
+
+We have prepared a snapshot of the AE reporting configuration of our
+current portfolio. For each study we have also measured a visit-specific
+AE rate which allows us to generate a synthetic portfolio with flexible
+AE rates across a study.
+
+df_config <- readr::read_csv("ae_conf_20240220.csv")
+df_ae_rates <- readr::read_csv("ae_rates_20240220.csv")
+
+df_config %>%
+ head(25) %>%
+ knitr::kable()
+
+
+
+
+study_id
+ae_per_visit_mean
+site_number
+max_visit_sd
+max_visit_mean
+n_pat
+
+
+
+0001
+0.2973806
+4746
+3.5355339
+16.50000
+2
+
+
+0001
+0.2973806
+4747
+0.7071068
+27.50000
+2
+
+
+0001
+0.2973806
+4750
+0.5000000
+19.25000
+4
+
+
+0001
+0.2973806
+4815
+8.3765546
+22.16667
+6
+
+
+0001
+0.2973806
+4816
+0.0000000
+27.00000
+1
+
+
+0001
+0.2973806
+4817
+0.0000000
+31.00000
+1
+
+
+0001
+0.2973806
+4818
+0.0000000
+30.00000
+1
+
+
+0001
+0.2973806
+4891
+2.3804761
+22.00000
+7
+
+
+0001
+0.2973806
+4893
+0.0000000
+20.00000
+2
+
+
+0001
+0.2973806
+4932
+1.0954451
+26.80000
+5
+
+
+0001
+0.2973806
+4941
+1.7320508
+27.00000
+3
+
+
+0001
+0.2973806
+4942
+1.4142136
+25.00000
+2
+
+
+0001
+0.2973806
+4943
+0.0000000
+25.00000
+1
+
+
+0001
+0.2973806
+4966
+0.0000000
+18.00000
+1
+
+
+0001
+0.2973806
+4967
+0.8164966
+18.00000
+4
+
+
+0001
+0.2973806
+4968
+0.0000000
+24.00000
+1
+
+
+0001
+0.2973806
+4969
+1.5275252
+21.33333
+3
+
+
+0001
+0.2973806
+4984
+1.7078251
+30.75000
+4
+
+
+0001
+0.2973806
+4985
+8.9087972
+21.16667
+6
+
+
+0001
+0.2973806
+4986
+0.5773503
+20.33333
+3
+
+
+0001
+0.2973806
+4988
+1.4142136
+19.00000
+2
+
+
+0001
+0.2973806
+5079
+2.1602469
+23.00000
+4
+
+
+0001
+0.2973806
+5080
+0.0000000
+23.00000
+1
+
+
+0001
+0.2973806
+5081
+0.0000000
+23.00000
+1
+
+
+0001
+0.2973806
+5082
+1.4142136
+22.00000
+2
+
+
+
+
+
+
+study_id
+cum_visit
+ae_rate
+n_pat
+
+
+
+0001
+1
+0.017699
+113
+
+
+0001
+2
+0.088496
+113
+
+
+0001
+3
+0.230088
+113
+
+
+0001
+4
+0.223214
+112
+
+
+0001
+5
+0.187500
+112
+
+
+0001
+6
+0.196429
+112
+
+
+0001
+7
+0.396396
+111
+
+
+0001
+8
+0.207207
+111
+
+
+0001
+9
+0.216216
+111
+
+
+0001
+10
+0.315315
+111
+
+
+0001
+11
+0.263636
+110
+
+
+0001
+12
+0.318182
+110
+
+
+0001
+13
+0.345455
+110
+
+
+0001
+14
+0.379630
+108
+
+
+0001
+15
+0.401869
+107
+
+
+0001
+16
+0.485981
+107
+
+
+0001
+17
+0.383178
+107
+
+
+0001
+18
+0.396226
+106
+
+
+0001
+19
+0.370000
+100
+
+
+0001
+20
+0.439560
+91
+
+
+0001
+21
+0.369863
+73
+
+
+0001
+22
+0.370968
+62
+
+
+0001
+23
+0.345455
+55
+
+
+0001
+24
+0.446809
+47
+
+
+0001
+25
+0.292683
+41
+
+
+
+
+
+Simulate Portfolio
+
+We generate two synthetic portfolios with no AE under-reporting
+sites. One portfolio with a fixed AE rate for all visits and another one
+with a flexible visit-specific AE rate.
+
+df_portf_fix <- sim_test_data_portfolio(df_config, parallel = TRUE, progress = TRUE)
+
+df_portf_fix %>%
+ head(25) %>%
+ knitr::kable()
+
+
+
+
+study_id
+ae_per_visit_mean
+site_number
+max_visit_sd
+max_visit_mean
+patnum
+visit
+n_ae
+
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+1
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+2
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+3
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+4
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+5
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+6
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+7
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+8
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+9
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+10
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0001
+11
+0
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+1
+1
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+2
+1
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+3
+1
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+4
+2
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+5
+2
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+6
+3
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+7
+4
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+8
+4
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+9
+5
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+10
+5
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+11
+6
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+12
+6
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+13
+6
+
+
+0001
+0.2973806
+4746
+3.535534
+16.5
+0002
+14
+6
+
+
+
+
+df_portf_flex <- sim_test_data_portfolio(df_config, df_ae_rates = df_ae_rates, parallel = TRUE, progress = TRUE)
+
+
+Compare AE rates
+
+Next we confirm the different AE rates in our two synthetic
+portfolios.
+
+df_rate_fix <- df_portf_fix %>%
+ mutate(ae_rate = coalesce(n_ae - lag(n_ae), n_ae), .by = c("study_id", "patnum")) %>%
+ summarise(ae_rate = mean(ae_rate), .by = c("study_id", "visit")) %>%
+ mutate(rate = "fix")
+
+
+df_rate_flex <- df_portf_flex %>%
+ mutate(ae_rate = coalesce(n_ae - lag(n_ae), n_ae), .by = c("study_id", "patnum")) %>%
+ summarise(ae_rate = mean(ae_rate), .by = c("study_id", "visit")) %>%
+ mutate(rate = "flex")
+
+bind_rows(df_rate_flex, df_rate_fix) %>%
+ ggplot(aes(visit, ae_rate)) +
+ geom_line(aes(group = study_id), alpha = 0.2) +
+ geom_smooth() +
+ facet_wrap(~ rate) +
+ labs(title = "Average AE rates per Study")
+
+
+bind_rows(df_rate_flex, df_rate_fix) %>%
+ filter(dense_rank(study_id) <= 16) %>%
+ ggplot(aes(visit, ae_rate)) +
+ geom_line(aes(group = rate, color = rate)) +
+ facet_wrap(~ study_id, scales = "free") +
+ labs(title = "Average AE rates for Selected Studies")
+
+We can confirm that the AE rates in the “flexible” portfolio are not
+constant. Moreover we see that the AE rate profile is very unique for
+each study.
+
+
+
+Apply {gsm}
+
+
+Example
+
+Here we demonstrate how to use the {gsm} package on our simulated
+portfolios.
+
+get_SUBJ <- function(df_portf) {
+ df_portf %>%
+ select(study_id, siteid = site_number, subjid = patnum, timeonstudy = visit) %>%
+ summarise(timeonstudy = max(timeonstudy), .by = c(study_id, siteid, subjid)) %>%
+ group_by(study_id) %>%
+ nest()
+}
+
+
+get_AE <- function(df_portf) {
+ df_portf_fix %>%
+ select(study_id, subjid = patnum, n_ae) %>%
+ summarise(n_ae = max(n_ae), .by = c(study_id, subjid)) %>%
+ filter(n_ae > 0) %>%
+ mutate(n_ae = map(n_ae, ~ tibble(n = seq(1, .)), .progress = TRUE)) %>%
+ unnest(n_ae) %>%
+ select(- n) %>%
+ group_by(study_id) %>%
+ nest()
+}
+
+dfSUBJ_fix <- get_SUBJ(df_portf_fix)
+dfAE_fix <- get_AE(df_portf_fix)
+
+dfInput <- gsm::AE_Map_Raw(list(dfSUBJ = dfSUBJ_fix$data[[1]], dfAE = dfAE_fix$data[[1]]))
+dfInput
+## # A tibble: 113 × 5
+## SubjectID SiteID Exposure Count Rate
+## <chr> <chr> <int> <int> <dbl>
+## 1 0001 4746 11 0 0
+## 2 0002 4746 17 6 0.353
+## 3 0003 4747 27 6 0.222
+## 4 0004 4747 27 5 0.185
+## 5 0005 4750 18 10 0.556
+## 6 0006 4750 19 8 0.421
+## 7 0007 4750 19 8 0.421
+## 8 0008 4750 18 7 0.389
+## 9 0009 4815 15 3 0.2
+## 10 0010 4815 31 2 0.0645
+## # ℹ 103 more rows
+
+dfTransformed <- gsm::Transform_Rate(
+ dfInput,
+ strNumeratorCol = "Count",
+ strDenominatorCol = "Exposure"
+ )
+dfTransformed
+## # A tibble: 44 × 4
+## GroupID Numerator Denominator Metric
+## <chr> <int> <int> <dbl>
+## 1 4746 6 28 0.214
+## 2 4747 11 54 0.204
+## 3 4750 33 74 0.446
+## 4 4815 33 145 0.228
+## 5 4816 8 27 0.296
+## 6 4817 5 31 0.161
+## 7 4818 5 30 0.167
+## 8 4891 45 144 0.312
+## 9 4893 13 40 0.325
+## 10 4932 36 132 0.273
+## # ℹ 34 more rows
+
+dfAnalyzed <- gsm::Analyze_NormalApprox(dfTransformed)
+dfAnalyzed
+## # A tibble: 44 × 7
+## GroupID Numerator Denominator Metric OverallMetric Factor Score
+## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
+## 1 4815 33 145 0.228 0.300 1.31 -1.65
+## 2 4941 17 80 0.212 0.300 1.31 -1.48
+## 3 4817 5 31 0.161 0.300 1.31 -1.47
+## 4 4818 5 30 0.167 0.300 1.31 -1.39
+## 5 5079 20 89 0.225 0.300 1.31 -1.35
+## 6 4747 11 54 0.204 0.300 1.31 -1.34
+## 7 5311 3 20 0.15 0.300 1.31 -1.27
+## 8 5084 4 21 0.190 0.300 1.31 -0.953
+## 9 4943 5 25 0.2 0.300 1.31 -0.949
+## 10 4746 6 28 0.214 0.300 1.31 -0.860
+## # ℹ 34 more rows
+
+dfFlagged <- gsm::Flag_NormalApprox(dfAnalyzed, vThreshold = c(-3, -2, 2, 3))
+dfFlagged
+## # A tibble: 44 × 8
+## GroupID Numerator Denominator Metric OverallMetric Factor Score Flag
+## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
+## 1 5168 9 15 0.6 0.300 1.31 2.22 1
+## 2 5194 76 197 0.386 0.300 1.31 2.31 1
+## 3 4750 33 74 0.446 0.300 1.31 2.40 1
+## 4 4815 33 145 0.228 0.300 1.31 -1.65 0
+## 5 4941 17 80 0.212 0.300 1.31 -1.48 0
+## 6 4817 5 31 0.161 0.300 1.31 -1.47 0
+## 7 4818 5 30 0.167 0.300 1.31 -1.39 0
+## 8 5079 20 89 0.225 0.300 1.31 -1.35 0
+## 9 4747 11 54 0.204 0.300 1.31 -1.34 0
+## 10 5311 3 20 0.15 0.300 1.31 -1.27 0
+## # ℹ 34 more rows
+
+dfSummary <- gsm::Summarize(dfFlagged)
+dfSummary
+## # A tibble: 44 × 6
+## GroupID Numerator Denominator Metric Score Flag
+## <chr> <int> <int> <dbl> <dbl> <dbl>
+## 1 5168 9 15 0.6 2.22 1
+## 2 4750 33 74 0.446 2.40 1
+## 3 5194 76 197 0.386 2.31 1
+## 4 5083 18 44 0.409 1.39 0
+## 5 5265 16 42 0.381 1.01 0
+## 6 5130 52 140 0.371 1.62 0
+## 7 4986 21 58 0.362 0.908 0
+## 8 5131 8 23 0.348 0.442 0
+## 9 4985 38 111 0.342 0.860 0
+## 10 5224 8 24 0.333 0.316 0
+## # ℹ 34 more rows
+
+dfBounds <- gsm::Analyze_NormalApprox_PredictBounds(dfTransformed, vThreshold = c(-3, -2, 2, 3))
+dfBounds
+## # A tibble: 1,246 × 5
+## Threshold Denominator LogDenominator Numerator Metric
+## <dbl> <dbl> <dbl> <dbl> <dbl>
+## 1 -3 28.1 3.34 0.0765 0.00272
+## 2 -3 28.8 3.36 0.187 0.00649
+## 3 -3 29.6 3.39 0.299 0.0101
+## 4 -3 30.3 3.41 0.413 0.0136
+## 5 -3 31.0 3.43 0.527 0.0170
+## 6 -3 31.7 3.46 0.643 0.0203
+## 7 -3 32.5 3.48 0.760 0.0234
+## 8 -3 33.2 3.50 0.878 0.0264
+## 9 -3 33.9 3.52 0.997 0.0294
+## 10 -3 34.7 3.55 1.12 0.0322
+## # ℹ 1,236 more rows
+
+chart <- gsm::Visualize_Scatter(dfFlagged, dfBounds)
+chart
+
+
+
+Simulate UR
+
+We write a function that removes a given ratio of AEs from one site
+in the data set and returns its z-score.
+
+sim_site_ur_gsm <- function(site, ur_rate, dfTransformed) {
+ dfTransformed <- dfTransformed %>%
+ mutate(
+ Numerator = ifelse(GroupID == site, Numerator * (1 - ur_rate), Numerator),
+ Metric = Numerator / Denominator
+ )
+
+ gsm::Analyze_NormalApprox(dfTransformed) %>%
+ filter(GroupID == site) %>%
+ pull(Score)
+}
+
+sim_site_ur_gsm("4747", ur_rate = 0.75, dfTransformed)
+## [1] -3.095796
+We write another function that systematically applies this
+sim_site_ur_gsm across all sites in all studies across a
+range of under-reporting ratios.
+
+sim_ur_gsm <- function(dfSUBJ, dfAE) {
+ dfSUBJ %>%
+ inner_join(dfAE, by = "study_id") %>%
+ ungroup() %>%
+ mutate(
+ trans = map2(data.x, data.y, ~ gsm::AE_Map_Raw(list(dfSUBJ = .x, dfAE = .y))),
+ trans = map(trans, ~ gsm::Transform_Rate(., strNumeratorCol = "Count", strDenominatorCol = "Exposure")),
+ sites = map(data.x, ~ distinct(., siteid))
+ ) %>%
+ select(- starts_with("data.")) %>%
+ unnest(sites) %>%
+ mutate(ur = list(tibble(ur_rate = c(0, 0.1, 0.25, 0.5, 0.75, 1)))) %>%
+ unnest(ur) %>%
+ mutate(
+ score = pmap_dbl(list(siteid, ur_rate, trans), sim_site_ur_gsm, .progress = TRUE)
+ )
+}
+
+df_sim_gsm_fix <- sim_ur_gsm(dfSUBJ_fix, dfAE_fix)
+## Warning: There were 3116 warnings in `mutate()`.
+## The first warning was:
+## ℹ In argument: `z_i = ifelse(...)`.
+## Caused by warning:
+## ! There was 1 warning in `mutate()`.
+## ℹ In argument: `z_0 = ifelse(...)`.
+## Caused by warning in `sqrt()`:
+## ! NaNs produced
+## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3115 remaining warnings.
+
+df_sim_gsm_fix
+## # A tibble: 129,966 × 5
+## study_id trans siteid ur_rate score
+## <chr> <list> <chr> <dbl> <dbl>
+## 1 0001 <tibble [44 × 4]> 4746 0 -0.860
+## 2 0001 <tibble [44 × 4]> 4746 0.1 -1.07
+## 3 0001 <tibble [44 × 4]> 4746 0.25 -1.38
+## 4 0001 <tibble [44 × 4]> 4746 0.5 -1.87
+## 5 0001 <tibble [44 × 4]> 4746 0.75 -2.33
+## 6 0001 <tibble [44 × 4]> 4746 1 -2.75
+## 7 0001 <tibble [44 × 4]> 4747 0 -1.34
+## 8 0001 <tibble [44 × 4]> 4747 0.1 -1.61
+## 9 0001 <tibble [44 × 4]> 4747 0.25 -1.99
+## 10 0001 <tibble [44 × 4]> 4747 0.5 -2.57
+## # ℹ 129,956 more rows
+We repeat the same steps for the portfolio with the flexible AE
+rates.
+
+dfSUBJ_flex <- get_SUBJ(df_portf_flex)
+dfAE_flex<- get_AE(df_portf_flex)
+
+
+df_sim_gsm_flex <- sim_ur_gsm(dfSUBJ_flex, dfAE_flex)
+## Warning: There were 3085 warnings in `mutate()`.
+## The first warning was:
+## ℹ In argument: `z_i = ifelse(...)`.
+## Caused by warning:
+## ! There was 1 warning in `mutate()`.
+## ℹ In argument: `z_0 = ifelse(...)`.
+## Caused by warning in `sqrt()`:
+## ! NaNs produced
+## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3084 remaining warnings.
+
+
+
+UR {simaerep}
+
+We simulate under-reporting for both portfolios using {simaerep}
+
+df_sim_simaerep_fix <- sim_ur_scenarios(
+ df_portf_fix,
+ extra_ur_sites = 0,
+ ur_rate = c(0, 0.1, 0.25, 0.5, 0.75, 1),
+ parallel = TRUE,
+ poisson = TRUE,
+ prob_lower = TRUE,
+ progress = TRUE
+)
+
+df_sim_simaerep_flex <- sim_ur_scenarios(
+ df_portf_flex,
+ extra_ur_sites = 0,
+ ur_rate = c(0, 0.1, 0.25, 0.5, 0.75, 1),
+ parallel = TRUE,
+ poisson = TRUE,
+ prob_lower = TRUE,
+ progress = TRUE
+)
+
+
+Evaluate
+
+
+Combine Results
+
+
+df_sim_gsm_fix$ae_rate <- "AE rate: fix"
+df_sim_gsm_flex$ae_rate <- "AE rate: flexible"
+df_sim_simaerep_fix$ae_rate <- "AE rate: fix"
+df_sim_simaerep_flex$ae_rate <- "AE rate: flexible"
+
+df_sim_thresh2 <- bind_rows(df_sim_gsm_fix, df_sim_gsm_flex) %>%
+ mutate(
+ is_ur = score <= -2,
+ type = "{gsm} - thresh: -2",
+ site_number = siteid
+ ) %>%
+ select(type, ae_rate, study_id, site_number, ur_rate, is_ur, score)
+
+
+df_sim_simaerep_threshp95 <- bind_rows(df_sim_simaerep_fix, df_sim_simaerep_flex) %>%
+ mutate(
+ is_ur = prob_low_prob_ur >= 0.95,
+ type = "{simaerep} - thresh: 0.95"
+ ) %>%
+ select(type, ae_rate, study_id, site_number, ur_rate, is_ur, score = prob_low_prob_ur)
+
+
+df_eval <- bind_rows(
+ df_sim_thresh2,
+ df_sim_simaerep_threshp95,
+)
+
+
+Aggregate
+
+
+get_prop_test_ci95 <- function(..., ix) {
+
+ stopifnot(ix %in% c(1, 2))
+
+ tryCatch(
+ prop.test(...)$conf.int[ix],
+ error = function(cnd) c(NA, NA)[ix]
+ )
+}
+
+aggr_results <- function(df_eval) {
+
+df_perf <- df_eval %>%
+ summarise(
+ n = n(),
+ .by = c(type, ae_rate, ur_rate, is_ur)
+ ) %>%
+ pivot_wider(
+ names_from = is_ur,
+ values_from = n,
+ names_prefix = "is_ur_",
+ values_fill = 0
+ ) %>%
+ mutate(
+ n_sites = is_ur_TRUE + is_ur_FALSE + is_ur_NA,
+ ratio = is_ur_TRUE / n_sites,
+ ratio_type = ifelse(ur_rate == 0, "fpr", "tpr"),
+ ci95_low = map2_dbl(is_ur_TRUE, n_sites, ~ get_prop_test_ci95(.x, .y, ix = 1)),
+ ci95_high = map2_dbl(is_ur_TRUE, n_sites, ~ get_prop_test_ci95(.x, .y, ix = 2))
+ )
+}
+
+df_perf <- aggr_results(df_eval)
+
+
+
+Results
+
+
+Table
+
+
+
+
+
+
+type
+ae_rate
+ur_rate
+is_ur_FALSE
+is_ur_TRUE
+is_ur_NA
+n_sites
+ratio
+ratio_type
+ci95_low
+ci95_high
+
+
+
+{gsm} - thresh: -2
+AE rate: fix
+0.00
+20830
+299
+532
+21661
+0.014
+fpr
+0.012
+0.015
+
+
+{gsm} - thresh: -2
+AE rate: fix
+0.10
+19851
+1278
+532
+21661
+0.059
+tpr
+0.056
+0.062
+
+
+{gsm} - thresh: -2
+AE rate: fix
+0.25
+15848
+5282
+531
+21661
+0.244
+tpr
+0.238
+0.250
+
+
+{gsm} - thresh: -2
+AE rate: fix
+0.50
+9154
+11987
+520
+21661
+0.553
+tpr
+0.547
+0.560
+
+
+{gsm} - thresh: -2
+AE rate: fix
+0.75
+5214
+15941
+506
+21661
+0.736
+tpr
+0.730
+0.742
+
+
+{gsm} - thresh: -2
+AE rate: fix
+1.00
+2915
+18251
+495
+21661
+0.843
+tpr
+0.838
+0.847
+
+
+{gsm} - thresh: -2
+AE rate: flexible
+0.00
+20899
+241
+521
+21661
+0.011
+fpr
+0.010
+0.013
+
+
+{gsm} - thresh: -2
+AE rate: flexible
+0.10
+20669
+471
+521
+21661
+0.022
+tpr
+0.020
+0.024
+
+
+{gsm} - thresh: -2
+AE rate: flexible
+0.25
+19742
+1399
+520
+21661
+0.065
+tpr
+0.061
+0.068
+
+
+{gsm} - thresh: -2
+AE rate: flexible
+0.50
+16448
+4698
+515
+21661
+0.217
+tpr
+0.211
+0.222
+
+
+{gsm} - thresh: -2
+AE rate: flexible
+0.75
+12320
+8830
+511
+21661
+0.408
+tpr
+0.401
+0.414
+
+
+{gsm} - thresh: -2
+AE rate: flexible
+1.00
+8727
+12437
+497
+21661
+0.574
+tpr
+0.568
+0.581
+
+
+{simaerep} - thresh: 0.95
+AE rate: fix
+0.00
+21367
+294
+0
+21661
+0.014
+fpr
+0.012
+0.015
+
+
+{simaerep} - thresh: 0.95
+AE rate: fix
+0.10
+20366
+1295
+0
+21661
+0.060
+tpr
+0.057
+0.063
+
+
+{simaerep} - thresh: 0.95
+AE rate: fix
+0.25
+17155
+4506
+0
+21661
+0.208
+tpr
+0.203
+0.214
+
+
+{simaerep} - thresh: 0.95
+AE rate: fix
+0.50
+10635
+11026
+0
+21661
+0.509
+tpr
+0.502
+0.516
+
+
+{simaerep} - thresh: 0.95
+AE rate: fix
+0.75
+5845
+15816
+0
+21661
+0.730
+tpr
+0.724
+0.736
+
+
+{simaerep} - thresh: 0.95
+AE rate: fix
+1.00
+4053
+17608
+0
+21661
+0.813
+tpr
+0.808
+0.818
+
+
+{simaerep} - thresh: 0.95
+AE rate: flexible
+0.00
+21340
+321
+0
+21661
+0.015
+fpr
+0.013
+0.017
+
+
+{simaerep} - thresh: 0.95
+AE rate: flexible
+0.10
+20234
+1427
+0
+21661
+0.066
+tpr
+0.063
+0.069
+
+
+{simaerep} - thresh: 0.95
+AE rate: flexible
+0.25
+16661
+5000
+0
+21661
+0.231
+tpr
+0.225
+0.237
+
+
+{simaerep} - thresh: 0.95
+AE rate: flexible
+0.50
+9878
+11783
+0
+21661
+0.544
+tpr
+0.537
+0.551
+
+
+{simaerep} - thresh: 0.95
+AE rate: flexible
+0.75
+5216
+16445
+0
+21661
+0.759
+tpr
+0.753
+0.765
+
+
+{simaerep} - thresh: 0.95
+AE rate: flexible
+1.00
+3595
+18066
+0
+21661
+0.834
+tpr
+0.829
+0.839
+
+
+
+
+
+Plot Performance Metrics
+
+
+- {gsm} has better performance than {simaerep} when the AE rate is
+fixed, while {simaerep} greatly outperforms {gsm} when the AE rate is
+flexible and mimics the AE rates encountered in real study data
+sets.
+
+
+plot_perf <- function(df_perf) {
+
+ df_perf %>%
+ mutate(ur_rate = paste0("under-reporting rate: ", ur_rate, " - ", ratio_type),
+ ur_rate = ifelse(str_detect(ur_rate, "fpr"), "fpr", ur_rate)) %>%
+ group_by(ur_rate) %>%
+ ggplot(aes(type, ratio)) +
+ geom_errorbar(aes(ymin = ci95_low, ymax = ci95_high, color = type), linewidth = 1) +
+ facet_grid(ur_rate ~ ae_rate) +
+ coord_flip() +
+ theme(legend.position = "bottom") +
+ labs(
+ x = "",
+ y = "CI95 Performance Ratio",
+ title = "{simaerep} vs {gsm} Performance"
+ ) +
+ scale_color_manual(values = c("#5491CC", "#F46626"))
+}
+
+plot_perf(df_perf)
+
+
+plan(sequential)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/articles/gsm_perf_files/figure-html/plot-1.png b/docs/articles/gsm_perf_files/figure-html/plot-1.png
new file mode 100644
index 0000000..2feedb9
Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/plot-1.png differ
diff --git a/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-10-1.png b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-10-1.png
new file mode 100644
index 0000000..5098a97
Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-10-1.png differ
diff --git a/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-3-1.png b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-3-1.png
new file mode 100644
index 0000000..91a7a28
Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-3-1.png differ
diff --git a/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-4-1.png b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-4-1.png
new file mode 100644
index 0000000..e27fc8a
Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-4-1.png differ
diff --git a/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-6-1.png b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-6-1.png
new file mode 100644
index 0000000..524a0f5
Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-6-1.png differ
diff --git a/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-9-1.png b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-9-1.png
new file mode 100644
index 0000000..7289510
Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-9-1.png differ
diff --git a/docs/articles/index.html b/docs/articles/index.html
index 8b6151d..d1f3485 100644
--- a/docs/articles/index.html
+++ b/docs/articles/index.html
@@ -23,7 +23,7 @@
+df_check_pval <- get_ecd_values(df_sim_studies, df_sim_sites, val_str = "pval") df_check_pval-## # A tibble: 300 × 11 -## study…¹ site_…² n_pat n_pat…³ visit…⁴ mean_…⁵ mean_…⁶ n_pat…⁷ pval prob_…⁸ -## <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <int> <dbl> <dbl> -## 1 ae_per… S0001 10 8 16 0 0.699 860 0.00834 0.004 -## 2 ae_per… S0002 10 9 16 0.333 0.696 859 0.307 0.132 -## 3 ae_per… S0003 10 10 15 0.3 0.643 898 0.230 0.106 -## 4 ae_per… S0004 10 9 16 0.111 0.698 859 0.0252 0.016 -## 5 ae_per… S0005 10 9 18 0.556 0.791 704 0.570 0.28 -## 6 ae_per… S0006 10 8 17 0 0.742 787 0.00559 0.007 -## 7 ae_per… S0007 10 10 16 0.6 0.693 858 1 0.44 -## 8 ae_per… S0008 10 10 15 0.5 0.640 898 0.840 0.358 -## 9 ae_per… S0009 10 9 17 0.111 0.742 786 0.0178 0.013 -## 10 ae_per… S0010 10 10 16 0.5 0.695 858 0.699 0.304 -## # … with 290 more rows, 1 more variable: pval_ecd <dbl>, and abbreviated -## # variable names ¹study_id, ²site_number, ³n_pat_with_med75, ⁴visit_med75, -## # ⁵mean_ae_site_med75, ⁶mean_ae_study_med75, ⁷n_pat_with_med75_study, -## # ⁸prob_lowdiff --git a/docs/articles/check_poisson_files/figure-html/unnamed-chunk-4-1.png b/docs/articles/check_poisson_files/figure-html/unnamed-chunk-4-1.png index 4eb24a0..e25f87d 100644 Binary files a/docs/articles/check_poisson_files/figure-html/unnamed-chunk-4-1.png and b/docs/articles/check_poisson_files/figure-html/unnamed-chunk-4-1.png differ diff --git a/docs/articles/check_poisson_files/figure-html/unnamed-chunk-5-1.png b/docs/articles/check_poisson_files/figure-html/unnamed-chunk-5-1.png index deeab82..8a3a9bb 100644 Binary files a/docs/articles/check_poisson_files/figure-html/unnamed-chunk-5-1.png and b/docs/articles/check_poisson_files/figure-html/unnamed-chunk-5-1.png differ diff --git a/docs/articles/gsm_perf.html b/docs/articles/gsm_perf.html new file mode 100644 index 0000000..f8ae437 --- /dev/null +++ b/docs/articles/gsm_perf.html @@ -0,0 +1,1643 @@ + + + + + + + ++## study_id site_number n_pat n_pat_with_med75 visit_med75 mean_ae_site_med75 +## <chr> <chr> <int> <int> <dbl> <dbl> +## 1 ae_per_vis… S0001 10 9 16 0.556 +## 2 ae_per_vis… S0002 10 10 14 0.3 +## 3 ae_per_vis… S0003 10 8 15 0.375 +## 4 ae_per_vis… S0004 10 8 14 0.375 +## 5 ae_per_vis… S0005 10 9 17 0.333 +## 6 ae_per_vis… S0006 10 7 15 0.143 +## 7 ae_per_vis… S0007 10 8 16 0.375 +## 8 ae_per_vis… S0008 10 7 19 0.714 +## 9 ae_per_vis… S0009 10 9 15 0.333 +## 10 ae_per_vis… S0010 10 7 15 0.714 +## # ℹ 290 more rows +## # ℹ 5 more variables: mean_ae_study_med75 <dbl>, n_pat_with_med75_study <int>, +## # pval <dbl>, prob_low <dbl>, pval_ecd <dbl>+@@ -408,7 +403,7 @@df_check_pval %>% ggplot(aes(log(pval, base = 10), log(pval_ecd, base = 10))) + geom_point(alpha = 0.5, size = 2) + @@ -342,28 +336,27 @@Perform Same Check on
Note: Dots on the edge of the graph that are cut off have zero value for the corresponding axis. For the tie simulations using 1000 repeats the smallest value greater zero we get for prob_low is 0.001 (1e-3).
-+df_check_prob <- get_ecd_values(df_sim_studies, df_sim_sites, val_str = "prob_low") df_check_prob-## # A tibble: 300 × 11 -## study…¹ site_…² n_pat n_pat…³ visit…⁴ mean_…⁵ mean_…⁶ n_pat…⁷ pval prob_…⁸ -## <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <int> <dbl> <dbl> -## 1 ae_per… S0001 10 8 16 0 0.699 860 0.00834 0.004 -## 2 ae_per… S0002 10 9 16 0.333 0.696 859 0.307 0.132 -## 3 ae_per… S0003 10 10 15 0.3 0.643 898 0.230 0.106 -## 4 ae_per… S0004 10 9 16 0.111 0.698 859 0.0252 0.016 -## 5 ae_per… S0005 10 9 18 0.556 0.791 704 0.570 0.28 -## 6 ae_per… S0006 10 8 17 0 0.742 787 0.00559 0.007 -## 7 ae_per… S0007 10 10 16 0.6 0.693 858 1 0.44 -## 8 ae_per… S0008 10 10 15 0.5 0.640 898 0.840 0.358 -## 9 ae_per… S0009 10 9 17 0.111 0.742 786 0.0178 0.013 -## 10 ae_per… S0010 10 10 16 0.5 0.695 858 0.699 0.304 -## # … with 290 more rows, 1 more variable: prob_low_ecd <dbl>, and abbreviated -## # variable names ¹study_id, ²site_number, ³n_pat_with_med75, ⁴visit_med75, -## # ⁵mean_ae_site_med75, ⁶mean_ae_study_med75, ⁷n_pat_with_med75_study, -## # ⁸prob_low+## study_id site_number n_pat n_pat_with_med75 visit_med75 mean_ae_site_med75 +## <chr> <chr> <int> <int> <dbl> <dbl> +## 1 ae_per_vis… S0001 10 9 16 0.556 +## 2 ae_per_vis… S0002 10 10 14 0.3 +## 3 ae_per_vis… S0003 10 8 15 0.375 +## 4 ae_per_vis… S0004 10 8 14 0.375 +## 5 ae_per_vis… S0005 10 9 17 0.333 +## 6 ae_per_vis… S0006 10 7 15 0.143 +## 7 ae_per_vis… S0007 10 8 16 0.375 +## 8 ae_per_vis… S0008 10 7 19 0.714 +## 9 ae_per_vis… S0009 10 9 15 0.333 +## 10 ae_per_vis… S0010 10 7 15 0.714 +## # ℹ 290 more rows +## # ℹ 5 more variables: mean_ae_study_med75 <dbl>, n_pat_with_med75_study <int>, +## # pval <dbl>, prob_low <dbl>, prob_low_ecd <dbl>+df_check_prob %>% ggplot(aes(log(prob_low, base = 10), log(prob_low_ecd, base = 10))) + geom_point(alpha = 0.5, size = 2) + @@ -387,6 +380,8 @@Conclusion
+plan(sequential)Conclusion -
Site built with pkgdown 2.0.5.
+Site built with pkgdown 2.0.7.
Comparing {simaerep} and {gsm} Performance • simaerep + + + + + + + + + + + + + + + + + + + +++ + + + + + + + diff --git a/docs/articles/gsm_perf_files/figure-html/plot-1.png b/docs/articles/gsm_perf_files/figure-html/plot-1.png new file mode 100644 index 0000000..2feedb9 Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/plot-1.png differ diff --git a/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-10-1.png b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-10-1.png new file mode 100644 index 0000000..5098a97 Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-10-1.png differ diff --git a/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-3-1.png b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-3-1.png new file mode 100644 index 0000000..91a7a28 Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-3-1.png differ diff --git a/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-4-1.png b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-4-1.png new file mode 100644 index 0000000..e27fc8a Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-4-1.png differ diff --git a/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-6-1.png b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-6-1.png new file mode 100644 index 0000000..524a0f5 Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-6-1.png differ diff --git a/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-9-1.png b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-9-1.png new file mode 100644 index 0000000..7289510 Binary files /dev/null and b/docs/articles/gsm_perf_files/figure-html/unnamed-chunk-9-1.png differ diff --git a/docs/articles/index.html b/docs/articles/index.html index 8b6151d..d1f3485 100644 --- a/docs/articles/index.html +++ b/docs/articles/index.html @@ -23,7 +23,7 @@+ + + + + ++ + + + +++ + + +++ + + +Comparing {simaerep} and {gsm} Performance
+ + + Source:vignettes/gsm_perf.Rmd++ +gsm_perf.Rmd++Load +
+++suppressPackageStartupMessages(library(tidyverse)) +suppressPackageStartupMessages(library(knitr)) +suppressPackageStartupMessages(library(furrr)) +suppressPackageStartupMessages(library(future)) +suppressPackageStartupMessages(library(simaerep)) + +plan(multisession, workers = 6)++Install {gsm} +
+++devtools::install_github("Gilead-BioStats/gsm", ref = "main")++Introduction +
+The {gsm} R +package provides a standardized Risk Based Quality Monitoring (RBQM) +framework for clinical trials that pairs a flexible data pipeline with +robust reports. It also uses Funnel Plots to flag outliers which provide +broader tolerance limits for sites with low exposure and narrower limits +for sites with higher exposure. This method is different to the event +rate based limits we have used in previous heuristics to measure +{simaerep} performance. Funnel plots are discussed in greater detail by +Zink et +al. 2018
+One of the draw backs of using funnel plots for flagging is that they +assume that the AE rate remains constant over the course of the +study.
+++Prepare Data +
+++Load Portfolio Configurations +
+We have prepared a snapshot of the AE reporting configuration of our +current portfolio. For each study we have also measured a visit-specific +AE rate which allows us to generate a synthetic portfolio with flexible +AE rates across a study.
+++df_config <- readr::read_csv("ae_conf_20240220.csv") +df_ae_rates <- readr::read_csv("ae_rates_20240220.csv") + +df_config %>% + head(25) %>% + knitr::kable()+
+ ++ ++ + + + + + + + +study_id +ae_per_visit_mean +site_number +max_visit_sd +max_visit_mean +n_pat ++ +0001 +0.2973806 +4746 +3.5355339 +16.50000 +2 ++ +0001 +0.2973806 +4747 +0.7071068 +27.50000 +2 ++ +0001 +0.2973806 +4750 +0.5000000 +19.25000 +4 ++ +0001 +0.2973806 +4815 +8.3765546 +22.16667 +6 ++ +0001 +0.2973806 +4816 +0.0000000 +27.00000 +1 ++ +0001 +0.2973806 +4817 +0.0000000 +31.00000 +1 ++ +0001 +0.2973806 +4818 +0.0000000 +30.00000 +1 ++ +0001 +0.2973806 +4891 +2.3804761 +22.00000 +7 ++ +0001 +0.2973806 +4893 +0.0000000 +20.00000 +2 ++ +0001 +0.2973806 +4932 +1.0954451 +26.80000 +5 ++ +0001 +0.2973806 +4941 +1.7320508 +27.00000 +3 ++ +0001 +0.2973806 +4942 +1.4142136 +25.00000 +2 ++ +0001 +0.2973806 +4943 +0.0000000 +25.00000 +1 ++ +0001 +0.2973806 +4966 +0.0000000 +18.00000 +1 ++ +0001 +0.2973806 +4967 +0.8164966 +18.00000 +4 ++ +0001 +0.2973806 +4968 +0.0000000 +24.00000 +1 ++ +0001 +0.2973806 +4969 +1.5275252 +21.33333 +3 ++ +0001 +0.2973806 +4984 +1.7078251 +30.75000 +4 ++ +0001 +0.2973806 +4985 +8.9087972 +21.16667 +6 ++ +0001 +0.2973806 +4986 +0.5773503 +20.33333 +3 ++ +0001 +0.2973806 +4988 +1.4142136 +19.00000 +2 ++ +0001 +0.2973806 +5079 +2.1602469 +23.00000 +4 ++ +0001 +0.2973806 +5080 +0.0000000 +23.00000 +1 ++ +0001 +0.2973806 +5081 +0.0000000 +23.00000 +1 ++ + +0001 +0.2973806 +5082 +1.4142136 +22.00000 +2 ++
++ + +study_id +cum_visit +ae_rate +n_pat ++ +0001 +1 +0.017699 +113 ++ +0001 +2 +0.088496 +113 ++ +0001 +3 +0.230088 +113 ++ +0001 +4 +0.223214 +112 ++ +0001 +5 +0.187500 +112 ++ +0001 +6 +0.196429 +112 ++ +0001 +7 +0.396396 +111 ++ +0001 +8 +0.207207 +111 ++ +0001 +9 +0.216216 +111 ++ +0001 +10 +0.315315 +111 ++ +0001 +11 +0.263636 +110 ++ +0001 +12 +0.318182 +110 ++ +0001 +13 +0.345455 +110 ++ +0001 +14 +0.379630 +108 ++ +0001 +15 +0.401869 +107 ++ +0001 +16 +0.485981 +107 ++ +0001 +17 +0.383178 +107 ++ +0001 +18 +0.396226 +106 ++ +0001 +19 +0.370000 +100 ++ +0001 +20 +0.439560 +91 ++ +0001 +21 +0.369863 +73 ++ +0001 +22 +0.370968 +62 ++ +0001 +23 +0.345455 +55 ++ +0001 +24 +0.446809 +47 ++ + +0001 +25 +0.292683 +41 +++Simulate Portfolio +
+We generate two synthetic portfolios with no AE under-reporting +sites. One portfolio with a fixed AE rate for all visits and another one +with a flexible visit-specific AE rate.
+++df_portf_fix <- sim_test_data_portfolio(df_config, parallel = TRUE, progress = TRUE) + +df_portf_fix %>% + head(25) %>% + knitr::kable()+
++ ++ + + + + + + + + + +study_id +ae_per_visit_mean +site_number +max_visit_sd +max_visit_mean +patnum +visit +n_ae ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0001 +1 +0 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0001 +2 +0 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0001 +3 +0 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0001 +4 +0 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0001 +5 +0 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0001 +6 +0 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0001 +7 +0 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0001 +8 +0 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0001 +9 +0 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0001 +10 +0 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0001 +11 +0 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0002 +1 +1 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0002 +2 +1 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0002 +3 +1 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0002 +4 +2 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0002 +5 +2 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0002 +6 +3 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0002 +7 +4 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0002 +8 +4 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0002 +9 +5 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0002 +10 +5 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0002 +11 +6 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0002 +12 +6 ++ +0001 +0.2973806 +4746 +3.535534 +16.5 +0002 +13 +6 ++ + +0001 +0.2973806 +4746 +3.535534 +16.5 +0002 +14 +6 +++df_portf_flex <- sim_test_data_portfolio(df_config, df_ae_rates = df_ae_rates, parallel = TRUE, progress = TRUE)++Compare AE rates +
+Next we confirm the different AE rates in our two synthetic +portfolios.
+++df_rate_fix <- df_portf_fix %>% + mutate(ae_rate = coalesce(n_ae - lag(n_ae), n_ae), .by = c("study_id", "patnum")) %>% + summarise(ae_rate = mean(ae_rate), .by = c("study_id", "visit")) %>% + mutate(rate = "fix") + + +df_rate_flex <- df_portf_flex %>% + mutate(ae_rate = coalesce(n_ae - lag(n_ae), n_ae), .by = c("study_id", "patnum")) %>% + summarise(ae_rate = mean(ae_rate), .by = c("study_id", "visit")) %>% + mutate(rate = "flex")++bind_rows(df_rate_flex, df_rate_fix) %>% + ggplot(aes(visit, ae_rate)) + + geom_line(aes(group = study_id), alpha = 0.2) + + geom_smooth() + + facet_wrap(~ rate) + + labs(title = "Average AE rates per Study")+
++bind_rows(df_rate_flex, df_rate_fix) %>% + filter(dense_rank(study_id) <= 16) %>% + ggplot(aes(visit, ae_rate)) + + geom_line(aes(group = rate, color = rate)) + + facet_wrap(~ study_id, scales = "free") + + labs(title = "Average AE rates for Selected Studies")+
We can confirm that the AE rates in the “flexible” portfolio are not +constant. Moreover we see that the AE rate profile is very unique for +each study.
+++Apply {gsm} +
+++Example +
+Here we demonstrate how to use the {gsm} package on our simulated +portfolios.
+++get_SUBJ <- function(df_portf) { + df_portf %>% + select(study_id, siteid = site_number, subjid = patnum, timeonstudy = visit) %>% + summarise(timeonstudy = max(timeonstudy), .by = c(study_id, siteid, subjid)) %>% + group_by(study_id) %>% + nest() +} + + +get_AE <- function(df_portf) { + df_portf_fix %>% + select(study_id, subjid = patnum, n_ae) %>% + summarise(n_ae = max(n_ae), .by = c(study_id, subjid)) %>% + filter(n_ae > 0) %>% + mutate(n_ae = map(n_ae, ~ tibble(n = seq(1, .)), .progress = TRUE)) %>% + unnest(n_ae) %>% + select(- n) %>% + group_by(study_id) %>% + nest() +} + +dfSUBJ_fix <- get_SUBJ(df_portf_fix) +dfAE_fix <- get_AE(df_portf_fix)++dfInput <- gsm::AE_Map_Raw(list(dfSUBJ = dfSUBJ_fix$data[[1]], dfAE = dfAE_fix$data[[1]])) +dfInput+## # A tibble: 113 × 5 +## SubjectID SiteID Exposure Count Rate +## <chr> <chr> <int> <int> <dbl> +## 1 0001 4746 11 0 0 +## 2 0002 4746 17 6 0.353 +## 3 0003 4747 27 6 0.222 +## 4 0004 4747 27 5 0.185 +## 5 0005 4750 18 10 0.556 +## 6 0006 4750 19 8 0.421 +## 7 0007 4750 19 8 0.421 +## 8 0008 4750 18 7 0.389 +## 9 0009 4815 15 3 0.2 +## 10 0010 4815 31 2 0.0645 +## # ℹ 103 more rows++dfTransformed <- gsm::Transform_Rate( + dfInput, + strNumeratorCol = "Count", + strDenominatorCol = "Exposure" + ) +dfTransformed+## # A tibble: 44 × 4 +## GroupID Numerator Denominator Metric +## <chr> <int> <int> <dbl> +## 1 4746 6 28 0.214 +## 2 4747 11 54 0.204 +## 3 4750 33 74 0.446 +## 4 4815 33 145 0.228 +## 5 4816 8 27 0.296 +## 6 4817 5 31 0.161 +## 7 4818 5 30 0.167 +## 8 4891 45 144 0.312 +## 9 4893 13 40 0.325 +## 10 4932 36 132 0.273 +## # ℹ 34 more rows++dfAnalyzed <- gsm::Analyze_NormalApprox(dfTransformed) +dfAnalyzed+## # A tibble: 44 × 7 +## GroupID Numerator Denominator Metric OverallMetric Factor Score +## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> +## 1 4815 33 145 0.228 0.300 1.31 -1.65 +## 2 4941 17 80 0.212 0.300 1.31 -1.48 +## 3 4817 5 31 0.161 0.300 1.31 -1.47 +## 4 4818 5 30 0.167 0.300 1.31 -1.39 +## 5 5079 20 89 0.225 0.300 1.31 -1.35 +## 6 4747 11 54 0.204 0.300 1.31 -1.34 +## 7 5311 3 20 0.15 0.300 1.31 -1.27 +## 8 5084 4 21 0.190 0.300 1.31 -0.953 +## 9 4943 5 25 0.2 0.300 1.31 -0.949 +## 10 4746 6 28 0.214 0.300 1.31 -0.860 +## # ℹ 34 more rows++dfFlagged <- gsm::Flag_NormalApprox(dfAnalyzed, vThreshold = c(-3, -2, 2, 3)) +dfFlagged+## # A tibble: 44 × 8 +## GroupID Numerator Denominator Metric OverallMetric Factor Score Flag +## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> +## 1 5168 9 15 0.6 0.300 1.31 2.22 1 +## 2 5194 76 197 0.386 0.300 1.31 2.31 1 +## 3 4750 33 74 0.446 0.300 1.31 2.40 1 +## 4 4815 33 145 0.228 0.300 1.31 -1.65 0 +## 5 4941 17 80 0.212 0.300 1.31 -1.48 0 +## 6 4817 5 31 0.161 0.300 1.31 -1.47 0 +## 7 4818 5 30 0.167 0.300 1.31 -1.39 0 +## 8 5079 20 89 0.225 0.300 1.31 -1.35 0 +## 9 4747 11 54 0.204 0.300 1.31 -1.34 0 +## 10 5311 3 20 0.15 0.300 1.31 -1.27 0 +## # ℹ 34 more rows++dfSummary <- gsm::Summarize(dfFlagged) +dfSummary+## # A tibble: 44 × 6 +## GroupID Numerator Denominator Metric Score Flag +## <chr> <int> <int> <dbl> <dbl> <dbl> +## 1 5168 9 15 0.6 2.22 1 +## 2 4750 33 74 0.446 2.40 1 +## 3 5194 76 197 0.386 2.31 1 +## 4 5083 18 44 0.409 1.39 0 +## 5 5265 16 42 0.381 1.01 0 +## 6 5130 52 140 0.371 1.62 0 +## 7 4986 21 58 0.362 0.908 0 +## 8 5131 8 23 0.348 0.442 0 +## 9 4985 38 111 0.342 0.860 0 +## 10 5224 8 24 0.333 0.316 0 +## # ℹ 34 more rows++dfBounds <- gsm::Analyze_NormalApprox_PredictBounds(dfTransformed, vThreshold = c(-3, -2, 2, 3)) +dfBounds+## # A tibble: 1,246 × 5 +## Threshold Denominator LogDenominator Numerator Metric +## <dbl> <dbl> <dbl> <dbl> <dbl> +## 1 -3 28.1 3.34 0.0765 0.00272 +## 2 -3 28.8 3.36 0.187 0.00649 +## 3 -3 29.6 3.39 0.299 0.0101 +## 4 -3 30.3 3.41 0.413 0.0136 +## 5 -3 31.0 3.43 0.527 0.0170 +## 6 -3 31.7 3.46 0.643 0.0203 +## 7 -3 32.5 3.48 0.760 0.0234 +## 8 -3 33.2 3.50 0.878 0.0264 +## 9 -3 33.9 3.52 0.997 0.0294 +## 10 -3 34.7 3.55 1.12 0.0322 +## # ℹ 1,236 more rows++chart <- gsm::Visualize_Scatter(dfFlagged, dfBounds) +chart+
++Simulate UR +
+We write a function that removes a given ratio of AEs from one site +in the data set and returns its z-score.
+++sim_site_ur_gsm <- function(site, ur_rate, dfTransformed) { + dfTransformed <- dfTransformed %>% + mutate( + Numerator = ifelse(GroupID == site, Numerator * (1 - ur_rate), Numerator), + Metric = Numerator / Denominator + ) + + gsm::Analyze_NormalApprox(dfTransformed) %>% + filter(GroupID == site) %>% + pull(Score) +} + +sim_site_ur_gsm("4747", ur_rate = 0.75, dfTransformed)+## [1] -3.095796We write another function that systematically applies this +
+sim_site_ur_gsmacross all sites in all studies across a +range of under-reporting ratios.++sim_ur_gsm <- function(dfSUBJ, dfAE) { + dfSUBJ %>% + inner_join(dfAE, by = "study_id") %>% + ungroup() %>% + mutate( + trans = map2(data.x, data.y, ~ gsm::AE_Map_Raw(list(dfSUBJ = .x, dfAE = .y))), + trans = map(trans, ~ gsm::Transform_Rate(., strNumeratorCol = "Count", strDenominatorCol = "Exposure")), + sites = map(data.x, ~ distinct(., siteid)) + ) %>% + select(- starts_with("data.")) %>% + unnest(sites) %>% + mutate(ur = list(tibble(ur_rate = c(0, 0.1, 0.25, 0.5, 0.75, 1)))) %>% + unnest(ur) %>% + mutate( + score = pmap_dbl(list(siteid, ur_rate, trans), sim_site_ur_gsm, .progress = TRUE) + ) +} + +df_sim_gsm_fix <- sim_ur_gsm(dfSUBJ_fix, dfAE_fix)+## Warning: There were 3116 warnings in `mutate()`. +## The first warning was: +## ℹ In argument: `z_i = ifelse(...)`. +## Caused by warning: +## ! There was 1 warning in `mutate()`. +## ℹ In argument: `z_0 = ifelse(...)`. +## Caused by warning in `sqrt()`: +## ! NaNs produced +## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3115 remaining warnings.++df_sim_gsm_fix+## # A tibble: 129,966 × 5 +## study_id trans siteid ur_rate score +## <chr> <list> <chr> <dbl> <dbl> +## 1 0001 <tibble [44 × 4]> 4746 0 -0.860 +## 2 0001 <tibble [44 × 4]> 4746 0.1 -1.07 +## 3 0001 <tibble [44 × 4]> 4746 0.25 -1.38 +## 4 0001 <tibble [44 × 4]> 4746 0.5 -1.87 +## 5 0001 <tibble [44 × 4]> 4746 0.75 -2.33 +## 6 0001 <tibble [44 × 4]> 4746 1 -2.75 +## 7 0001 <tibble [44 × 4]> 4747 0 -1.34 +## 8 0001 <tibble [44 × 4]> 4747 0.1 -1.61 +## 9 0001 <tibble [44 × 4]> 4747 0.25 -1.99 +## 10 0001 <tibble [44 × 4]> 4747 0.5 -2.57 +## # ℹ 129,956 more rowsWe repeat the same steps for the portfolio with the flexible AE +rates.
+++dfSUBJ_flex <- get_SUBJ(df_portf_flex) +dfAE_flex<- get_AE(df_portf_flex) + + +df_sim_gsm_flex <- sim_ur_gsm(dfSUBJ_flex, dfAE_flex)+## Warning: There were 3085 warnings in `mutate()`. +## The first warning was: +## ℹ In argument: `z_i = ifelse(...)`. +## Caused by warning: +## ! There was 1 warning in `mutate()`. +## ℹ In argument: `z_0 = ifelse(...)`. +## Caused by warning in `sqrt()`: +## ! NaNs produced +## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3084 remaining warnings.++UR {simaerep} +
+We simulate under-reporting for both portfolios using {simaerep}
+++df_sim_simaerep_fix <- sim_ur_scenarios( + df_portf_fix, + extra_ur_sites = 0, + ur_rate = c(0, 0.1, 0.25, 0.5, 0.75, 1), + parallel = TRUE, + poisson = TRUE, + prob_lower = TRUE, + progress = TRUE +)++df_sim_simaerep_flex <- sim_ur_scenarios( + df_portf_flex, + extra_ur_sites = 0, + ur_rate = c(0, 0.1, 0.25, 0.5, 0.75, 1), + parallel = TRUE, + poisson = TRUE, + prob_lower = TRUE, + progress = TRUE +)++Evaluate +
+++Combine Results +
+++df_sim_gsm_fix$ae_rate <- "AE rate: fix" +df_sim_gsm_flex$ae_rate <- "AE rate: flexible" +df_sim_simaerep_fix$ae_rate <- "AE rate: fix" +df_sim_simaerep_flex$ae_rate <- "AE rate: flexible" + +df_sim_thresh2 <- bind_rows(df_sim_gsm_fix, df_sim_gsm_flex) %>% + mutate( + is_ur = score <= -2, + type = "{gsm} - thresh: -2", + site_number = siteid + ) %>% + select(type, ae_rate, study_id, site_number, ur_rate, is_ur, score) + + +df_sim_simaerep_threshp95 <- bind_rows(df_sim_simaerep_fix, df_sim_simaerep_flex) %>% + mutate( + is_ur = prob_low_prob_ur >= 0.95, + type = "{simaerep} - thresh: 0.95" + ) %>% + select(type, ae_rate, study_id, site_number, ur_rate, is_ur, score = prob_low_prob_ur) + + +df_eval <- bind_rows( + df_sim_thresh2, + df_sim_simaerep_threshp95, +)++Aggregate +
+++get_prop_test_ci95 <- function(..., ix) { + + stopifnot(ix %in% c(1, 2)) + + tryCatch( + prop.test(...)$conf.int[ix], + error = function(cnd) c(NA, NA)[ix] + ) +} + +aggr_results <- function(df_eval) { + +df_perf <- df_eval %>% + summarise( + n = n(), + .by = c(type, ae_rate, ur_rate, is_ur) + ) %>% + pivot_wider( + names_from = is_ur, + values_from = n, + names_prefix = "is_ur_", + values_fill = 0 + ) %>% + mutate( + n_sites = is_ur_TRUE + is_ur_FALSE + is_ur_NA, + ratio = is_ur_TRUE / n_sites, + ratio_type = ifelse(ur_rate == 0, "fpr", "tpr"), + ci95_low = map2_dbl(is_ur_TRUE, n_sites, ~ get_prop_test_ci95(.x, .y, ix = 1)), + ci95_high = map2_dbl(is_ur_TRUE, n_sites, ~ get_prop_test_ci95(.x, .y, ix = 2)) + ) +} + +df_perf <- aggr_results(df_eval)++Results +
+++Table +
+ ++
++ ++ + + + + + + + + + + + + +type +ae_rate +ur_rate +is_ur_FALSE +is_ur_TRUE +is_ur_NA +n_sites +ratio +ratio_type +ci95_low +ci95_high ++ +{gsm} - thresh: -2 +AE rate: fix +0.00 +20830 +299 +532 +21661 +0.014 +fpr +0.012 +0.015 ++ +{gsm} - thresh: -2 +AE rate: fix +0.10 +19851 +1278 +532 +21661 +0.059 +tpr +0.056 +0.062 ++ +{gsm} - thresh: -2 +AE rate: fix +0.25 +15848 +5282 +531 +21661 +0.244 +tpr +0.238 +0.250 ++ +{gsm} - thresh: -2 +AE rate: fix +0.50 +9154 +11987 +520 +21661 +0.553 +tpr +0.547 +0.560 ++ +{gsm} - thresh: -2 +AE rate: fix +0.75 +5214 +15941 +506 +21661 +0.736 +tpr +0.730 +0.742 ++ +{gsm} - thresh: -2 +AE rate: fix +1.00 +2915 +18251 +495 +21661 +0.843 +tpr +0.838 +0.847 ++ +{gsm} - thresh: -2 +AE rate: flexible +0.00 +20899 +241 +521 +21661 +0.011 +fpr +0.010 +0.013 ++ +{gsm} - thresh: -2 +AE rate: flexible +0.10 +20669 +471 +521 +21661 +0.022 +tpr +0.020 +0.024 ++ +{gsm} - thresh: -2 +AE rate: flexible +0.25 +19742 +1399 +520 +21661 +0.065 +tpr +0.061 +0.068 ++ +{gsm} - thresh: -2 +AE rate: flexible +0.50 +16448 +4698 +515 +21661 +0.217 +tpr +0.211 +0.222 ++ +{gsm} - thresh: -2 +AE rate: flexible +0.75 +12320 +8830 +511 +21661 +0.408 +tpr +0.401 +0.414 ++ +{gsm} - thresh: -2 +AE rate: flexible +1.00 +8727 +12437 +497 +21661 +0.574 +tpr +0.568 +0.581 ++ +{simaerep} - thresh: 0.95 +AE rate: fix +0.00 +21367 +294 +0 +21661 +0.014 +fpr +0.012 +0.015 ++ +{simaerep} - thresh: 0.95 +AE rate: fix +0.10 +20366 +1295 +0 +21661 +0.060 +tpr +0.057 +0.063 ++ +{simaerep} - thresh: 0.95 +AE rate: fix +0.25 +17155 +4506 +0 +21661 +0.208 +tpr +0.203 +0.214 ++ +{simaerep} - thresh: 0.95 +AE rate: fix +0.50 +10635 +11026 +0 +21661 +0.509 +tpr +0.502 +0.516 ++ +{simaerep} - thresh: 0.95 +AE rate: fix +0.75 +5845 +15816 +0 +21661 +0.730 +tpr +0.724 +0.736 ++ +{simaerep} - thresh: 0.95 +AE rate: fix +1.00 +4053 +17608 +0 +21661 +0.813 +tpr +0.808 +0.818 ++ +{simaerep} - thresh: 0.95 +AE rate: flexible +0.00 +21340 +321 +0 +21661 +0.015 +fpr +0.013 +0.017 ++ +{simaerep} - thresh: 0.95 +AE rate: flexible +0.10 +20234 +1427 +0 +21661 +0.066 +tpr +0.063 +0.069 ++ +{simaerep} - thresh: 0.95 +AE rate: flexible +0.25 +16661 +5000 +0 +21661 +0.231 +tpr +0.225 +0.237 ++ +{simaerep} - thresh: 0.95 +AE rate: flexible +0.50 +9878 +11783 +0 +21661 +0.544 +tpr +0.537 +0.551 ++ +{simaerep} - thresh: 0.95 +AE rate: flexible +0.75 +5216 +16445 +0 +21661 +0.759 +tpr +0.753 +0.765 ++ + +{simaerep} - thresh: 0.95 +AE rate: flexible +1.00 +3595 +18066 +0 +21661 +0.834 +tpr +0.829 +0.839 +++Plot Performance Metrics +
++
+- {gsm} has better performance than {simaerep} when the AE rate is +fixed, while {simaerep} greatly outperforms {gsm} when the AE rate is +flexible and mimics the AE rates encountered in real study data +sets.
+++plot_perf <- function(df_perf) { + + df_perf %>% + mutate(ur_rate = paste0("under-reporting rate: ", ur_rate, " - ", ratio_type), + ur_rate = ifelse(str_detect(ur_rate, "fpr"), "fpr", ur_rate)) %>% + group_by(ur_rate) %>% + ggplot(aes(type, ratio)) + + geom_errorbar(aes(ymin = ci95_low, ymax = ci95_high, color = type), linewidth = 1) + + facet_grid(ur_rate ~ ae_rate) + + coord_flip() + + theme(legend.position = "bottom") + + labs( + x = "", + y = "CI95 Performance Ratio", + title = "{simaerep} vs {gsm} Performance" + ) + + scale_color_manual(values = c("#5491CC", "#F46626")) +} + +plot_perf(df_perf)+
++plan(sequential)
