25 remove hardcoded intervals

R/summarize_rtestimate.R

@@ -6,6 +6,7 @@
 #' @param median Double vector. vector of median values.
 #' @param lb Double vector. vector of lower bounds.
 #' @param ub Double vector. vector of upper bounds.
+#' @param level Double scalar. the confidence level associated with `lb`/`ub`
 #' @param package String. Name of the package.
 #' @param notes String. Notes about the summary.
 #' @export
@@ -15,17 +16,19 @@
 #'   - `median`: Double vector. vector of median values.
 #'   - `lb`: Double vector. vector of lower bounds.
 #'   - `ub`: Double vector. vector of upper bounds.
+#' - `level`: Double scalar.
 #' - `package`: String. Name of the package.
 #' - `notes`: String. Notes about the summary.
 new_summrt <- function(
-    date, median, lb, ub, package, notes
+    date, median, lb, ub, level, package, notes
 ) {
 
   # Asserting the types
   checkmate::assert_integer(date)
   checkmate::assert_double(median)
   checkmate::assert_double(lb)
   checkmate::assert_double(ub)
+  checkmate::assert_number(level, lower = 0, upper = 1)
   checkmate::assert_string(package)
   checkmate::assert_string(notes)
 
@@ -46,6 +49,7 @@ new_summrt <- function(
         lb = lb,
         ub = ub
       ),
+      level = level,
       package = package,
       notes = notes
     ),
@@ -59,25 +63,29 @@ new_summrt <- function(
 #' @param ... Additional arguments passed to methods.
 print.summrt_summary <- function(x, ...) {
   cat("Summary of Rt estimation\n")
+  cat("Confidence level : ", x$level, "\n")
   cat("Package : ", x$package, "\n")
-  cat("Notes   : ", x$notes, "\n")
+  if (nchar(x$notes) > 0L) cat("Notes   : ", x$notes, "\n")
   print(x$estimates)
   invisible(x)
 }
 
 #' Extract Rt estimation from a model fit
 #' @param x Object to extract Rt from.
+#' @param level Desired confidence level for the estimate
 #' @param ... Additional arguments passed to methods.
 #' @param notes String. Optional notes to add to the summary.
 #' @export
-summarize_rtestimate <- function(x, ..., notes = "") {
+summarize_rtestimate <- function(x, level = 0.95, ..., notes = "") {
+  checkmate::assert_number(level, lower = 0, upper = 1)
+  checkmate::assert_string(notes)
   UseMethod("summarize_rtestimate")
 }
 
 #' @rdname summarize_rtestimate
 #' @importFrom cli cli_abort
 #' @export
-summarize_rtestimate.default <- function(x, ..., notes = "") {
+summarize_rtestimate.default <- function(x, level = 0.95, ..., notes = "") {
   cli::cli_abort("Your Rt method isn't supported yet. You should create a method.")
 }
 
@@ -100,14 +108,14 @@ summarize_rtestimate.cv_poisson_rt <- function(
   } else {
     checkmate::assert_number(lambda, lower = 0)
   }
-  checkmate::assert_number(level, lower = 0, upper = 1)
   cb <- rtestim::confband(x, lambda = lambda, level = level, ...)
 
   new_summrt(
     date = as.integer(x$full_fit$x),
     median = cb$fit,
     lb = cb[[2]], # danger
     ub = cb[[3]],
+    level = level,
     package = "rtestim",
     notes = notes
   )
@@ -126,14 +134,14 @@ summarize_rtestimate.poisson_rt <- function(x, level = 0.95, lambda = NULL, ...,
     lambda <- 10^stats::median(log10(x$lambda))
   }
   checkmate::assert_number(lambda, lower = 0)
-  checkmate::assert_number(level, lower = 0, upper = 1)
   cb <- rtestim::confband(x, lambda = lambda, level = level, ...)
 
   new_summrt(
     date = as.integer(x$x),
     median = cb$fit,
     lb = cb[[2]],
     ub = cb[[3]],
+    level = level,
     package = "rtestim",
     notes = notes
   )
@@ -154,30 +162,45 @@ summarize_rtestimate.epinow <- function(x, level = 0.95, ..., notes = "") {
   t_length <- as.integer(t_max - t_min)
   t_length_forecast <- ncol(rstan::extract(x$estimates$fit)$R) - nrow(x$estimates$observations)
 
+  probs <- level_to_probs(level)
+
   return(new_summrt(
     date = c(0:t_length, (t_length + 1):(t_length + t_length_forecast)),
     median = apply(y_extract, 2, stats::quantile, probs = 0.5),
-    lb = apply(y_extract, 2, stats::quantile, probs = 0.025),
-    ub = apply(y_extract, 2, stats::quantile, probs = 0.975),
+    lb = apply(y_extract, 2, stats::quantile, probs = probs$lower),
+    ub = apply(y_extract, 2, stats::quantile, probs = probs$upper),
+    level = level,
     package = "EpiNow2",
     notes = notes
   ))
 
 }
 
 #' @export
-#' @details The `estimate_R` method is for the `EpiEstim` package.
+#' @details The `estimate_R` method is for the `EpiEstim` package. Currently,
+#'   only levels in 50%, 90% and 95% confidence levels are allowed.
 #' @rdname summarize_rtestimate
-summarize_rtestimate.estimate_R <- function(x, ..., notes = "") {
+summarize_rtestimate.estimate_R <- function(x, level = 0.95, ..., notes = "") {
   if (!requireNamespace("EpiEstim", quietly = TRUE)) {
     cli::cli_abort("You must install the {.pkg EpiEstim} package for this functionality.")
   }
-
+  allowed_levels <- c(0.95, 0.9, 0.5)
+  if (min(abs(level - allowed_levels)) > sqrt(.Machine$double.eps)) {
+    cli::cli_abort(paste(
+      "For {.pkg EpiEstim}, allowable confidence levels are",
+      "{.val {allowed_levels}}, not {.val {level}}."
+    ))
+  }
+  probs <- level_to_probs(level)
+  probs <- lapply(probs, probs_to_char)
+  lb_name <- paste0("Quantile.", probs$lower, "(R)")
+  ub_name <- paste0("Quantile.", probs$upper, "(R)")
   new_summrt(
     date    = as.integer((x$R$t_end + x$R$t_start) / 2),
-    median  = x$R$`Median(R)`,
-    lb      = x$R$`Quantile.0.025(R)`,
-    ub      = x$R$`Quantile.0.975(R)`,
+    median  = x$R[["Median(R)"]],
+    lb      = x$R[[lb_name]],
+    ub      = x$R[[ub_name]],
+    level   = level,
     package = "EpiEstim",
     notes   = notes
   )
@@ -186,16 +209,29 @@ summarize_rtestimate.estimate_R <- function(x, ..., notes = "") {
 #' @export
 #' @details The `Rt` method is for the `EpiLPS` package.
 #' @rdname summarize_rtestimate
-summarize_rtestimate.Rt <- function(x, ..., notes = "") {
+summarize_rtestimate.Rt <- function(x, level = 0.95, ..., notes = "") {
   if (!requireNamespace("EpiLPS", quietly = TRUE)) {
     cli::cli_abort("You must install the {.pkg EpiLPS} package for this functionality.")
   }
 
+  allowed_levels <- c(0.95, 0.9, 0.5)
+  if (min(abs(level - allowed_levels)) > sqrt(.Machine$double.eps)) {
+    cli::cli_abort(paste(
+      "For {.pkg EpiLPS}, allowable confidence levels are",
+      "{.val {allowed_levels}}, not {.val {level}}."
+    ))
+  }
+  probs <- level_to_probs(level)
+  probs <- lapply(probs, probs_to_char)
+  lb_name <- paste0("Rq", probs$lower)
+  ub_name <- paste0("Rq", probs$upper)
+
   new_summrt(
     date    = x$RLPS$Time,
     median  = x$RLPS$Rq0.50,
-    lb      = x$RLPS$Rq0.025,
-    ub      = x$RLPS$Rq0.975,
+    lb      = x$RLPS[[lb_name]],
+    ub      = x$RLPS[[ub_name]],
+    level   = level,
     package = "EpiLPS",
     notes   = notes
   )

R/utils.R

@@ -0,0 +1,19 @@
+level_to_probs <- function(level) {
+  checkmate::assert_number(level, lower = 0, upper = 1)
+  lower <- (1 - level) / 2
+  upper <- 1 - lower
+  list(lower = lower, upper = upper)
+}
+
+probs_to_char <- function(prob, digits = 3L,
+                          remove_trailing_zero = TRUE,
+                          remove_leading_zero = FALSE) {
+  checkmate::assert_number(prob, lower = 0, upper = 1)
+  checkmate::assert_int(digits, lower = 0)
+  checkmate::assert_logical(remove_trailing_zero, len = 1L)
+  checkmate::assert_logical(remove_leading_zero, len = 1L)
+  .prob <- sprintf(paste0("%.", digits, "f"), prob)
+  if (remove_trailing_zero) .prob <- sub("(0+)$", "", .prob)
+  if (remove_leading_zero) .prob <- sub("^0", "", .prob)
+  .prob
+}

inst/tinytest/test_EpiEstim.R

@@ -1,4 +1,4 @@
-# Example data 
+# Example data
 EpiEstim_obj <- readRDS(
   system.file("extdata", "EpiEstim_example.rds",
               package = "summrt")
@@ -7,9 +7,9 @@ EpiEstim_obj <- readRDS(
 std_epiestim <- summarize_rtestimate(EpiEstim_obj)
 
 message("Check that names of EpiEstim are correct")
-checkmate::expect_names( names(std_epiestim), 
-                         must.include = c("estimates", 
-                                          "package", 
+checkmate::expect_names( names(std_epiestim),
+                         must.include = c("estimates",
+                                          "package",
                                           "notes"))
 
 message("Check that the date column is actually an integer")
@@ -18,9 +18,15 @@ expect_true(is.integer(std_epiestim$estimates$date))
 message("Check that the package name is correct for EpiEstim")
 expect_equal(std_epiestim$package, "EpiEstim")
 
+expect_error(summarize_rtestimate(std_epiestim, level = 0.8))
+
 message("Check that there are no NAs in median, lbs, ubs")
 expect_true(all(!is.na(std_epiestim$estimates$median)))
 expect_true(all(!is.na(std_epiestim$estimates$lb)))
 expect_true(all(!is.na(std_epiestim$estimates$ub)))
 
-
+message("Check that different levels can be used")
+expect_error(summarize_rtestimate(EpiEstim_obj, level = 0.8))
+expect_silent(summarize_rtestimate(EpiEstim_obj))
+expect_silent(summarize_rtestimate(EpiEstim_obj, level = 0.5))
+expect_silent(summarize_rtestimate(EpiEstim_obj, level = 0.9))

inst/tinytest/test_EpiLPS.R

@@ -1,4 +1,4 @@
-# Example data 
+# Example data
 EpiLPS_obj <- readRDS(
   system.file("extdata", "EpiLPS_example.rds",
               package = "summrt")
@@ -7,9 +7,9 @@ EpiLPS_obj <- readRDS(
 std_epilps <- summarize_rtestimate(EpiLPS_obj)
 
 message("Check that names of EpiLPS are correct")
-checkmate::expect_names( names(std_epilps), 
-                         must.include = c("estimates", 
-                                          "package", 
+checkmate::expect_names( names(std_epilps),
+                         must.include = c("estimates",
+                                          "package",
                                           "notes"))
 
 message("Check that the date column is actually an integer")
@@ -21,4 +21,10 @@ expect_equal(std_epilps$package, "EpiLPS")
 message("Check that there are no NAs in median, lbs, ubs")
 expect_true(all(!is.na(std_epilps$estimates$median)))
 expect_true(all(!is.na(std_epilps$estimates$lb)))
-expect_true(all(!is.na(std_epilps$estimates$ub)))
+expect_true(all(!is.na(std_epilps$estimates$ub)))
+
+message("Check that different levels can be used")
+expect_error(summarize_rtestimate(EpiLPS_obj, level = 0.8))
+expect_silent(summarize_rtestimate(EpiLPS_obj))
+expect_silent(summarize_rtestimate(EpiLPS_obj, level = 0.5))
+expect_silent(summarize_rtestimate(EpiLPS_obj, level = 0.9))

inst/tinytest/test_utils.R

@@ -0,0 +1,19 @@
+message("level_to_probs works")
+expect_error(summrt:::level_to_probs(c(.1, .9)))
+expect_error(summrt:::level_to_probs(2))
+expect_equal(summrt:::level_to_probs(.9)$lower, .05)
+expect_equal(summrt:::level_to_probs(.9)$upper, .95)
+
+message("probs_to_char works")
+expect_error(summrt:::probs_to_char(c(.1, .9)))
+expect_error(summrt:::probs_to_char(2))
+expect_error(summrt:::probs_to_char(.5, digits = -1))
+expect_error(summrt:::probs_to_char(.5, remove_trailing_zero = 2))
+expect_error(summrt:::probs_to_char(.5, remove_trailing_zero = c(TRUE, FALSE)))
+expect_error(summrt:::probs_to_char(.5, remove_leading_zero = 2))
+expect_error(summrt:::probs_to_char(.5, remove_leading_zero = c(TRUE, FALSE)))
+expect_identical(summrt:::probs_to_char(0.05), "0.05")
+expect_identical(summrt:::probs_to_char(0.05, digits = 4L), "0.05")
+expect_identical(summrt:::probs_to_char(0.005, digits = 4L), "0.005")
+
+