write and test binprop_ci(), not integrated into calcualte_2d_metric(…

…) yet
lgessl · Jan 22, 2024 · 5ff2411 · 5ff2411
1 parent 6dbb6a5
commit 5ff2411
Show file tree

Hide file tree

Showing 3 changed files with 86 additions and 13 deletions.
diff --git a/R/calculate_metric.R b/R/calculate_metric.R
@@ -40,6 +40,11 @@ calculate_2d_metric <- function(
                     y_metric = perf_plot_spec$y_metric
                 )
             }
+            names(tbl) <- c(
+                perf_plot_spec$x_metric, 
+                perf_plot_spec$y_metric, 
+                "cutoff"
+            )
             tbl[["split"]] <- i
             tbl[["model"]] <- estimate_name       
             tbl_list <- c(tbl_list, list(tbl))
@@ -55,6 +60,50 @@ calculate_2d_metric <- function(
 }
 
 
+binprop_ci <- function(
+    estimate,
+    actual,
+    confidence_level = 0.95,
+    y_metric = "ci_boundary",
+    x_metric = "prevalence",
+    lower_boundary = TRUE
+){
+    estimate_actual <- intersect_by_names(
+        estimate, 
+        actual, 
+        rm_na = TRUE
+    )
+    estimate <- estimate_actual[[1]]
+    actual <- estimate_actual[[2]]
+    cutoffs <- estimate[estimate > min(estimate)] |> unique() |> sort()
+    prevalence <- numeric(length(cutoffs))
+    ci_boundary <- numeric(length(cutoffs))
+
+    binprop_ci_core <- function(
+        cutoff
+    ){
+        positive <- ifelse(estimate >= cutoff, 1, 0)
+        prevalence <- mean(positive)
+        htest <- stats::binom.test(
+            x = sum(positive * actual),
+            n = sum(positive),
+            conf.level = confidence_level
+        )
+        ci_boundary <- ifelse(lower_boundary, htest$conf.int[1], htest$conf.int[2])
+        c("prevalence" = prevalence, "ci_boundary" = ci_boundary)
+    }
+    mat <- sapply(cutoffs, binprop_ci_core)
+
+    # Store them in a tibble
+    tbl <- tibble::tibble(
+        mat["prevalence", ],
+        mat["ci_boundary", ],
+        cutoffs
+    )
+    return(tbl)
+}
+
+
 metric_with_rocr <- function(
     estimate,
     actual,
@@ -82,11 +131,6 @@ metric_with_rocr <- function(
         rocr_perf@y.values[[1]],
         rocr_perf@alpha.values[[1]]
     )
-    names(tbl) <- c(
-        x_metric, 
-        y_metric,
-        "cutoff"
-    )
     return(tbl)
 }
 
@@ -134,10 +178,5 @@ logrank_metric <- function(
         logrank_p,
         cutoffs
     )
-    names(tbl) <- c(
-        x_metric, 
-        y_metric,
-        "cutoff"
-    )
     return(tbl)
 }
diff --git a/tests/testthat/test-calculate_metric.R b/tests/testthat/test-calculate_metric.R
@@ -63,6 +63,40 @@ test_that("calculate_2d_metric() works", {
 })
 
 
+test_that("binprop_ci() works", {
+
+  set.seed(143)
+
+  n_samples <- 10
+  n_extra <- 2
+  gamma <- .85
+
+  estimate <- rnorm(n_samples)
+  actual <- sample(c(0, 1), n_samples + n_extra, replace = TRUE)
+  names(estimate) <- paste0("s", sample(n_samples))
+  names(actual) <- paste0("s", sample(length(actual)))
+  actual[1] <- NA
+
+  tbl_low <- binprop_ci(
+    estimate = estimate,
+    actual = actual,
+    confidence_level = gamma,
+    y_metric = "ci_boundary",
+    x_metric = "prevalence",
+    lower_boundary = TRUE
+  )
+  tbl_high <- binprop_ci(
+    estimate = estimate,
+    actual = actual,
+    confidence_level = gamma,
+    y_metric = "ci_boundary",
+    x_metric = "prevalence",
+    lower_boundary = FALSE
+  )
+  expect_true(all(tbl_low[, 1] <= tbl_high[, 1]))
+})
+
+
 test_that("metric_with_rocr() works", {
 
   set.seed(354)
@@ -83,7 +117,7 @@ test_that("metric_with_rocr() works", {
     y_metric = "prec"
   )
   expect_s3_class(tbl, "tbl_df")
-  expect_equal(names(tbl), c("rpp", "prec", "cutoff"))
+  expect_equal(ncol(tbl), 3)
 })
 
 
@@ -118,6 +152,6 @@ test_that("logrank_metric() works", {
     x_metric = "my_x"
   )
   expect_s3_class(tbl, "tbl_df")
-  expect_equal(names(tbl), c("my_x", "my_y", "cutoff"))
+  expect_equal(ncol(tbl), 3)
   expect_equal(nrow(tbl), length(estimate) - 1 - n_na_in_estimate)
 })
diff --git a/tests/testthat/test-plot.R b/tests/testthat/test-plot.R
@@ -12,7 +12,7 @@ test_that("plot_2d_metric() works", {
     y_metric = "prec",
     pivot_time_cutoff = 2.,
     benchmark = "bm",
-    show_plots = TRUE,
+    show_plots = FALSE,
     title = "this title",
     x_lab = "this x lab",
     y_lab = "that y lab",