Quick start
In the following section we give an overview of the simple use case
for epinow() and regional_epinow().
The first step to using the package is to load it as follows.
Reporting delays, incubation period and generation time
Distributions can be supplied in two ways. First, one can supplying
delay data to estimate_delay(), where a subsampled
bootstrapped lognormal will be fit to account for uncertainty in the
observed data without being biased by changes in incidence (see
?EpiNow2::estimate_delay()).
Second, one can specify predetermined delays with uncertainty using
the distribution functions such as Gamma or
Lognormal. An arbitrary number of delay distributions are
supported in dist_spec() with a common use case being an
incubation period followed by a reporting delay. For more information on
specifying distributions see (see
?EpiNow2::Distributions).
For example if data on the delay between onset and infection was
available we could fit a distribution to it, using
estimate_delay(), with appropriate uncertainty as follows
(note this is a synthetic example),
reporting_delay <- estimate_delay(
rlnorm(1000, log(2), 1),
max_value = 14, bootstraps = 1
)If data was not available we could instead specify an informed
estimate of the likely delay using the distribution functions
Gamma or LogNormal. To demonstrate, we choose
a lognormal distribution with mean 2, standard deviation 1 and a maximum
of 10. This is just an example and unlikely to apply in any
particular use case.
reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10)
reporting_delay## - lognormal distribution (max: 10):
## meanlog:
## 0.58
## sdlog:
## 0.47For the rest of this vignette, we will use inbuilt example literature estimates for the incubation period and generation time of Covid-19 (see here for the code that generates these estimates). These distributions are unlikely to be applicable for your use case. We strongly recommend investigating what might be the best distributions to use in any given use case.
example_generation_time## - gamma distribution (max: 14):
## shape:
## - normal distribution:
## mean:
## 1.4
## sd:
## 0.48
## rate:
## - normal distribution:
## mean:
## 0.38
## sd:
## 0.25
example_incubation_period## - lognormal distribution (max: 14):
## meanlog:
## - normal distribution:
## mean:
## 1.6
## sd:
## 0.064
## sdlog:
## - normal distribution:
## mean:
## 0.42
## sd:
## 0.069Now, to the functions.
epinow()
This function represents the core functionality of the package and includes results reporting, plotting, and optional saving. It requires a data frame of cases by date of report and the distributions defined above.
Load example case data from EpiNow2.
reported_cases <- example_confirmed[1:60]
head(reported_cases)## date confirm
## <Date> <num>
## 1: 2020-02-22 14
## 2: 2020-02-23 62
## 3: 2020-02-24 53
## 4: 2020-02-25 97
## 5: 2020-02-26 93
## 6: 2020-02-27 78Estimate cases by date of infection, the time-varying reproduction
number, the rate of growth, and forecast these estimates into the future
by 7 days. Summarise the posterior and return a summary table and plots
for reporting purposes. If a target_folder is supplied
results can be internally saved (with the option to also turn off
explicit returning of results). Here we use the default model
parameterisation that prioritises real-time performance over run-time or
other considerations. For other formulations see the documentation for
estimate_infections().
estimates <- epinow(
data = reported_cases,
generation_time = generation_time_opts(example_generation_time),
delays = delay_opts(example_incubation_period + reporting_delay),
rt = rt_opts(prior = list(mean = 2, sd = 0.2)),
stan = stan_opts(cores = 4, control = list(adapt_delta = 0.99)),
verbose = interactive()
)## DEBUG [2024-03-06 10:56:48] epinow: Running in exact mode for 2000 samples (across 4 chains each with a warm up of 250 iterations each) and 81 time steps of which 7 are a forecast
## WARN [2024-04-23 16:58:18] epinow: There were 1 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them. -
## WARN [2024-04-23 16:58:18] epinow: Examine the pairs() plot to diagnose sampling problems
## -
names(estimates)## [1] "estimates" "estimated_reported_cases" "summary"
## [4] "plots" "timing"Both summary measures and posterior samples are returned for all
parameters in an easily explored format which can be accessed using
summary. The default is to return a summary table of
estimates for key parameters at the latest date partially supported by
data.
| measure | estimate |
|---|---|
| New infections per day | 2200 (951 – 4666) |
| Expected change in daily reports | Likely decreasing |
| Effective reproduction no. | 0.87 (0.57 – 1.2) |
| Rate of growth | -0.035 (-0.16 – 0.067) |
| Doubling/halving time (days) | -20 (10 – -4.4) |
Summarised parameter estimates can also easily be returned, either filtered for a single parameter or for all parameters.
## date variable strat type median mean sd lower_90 lower_50 lower_20 upper_20
## <Date> <char> <char> <char> <num> <num> <num> <num> <num> <num> <num>
## 1: 2020-02-22 R <NA> estimate 2.291017 2.295305 0.1622730 2.038582 2.185297 2.246205 2.333835
## 2: 2020-02-23 R <NA> estimate 2.257591 2.260566 0.1377286 2.045927 2.162193 2.218197 2.296625
## 3: 2020-02-24 R <NA> estimate 2.221472 2.223189 0.1196990 2.038071 2.137169 2.187203 2.250281
## 4: 2020-02-25 R <NA> estimate 2.177615 2.183177 0.1076970 2.014396 2.106603 2.152258 2.206708
## 5: 2020-02-26 R <NA> estimate 2.135769 2.140651 0.1005877 1.983849 2.069109 2.110227 2.160111
## 6: 2020-02-27 R <NA> estimate 2.091666 2.095831 0.0968505 1.942662 2.027977 2.068433 2.115299
## upper_50 upper_90
## <num> <num>
## 1: 2.409435 2.570970
## 2: 2.356789 2.490317
## 3: 2.305457 2.422693
## 4: 2.256158 2.361948
## 5: 2.206701 2.312967
## 6: 2.158071 2.260933Reported cases are returned in a separate data frame in order to streamline the reporting of forecasts and for model evaluation.
## date type median mean sd lower_90 lower_50 lower_20 upper_20 upper_50 upper_90
## <Date> <char> <num> <num> <num> <num> <num> <num> <num> <num> <num>
## 1: 2020-02-22 gp_rt 77 77.9745 22.03483 46 61 71 82 92 116.00
## 2: 2020-02-23 gp_rt 87 90.1460 25.20663 54 73 82 93 105 135.00
## 3: 2020-02-24 gp_rt 87 89.2040 25.05848 53 72 81 93 104 135.05
## 4: 2020-02-25 gp_rt 80 82.0715 22.13329 49 67 75 86 95 120.00
## 5: 2020-02-26 gp_rt 79 82.0965 22.81924 50 66 74 85 96 123.00
## 6: 2020-02-27 gp_rt 106 109.0665 30.32398 66 88 99 114 127 164.00A range of plots are returned (with the single summary plot shown
below). These plots can also be generated using the following
plot method.
plot(estimates)
regional_epinow()
The regional_epinow() function runs the
epinow() function across multiple regions in an efficient
manner.
Define cases in multiple regions delineated by the region variable.
reported_cases <- data.table::rbindlist(list(
data.table::copy(reported_cases)[, region := "testland"],
reported_cases[, region := "realland"]
))
head(reported_cases)## date confirm region
## <Date> <num> <char>
## 1: 2020-02-22 14 testland
## 2: 2020-02-23 62 testland
## 3: 2020-02-24 53 testland
## 4: 2020-02-25 97 testland
## 5: 2020-02-26 93 testland
## 6: 2020-02-27 78 testlandCalling regional_epinow() runs the epinow()
on each region in turn (or in parallel depending on the settings used).
Here we switch to using a weekly random walk rather than the full
Gaussian process model giving us piecewise constant estimates by
week.
estimates <- regional_epinow(
data = reported_cases,
generation_time = generation_time_opts(example_generation_time),
delays = delay_opts(example_incubation_period + reporting_delay),
rt = rt_opts(prior = list(mean = 2, sd = 0.2), rw = 7),
gp = NULL,
stan = stan_opts(cores = 4, warmup = 250, samples = 1000)
)## INFO [2024-04-23 16:58:21] Producing following optional outputs: regions, summary, samples, plots, latest
## INFO [2024-04-23 16:58:21] Reporting estimates using data up to: 2020-04-21
## INFO [2024-04-23 16:58:21] No target directory specified so returning output
## INFO [2024-04-23 16:58:21] Producing estimates for: testland, realland
## INFO [2024-04-23 16:58:21] Regions excluded: none
## DEBUG [2024-04-23 16:58:21] testland: Running in exact mode for 1000 samples (across 4 chains each with a warm up of 250 iterations each) and 81 time steps of which 7 are a forecast
## INFO [2024-04-23 16:58:33] Completed estimates for: testland
## DEBUG [2024-04-23 16:58:33] realland: Running in exact mode for 1000 samples (across 4 chains each with a warm up of 250 iterations each) and 81 time steps of which 7 are a forecast
## INFO [2024-04-23 16:58:44] Completed estimates for: realland
## INFO [2024-04-23 16:58:44] Completed regional estimates
## INFO [2024-04-23 16:58:44] Regions with estimates: 2
## INFO [2024-04-23 16:58:44] Regions with runtime errors: 0
## INFO [2024-04-23 16:58:44] Producing summary
## INFO [2024-04-23 16:58:44] No summary directory specified so returning summary output
## INFO [2024-04-23 16:58:44] No target directory specified so returning timingsResults from each region are stored in a regional list
with across region summary measures and plots stored in a
summary list. All results can be set to be internally saved
by setting the target_folder and summary_dir
arguments. Each region can be estimated in parallel using the
future package (when in most scenarios cores
should be set to 1). For routine use each MCMC chain can also be run in
parallel (with future = TRUE) with a time out
(max_execution_time) allowing for partial results to be
returned if a subset of chains is running longer than expected. See the
documentation for the future package for details on
nested futures.
Summary measures that are returned include a table formatted for reporting (along with raw results for further processing). Futures updated will extend the S3 methods used above to smooth access to this output.
knitr::kable(estimates$summary$summarised_results$table)| Region | New infections per day | Expected change in daily reports | Effective reproduction no. | Rate of growth | Doubling/halving time (days) |
|---|---|---|---|---|---|
| realland | 2107 (1055 – 4430) | Likely decreasing | 0.86 (0.62 – 1.2) | -0.038 (-0.11 – 0.051) | -18 (14 – -6.3) |
| testland | 2121 (1035 – 4452) | Likely decreasing | 0.87 (0.61 – 1.2) | -0.037 (-0.11 – 0.05) | -19 (14 – -6.2) |
A range of plots are again returned (with the single summary plot shown below).
estimates$summary$summary_plot