Compare timeseries and forecast models

Usage

compare_timeseries(
  obs_rts = NULL,
  obs_cases = NULL,
  models = NULL,
  horizon = 7,
  samples = 1000,
  bound_rt = TRUE,
  min_points = 3,
  timeout = 30,
  serial_interval = NULL,
  rdist = NULL,
  return_raw = FALSE
)

Arguments

obs_rts: A dataframe of observed Rts including a timeseries variable to denote each timeseris and a rt vector (to forecast) and a date vector (to denote time). Optionally this dataframe can contain samples for each timeseries in which case this should be denoted using a sample variable.
obs_cases: A dataframe of observed cases including a timeseries variable to denote each timeseris and a cases vector (to forecast) and a date vector (to denote time). Optionally this dataframe can contain samples for each timeseries in which case this should be denoted using a sample variable.
models: A list of models. A configuration is given in the examples. Each model needs to be wrapped in a function that takes a ... argument and returns a dataframe of samples with each column representing a time horizon. Example: function(...) {EpiSoon::bsts_model(model = function(ss, y){bsts::AddAr(ss, y = y, lags = 3)}, ...)}.
horizon: Numeric, the time horizon over which to predict.
samples: Numeric, number of samples to take.
bound_rt: Logical, defaults to TRUE. Should Rt values be bounded to be greater than or equal to 0.
min_points: Numeric, defaults to 3. The minimum number of time points at which to begin iteratively evaluating the forecast.
timeout: Numeric, timeout of model fitting in seconds. Defaults to 30 seconds.
serial_interval: A numeric vector describing the probability distribution the serial interval. See EpiNow::covid_serial_interval for an example of the format.
rdist: A function to be used to sample the number of cases. Must take two arguments with the first specfying the number of samples and the second the mean. Defaults to rpois if not supplied
return_raw: Logical, should raw cases and rt forecasts be returned. Defaults to FALSE.

Value

A list of dataframes as produced by evaluate model but with an additional model column.

Examples

if (FALSE) {
## Example data
obs_rts <- EpiSoon::example_obs_rts %>%
  dplyr::mutate(timeseries = "Region 1") %>%
  dplyr::bind_rows(EpiSoon::example_obs_rts %>%
    dplyr::mutate(timeseries = "Region 2"))

obs_cases <- EpiSoon::example_obs_cases %>%
  dplyr::mutate(timeseries = "Region 1") %>%
  dplyr::bind_rows(EpiSoon::example_obs_cases %>%
    dplyr::mutate(timeseries = "Region 2"))

## List of forecasting bsts models wrapped in functions.
models <- list(
  "AR 3" =
    function(...) {
      EpiSoon::bsts_model(
        model =
          function(ss, y) {
            bsts::AddAr(ss, y = y, lags = 3)
          }, ...
      )
    },
  "Semi-local linear trend" =
    function(...) {
      EpiSoon::bsts_model(
        model =
          function(ss, y) {
            bsts::AddSemilocalLinearTrend(ss, y = y)
          }, ...
      )
    },
  "ARIMA" =
    function(...) {
      fable_model(model = fable::ARIMA(y ~ time), ...)
    }
)


## Compare models
evaluations <- compare_timeseries(obs_rts, obs_cases, models,
  horizon = 7, samples = 10,
  serial_interval = EpiSoon::example_serial_interval
)

evaluations

## Example evaluation plot for comparing forecasts
## with actuals for a range of models and timeseries.
plot_forecast_evaluation(evaluations$forecast_rts, obs_rts, c(7)) +
  ggplot2::facet_grid(model ~ timeseries) +
  cowplot::panel_border()

## Hack to plot observed cases vs predicted
plot_forecast_evaluation(
  evaluations$forecast_cases,
  obs_cases, c(7)
) +
  ggplot2::facet_grid(model ~ timeseries, scales = "free") +
  cowplot::panel_border()
}