Need to first implement pseudo-code and then actual code for how to predict incidences from a single sample of Rts (with dates in a dataframe).

See scaffold here: https://github.com/epiforecasts/EpiSoon/blob/master/R/predict_incidences.R

a wrapper for infectiousness should be used (https://github.com/epiforecasts/EpiSoon/blob/master/R/infectiousness.R)

Please include a working example.

The code style is currently not consistent. I suggest we run styler over it and then make sure to lint etc as needed to keep consistent.

plot_forecast_evaluation(evaluations$forecast, observations, c(7)) +
  ggplot2::facet_grid(model ~ region) +
  cowplot::panel_border()

in the example, the function only takes one observation in.
instead it should be able to handle a data.frame (or ideally a list, since this is the format that was assigned to timeseries = list(observations, observations) with the observed values for the regions.

The newly added brms support is a good template for how we should be looking to add modelling packages. It would be great to lighten the dependencies by doing the same for fable_model and bsts_model.

There are a few different uses/definitions of the serial interval/generation time. If sticking with the currently implemented generative model this needs to be reviewed and discussed.

fable has been shown to have a memory leak and so is not suitable for large scale applications. forecastHybrid looks like a good drop-in replacement and appears to be a robust package with a good feature set. A potential issue is that it does not support sampling of the forecast out of the box and so an approximation is needed (potentially

A small example of a simple implementation:

library(forecastHybrid)


rts <- EpiSoon::example_obs_rts

y <- rts$rt

model <- forecastHybrid::hybridModel(y, weights = "cv.errors", 
                                     rolling = TRUE, cvHorizon = 7,
                                     horizonAverage = TRUE, 
                                     windowSize = 7)


preds <- forecast(model)


plot(preds)

Many of the functions may be better provided as S3 methods. It would be great to get a discussion going on whether we want to do this.

Hi,

As we are working on developing unit tests for the package functions (and considering I'm quite new to epidemiology), I could use some clarification regarding the behavior expected from forecast_cases().

I start here from the examples in forecast_cases() manual :

forecast <- forecast_rt(
	EpiSoon::example_obs_rts[1:10,],
	model = function(...) {
		EpiSoon::bsts_model(
			model = function(ss, y){
				bsts::AddAutoAr(ss, y=y, lags=10)
			},
			...
		)
	},
	horizon = 7,
	samples = 10
)

set.seed(42)
forecast_cases(
	EpiSoon::example_obs_cases,
	fit_samples = forecast,
	serial_interval = EpiSoon::example_serial_interval,
	horizon = 1
)

set.seed(42)
forecast_cases(
	EpiSoon::example_obs_cases,
	fit_samples = forecast,
	serial_interval = EpiSoon::example_serial_interval,
	horizon = 5
)

set.seed(42)
forecast_cases(
	EpiSoon::example_obs_cases,
	fit_samples = forecast,
	serial_interval = EpiSoon::example_serial_interval,
	horizon = 10
)

Isn't the output expected to be limited to rows with relevant "horizon" values ? How should forecast_cases() behave when horizon exceed the value provided to fit forecast ?

I am working with a version pulled from Github last week (8eddc07).

Best regards,
Sylvain

To do this we need to map back to reported cases and the design a measure of fit quality (could be current scores).

models are currently specified as functions in functions.
Maybe there is a better / more easily understandable alternative to that

The output of compare_timeseries is complex and hard to use. Sensible visualizations can make understanding the scores easier.

I have started (note this is a public-facing version of some work based on confidential data so may be slightly out of step) to think about how this can be visualised here: https://github.com/nikosbosse/covidevaluation/blob/master/R/evaluate_models.R

Wrapping some of these plots into package functions would be a great start to making our approach more easily adoptable by others.

## 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), ...)})

returns

Error: Element 3 must be a vector, not NULL
Run `rlang::last_error()` to see where the error occurred

The forecastHybrid package can produce some robust cross-validation weighted ensembles. Unfortunately, it handles confidence intervals by taking an unweighted mean of all candidate models or using the max/min for the candidate models.

As we only use samples from models (except when samples = 1 and a point forecast is used) this behaviour is a problem. In the current wrapper forecastHybrid_model we take samples by assuming the forecast is normally distributed and reconstructing a mean and standard deviation from the upper and lower bounds. This means that we see no benefit of any kind from a weighted ensemble. In an ideal scenario, the underlying package would produce better samples for its forecasts but failing that a more elegant workaround would be ideal.

The renewal code can be sped up using a dot product vs its current implementation. This a core part of the simulation code and so would give a good speed up. However, if planning on using thee EpiNow2 simulation model then this change is not needed as this part of the code base will be depreciated I would imagine.

Example:

 rts <- data.frame(rt = 1:10, date = as.Date("2020-01-01") + lubridate::days(1:10))


fit_model(rts, model = function(ss, y){bsts::AddAutoAr(ss, y = y, lags = 1000)},
         horizon = 7, samples = 10)

Whilst fable is the replacement for the forecast package it looks like much of the functionality has yet to be ported across. It would make sense to provide a wrapper for forecast.

I was experimenting with higher horizon number like 30

In this cases I think the Rt factor when going close to 0 for the forecast is overestimating the decline

Maybe one can fine a more conservative forecasting method or do you think it can be that for any european country Rt goes much below 1?

As the codebase matures (thanks to all contributors so far!) it might make sense to start thinking about having it peer-reviewed (this is linked to #56 ). I have had great experiences with JOSS so this would be my first thought.

To keep up with current practice. This should be done using the tooling in usethis.

The example:

observations <- data.frame(rt = 1:20,
                           date = as.Date("2020-01-01")
                                   + lubridate::days(1:20))

 ## Fit a model
 samples <- fit_model(observations[1:10],
                      model = function(ss, y){bsts::AddSemilocalLinearTrend(ss, y = y)},
                      horizon = 7, samples = 10)

 ## Score the model fit
 score_model(samples, observations)

doesn't work anymore if observations is used instead of observations[1:10]. In this case the function score_model throws an error

Can we use current plot functionality in a wrapper for plotting

When running compare_timeseries on a large grid I see this warning. Reprex needed for further debug.

Warning message:
In rdist(1, predictions[i, ]$infectiousness * predictions[i, ]$rt) :
NAs produced

We have been developing this package very much for our own use case and tightly integrated with EpiNow (our Rt estimation tooling). It would be great to hear if others are interested in EpiSoon as a standalone tool and if so would a CRAN release be helpful.

When entering a ts that is too short, current error is "Error: Column forecast_date is unknown"

We currently are using both - it would be great to standardise around one. I have not personally used GitHub's actions but hear it may be the new kid on the block.

Ideally, we would also automate documentation generation (i.e the pkgdown site) on new pull requests.

As we move to finalise the package it makes sense to remove the warning and improve the documentation. A sensible starting point would be to add a first paragraph that contains the information required by the JOSS guidelines.

This package focusses on providing case forecasts via first forecasting Rt values. An outstanding question is the value this adds over a simple case only forecast. It would be interesting to add tooling to enable case only forecasts to be easily compared to Rt based forecasts but this may require a complex rewrite. Opening this as a starting point for thinking about this.

As title, this likely indicates a depreciation issue or other bug from long term lack of development.

In general, the package has quite a few dependencies it would be great to lighten this where possible. It would also be great to try and improve performance. For example compare_timeseries can be very slow to run when considering many samples + models.

When usingforecastHybrid more than one core is used despite parallelization being turned off explicitly. This means that large runs are slow as cores conflict. Identifying someway to internally turn this off or externally prevent more than one core being accessed would both be good solutions.

This could very easily just be user error rather than a bug!

Issue opened on the package repo (ellisp/forecastHybrid#92)

obs_rts_aus <- readRDS("obs_rts_aus.rds")
obs_cases_aus <- readRDS("obs_cases_aus.rds")
models <- list("AR 1" = function(ss, y){bsts::AddAr(ss, y = y, lags = 1)})

forecast_eval <- compare_timeseries(obs_rts = obs_rts_aus,
obs_cases = obs_cases_aus,
serial_interval = EpiNow::covid_serial_intervals,
models = models,
horizon = 7, samples = 20)

Package needs to be tested again real Rt traces and the limitations explored.

For few regions, it makes sense to run in parallel over both regions and models. To do this need to strip out the reliance on compare_models and roll a new version that makes a data frame/list of all possible combinations and then maps over them.

It also makes sense to make compare_models potentially be parallel when used in standalone mode as well as evaluate_model potentially with a pass-through.

• Add integer and continuous scores
• Add model calibration check
• Score within samples