As the latest estimates are essentially forecasts it makes sense in many contexts to report estimates that are based on more complete data. This can either be via the internal categorical switch (estimate, partial estimate etc.) which is based on the mean of all input delays, or be based on a numeric number of days on which to lag.

This feature needs to be added to the following

• Summary tables produced within individual epinow calls
• Summary tables produced using regional_epinow
• Captions in the Rmarkdown templates (in inst/templates).

Sensible to explicitly model imports. Ideally, this should vary with time but a sensible first approximation is to provide a background constant rate.

Related to #30.

We ran regional Epinow2 on a full time-series (Covid-19 in the UK from March to August), but we had a couple of problems. Below is a description of what we tried, and the model specifications that produced sensible estimates for all regions.
The actual work on this is here. Some of the problems can be seen in the various iterations of the plots

Model inputs:

• date = time-series of ~130 days
• region = 11 regions (UK data: 7 NHS regions + 4 nations)
• confirm = sequentially ran each of
  • test-positive cases
  • hospital admissions
  • deaths

Model spec:

• In all runs we set the following parameters:
  • samples = 2000
  • burn_in = 0
  • chains = determined by number of cores
• Initial run, which produced strange outputs:
  • adapt_delta = 0.98
  • warmup = 1,000
• Final run, which produced sensible outputs:
  • adapt_delta = 0.99
  • warmup = 2,000
• Model family: we tried specifying the model family as Poisson, but this produced even stranger results. Switched back to default negative binomial.

Issues:

• It looked like estimation got stuck on some regions for very long periods, eg 4 hours. Overall runtime took up to 10 hours (this is for three runs, one for each data source). Running about half the time-series (the most recent 12 weeks) took <2hrs.
• Hugely inflated uncertainty for some regions in only one data source (which changed depending on the run). Unclear why this was happening. Could be due to missing and low counts in data.
• Example - here

Solution: Increasing warm-up seemed to be the parameter that produced the "most sensible" results.

Discovered this while working on #29 and I am happy to fix at the same time. Here is a reprex of the problem:

library(rnaturalearth)
library(EpiNow2)

df <- data.table::data.table(variable = c("Increasing", "Increasing", "Increasing"),
                             country = c("France", "Germany", "United Kingdom") )

global_map(df, variable = "variable")

Which produces this:

The problem relates to missing iso_a3 values in the natural earth data (iso_a3 is used in global_map to join countries' data to the map):

> world <- rnaturalearth::ne_countries(scale='medium',returnclass = 'sf')
> dplyr::filter(world, is.na(iso_a3))$name_en
[1] "Somaliland"                          "Norway"                             
[3] "Kosovo"                              "France"                             
[5] "Turkish Republic of Northern Cyprus" "Australian Indian Ocean Territories"
[7] "Ashmore and Cartier Islands"         "Siachen Glacier"

This seems to be a known issue for some time. It might relate to only some versions. I will find a more robust way to get around the problem.

EpiNow2 contains multiple long running large functions that each have multiple features. These should ideally be tested using unit tests. The complication is that some of this functionality is very long running and stochastic and hence hard to test. Unit tests should ideally be build using testthat and coverage functionality that can't be tested using lower level (#53) unit tests.

By default EpiNow2 is now noisy and produces more messages than previously. It would be sensible to give users more control over this.

The package standard has been to use missing arguments when no default is given. The issue with doing this is that a missing argument from a top-level function cannot be passed into a lower-level function. This means converting to a NULL is required and adding handling of missing -> NULL. It might be cleaner to instead default arguments to NULL.

From review comments by @joeHickson on #74

Use: github.com/epiforecasts/backcalc as a starting point

Provide a (crude) mechanism for ending individual location processing after a timeout.
https://www.rdocumentation.org/packages/R.utils/versions/2.9.2/topics/withTimeout

Note (re: stan code keeping running):
"More precisely, if a function is implemented in native code (e.g. C) and the developer of that function does not check for user interrupts, then you cannot interrupt that function neither via a user interrupt (e.g. Ctrl-C) nor via the built-in time out mechanism. To change this, you need to contact the developer of that piece of code and ask them to check for R user interrupts in their native code."
I have found reference to Rstan checking for interrupts every iteration so we'll have to see what happens.

Unit tests for low level functionality (i.e plotting functions, data manipulation functions etc) would greatly improve the robustness of the package. These should be implemented using testthat and each test a single feature

The README is very helpful but more could be added on the detail of the estimation method. Attached below a first quick summary written by @seabbs. This could be put in the README and/or function docs.

Not finalised: first draft version

EpiNow2 implements a fully bayesian latent variable approach, which works as follows. Initial seed infections are combined with time-varying estimates of Rt, via an uncertain estimate of the generation time, to produce a time-series of imputed infections using the Cori et al approach [1]. Infections are then mapped to reported case counts by convolving over an uncertain incubation period and report delay, and then applying a negative binomial observation model combined with a day of the week effect.

Temporal variation is controlled using a Gaussian process with a squared exponential kernel, with the length scale and magnitude estimated during the model fitting process.

We used a

• gamma distributed generation time, based on [2], with
  • mean 3.6 (standard deviation 0.7), and
  • standard deviation 3.1 (SD 0.8)
• a log-normal incubation period distribution, based on [3], with
  • mean of 5.2 (SD 1.1) and
  • SD 1.52 (SD 1.1)

Report delays included

• a report delay for cases with a
  • mean of 4 (SD 1.4) and
  • SD 2.3 (SD 1.1)
• a report delay for deaths with a
  • mean of 9.9 (SD 1.1)
  • SD 2.1 (SD 1.1)

References:
[1] Cori A, Ferguson NM, Fraser C, Cauchemez S. A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics. Am J Epidemiol. 2013;178(9):1505–12.
[2] Ganyani T, Kremer C, Chen D, Torneri A, Faes C, Wallinga J, et al. Estimating the generation interval for coronavirus disease (COVID-19) based on symptom onset data, March 2020. Eurosurveillance. 2020 Apr 30;25(17):2000257.
[3] Lauer SA, Grantz KH, Bi Q, et al.: The incubation period of coronavirus disease 2019 (covid-19) from publicly reported confirmed cases: Estimation and application. Ann Intern Med. 2020; 172(9): 577–82.

The package currently provides some undocumented templates (inst/templates) that we use internally for reporting purposes and as a backend for our website. It would be great to expand on these to provide easy reporting for package users. Ideally, this would support both interactive html and static pdf reports.

An interesting alternative might be to provided a packaged shiny dashboard that can be used with results saved in the standard package format. An example of this idea is here. The same repo also contains an example of a basic automated report.

Interventions coded as breakpoints may be explored using the current model. This should be done in a model comparison framework scoring models with LFO.

Scoping from @sbfnk indicates that variational bayes gives sensible results very rapidly (though at the cost of some reduction in uncertainty and occasional fitting failures).

Supporting this would be useful for applications where little compute is available and reduced robustness in results is acceptable.

It would make sense to add this via fit_model with an argument passed in to estimate_infections to control the fitting used. Defaults supplied to stan can then be altered. An alternative approach that led to fewer arguments being surfaced might be sensible as estimate_infections has some quite severe argument bloat. Key arguments that are needed would be iter (iterations to use) and something to indicate the number of trials (in case of variational Bayes stochastically failing). Suggest combining resolving this issue with #73.

The suggestion from @sbfnk is to surface this as an "exact"/"approximate" option with the technical changes handled internally (but ultimately able to be changed by the user).

Test from @sbfnk

vi_test.pdf

using:

vb_args <- stan_args
vb_args$iter <- 10000
vb_args$refresh <- 1
vb_args$warmup <- NULL
vb_args$cores <- NULL
vb_args$chains <- NULL
vb_args$control <- NULL
vb_args$save_warmup <- NULL
fit <- do.call(rstan::vb, vb_args)

In order for this to work the model had to changed to explicitly stop negative reports and infections.

Method paper for reference: https://arxiv.org/pdf/1506.03431.pdf

The current approach uses a parameter for each day of the week making it hard to fit. Approximation approaches using Fourier series etc exist and maybe a sensible alternative. Ideally if implemented in a robust fashion these alternative approaches should be scalable so that, in the future, they can be made to vary slowly over time.

Test a model where the distribution of Rt remains at the last bit of orange; i.e. be a bit more conservative with uncertainty in the forecast whilst using the maximal amount of data.

Currently, we use a crude approach that timeouts whole regions when running over a time limit. This approach is problematic because most of the time only a single chain is taking a long time to finish whilst others converge rapidly.

It would be sensible to provide chain by chain timing out with failure only occurring if less than two chains (in order to assess convergence) are not timed out.

An example of one implementation of this approach is here: https://github.com/mike-lawrence/ezStan/blob/master/R/collect_stan.R

There is another example of this in brms: https://github.com/paul-buerkner/brms/blob/87afed005e67f7685fe1c53584084775406066d3/R/backends.R#L112

The stationary GP formulation may be more robust. Using the currently implemented kernals it does not work for sparse. Using a kernel that combines a long and short lengthscale may solve this.

At the moment the default setting is to assume that all NAs in timeseries are in fact 0 case counts. This can lead to some very odd results when NAs are sporadic but not extremely rare.

It would be sensible to resolve this by adding a heuristic (ideally user settable) that treats NAs as either NAs or true 0's.

This requires:

• Updating create_clean_reported_cases to keep some NAs. This likely needs some kind of detection rule in order to either allocate NAs to be 0 or NA (some regions never report cases on Saturday for example).
• Adding support to estimate_infections.stan so that NAs are skipped in the likelihood. using an ifelse or similar).

Suggested by @sbfnk.

Just a general question. I built a fork using cmdstanr which should be extensible and avoid some of the rstan issues in the future. Downside is that cmdstanr isn't (yet) on CRAN and it adds an external system requirement for cmdstan. I don't know if there is interest, but I could update my code and include cmstanr as an option. It could be an argument to estimate_infections

 use_cmdstan = FALSE # in arguments

Then in body of function:

if (use_cmdstan){
    
    if(is.null(model)){
      model <- system.file("stan", "estimate_infections.stan", package = "EpiNow2")
    }
    
    #stan_file <- here::here("inst", "stan", "estimate_infections.stan")
    
    mod <- cmdstanr::cmdstan_model(model, compile = TRUE)
    
    if (verbose) {
      message(paste0("Running for ", samples," samples (across ", chains, 
                     " chains each with a warm up of ", warmup, " iterations each) and ",
                     data$t," time steps of which ", horizon, " are a forecast"))
    }
    
    if(is.null(threads_per_chain)){
      threads_per_chain <- 1L
    }
    
    cmdfit <- mod$sample(data = data, chains = chains,
                      parallel_chains = cores,
                      init = init_fun, iter_warmup = warmup,
                      iter_sampling =  ceiling(samples / chains), 
                      max_treedepth = max_treedepth,
                      threads_per_chain = threads_per_chain,
                      adapt_delta = adapt_delta,
                      refresh = ifelse(verbose, 50, 0),
                      save_warmup = debug)
    
    fit <- rstan::read_stan_csv(cmdfit$output_files())

I don't know how much IO time could add, but I have been finding cmdstanr is more stable.

As an optional module it makes sense to fit to deaths and cases (in a similar way to the implementation of breakpoints). There are then choices to make regarding how deaths and cases are linked. An obvious choice is to assume some kind of underreporting parameter (that also varies over time using a GP). This could be parameterized using testing data as a prior.

This needs #10 to be done before it can be implemented.

This should be a simple case of updating the nowcast scripts in each results repo. The methods section also needs to be updated. I would suggest keeping this brief so we can move the methods to a methods paper that also have validation (and make the website paper a our first version of the methods + its a website paper)

This package's predecessor EpiNow could distinguish the import_status of cases between local and imported. Is this something that is already feasible with EpiNow2 and I have just missed how to do it? If not, is it a possible enhancement for the future?

Feature creep has led to increasingly complex and difficult to maintain stan and R interface code. A refactor using stan functions, stan snippets and R functions for the R interface code would both improve maintainability and open the door for feature features/uses of the current code base.

Estimates may appear spuriously precise if given to greater levels of precision.

It makes sense to highlight downstream uses cases both for other users and in order to advertise features of EpiNow2.

In order for this to exist an initial list of downstream uses is required.

1. Global and subnational estimates of Rt for Covid-19 in real-time: https://github.com/epiforecasts/covid-rt-estimates
2. Estimates of R0 and Rt across 500 cities globally (used as data in a study to determine evidence for seasonality in Covid-19 tranmission): Link to follow
3. Rt estimates pre and post-intervention in LMIC settings globally. Developed by @pearsonca for a HPC setting. Feeds in to regularly updated LMIC reports supplied to policy makers. Link: https://github.com/cmmid/epinow2_hpcarray

For interest: @sbfnk @joeHickson

As apart of the next release it might make sense to deal with some of the current argument bloat.

One major contributor to the number of arguments are the stan controls. As stan_args allows any arguments to be passed to stan the number of surfaced arguments could potentially be reduced (i.e removing warmup, chains, adapt delta etc) and reducing these to internal defaults. This would however be a fairly substantial breaking change for any downstream users.

At the moment these helpfiles are minimal. I'd be happy to do a more illustrative example if you're up for pull requests.

A call to epinow with the data in a data frame that works fine if you aren't saving results will cause an error if you try to save the results by specifying target_folder. Reproducible example:

library(EpiNow2)
## Construct example distributions
generation_time <- list(mean = EpiNow2::covid_generation_times[1, ]$mean,
                        mean_sd = EpiNow2::covid_generation_times[1, ]$mean_sd,
                        sd = EpiNow2::covid_generation_times[1, ]$sd,
                        sd_sd = EpiNow2::covid_generation_times[1, ]$sd_sd,
                        max = 30)

incubation_period <- list(mean = EpiNow2::covid_incubation_period[1, ]$mean,
                          mean_sd = EpiNow2::covid_incubation_period[1, ]$mean_sd,
                          sd = EpiNow2::covid_incubation_period[1, ]$sd,
                          sd_sd = EpiNow2::covid_incubation_period[1, ]$sd_sd,
                          max = 30)

reporting_delay <- EpiNow2::bootstrapped_dist_fit(rlnorm(100, log(6), 1))
## Set max allowed delay to 60 days to truncate computation
reporting_delay$max <- 60

## Example case data
reported_cases <- as.data.frame(EpiNow2::example_confirmed[1:40] )

## This works:
out <- epinow(reported_cases = reported_cases, generation_time = generation_time,
              delays = list(incubation_period, reporting_delay),
              rt_prior = list(mean = 1, sd = 1),
              samples = 1000, warmup = 200, cores = 1, chains = 1,
              verbose = TRUE, return_fit = TRUE)

## This doesn't:
out <- epinow(reported_cases = reported_cases, generation_time = generation_time,
              delays = list(incubation_period, reporting_delay),
              rt_prior = list(mean = 1, sd = 1),
              samples = 1000, warmup = 200, cores = 1, chains = 1,
              verbose = TRUE, return_fit = TRUE, target_folder = 'saved-results')

## This works too:
out <- epinow(reported_cases = as.data.table(reported_cases), generation_time = generation_time,
              delays = list(incubation_period, reporting_delay),
              rt_prior = list(mean = 1, sd = 1),
              samples = 1000, warmup = 200, cores = 1, chains = 1,
              verbose = TRUE, return_fit = TRUE, target_folder = 'saved-results')

This is the error:

Error in `[.data.frame`(reported_cases, confirm > 0) : 
  object 'confirm' not found

> traceback()
7: `[.data.frame`(reported_cases, confirm > 0)
6: reported_cases[confirm > 0]
5: epinow(reported_cases = reported_cases, generation_time = generation_time, 
       delays = list(incubation_period, reporting_delay), rt_prior = list(mean = 1, 
           sd = 1), samples = 1000, warmup = 200, cores = 1, chains = 1, 
       verbose = TRUE, return_fit = TRUE, target_folder = "saved-results") at .active-rstudio-document#23
4: eval(ei, envir)
3: eval(ei, envir)
2: withVisible(eval(ei, envir))
1: source("~/.active-rstudio-document", echo = TRUE)

This could probably be fixed with an appropriate as.data.table() somewhere but I will leave it to you whether you want to do that, or just specify that reported_cases must be a data.table (I would vote you allow data frames and tibbles but that's your call).

I think the StanHeaders dependency in DESCRIPTION needs to be updated. I had StanHeaders version 2.21.0-3 installed, and I got an error about a missing header file. It worked when I upgraded StanHeaders to the latest version (2.21.0-5).

Currently this method is a little muddled - half the optional data arrives as parameters, the rest in a list. Consider revising so they are all exposed (and therefore validated - e.g. stan_args = list(iterations=1000) wouldn't pick up the misnamed variable).

This might be a good place for an object. I feel like it should be possible in R4 but the following is in R6:

AbstractStanArgs <- R6Class("AbstractStanArgs", list(
  object = NA, # consider renaming this model or renaming the arg object
  data = NA,
  init = NA,
  refresh = 0,
  initialize = function(model = stanmodels$estimate_infections,
                        data = NULL,
                        init = "random",
                        verbose = FALSE) {
    self$object <- model
    self$data <- data
    self$init <- init
    self$refresh <- ifelse(verbose, 50, 0)
  }
))

StanArgsExact <- R6Class("StanArgsExact",
                         inherit = AbstractStanArgs,
                         public = list(cores = NA,
                                       warmup = NA,
                                       save_warmup = FALSE,
                                       chains = NA,
                                       control.adapt_delta = NA,
                                       control.max_treedepth = NA,
                                       iter = NA,
                                       initialize = function(...,
                                                             cores = 1,
                                                             warmup = 500,
                                                             save_warmup = FALSE,
                                                             chains = 4,
                                                             control.adapt_delta = 0.99,
                                                             control.max_treedepth = 15,
                                                             samples = 1000
                                       ) {
                                         super$initialize(...)
                                         self$cores <- cores
                                         self$warmup <- warmup
                                         self$save_warmup <- save_warmup
                                         self$chains <- chains
                                         self$control.adapt_delta <- control.adapt_delta
                                         self$control.max_treedepth <- control.max_treedepth
                                         self$iter <- ceiling(samples / chains) + warmup
                                       }))

StanArgsApproximate <- R6Class("StanArgsApproximate",
                               inherit = AbstractStanArgs,
                               public = list(trials = NA,
                                             iter = NA,
                                             output_samples = NA,
                                             initialize = function(...,
                                                                   trials = 5,
                                                                   iter = 10000,
                                                                   output_samples = 1000
                                             ) {
                                               super$initialize(...)
                                               self$trials <- trials
                                               self$iter <- iter
                                               self$output_samples <- output_samples
                                             }))

This is slightly longer code but produces a strongly typed object. It does require a different call depending on the type of object you require but that gives a cleaner method signature.

e.g.

create_stan_args(method = "approximate")  -> StanArgsApproximate$new()
create_stan_args(stan_args = list(warmup = 1000)) -> StanArgsExact$new(warmup=1000)

This object could replace the list of stan_args that is being passed around.

e.g.

regional_epinow(..... , stan_args = StanArgsExact$new(warmup = 100, cores = 1, chains = 2), .......)

It can then be further modified by inner calls replacing defaults.

# estimate_infections.R::estimate_infections (L305 and similarly for data)
stan_args$init <- create_initial_conditions(data, delays, rt_prior, 
                                                                  generation_time, mean_shift)

If the init / data are only ever generated at lower level it might be that they shouldn't be part of the stan_args object or a simpler object is created that the user interacts with that is included in the lower level object that is passed to stan.

Originally posted by @joeHickson in #74 (comment)

When estimates are very skewed central high-density regions frequently do not contain the mean. This may mislead users/consumers of results. Quantiles are likely to be more generally understood and user friendly.

Suggested by @sbfnk

Hi. I am trying to install EpiNow2 on RStudio cloud which uses Ubuntu software. I have tried via drat and remotes. Same issue occurred in a Windows laptop. I got this error:

g++: internal compiler error: Killed (program cc1plus)
Please submit a full bug report,
with preprocessed source if appropriate.
See <file:///usr/share/doc/gcc-5/README.Bugs> for instructions.
/opt/R/3.6.3/lib/R/etc/Makeconf:177: recipe for target 'RcppExports.o' failed
make: *** [RcppExports.o] Error 4
ERROR: compilation failed for package ‘EpiNow2’
* removing ‘/home/rstudio-user/R/x86_64-pc-linux-gnu-library/3.6/EpiNow2’
Warning in install.packages :
  installation of package ‘EpiNow2’ had non-zero exit status

The downloaded source packages are in
	‘/tmp/RtmpRQhR0s/downloaded_packages’

Seeing -NANs during convolutions and blow ups in the Poisson rate paramter. Issue only occasionally happens on some chains indiciating perhaps its an issue in some of the parameter space but not all.

Due to the speed of implementation and code from multiple authors package switches styles. Sensible to implement a standard style and adding linting to check new additions.

Based on comments by @joeHickson in #74

I think there might be a typo error with a repeated incubation_def causing the pipeline to fail

https://github.com/epiforecasts/EpiNow2/runs/865961150#step:10:288

Code is largely unreviewed. It would be good to change this.

Question: Can I use onset date instead of reporting date in my analysis? Any recommendations on the interpretation of the results obtained when using onset date?

Underlying latent infections are reported at some rate (in reality varying over time but as a first pass can assume to be static). This rate should be added as an optional parameter to estimate_infections so that the true underlying infections can be estimated vs infections who are then reported (which is the current default).

Suggested by @sbfnk

Model fitting can take anywhere from a few minutes to an hour on a short timeseries. The data used in covid-rt-estimates highlight this.

This could be because the model is misspecified in some outbreaks (i.e due to high levels of imports) or because the data is not suitable to be used for estimation (i.e very few to no cases or extremely high reporting variance). Adding in imports (#39) may help account for the first of these issues and increasing data checks may help with the second.

I think documentation of epinow() function must be updated to state that the reported_cases parameter must be a data.table and not a data.frame as it is currently stated.
A check of that condition would be worth.
It took me hours to find where the

Error in `>=.default`(date, min(estimate$date, na.rm = TRUE)) : 
  comparison (5) is possible only for atomic and list types

comes from.

Thanks.

John Aponte

A SEIR modelling framework has not been able to repoduce the performance we saw using a branching process. This may point to a bug. A methodological issue or a input issue.

SEIR investigation is here: https://github.com/kgostic/epinow2_tests/tree/master/07-09-2020_constRt

Branching process investigation is here (out of date): github.com/epiforecasts/backcalc

I recently downloaded the development version and get a
'list' object cannot be coerced to type 'double' error from regional_epinow.

I believe this arises from the following test on regional_errors (regional_errors is a list of NULLs if there are no errors).

  if (length(regional_errors != 0)) {
    message("Runtime errors caught: ")
    print(regional_errors)
  }

Target completion date: 2nd October

• Resolve:

• #77,
• #76
• #68
• #67

• As CRAN check
• Review docs and news
• @seabbs -> devtools::release()