add error message for missing strata

NEWS.md

@@ -1,7 +1,9 @@
-# serocalculator (development version)
+# serocalculator 1.3.0
 
 ## New features
 
+* Add test for missing strata in `est.incidence.by` (#227)
+
 * Added example datasets with documentation for examples and testing (#314)
 
 * Improved error messaging for `autoplot.pop_data()` (#234).

R/est.incidence.by.R

@@ -1,14 +1,22 @@
 #' Estimate Seroincidence
-#'
 #' @description
-#' Function to estimate seroincidences based on cross-section serology data and longitudinal
+#' Function to estimate seroincidences based on cross-sectional
+#' serology data and longitudinal
 #' response model.
 #'
-#' @param pop_data a [data.frame] with cross-sectional serology data per antibody and age, and additional columns corresponding to each element of the `strata` input
-#' @param strata a [character] vector of stratum-defining variables. Values must be variable names in `pop_data`.
-#' @param curve_strata_varnames A subset of `strata`. Values must be variable names in `curve_params`. Default = "".
-#' @param noise_strata_varnames A subset of `strata`. Values must be variable names in `noise_params`. Default = "".
-#' @param num_cores Number of processor cores to use for calculations when computing by strata. If set to more than 1 and package \pkg{parallel} is available, then the computations are executed in parallel. Default = 1L.
+#' @param pop_data a [data.frame] with cross-sectional serology data per
+#' antibody and age, and additional columns corresponding to
+#' each element of the `strata` input
+#' @param strata a [character] vector of stratum-defining variables.
+#' Values must be variable names in `pop_data`.
+#' @param curve_strata_varnames A subset of `strata`.
+#' Values must be variable names in `curve_params`. Default = "".
+#' @param noise_strata_varnames A subset of `strata`.
+#' Values must be variable names in `noise_params`. Default = "".
+#' @param num_cores Number of processor cores to use for
+#' calculations when computing by strata. If set to
+#' more than 1 and package \pkg{parallel} is available,
+#' then the computations are executed in parallel. Default = 1L.
 
 #' @details
 #'
@@ -17,7 +25,8 @@
 #' and then the data will be passed to [est.incidence()].
 #' If for some reason you want to use [est.incidence.by()]
 #' with no strata instead of calling [est.incidence()],
-#' you may use `NA`, `NULL`, or `""` as the `strata` argument to avoid that warning.
+#' you may use `NA`, `NULL`, or `""` as the `strata`
+#' argument to avoid that warning.
 #'
 #'
 #' @inheritParams est.incidence
@@ -26,7 +35,9 @@
 #'
 #' @return
 #' * if `strata` has meaningful inputs:
-#' An object of class `"seroincidence.by"`; i.e., a list of `"seroincidence"` objects from [est.incidence()], one for each stratum, with some meta-data attributes.
+#' An object of class `"seroincidence.by"`; i.e., a list of `
+#' "seroincidence"` objects from [est.incidence()], one for each stratum,
+#' with some meta-data attributes.
 #' * if `strata` is missing, `NULL`, `NA`, or `""`:
 #' An object of class `"seroincidence"`.
 #'
@@ -39,7 +50,8 @@
 #'
 #' curve <- load_curve_params("https://osf.io/download/rtw5k/") %>%
 #'   filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG")) %>%
-#'   slice(1:100, .by = antigen_iso) # Reduce dataset for the purposes of this example
+#'  # Reduce dataset for the purposes of this example
+#'   slice(1:100, .by = antigen_iso)
 #'
 #' noise <- load_noise_params("https://osf.io/download//hqy4v/")
 #'
@@ -49,13 +61,13 @@
 #'   curve_params = curve,
 #'   noise_params = noise %>% filter(Country == "Pakistan"),
 #'   antigen_isos = c("HlyE_IgG", "HlyE_IgA"),
-#'   #num_cores = 8 # Allow for parallel processing to decrease run time
+#'   # num_cores = 8 # Allow for parallel processing to decrease run time
 #'   iterlim = 5 # limit iterations for the purpose of this example
 #' )
 #'
 #' summary(est2)
 #'
 est.incidence.by <- function(
     pop_data,
     curve_params,
     noise_params,
@@ -75,16 +87,19 @@
     warning(
       "The `strata` argument to `est.incidence.by()` is missing.",
       "\n\n  If you do not want to stratify your data, ",
-      "consider using the `est.incidence()` function to simplify your code and avoid this warning.",
-      "\n\n Since the `strata` argument is empty, `est.incidence.by()` will return a `seroincidence` object, instead of a `seroincidence.by` object.\n"
+      "consider using the `est.incidence()` function to
+      simplify your code and avoid this warning.",
+      "\n\n Since the `strata` argument is empty, `est.incidence.by()`
+      will return a `seroincidence` object, instead of a
+      `seroincidence.by` object.\n"
     )
   }
 
   strata_is_empty <-
     missing(strata) ||
-      is.null(strata) ||
-      setequal(strata, NA) ||
-      setequal(strata, "")
+    is.null(strata) ||
+    setequal(strata, NA) ||
+    setequal(strata, "")
 
   if (strata_is_empty) {
     to_return <-
@@ -110,7 +125,7 @@
   )
 
   # Split data per stratum
   stratumDataList <- stratify_data(
     antigen_isos = antigen_isos,
     data = pop_data %>% filter(.data$antigen_iso %in% antigen_isos),
     curve_params = curve_params %>% filter(.data$antigen_iso %in% antigen_isos),
@@ -123,15 +138,17 @@
   strata_table <- stratumDataList %>% attr("strata")
 
   if (verbose) {
-    message("Data has been stratified.")
-    message("Here are the strata that will be analyzed:")
+    cli::cli_inform("Data has been stratified.")
+    cli::cli_inform("Here are the strata that will be analyzed:")
     print(strata_table)
   }
 
   if (num_cores > 1L && !requireNamespace("parallel", quietly = TRUE)) {
-    warning(
-      "The `parallel` package is not installed, so `num_cores > 1` has no effect.",
-      "To install `parallel`, run `install.packages('parallel')` in the console."
+    cli::cli_warn(
+      "The `parallel` package is not installed,
+      so `num_cores > 1` has no effect.",
+      "To install `parallel`, run `install.packages('parallel')`
+      in the console."
     )
   }
 
@@ -142,11 +159,11 @@
     num_cores <- num_cores %>% check_parallel_cores()
 
     if (verbose) {
-      message("Setting up parallel processing with `num_cores` = ", num_cores, ".")
+      message("Setting up parallel processing with
+              `num_cores` = ", num_cores, ".")
     }
 
-
     libPaths <- .libPaths()
     cl <-
       num_cores %>%
       parallel::makeCluster() %>%
@@ -155,14 +172,19 @@
       parallel::stopCluster(cl)
     })
 
+    # Export library paths to the cluster
     parallel::clusterExport(cl, c("libPaths"), envir = environment())
+
+    # Evaluate library loading on the cluster
     parallel::clusterEvalQ(cl, {
       .libPaths(libPaths)
-      require(serocalculator) # note - this gets out of sync when using load_all() in development
+      # note - this gets out of sync when using load_all() in development
+      require(serocalculator)
       require(dplyr)
     })
 
-    {
+    # Perform parallel computation and record execution time
+    time <- system.time({
       fits <- parallel::parLapplyLB(
         cl = cl,
         X = stratumDataList,
@@ -183,49 +205,41 @@
           )
         }
       )
-    } %>% system.time() -> time
+    })
 
     if (verbose) {
       message("Elapsed time for parallelized code: ")
       print(time)
     }
   } else {
-    # fits <- lapply(
-    #   X = stratumDataList,
-    #   FUN = function(x) est.incidence(dataList = x, verbose = verbose, ...))
-
-    fits <- list()
+    fits <- list()  # Initialize an empty list for fits
 
-    { # time progress
+    # Time progress
+    for (cur_stratum in names(stratumDataList)) {
+      cur_stratum_vars <- strata_table %>%
+        dplyr::filter(.data$Stratum == cur_stratum)
 
-      for (cur_stratum in names(stratumDataList))
-      {
-        cur_stratum_vars <-
-          strata_table %>%
-          dplyr::filter(.data$Stratum == cur_stratum)
-
-        if (verbose) {
-          message("starting new stratum: ", cur_stratum)
-          print(cur_stratum_vars)
-        }
+      if (verbose) {
+        message("starting new stratum: ", cur_stratum)
+        print(cur_stratum_vars)
+      }
 
-        fits[[cur_stratum]] <-
-          do.call(
-            what = est.incidence,
-            args = c(
-              stratumDataList[[cur_stratum]],
-              list(
-                lambda_start = lambda_start,
-                antigen_isos = antigen_isos,
-                build_graph = build_graph,
-                print_graph = print_graph,
-                verbose = verbose,
-                ...
-              )
-            )
+      fits[[cur_stratum]] <- do.call(
+        what = est.incidence,
+        args = c(
+          stratumDataList[[cur_stratum]],
+          list(
+            lambda_start = lambda_start,
+            antigen_isos = antigen_isos,
+            build_graph = build_graph,
+            print_graph = print_graph,
+            verbose = verbose,
+            ...
           )
-      }
-    } %>% system.time() -> time
+        )
+      )
+    }
+
 
     if (verbose) {
       message("Elapsed time for loop over strata: ")
@@ -233,7 +247,7 @@
     }
   }
 
   incidenceData <- structure(
     fits,
     antigen_isos = antigen_isos,
     Strata = strata_table,

tests/testthat/_snaps/est.incidence.by.md

@@ -0,0 +1,30 @@
+# est.incidence.by() produces expected results for typhoid data
+
+    Code
+      typhoid_results
+    Output
+      `seroincidence.by` object estimated given the following setup:
+      a) Antigen isotypes   : HlyE_IgG, HlyE_IgA 
+      b) Strata       : catchment 
+
+      This object is a list of `seroincidence` objects, with added meta-data attributes:`antigen_isos`   - Character vector of antigen isotypes used in analysis.
+      `Strata`       - Input parameter strata of function `est.incidence.by()`
+
+      Call the `summary()` function to obtain output results.
+
+---
+
+    structure(list("Stratum 1" = structure(list(minimum = 111.463434776159, 
+        estimate = -2.41983536033088, gradient = -3.74675296382939e-06, 
+        hessian = structure(15.251407603436, dim = c(1L, 1L)), code = 1L, 
+        iterations = 4L), class = c("seroincidence", "list"), lambda_start = 0.1, antigen_isos = c("HlyE_IgG", 
+    "HlyE_IgA")), "Stratum 2" = structure(list(minimum = 146.290052920951, 
+        estimate = -1.640278556419, gradient = -8.66368377467952e-08, 
+        hessian = structure(21.109488557163, dim = c(1L, 1L)), code = 1L, 
+        iterations = 5L), class = c("seroincidence", "list"), lambda_start = 0.1, antigen_isos = c("HlyE_IgG", 
+    "HlyE_IgA"))), antigen_isos = c("HlyE_IgG", "HlyE_IgA"), Strata = structure(list(
+        Stratum = c("Stratum 1", "Stratum 2"), catchment = c("aku", 
+        "kgh"), n = c(25L, 25L)), row.names = c(NA, -2L), class = c("tbl_df", 
+    "tbl", "data.frame"), strata_vars = "catchment"), graphs_included = FALSE, class = c("seroincidence.by", 
+    "list"))
+