add estimate_secondarry and estimate_truncation simple tests

epiforecasts · seabbs · Oct 14, 2022 · Oct 14, 2022 · Oct 14, 2022 · Oct 14, 2022
commit 7b41bcce96dafcededc14e652e274a00b55c9097
diff --git a/NEWS.md b/NEWS.md
@@ -1,7 +1,6 @@
 # EpiNow2 1.3.3
 
-This maintenance release adds a range of new minor features, squashes bugs, 
-and removes some obsolete features.
+This maintenance release adds a range of new minor features, squashes bugs, slightly expands unit testing, and removes some obsolete features.
 
 Thanks to @Bisaloo, @hsbadr, @LloydChapman, @medewitt, and @sbfnk.
 
@@ -35,6 +34,7 @@ Thanks to @Bisaloo, @hsbadr, @LloydChapman, @medewitt, and @sbfnk.
 * Better test skipping thanks to @Bisaloo.
 * Switched from `cowplot::theme_cowplot()` to `ggplot2::theme_bw()`. This allows the removal of `cowplot` as a dependency as well making plots visuable for users saving as pngs and using a dark theme. By @seabbs.
 * By default `epinow` and downstream functions remove leading zeros. Now this is optional with the new `filter_leading_zeros` option. Thanks to @LloydChapman in #285.
+* Basic tests have been added to cover `estimate_secondary()`, `forecast_secondary()`, and `estimate_truncation()`. By @seabbs in 
 
 ## Other changes
 

diff --git a/R/estimate_truncation.R b/R/estimate_truncation.R
@@ -184,7 +184,13 @@ estimate_truncation <- function(obs, max_truncation = 10,
     estimates <- recon_obs[dataset == index][, c("id", "dataset") := NULL]
     estimates <- estimates[, lapply(.SD, as.integer)]
     estimates <- estimates[, index := .N - 0:(.N - 1)]
-    estimates[, c("n_eff", "Rhat") := NULL]
+    if (!is.null(estimates$n_eff)) {
+      estimates[, c("n_eff") := NULL]
+    }
+    if (!is.null(estimates$Rhat)) {
+      estimates[, c("Rhat") := NULL]
+    }
+
     target_obs <-
       data.table::merge.data.table(
         target_obs, last_obs,

diff --git a/man/setup_future.Rd b/man/setup_future.Rd
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
@@ -0,0 +1,41 @@
+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
+// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#include <RcppEigen.h>
+#include <Rcpp.h>
+
+using namespace Rcpp;
+
+#ifdef RCPP_USE_GLOBAL_ROSTREAM
+Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
+Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
+#endif
+
+
+RcppExport SEXP _rcpp_module_boot_stan_fit4estimate_infections_mod();
+RcppExport SEXP _rcpp_module_boot_stan_fit4estimate_secondary_mod();
+RcppExport SEXP _rcpp_module_boot_stan_fit4estimate_truncation_mod();
+RcppExport SEXP _rcpp_module_boot_stan_fit4exp_mod();
+RcppExport SEXP _rcpp_module_boot_stan_fit4gamma_mod();
+RcppExport SEXP _rcpp_module_boot_stan_fit4lnorm_mod();
+RcppExport SEXP _rcpp_module_boot_stan_fit4simulate_infections_mod();
+RcppExport SEXP _rcpp_module_boot_stan_fit4simulate_secondary_mod();
+RcppExport SEXP _rcpp_module_boot_stan_fit4tune_inv_gamma_mod();
+
+static const R_CallMethodDef CallEntries[] = {
+    {"_rcpp_module_boot_stan_fit4estimate_infections_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4estimate_infections_mod, 0},
+    {"_rcpp_module_boot_stan_fit4estimate_secondary_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4estimate_secondary_mod, 0},
+    {"_rcpp_module_boot_stan_fit4estimate_truncation_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4estimate_truncation_mod, 0},
+    {"_rcpp_module_boot_stan_fit4exp_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4exp_mod, 0},
+    {"_rcpp_module_boot_stan_fit4gamma_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4gamma_mod, 0},
+    {"_rcpp_module_boot_stan_fit4lnorm_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4lnorm_mod, 0},
+    {"_rcpp_module_boot_stan_fit4simulate_infections_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4simulate_infections_mod, 0},
+    {"_rcpp_module_boot_stan_fit4simulate_secondary_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4simulate_secondary_mod, 0},
+    {"_rcpp_module_boot_stan_fit4tune_inv_gamma_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4tune_inv_gamma_mod, 0},
+    {NULL, NULL, 0}
+};
+
+RcppExport void R_init_EpiNow2(DllInfo *dll) {
+    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+    R_useDynamicSymbols(dll, FALSE);
+}
diff --git a/tests/testthat/test-estimate_secondary.R b/tests/testthat/test-estimate_secondary.R
@@ -0,0 +1,62 @@
+ skip_on_cran()
+ library(data.table)
+
+ # make some example secondary incidence data
+ cases <- example_confirmed
+ cases <- as.data.table(cases)[, primary := confirm]
+ # Assume that only 40 percent of cases are reported
+ cases[, scaling := 0.4]
+ # Parameters of the assumed log normal delay distribution
+ cases[, meanlog := 1.8][, sdlog := 0.5]
+
+ # apply a convolution of a log normal to a vector of observations
+ weight_cmf <- function(x, ...) {
+   set.seed(x[1])
+   meanlog <- rnorm(1, 1.6, 0.2)
+   sdlog <- rnorm(1, 0.8, 0.1)
+   cmf <- cumsum(dlnorm(1:length(x), meanlog, sdlog)) -
+     cumsum(dlnorm(0:(length(x) - 1), meanlog, sdlog))
+   cmf <- cmf / plnorm(length(x), meanlog, sdlog)
+   conv <- sum(x * rev(cmf), na.rm = TRUE)
+   conv <- round(conv, 0)
+   return(conv)
+ }
+ # roll over observed cases to produce a convolution
+ cases <- cases[, .(date, primary = confirm, secondary = confirm)]
+ cases <- cases[, secondary := frollapply(secondary, 15, weight_cmf, align = "right")]
+ cases <- cases[!is.na(secondary)]
+ # add a day of the week effect and scale secondary observations at 40\% of primary
+ cases <- cases[lubridate::wday(date) == 1, secondary := round(0.5 * secondary, 0)]
+ cases <- cases[, secondary := round(secondary * rnorm(.N, 0.4, 0.025), 0)]
+ cases <- cases[secondary < 0, secondary := 0]
+ cases <- cases[, secondary := map_dbl(secondary, ~ rpois(1, .))]
+
+ # fit model to example data assuming only a given fraction of primary observations
+ # become secondary observations
+ inc <- estimate_secondary(cases[1:60],
+   obs = obs_opts(scale = list(mean = 0.2, sd = 0.2)),
+   chains = 2, warmup = 250, iter = 750, cores = 2,
+   refresh = 0
+ )
+
+test_that("estimate_secondary can return values from simulated data and plot
+           them", {
+  expect_equal(names(inc), c("predictions", "data", "fit"))
+  expect_equal(
+    names(inc$predictions),
+    c("date", "primary", "secondary", "mean", "se_mean", "sd", "lower_90",
+      "lower_50", "lower_20", "median", "upper_20", "upper_50", "upper_90"
+    )
+  )
+  expect_true(is.list(inc$data))
+  # validation plot of observations vs estimates
+  expect_error(plot(inc, primary = TRUE), NA)
+})
+
+test_that("forecast_secondary can return values from simulated data and plot
+           them", {
+  inc_preds <- forecast_secondary(inc, cases[61:.N][, value := primary])
+  expect_equal(names(inc_preds), c("samples", "forecast", "predictions"))
+  # validation plot of observations vs estimates
+  expect_error(plot(inc_preds, new_obs = cases, from = "2020-05-01"), NA)
+})
diff --git a/tests/testthat/test-estimate_truncation.R b/tests/testthat/test-estimate_truncation.R
@@ -0,0 +1,57 @@
+# Setup for testing -------------------------------------------------------
+skip_on_cran()
+futile.logger::flog.threshold("FATAL")
+
+ # set number of cores to use
+ options(mc.cores = ifelse(interactive(), 4, 1))
+ # get example case counts
+ reported_cases <- example_confirmed[1:60]
+
+ # define example truncation distribution (note not integer adjusted)
+ trunc_dist <- list(
+   mean = convert_to_logmean(3, 2),
+   mean_sd = 0.1,
+   sd = convert_to_logsd(3, 2),
+   sd_sd = 0.1,
+   max = 10
+ )
+
+ # apply truncation to example data
+ construct_truncation <- function(index, cases, dist) {
+   set.seed(index)
+   cmf <- cumsum(
+     dlnorm(
+       1:(dist$max + 1),
+       rnorm(1, dist$mean, dist$mean_sd),
+       rnorm(1, dist$sd, dist$sd_sd)
+     )
+   )
+   cmf <- cmf / cmf[dist$max + 1]
+   cmf <- rev(cmf)[-1]
+   trunc_cases <- data.table::copy(cases)[1:(.N - index)]
+   trunc_cases[(.N - length(cmf) + 1):.N, confirm := as.integer(confirm * cmf)]
+   return(trunc_cases)
+ }
+ example_data <- purrr::map(c(20, 15, 10, 0),
+   construct_truncation,
+   cases = reported_cases,
+   dist = trunc_dist
+ )
+
+test_that("estimate_truncation can return values from simulated data and plot
+           them", {
+  # fit model to example data
+  est <- estimate_truncation(example_data,
+    verbose = interactive(), refresh = 0,
+    chains = 2, iter = 1000, warmup = 250
+  )
+  expect_equal(
+    names(est),
+    c("dist", "obs", "last_obs", "cmf", "data", "fit")
+  )
+  expect_equal(
+    names(est$dist),
+    c("mean", "mean_sd", "sd", "sd_sd", "max")
+  )
+  expect_error(plot(est), NA)
+})