Add check and warning for PMFs longer than the data

NEWS.md

@@ -1,5 +1,9 @@
 # EpiNow2 (development version)
 
+## Documentation
+
+- If users supply PMFs that are longer than the data, they are now informed that this will be trimmed to match the length of the data. By @jamesmbaazam in # and reviewed by <REVIEWER>.
+
 # EpiNow2 1.7.1
 
 This is a patch release in response to an upstream issue in `rstan`, as flagged in CRAN checks.

R/checks.R

@@ -109,7 +109,7 @@
   }
   assert_numeric(attr(dist, "cdf_cutoff"), lower = 0, upper = 1)
   # Check that `dist` has a finite maximum
   if (any(is.infinite(max(dist))) && !(attr(dist, "cdf_cutoff") > 0)) {
     cli_abort(
       c(
         "i" = "All distributions passed to the model need to have a
@@ -180,3 +180,54 @@
     )
   }
 }
+
+
+#' Check that supplied PMFs are not longer than the data
+#'
+#' @param ... Delay distributions
+#' @inheritParams estimate_infections
+#' @importFrom cli cli_warn col_red
+#'
+#' @returns Called for its side effects
+#' @keywords internal
+check_pmf_length_against_data <- function(..., data) {
+  delays <- list(...)
+  flat_delays <- do.call(c, delays)
+  # Track which component each delay came from
+  delay_names <- rep(names(delays), vapply(delays, EpiNow2:::ndist, numeric(1)))
+  names(flat_delays) <- delay_names
+  # Find the non-parametric distributions
+  np_delays <- which(unname(vapply(
+    flat_delays, function(x) {
+      get_distribution(x) == "nonparametric"
+    }, logical(1)
+  )))
+
+  if (length(np_delays) == 0) return(invisible())
+
+  # Check lengths and collect info about exceeding PMFs
+  pmf_longer_than_data <- vapply(flat_delays[np_delays], function(x) {
+    length(x$pmf) > nrow(data)
+  }, logical(1))
+
+  if (any(pmf_longer_than_data)) {
+    # Get details for each long PMF
+    long_pmf_lengths <- vapply(
+      flat_delays[np_delays][pmf_longer_than_data], function(x) {
+        length(x$pmf)
+      }, numeric(1)
+    )
+  }
+
+    cli::cli_warn(
+      c(
+        "!" = "You have supplied PMFs that are longer than the data. ",
+        "{names(long_pmf_lengths)} {?has/have} length{?s}
+        {.val {long_pmf_lengths}} but data has
+        {.val {nrow(data)}} rows.",
+        "i" = "{cli::col_red('These will be trimmed to match the rows in the
+        data. To remove this message, make sure the PMFs have the same length
+        as the data')}"
+      )
+    )
+}

R/estimate_infections.R

@@ -190,6 +190,14 @@
   # store dirty reported case data
   dirty_reported_cases <- data.table::copy(data)
 
+  # Check that no PMF is longer than the data
+  check_pmf_length_against_data(
+    generation_time = generation_time,
+    truncation = truncation,
+    delays = delays,
+    data = data
+  )
+
   if (!is.null(rt) && !rt$use_rt) {
     rt <- NULL
   }
@@ -265,7 +273,7 @@
   ))
 
   # Set up default settings
   args <- create_stan_args(
     stan = stan,
     data = stan_data,
     init = create_initial_conditions(stan_data),
@@ -378,5 +386,5 @@
     order_by = c("variable", "date"),
     CrIs = CrIs
   )
   return(format_out)
 }

R/estimate_secondary.R

@@ -187,6 +187,13 @@
   assert_logical(weigh_delay_priors)
   assert_logical(verbose)
 
+  # Check that no PMF is longer than the data
+  check_pmf_length_against_data(
+    truncation = truncation,
+    delays = delays,
+    data = data
+  )
+
   reports <- data.table::as.data.table(data)
 
   reports <- default_fill_missing_obs(reports, obs, "secondary")
@@ -260,7 +267,7 @@
     c(stan_data, list(estimate_r = 0, fixed = 1, bp_n = 0))
   )
   # fit
   args <- create_stan_args(
     stan = stan, data = stan_data, init = inits, model = "estimate_secondary"
   )
   fit <- fit_model(args, id = "estimate_secondary")
@@ -320,7 +327,7 @@
       )
     }
     # replace scaling if present in the prior
     scale <- priors[grepl("frac_obs", variable, fixed = TRUE)]
     if (nrow(scale) > 0) {
       data$obs_scale_mean <- as.array(signif(scale$mean, 3))
       data$obs_scale_sd <- as.array(signif(scale$sd, 3))
@@ -398,14 +405,14 @@
     predictions <- predictions[date <= to]
   }
 
   plot <- ggplot2::ggplot(predictions, ggplot2::aes(x = date, y = secondary)) +
     ggplot2::geom_col(
       fill = "grey", col = "white",
       show.legend = FALSE, na.rm = TRUE
     )
 
   if (primary) {
     plot <- plot +
       ggplot2::geom_point(
         data = predictions,
         ggplot2::aes(y = primary),
@@ -416,7 +423,7 @@
         ggplot2::aes(y = primary), alpha = 0.4
       )
   }
   plot <- plot_CrIs(plot, extract_CrIs(predictions),
     alpha = 0.6, linewidth = 1
   )
   plot <- plot +