From 47e7ec71239394fd614bc1ed2f9e2f5e064c59f9 Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Tue, 28 Jan 2025 13:10:19 +0000 Subject: [PATCH 1/3] fix detection of gaps for forecasts --- R/create.R | 37 ++-------------------- R/estimate_infections.R | 9 ++++-- R/opts.R | 8 ++--- R/preprocessing.R | 18 +++++++++-- tests/testthat/test-create_forecast_data.R | 27 ---------------- tests/testthat/test-forecast-opts.R | 5 ++- tests/testthat/test-preprocessing.R | 20 ++++++++++++ 7 files changed, 50 insertions(+), 74 deletions(-) delete mode 100644 tests/testthat/test-create_forecast_data.R diff --git a/R/create.R b/R/create.R index 5eac4a5dc..b32da40c6 100644 --- a/R/create.R +++ b/R/create.R @@ -411,35 +411,6 @@ create_obs_model <- function(obs = obs_opts(), dates) { return(data) } -##' Create forecast settings -##' -##' @param forecast A list of options as generated by [forecast_opts()] defining -##' the forecast opitions. Defaults to [forecast_opts()]. If NULL then no -##' forecasting will be done. -##' @inheritParams create_stan_data -##' @return A list of settings ready to be passed to stan defining -##' the Observation Model -##' @keywords internal -create_forecast_data <- function(forecast = forecast_opts(), data) { - if (forecast$infer_accumulate && any(data$accumulate)) { - accumulation_times <- which(!data$accumulate) - gaps <- unique(diff(accumulation_times)) - if (length(gaps) == 1 && gaps > 1) { ## all gaps are the same - forecast$accumulate <- gaps - cli_inform(c( - "i" = "Forecasts accumulated every {gaps} days, same as accumulation - used in the likelihood. To change this behaviour or silence this - message set {.var accumulate} explicitly in {.fn forecast_opts}." - )) - } - } - data <- list( - horizon = forecast$horizon, - future_accumulate = forecast$accumulate - ) - return(data) -} - #' Create Stan Data Required for estimate_infections #' #' @description`r lifecycle::badge("stable")` @@ -490,12 +461,8 @@ create_stan_data <- function(data, seeding_time, rt, gp, obs, backcalc, shifted_cases = shifted_cases, t = length(data$date), burn_in = 0, - seeding_time = seeding_time - ) - # add forecast data - stan_data <- c( - stan_data, - create_forecast_data(forecast, cases) + seeding_time = seeding_time, + horizon = forecast$horizon ) # add Rt data stan_data <- c( diff --git a/R/estimate_infections.R b/R/estimate_infections.R index 515717fd4..bb0e2a8f6 100644 --- a/R/estimate_infections.R +++ b/R/estimate_infections.R @@ -225,9 +225,14 @@ estimate_infections <- function(data, ## add forecast horizon if forecasting is required if (forecast$horizon > 0) { - reported_cases <- add_horizon( - reported_cases, forecast$horizon, forecast$accumulate + horizon_args <- list( + data = reported_cases, + horizon = forecast$horizon ) + if (!is.null(forecast$accumulate)) { + horizon_args$accumulate <- forecast$accumulate + } + reported_cases <- do.call(add_horizon, horizon_args) } # Create clean and complete cases diff --git a/R/opts.R b/R/opts.R index ca5c11dda..1eb4b3f0f 100644 --- a/R/opts.R +++ b/R/opts.R @@ -1113,13 +1113,11 @@ stan_opts <- function(object = NULL, #' forecast_opts(horizon = 28, accumulate = 7) forecast_opts <- function(horizon = 7, accumulate) { opts <- list( - horizon = horizon, - infer_accumulate = missing(accumulate) + horizon = horizon ) - if (missing(accumulate)) { - accumulate <- 1 + if (!missing(accumulate)) { + opts$accumulate <- accumulate } - opts$accumulate <- accumulate attr(opts, "class") <- c("forecast_opts", class(opts)) return(opts) } diff --git a/R/preprocessing.R b/R/preprocessing.R index a7aa7ac71..6056f8022 100644 --- a/R/preprocessing.R +++ b/R/preprocessing.R @@ -185,7 +185,8 @@ default_fill_missing_obs <- function(data, obs, obs_column) { ##' [estimate_secondary()]. See the documentation there for the expected ##' format. ##' @param accumulate The number of days to accumulate when generating posterior -##' prediction, e.g. 7 for weekly accumulated forecasts. +##' prediction, e.g. 7 for weekly accumulated forecasts. If this is not set an +##' attempt will be made to detect the accumulation frequency in the data. ##' @inheritParams fill_missing ##' @inheritParams estimate_infections ##' @importFrom data.table copy merge.data.table setDT @@ -210,8 +211,21 @@ add_horizon <- function(data, horizon, accumulate = 1L, .(date = seq(max(date) + 1, max(date) + horizon, by = "days")), by = by ] - ## if we accumulate add the column + ## detect accumulation + if (missing(accumulate) && "accumulate" %in% colnames(data)) { + accumulation_times <- which(!data$accumulate) + gaps <- unique(diff(accumulation_times)) + if (length(gaps) == 1 && gaps > 1) { ## all gaps are the same + accumulate <- gaps + cli_inform(c( + "i" = "Forecasts accumulated every {gaps} days, same as accumulation + used in the likelihood. To change this behaviour or silence this + message set {.var accumulate} explicitly in {.fn forecast_opts}." + )) + } + } if (accumulate > 1 || "accumulate" %in% colnames(data)) { + ## if we accumulate add the column initial_future_accumulate <- sum(cumsum(rev(!data$accumulate)) == 0) reported_cases_future[, accumulate := TRUE] ## set accumulation to FALSE where appropriate diff --git a/tests/testthat/test-create_forecast_data.R b/tests/testthat/test-create_forecast_data.R deleted file mode 100644 index 31acda247..000000000 --- a/tests/testthat/test-create_forecast_data.R +++ /dev/null @@ -1,27 +0,0 @@ -test_that("create_forecast_data returns expected default values", { - result <- create_forecast_data(data = example_confirmed) - - expect_type(result, "list") - expect_equal(result$horizon, 7) - expect_equal(result$future_accumulate, 1) -}) - -test_that("create_rt_data identifies gaps correctly", { - reported_weekly <- suppressWarnings(fill_missing( - example_confirmed[seq(1, 60, by = 7)], - missing_dates = "accumulate" - )) - expect_message( - result <- create_forecast_data(data = reported_weekly), - "same as accumulation used in the likelihood" - ) - expect_equal(result$future_accumulate, 7) -}) - -test_that("create_rt_data doesn't try to identify non-equally spaced gaps", { - reported_irregular <- example_confirmed[c(seq(1, 43, by = 7), 45)] - expect_no_message( - result <- create_forecast_data(data = reported_irregular) - ) - expect_equal(result$future_accumulate, 1) -}) diff --git a/tests/testthat/test-forecast-opts.R b/tests/testthat/test-forecast-opts.R index 6e5ac5ac5..005c3cb70 100644 --- a/tests/testthat/test-forecast-opts.R +++ b/tests/testthat/test-forecast-opts.R @@ -1,11 +1,10 @@ test_that("forecast_opts returns correct default values", { forecast <- forecast_opts() expect_equal(forecast$horizon, 7) - expect_equal(forecast$accumulate, 1) - expect_equal(forecast$infer_accumulate, TRUE) + expect_null(forecast$accumulate) }) test_that("forecast_opts sets infer_accumulate to FALSE if accumulate is given", { forecast <- forecast_opts(accumulate = 7) - expect_equal(forecast$infer_accumulate, FALSE) + expect_equal(forecast$accumulate, 7) }) diff --git a/tests/testthat/test-preprocessing.R b/tests/testthat/test-preprocessing.R index 02b1621c9..667ae5f6f 100644 --- a/tests/testthat/test-preprocessing.R +++ b/tests/testthat/test-preprocessing.R @@ -40,3 +40,23 @@ test_that("add_horizon works", { "Initial data point not marked as accumulated" ) }) + +test_that("add_horizon identifies gaps correctly", { + filled <- fill_missing(cases, missing_dates = "accumulate", initial_accumulate = 7) + expect_message( + result <- add_horizon(filled, horizon = 7), + "Forecasts accumulated every 7 days" + ) + result <- add_horizon(filled, horizon = 7, accumulate = 7) + expect_true(all(result[seq(.N - 6, .N - 1), accumulate])) + expect_false(result[.N, accumulate]) +}) + +test_that("add_horizon doesn't try to identify non-equally spaced gaps", { + reported_irregular <- example_confirmed[c(seq(1, 43, by = 7), 45)] + filled <- suppressWarnings( + fill_missing(reported_irregular, missing_dates = "accumulate") + ) + result <- add_horizon(filled, horizon = 7) + expect_false(any(result[seq(.N - 6, .N), accumulate])) +}) From d3253f1e241c01c4bfa5f835505fb113a89a199a Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Tue, 28 Jan 2025 14:41:19 +0000 Subject: [PATCH 2/3] update docs --- man/add_horizon.Rd | 3 ++- man/create_forecast_data.Rd | 30 ------------------------------ man/create_stan_data.Rd | 4 ---- man/epinow.Rd | 4 ---- man/estimate_infections.Rd | 4 ---- man/forecast_opts.Rd | 1 - man/regional_epinow.Rd | 4 ---- man/setup_forecast.Rd | 4 ---- 8 files changed, 2 insertions(+), 52 deletions(-) delete mode 100644 man/create_forecast_data.Rd diff --git a/man/add_horizon.Rd b/man/add_horizon.Rd index 72b5ad18d..f3590cdef 100644 --- a/man/add_horizon.Rd +++ b/man/add_horizon.Rd @@ -16,7 +16,8 @@ format.} horizon} \item{accumulate}{The number of days to accumulate when generating posterior -prediction, e.g. 7 for weekly accumulated forecasts.} +prediction, e.g. 7 for weekly accumulated forecasts. If this is not set an +attempt will be made to detect the accumulation frequency in the data.} \item{obs_column}{Character (default: "confirm"). If given, only the column specified here will be used for checking missingness. This is useful if diff --git a/man/create_forecast_data.Rd b/man/create_forecast_data.Rd deleted file mode 100644 index 3ca7b0ccb..000000000 --- a/man/create_forecast_data.Rd +++ /dev/null @@ -1,30 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/create.R -\name{create_forecast_data} -\alias{create_forecast_data} -\title{Create forecast settings} -\usage{ -create_forecast_data(forecast = forecast_opts(), data) -} -\arguments{ -\item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining -the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no -forecasting will be one.} - -\item{data}{A \verb{} of confirmed cases (confirm) by date (date). -\code{confirm} must be numeric and \code{date} must be in date format. Optionally -this can also have a logical \code{accumulate} column which indicates whether -data should be added to the next data point. This is useful when modelling -e.g. weekly incidence data. See also the \code{\link[=fill_missing]{fill_missing()}} function which -helps add the \code{accumulate} column with the desired properties when dealing -with non-daily data. If any accumulation is done this happens after -truncation as specified by the \code{truncation} argument.} -} -\value{ -A list of settings ready to be passed to stan defining -the Observation Model -} -\description{ -Create forecast settings -} -\keyword{internal} diff --git a/man/create_stan_data.Rd b/man/create_stan_data.Rd index e0f05ff26..9bcc76a1f 100644 --- a/man/create_stan_data.Rd +++ b/man/create_stan_data.Rd @@ -45,10 +45,6 @@ observation model. Defaults to \code{\link[=obs_opts]{obs_opts()}}.} define the back calculation. Defaults to \code{\link[=backcalc_opts]{backcalc_opts()}}.} \item{shifted_cases}{A \verb{} of delay shifted cases} - -\item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining -the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no -forecasting will be one.} } \value{ A list of stan data diff --git a/man/epinow.Rd b/man/epinow.Rd index 8a62820a5..9c5761e5b 100644 --- a/man/epinow.Rd +++ b/man/epinow.Rd @@ -71,10 +71,6 @@ Gaussian process.} \item{obs}{A list of options as generated by \code{\link[=obs_opts]{obs_opts()}} defining the observation model. Defaults to \code{\link[=obs_opts]{obs_opts()}}.} -\item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining -the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no -forecasting will be one.} - \item{stan}{A list of stan options as generated by \code{\link[=stan_opts]{stan_opts()}}. Defaults to \code{\link[=stan_opts]{stan_opts()}}. Can be used to override \code{data}, \code{init}, and \code{verbose} settings if desired.} diff --git a/man/estimate_infections.Rd b/man/estimate_infections.Rd index 992b65d46..3a4144a95 100644 --- a/man/estimate_infections.Rd +++ b/man/estimate_infections.Rd @@ -67,10 +67,6 @@ Gaussian process.} \item{obs}{A list of options as generated by \code{\link[=obs_opts]{obs_opts()}} defining the observation model. Defaults to \code{\link[=obs_opts]{obs_opts()}}.} -\item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining -the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no -forecasting will be one.} - \item{stan}{A list of stan options as generated by \code{\link[=stan_opts]{stan_opts()}}. Defaults to \code{\link[=stan_opts]{stan_opts()}}. Can be used to override \code{data}, \code{init}, and \code{verbose} settings if desired.} diff --git a/man/forecast_opts.Rd b/man/forecast_opts.Rd index 2bd02766b..cadd823aa 100644 --- a/man/forecast_opts.Rd +++ b/man/forecast_opts.Rd @@ -17,7 +17,6 @@ forecasts unless set explicitly here.} } \value{ A \verb{} object of forecast setting. -rstan functions. } \description{ \ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} diff --git a/man/regional_epinow.Rd b/man/regional_epinow.Rd index efaa0ee6b..c95586bbf 100644 --- a/man/regional_epinow.Rd +++ b/man/regional_epinow.Rd @@ -64,10 +64,6 @@ Gaussian process.} \item{obs}{A list of options as generated by \code{\link[=obs_opts]{obs_opts()}} defining the observation model. Defaults to \code{\link[=obs_opts]{obs_opts()}}.} -\item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining -the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no -forecasting will be one.} - \item{stan}{A list of stan options as generated by \code{\link[=stan_opts]{stan_opts()}}. Defaults to \code{\link[=stan_opts]{stan_opts()}}. Can be used to override \code{data}, \code{init}, and \code{verbose} settings if desired.} diff --git a/man/setup_forecast.Rd b/man/setup_forecast.Rd index 0940dabfc..dea9b610d 100644 --- a/man/setup_forecast.Rd +++ b/man/setup_forecast.Rd @@ -7,10 +7,6 @@ setup_forecast(forecast, horizon = NULL) } \arguments{ -\item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining -the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no -forecasting will be done.} - \item{horizon}{Deprecated; use \code{forecast} instead to specify the predictive horizon} } From 25e58f6d4791e062903f99a2d0ba90d698fecfc3 Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Tue, 28 Jan 2025 15:14:04 +0000 Subject: [PATCH 3/3] update docs --- R/create.R | 1 - R/estimate_infections.R | 5 ++++- man/create_stan_data.Rd | 4 ++++ man/epinow.Rd | 4 ++++ man/estimate_infections.Rd | 4 ++++ man/regional_epinow.Rd | 4 ++++ man/setup_forecast.Rd | 4 ++++ 7 files changed, 24 insertions(+), 2 deletions(-) diff --git a/R/create.R b/R/create.R index b32da40c6..07dda74eb 100644 --- a/R/create.R +++ b/R/create.R @@ -429,7 +429,6 @@ create_obs_model <- function(obs = obs_opts(), dates) { #' @inheritParams create_obs_model #' @inheritParams create_rt_data #' @inheritParams create_backcalc_data -#' @inheritParams create_forecast_data #' @importFrom stats lm #' @importFrom purrr safely #' @return A list of stan data diff --git a/R/estimate_infections.R b/R/estimate_infections.R index bb0e2a8f6..99d308e14 100644 --- a/R/estimate_infections.R +++ b/R/estimate_infections.R @@ -43,6 +43,10 @@ #' used as the `truncation` argument here, thereby propagating the uncertainty #' in the estimate. #' +#' @param forecast A list of options as generated by [forecast_opts()] defining +#' the forecast opitions. Defaults to [forecast_opts()]. If NULL then no +#' forecasting will be done. +#' #' @param horizon Deprecated; use `forecast` instead to specify the predictive #' horizon #' @@ -65,7 +69,6 @@ #' [estimate_truncation()] #' @inheritParams create_stan_args #' @inheritParams create_stan_data -#' @inheritParams create_forecast_data #' @inheritParams create_gp_data #' @inheritParams fit_model_with_nuts #' @inheritParams create_clean_reported_cases diff --git a/man/create_stan_data.Rd b/man/create_stan_data.Rd index 9bcc76a1f..e5680e01a 100644 --- a/man/create_stan_data.Rd +++ b/man/create_stan_data.Rd @@ -45,6 +45,10 @@ observation model. Defaults to \code{\link[=obs_opts]{obs_opts()}}.} define the back calculation. Defaults to \code{\link[=backcalc_opts]{backcalc_opts()}}.} \item{shifted_cases}{A \verb{} of delay shifted cases} + +\item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining +the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no +forecasting will be done.} } \value{ A list of stan data diff --git a/man/epinow.Rd b/man/epinow.Rd index 9c5761e5b..008759969 100644 --- a/man/epinow.Rd +++ b/man/epinow.Rd @@ -71,6 +71,10 @@ Gaussian process.} \item{obs}{A list of options as generated by \code{\link[=obs_opts]{obs_opts()}} defining the observation model. Defaults to \code{\link[=obs_opts]{obs_opts()}}.} +\item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining +the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no +forecasting will be done.} + \item{stan}{A list of stan options as generated by \code{\link[=stan_opts]{stan_opts()}}. Defaults to \code{\link[=stan_opts]{stan_opts()}}. Can be used to override \code{data}, \code{init}, and \code{verbose} settings if desired.} diff --git a/man/estimate_infections.Rd b/man/estimate_infections.Rd index 3a4144a95..4f5a97dde 100644 --- a/man/estimate_infections.Rd +++ b/man/estimate_infections.Rd @@ -67,6 +67,10 @@ Gaussian process.} \item{obs}{A list of options as generated by \code{\link[=obs_opts]{obs_opts()}} defining the observation model. Defaults to \code{\link[=obs_opts]{obs_opts()}}.} +\item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining +the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no +forecasting will be done.} + \item{stan}{A list of stan options as generated by \code{\link[=stan_opts]{stan_opts()}}. Defaults to \code{\link[=stan_opts]{stan_opts()}}. Can be used to override \code{data}, \code{init}, and \code{verbose} settings if desired.} diff --git a/man/regional_epinow.Rd b/man/regional_epinow.Rd index c95586bbf..1dfd319df 100644 --- a/man/regional_epinow.Rd +++ b/man/regional_epinow.Rd @@ -64,6 +64,10 @@ Gaussian process.} \item{obs}{A list of options as generated by \code{\link[=obs_opts]{obs_opts()}} defining the observation model. Defaults to \code{\link[=obs_opts]{obs_opts()}}.} +\item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining +the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no +forecasting will be done.} + \item{stan}{A list of stan options as generated by \code{\link[=stan_opts]{stan_opts()}}. Defaults to \code{\link[=stan_opts]{stan_opts()}}. Can be used to override \code{data}, \code{init}, and \code{verbose} settings if desired.} diff --git a/man/setup_forecast.Rd b/man/setup_forecast.Rd index dea9b610d..0940dabfc 100644 --- a/man/setup_forecast.Rd +++ b/man/setup_forecast.Rd @@ -7,6 +7,10 @@ setup_forecast(forecast, horizon = NULL) } \arguments{ +\item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining +the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no +forecasting will be done.} + \item{horizon}{Deprecated; use \code{forecast} instead to specify the predictive horizon} }