Robust fxns (#121)

* Add assertions and clean up comments

* Update ordinal test vector names

* Write tests for ordinal_tests fxn

* Use tryCatches on ordinal tests

* Use tryCatch on calculating T1 error

* Minor formatting

R/assign_groups.R

@@ -19,6 +19,23 @@
 #'
 assign_groups <- function(sample_size, sample_prob, prob0, prob1, seed,
                           .rng_kind = NULL, .rng_normal_kind = NULL, .rng_sample_kind = NULL) {
+
+  assertthat::assert_that(
+    length(prob0) == length(prob1),
+    msg = "prob0 and prob1 must have the same length"
+  )
+
+  assertthat::assert_that(
+    near(sum(prob0), 1),
+    msg = "prob0 must sum to 1"
+  )
+
+  assertthat::assert_that(
+    near(sum(prob1), 1),
+    msg = "prob1 must sum to 1"
+  )
+
+
   withr::with_seed(seed,
     {
       y <- factor(sample(x = 0:1, size = sample_size, replace = TRUE, prob = sample_prob))
@@ -29,13 +46,6 @@ assign_groups <- function(sample_size, sample_prob, prob0, prob1, seed,
 
       x[y == 0] <- sample(1:K, n_null, replace = TRUE, prob = prob0)
       x[y == 1] <- sample(1:K, n_intervene, replace = TRUE, prob = prob1)
-      #
-      #   if( length(unique(x[y==0])) < K ) {
-      #     warning("Not all possible outcomes observed in the null group.")
-      #   }
-      #   if( length(unique(x[y==1])) < K ) {
-      #     warning("Not all possible outcomes observed in the intervention group.")
-      #   }
 
       list(
         y = y, x = x, n_null = n_null, n_intervene = n_intervene, sample_size = sample_size, K = K,

R/mod_plot_distributions.R

@@ -91,7 +91,7 @@ mod_plot_distributions_server <- function(id, p_value_table, n) {
     # !!!plot!!!
     distribution_plot <- reactive({
       p_value_reactive_table() %>%
-        dplyr::select(.data$wilcox:.data$coinasymp, .data$sample_size) %>%
+        dplyr::select(.data$Wilcoxon:.data$`Coin Indep. Test`, .data$sample_size) %>%
         plot_distribution_results(
           outlier_removal = outlier_percent_removal(),
           alpha = p_val_threshold()
@@ -106,7 +106,7 @@ mod_plot_distributions_server <- function(id, p_value_table, n) {
     # !!!statistics!!!
     distribution_statistics <- reactive({
       p_value_reactive_table() %>%
-        select(.data$wilcox:.data$coinasymp, .data$sample_size) %>%
+        select(.data$Wilcoxon:.data$`Coin Indep. Test`, .data$sample_size) %>%
         calculate_power_t2error(
           alpha = p_val_threshold(),
           n = n(),
@@ -141,7 +141,7 @@ mod_plot_distributions_server <- function(id, p_value_table, n) {
     group1_t1_reactive_table <- reactive({
       p_value_table$group1_results() %>%
         bind_rows() %>%
-        dplyr::select(.data$wilcox:.data$coinasymp, .data$sample_size) %>%
+        dplyr::select(.data$Wilcoxon:.data$`Coin Indep. Test`, .data$sample_size) %>%
         group_by(.data$sample_size) %>%
         calculate_t1_error(t1_error_confidence_int = input$t1_error_group1_confidence_int)
     })
@@ -163,7 +163,7 @@ mod_plot_distributions_server <- function(id, p_value_table, n) {
     group2_t1_reactive_table <- reactive({
       p_value_table$group2_results() %>%
         bind_rows() %>%
-        dplyr::select(.data$wilcox:.data$coinasymp, .data$sample_size) %>%
+        dplyr::select(.data$Wilcoxon:.data$`Coin Indep. Test`, .data$sample_size) %>%
         group_by(.data$sample_size) %>%
         calculate_t1_error(t1_error_confidence_int = input$t1_error_group2_confidence_int)
     })

R/mod_stats_calculations.R

@@ -115,7 +115,7 @@ mod_stats_calculations_server <- function(id, probability_data, sample_prob, ite
       # browser()
       comparison_results() %>%
         bind_rows() %>%
-        dplyr::select(.data$sample_size, .data$wilcox:.data$coinasymp) %>%
+        dplyr::select(.data$sample_size, .data$Wilcoxon:.data$`Coin Indep. Test`) %>%
         DT::datatable(options = list(scrollX = TRUE)) %>%
         DT::formatRound(2:7, digits = 5)
     })
@@ -126,7 +126,7 @@ mod_stats_calculations_server <- function(id, probability_data, sample_prob, ite
     output$group1_pvalues <- DT::renderDataTable(
       group1_results() %>%
         bind_rows() %>%
-        dplyr::select(.data$sample_size, .data$wilcox:.data$coinasymp) %>%
+        dplyr::select(.data$sample_size, .data$Wilcoxon:.data$`Coin Indep. Test`) %>%
         DT::datatable(options = list(scrollX = TRUE)) %>%
         DT::formatRound(2:7, digits = 5)
     )
@@ -135,7 +135,7 @@ mod_stats_calculations_server <- function(id, probability_data, sample_prob, ite
     output$group2_pvalues <- DT::renderDataTable(
       group2_results() %>%
         bind_rows() %>%
-        dplyr::select(.data$sample_size, .data$wilcox:.data$coinasymp) %>%
+        dplyr::select(.data$sample_size, .data$Wilcoxon:.data$`Coin Indep. Test`) %>%
         DT::datatable(options = list(scrollX = TRUE)) %>%
         DT::formatRound(2:7, digits = 5)
     )

R/ordinal_tests.R

@@ -21,11 +21,31 @@
 #'
 ordinal_tests <- function(x, y, ...) {
   c(
-    wilcox = stats::wilcox.test(x[y == 0], x[y == 1])[["p.value"]],
-    fisher = stats::fisher.test(x, y, simulate.p.value = FALSE, workspace = 2e7)[["p.value"]],
-    chi_sq_false = stats::chisq.test(x, y, correct = FALSE)[["p.value"]],
-    chi_sq_true = stats::chisq.test(x, y, correct = TRUE)[["p.value"]],
-    lrm = rms::lrm(x ~ y)$stats[["P"]],
-    coinasymp = coin::pvalue(coin::independence_test(x ~ y, ytrafo = coin::rank_trafo))
+    Wilcoxon = tryCatch(
+      {stats::wilcox.test(x[y == 0], x[y == 1])[["p.value"]]},
+      error = function(e) {NA_real_}
+      ),
+    Fisher = tryCatch(
+      {stats::fisher.test(x, y, simulate.p.value = FALSE, workspace = 2e7)[["p.value"]]},
+      error = function(e) {NA_real_}
+    ),
+    "Chi Squared\n(No Correction)" = tryCatch(
+      {stats::chisq.test(x, y, correct = FALSE)[["p.value"]]},
+      error = function(e) {NA_real_}
+    ),
+    "Chi Squared\n(Correction)" = tryCatch(
+      {stats::chisq.test(x, y, correct = TRUE)[["p.value"]]},
+      error = function(e) {NA_real_}
+    ),
+    "Prop. Odds" = tryCatch(
+      {rms::lrm(x ~ y)$stats[["P"]]},
+      error = function(e) {NA_real_}
+    ),
+    "Coin Indep. Test" = tryCatch(
+      {coin::pvalue(coin::independence_test(x ~ y, ytrafo = coin::rank_trafo))},
+      error = function(e) {NA_real_}
+    )
   )
 }
+
+

R/run_simulations.R

@@ -27,10 +27,19 @@
 run_simulations <- function(sample_size, sample_prob, prob0, prob1, niter,
                             .rng_kind = NULL, .rng_normal_kind = NULL, .rng_sample_kind = NULL) {
   # Check equal vector lengths
-  assert_that(length(prob0) == length(prob1))
+  assert_that(
+    length(prob0) == length(prob1),
+    msg = "prob0 and prob1 must have the same length"
+    )
   # Check probabilities for both groups sum to 1
-  assert_that(dplyr::near(sum(prob0), 1), msg = "Probability for Group 1 does not sum to 1.")
-  assert_that(dplyr::near(sum(prob1), 1), msg = "Probability for Group 2 does not sum to 1")
+  assertthat::assert_that(
+    near(sum(prob0), 1),
+    msg = "prob0 must sum to 1"
+  )
+  assertthat::assert_that(
+    near(sum(prob1), 1),
+    msg = "prob0 must sum to 1"
+  )
 
 
   K <- length(prob0)
@@ -44,15 +53,18 @@ run_simulations <- function(sample_size, sample_prob, prob0, prob1, niter,
   purrr::map(sample_size,
     ~ {
       sample_size_nested <- .x
-      initial_groups <- purrr::map(1:niter, ~ assign_groups(
-        sample_size = sample_size_nested,
-        sample_prob = sample_prob,
-        prob0 = prob0, prob1 = prob1,
-        seed = .x,
-        .rng_kind = .rng_kind,
-        .rng_normal_kind = .rng_normal_kind,
-        .rng_sample_kind = .rng_sample_kind
-      ))
+      initial_groups <- purrr::map(1:niter, ~{
+        assign_groups(
+          sample_size = sample_size_nested,
+          sample_prob = sample_prob,
+          prob0 = prob0, prob1 = prob1,
+          seed = .x,
+          .rng_kind = .rng_kind,
+          .rng_normal_kind = .rng_normal_kind,
+          .rng_sample_kind = .rng_sample_kind
+        )
+      }
+      )
 
       p_values <- initial_groups %>%
         sapply(., function(x) ordinal_tests(x[["x"]], x[["y"]])) %>%

R/utils.R

@@ -95,13 +95,30 @@ calculate_t1_error <- function(df, alpha = 0.05, t1_error_confidence_int = 95, n
     dplyr::group_modify(
       ~ {
         purrr::map(.x, ~ {
-          binom_results <- binom.test(sum(.x < alpha), length(.x), conf.level = t1_error_confidence_int / 100)
-          tibble(
-            lower_t1_bound = binom_results$conf.int[[1]],
-            upper_t1_bound = binom_results$conf.int[[2]],
-            t1_error = binom_results$estimate,
-            !!ci_label := paste0("[",round(lower_t1_bound, 4), ", ", round(upper_t1_bound, 4), "]")
-            )
+
+          tryCatch({
+
+            binom_results <- binom.test(sum(.x < alpha, na.rm = T),
+                                        sum(!is.na(.x)),
+                                        conf.level = t1_error_confidence_int / 100)
+
+              tibble(
+                lower_t1_bound = binom_results$conf.int[[1]],
+                upper_t1_bound = binom_results$conf.int[[2]],
+                t1_error = binom_results$estimate,
+                !!ci_label := paste0("[",round(lower_t1_bound, 4), ", ", round(upper_t1_bound, 4), "]")
+              )
+
+              },
+            error = function(e) {
+              tibble(
+                lower_t1_bound = NA_real_,
+                upper_t1_bound = NA_real_,
+                t1_error = NA_real_,
+                !!ci_label := NA_character_
+              )
+            })
+
         }) %>%
           purrr::list_rbind(names_to = "test")
       }

tests/testthat/test-ordinal_tests.R

@@ -0,0 +1,38 @@
+
+# write tests for the ordinal_tests function
+
+sample_size = 30
+prob0 <- c(.4, .3, .3)
+prob1 <- c(.9, .05, .05)
+seed <- 10
+sample_prob <- c(0.50, 0.50)
+
+groups <- assign_groups(sample_size = sample_size,
+              prob0 = prob0,
+              prob1 = prob1,
+              seed = 10,
+              sample_prob = sample_prob)
+
+
+
+test_that("ordinal_tests returns a vector", {
+  suppressWarnings(
+    expect_vector(ordinal_tests(groups[["x"]], groups[["y"]]))
+  )
+})
+
+test_that("ordinal_tests returns a vector of length 3", {
+  suppressWarnings(
+    expect_length(ordinal_tests(groups[["x"]], groups[["y"]]), 6)
+  )
+})
+
+test_that("ordinal_tests names the vector correctly", {
+  suppressWarnings(
+    expect_named(
+      ordinal_tests(groups[["x"]], groups[["y"]]),
+    c("Wilcoxon", "Fisher", "Chi Squared\n(No Correction)", "Chi Squared\n(Correction)", "Prop. Odds", "Coin Indep. Test")
+  )
+  )
+})
+