Add season argument and refactor epi_calendar

DESCRIPTION

@@ -19,7 +19,7 @@ Description: A powerful tool for automating the early detection of seasonal epid
 License: MIT + file LICENSE
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.2
 URL: https://github.com/ssi-dk/aedseo, https://ssi-dk.github.io/aedseo/
 BugReports: https://github.com/ssi-dk/aedseo/issues
 Suggests:

R/aedseo.R

@@ -21,10 +21,14 @@
 #' Choose between "poisson," or "quasipoisson".
 #' @param na_fraction_allowed Numeric value between 0 and 1 specifying the
 #' fraction of observables in the window of size k that are allowed to be NA.
+#' @param season A numeric vector of length 2 c(start,end) with the start and
+#' end weeks of the seasons to stratify the observations by.
+#' Ex: season = c(21,20). Default is NULL.
 #'
 #' @return A `aedseo` object containing:
 #'   - 'reference_time': The time point for which the growth rate is estimated.
 #'   - 'observed': The observed value in the reference time point.
+#'   - 'season': A stratification of observables in corresponding seasons.
 #'   - 'growth_rate': The estimated growth rate.
 #'   - 'lower_growth_rate': The lower bound of the growth rate's confidence
 #'   interval.
@@ -36,8 +40,8 @@
 #'   - 'sum_of_cases_warning': Logical. Does the Sum of Cases exceed the
 #'   disease threshold?
 #'   - 'seasonal_onset_alarm': Logical. Is there a seasonal onset alarm?
-#'   - 'converged': Logical. Was the IWLS judged to have converged?
 #'   - 'skipped_window': Logical. Was the window skipped due to missing?
+#'   - 'converged': Logical. Was the IWLS judged to have converged?
 #'
 #' @export
 #'
@@ -64,6 +68,7 @@
 #'   disease_threshold = 200,
 #'   family = "poisson",
 #'   na_fraction_allowed = 0.4,
+#'   season = NULL
 #' )
 #'
 #' # Print the AEDSEO results
@@ -78,7 +83,8 @@ aedseo <- function(
       "quasipoisson"
       # TODO: #10 Include negative.binomial regressions. @telkamp7
     ),
-    na_fraction_allowed = 0.4) {
+    na_fraction_allowed = 0.4,
+    season = NULL) {
   # Check input arguments
   coll <- checkmate::makeAssertCollection()
   checkmate::assert_data_frame(tsd)
@@ -101,6 +107,13 @@ aedseo <- function(
   res <- tibble::tibble()
   skipped_window <- base::rep(FALSE, base::nrow(tsd))
 
+  # Add the seasons to tsd if available
+  if (!is.null(season)) {
+    tsd <- tsd |> dplyr::mutate(Season = epi_calendar(time))
+  } else {
+    tsd <- tsd |> dplyr::mutate(Season = "Not defined")
+  }
+
   for (i in k:n) {
     # Index observations for this iteration
     obs_iter <- tsd[(i - k + 1):i, ]
@@ -138,6 +151,7 @@ aedseo <- function(
       tibble::tibble(
         reference_time = tsd$time[i],
         observed = tsd$observed[i],
+        season = tsd$Season[i],
         growth_rate = growth_rates$estimate[1],
         lower_growth_rate = growth_rates$estimate[2],
         upper_growth_rate = growth_rates$estimate[3],

R/epi_calendar.R

@@ -23,34 +23,32 @@
 #'
 #' @examples
 #' # Check if a date is within the epidemiological season
-#' epi_calendar(as.Date("2023-09-15"), start = 40, end = 20)
-#' # Expected output: "2022/2023"
+#' epi_calendar(as.Date("2023-09-15"), start = 21, end = 20)
+#' # Expected output: "2023/2024"
 #'
-#' epi_calendar(as.Date("2023-05-01"), start = 40, end = 20)
+#' epi_calendar(as.Date("2023-05-30"), start = 40, end = 20)
 #' # Expected output: "out_of_season"
 #'
-#' epi_calendar(as.Date("2023-01-15"), start = 40, end = 20)
+#' epi_calendar(as.Date("2023-01-15"), start = 21, end = 20)
 #' # Expected output: "2022/2023"
 #'
-#' epi_calendar(as.Date("2023-12-01"), start = 40, end = 20)
+#' epi_calendar(as.Date("2023-12-01"), start = 21, end = 20)
 #' # Expected output: "2023/2024"
-epi_calendar <- Vectorize(function(date, start = 40, end = 20) {
+epi_calendar <- Vectorize(function(date, start = 21, end = 20) {
   # Compute the current week
-  current_week <- as.integer(format(x = date, "%U"))
+  current_week <- as.integer(format(x = date, "%V"))
 
-  if (current_week <= start & end <= current_week) {
+  if (!(current_week >= start | current_week <= end)) {
     return("out_of_season")
   }
 
   # Compute the current year
-  current_year <- format(date, "%Y")
-  # ... and turn into integer
-  current_year_integer <- as.integer(current_year)
+  current_year <- as.integer(strftime(date, format = "%G"))
 
   if (current_week <= end) {
-    ans <- paste0(current_year_integer - 1, "/", current_year_integer)
+    ans <- paste0(current_year - 1, "/", current_year)
   } else {
-    ans <- paste0(current_year_integer, "/", current_year_integer + 1)
+    ans <- paste0(current_year, "/", current_year + 1)
   }
 
   return(ans)

R/summary.R

@@ -46,6 +46,13 @@ summary.aedseo <- function(object, ...) {
   # Extract the reference time
   reference_time <- last_observation$reference_time
 
+  # Extract the time_interval
+  time_interval <- attr(object, "time_interval")
+
+  # Extract the seasons
+  seasons <- unique(object$season)
+  seasons <- paste(seasons, collapse = ", ")
+
   # Latest sum of cases
   latest_sum_of_cases <- object %>%
     dplyr::filter(dplyr::row_number() == dplyr::n()) %>%
@@ -102,6 +109,9 @@ summary.aedseo <- function(object, ...) {
     calculation of sum of cases:
       %d
 
+    The time interval for the observations:
+      %s
+
     Disease specific threshold:
       %d
 
@@ -123,9 +133,13 @@ summary.aedseo <- function(object, ...) {
       %s
 
     Latest seasonal onset alarm:
+      %s
+
+    The seasons defined in the series:
       %s",
     family,
     k,
+    time_interval,
     disease_threshold,
     as.character(reference_time),
     latest_sum_of_cases,
@@ -137,7 +151,8 @@ summary.aedseo <- function(object, ...) {
     last_observation$upper_growth_rate,
     sum_of_growth_warnings,
     as.character(latest_growth_warning),
-    as.character(latest_seasonal_onset_alarm)
+    as.character(latest_seasonal_onset_alarm),
+    seasons
   )
 
   # Print the summary message

tests/testthat/test-epi_calendar.R

@@ -1,16 +1,16 @@
 # Test if a date within the season returns the correct season
 test_that("Within the season, correct season is returned", {
   expect_equal(
-    epi_calendar(as.Date("2023-03-15"), start = 40, end = 20), "2022/2023"
+    epi_calendar(as.Date("2023-03-15"), start = 21, end = 20), "2022/2023"
   )
   expect_equal(
-    epi_calendar(as.Date("2023-05-01"), start = 40, end = 20), "2022/2023"
+    epi_calendar(as.Date("2023-05-01"), start = 21, end = 20), "2022/2023"
   )
   expect_equal(
-    epi_calendar(as.Date("2023-01-15"), start = 40, end = 20), "2022/2023"
+    epi_calendar(as.Date("2023-01-15"), start = 21, end = 20), "2022/2023"
   )
   expect_equal(
-    epi_calendar(as.Date("2023-12-01"), start = 40, end = 20), "2023/2024"
+    epi_calendar(as.Date("2023-12-01"), start = 21, end = 20), "2023/2024"
   )
 })
 
@@ -26,3 +26,15 @@ test_that("Outside the season, 'out_of_season' is returned", {
     epi_calendar(as.Date("2023-06-30"), start = 40, end = 20), "out_of_season"
   )
 })
+
+# Test that all dates from week 53 belongs to correct season 2015/2016
+test_that("All dates from week 53 belongs to correct season 2015/2016", {
+  week_53_season_15_16 <- c("2015-12-28", "2015-12-29", "2015-12-30",
+                            "2015-12-31", "2016-01-01", "2016-01-02",
+                            "2016-01-03")
+  results <- purrr::map(week_53_season_15_16,
+                        ~ epi_calendar(as.Date(.x), start = 21, end = 20))
+  results_vector <- unlist(results)
+
+  expect_equal(results_vector, rep("2015/2016", length(week_53_season_15_16)))
+})