Real-time Rt Estimation, Forecasting and Reporting

Estimate Rt and cases by date of infection, forecast into the future, transform to date of report and then save summary measures and plots.

epinow(
  reported_cases,
  model,
  samples = 1000,
  stan_args,
  method = "exact",
  family = "negbin",
  generation_time,
  delays,
  gp = list(basis_prop = 0.3, boundary_scale = 2, lengthscale_mean = 0, lengthscale_sd
    = 2),
  rt_prior = list(mean = 1, sd = 1),
  horizon = 7,
  estimate_rt = TRUE,
  estimate_week_eff = TRUE,
  estimate_breakpoints = FALSE,
  burn_in = 0,
  stationary = FALSE,
  fixed = FALSE,
  fixed_future_rt = FALSE,
  future = FALSE,
  max_execution_time = Inf,
  prior_smoothing_window = 7,
  return_fit = FALSE,
  forecast_model,
  ensemble_type = "mean",
  return_estimates = TRUE,
  keep_samples = TRUE,
  make_plots = TRUE,
  target_folder,
  target_date,
  verbose = FALSE
)

Arguments

reported_cases	A data frame of confirmed cases (confirm) by date (date). confirm must be integer and date must be in date format.
model	A compiled stan model. By default uses the internal package model.
samples	Numeric, defaults to 1000. Number of samples post warmup.
stan_args	A list of stan arguments to be passed to `rstan::sampling` or `rstan::vb` (when using the "exact" or "approximate" method). For `method = approximate` an additional argument `trials` indicates the number of attempts to make using variational inference before returning an error (as stochastic failure is possible). The default for this is 5.
method	A character string defaults to "exact". Also accepts "approximate". Indicates the fitting method to be used this can either be "exact" (NUTs sampling) or "approximate" (variational inference). The exact approach returns samples from the posterior whilst the approximate method returns approximate samples. The approximate method is likely to return results several order of magnitudes faster than the exact method.
family	A character string indicating the reporting model to use. Defaults to negative binomial ("negbin") with poisson ("poisson") also supported.
generation_time	A list containing the mean, standard deviation of the mean (mean_sd), standard deviation (sd), standard deviation of the standard deviation and the maximum allowed value for the generation time (assuming a gamma distribution).
delays	A list of delays (i.e incubation period/reporting delay) between infection and report. Each list entry must also be a list containing the mean, standard deviation of the mean (mean_sd), standard deviation (sd), standard deviation of the standard deviation and the maximum allowed value for the that delay (assuming a lognormal distribution with all parameters excepting the max allowed value on the log scale).
gp	List controlling the Gaussian process approximation. Must contain the `basis_prop` (number of basis functions based on scaling the time points) which defaults to 0.3 and must be between 0 and 1 (increasing this increases the accuracy of the approximation and the cost of additional compute. Must also contain the `boundary_scale` (multiplied by half the range of the input time series). Increasing this increases the accuracy of the approximation at the cost of additional compute. See here: https://arxiv.org/abs/2004.11408 for more information on setting these parameters. Can optionally also contain the `lengthscale_mean` and `lengthscale_sd`. If these are specified this will override the defaults of 0 and 2 (normal distributed truncated at zero).
rt_prior	A list contain the mean and standard deviation (sd) of the gamma distributed prior for Rt. By default this is assumed to be mean 1 with a standard deviation of 1.
horizon	Numeric, defaults to 7. Number of days into the future to forecast.
estimate_rt	Logical, defaults TRUE. Should Rt be estimated when imputing infections.
estimate_week_eff	Logical, defaults TRUE. Should weekly reporting effects be estimated.
estimate_breakpoints	Logical, defaults to FALSE. Should breakpoints in Rt be estimated. If true then `reported_cases` must contain a `breakpoint` variable that is 1 on the dates with breakpoints and otherwise 0. Breakpoints are fit jointly with a global non-parametric effect and so represent a conservative estimate of breakpoint changes.
burn_in	Numeric, defaults to 0. The number of initial estimates to discard. This argument may be used to reduce spurious findings when running `estimate_infections` on a partial timeseries (as the earliest estimates will not be informed by all cases that occurred only those supplied to `estimate_infections`). The combined delays used will inform the appropriate length of this burn in but 7 days is likely a sensible starting point.
stationary	Logical, defaults to FALSE. Should Rt be estimated with a global mean. When estimating Rt this should substantially improve run times but will revert to the global average for real time and forecasted estimates. This setting is most appropriate when estimating historic Rt or when combined with breakpoints.
fixed	Logical, defaults to FALSE. If TRUE then a Gaussian process is not used and Rt is assumed to be constant over time (apart from any manually included breakpoints). If `estimate_rt` is FALSE then this reduces the backcalculation to a simple mean shift. This option can be used to produce a null model estimate, to produce a single Rt estimate for a short timeseries or as part of a wider analysis on the impact of interventions.
fixed_future_rt	Logical, defaults to FALSE. IF TRUE then the estimated Rt from the last time point with data is used for all future time points without data.
future	Logical, defaults to `FALSE`. Should stan chains be run in parallel using `future`. This allows users to have chains fail gracefully (i.e when combined with `max_execution_time`). Should be combined with a call to `future::plan`
max_execution_time	Numeric, defaults to Inf. If set will kill off processing of each chain if not finished within the specified timeout. When more than 2 chains finish successfully estimates will still be returned. If less than 2 chains return within the allowed time then estimation will fail with an informative error.
prior_smoothing_window	Numeric defaults to 7. The number of days over which to take a rolling average for the prior based on reported cases.
return_fit	Logical, defaults to FALSE. Should the fitted stan model be returned.
forecast_model	An uninitialised forecast model function to be passed to `EpiSoon::forecast_rt`. Used for forecasting future Rt and case co An example of the required structure is: `function(ss, y){bsts::AddSemilocalLinearTrend(ss, y = y)}`.
ensemble_type	Character string indicating the type of ensemble to use. By default this is an unweighted ensemble ("mean") with no other types currently supported.
return_estimates	Logical, defaults to TRUE. Should estimates be returned.
keep_samples	Logical, defaults to TRUE. Should samples be kept (either returned or saved or both depending on other settings).
make_plots	Logical, defaults to TRUE. Should plots be made and returned (or saved).
target_folder	Character string specifying where to save results (will create if not present).
target_date	Date, defaults to maximum found in the data if not specified.
verbose	Logical, defaults to `FALSE`. Should verbose progress messages be printed.

Value

A list of output from estimate_infections, forecast_infections, report_cases, and report_summary.

Examples

# \donttest{
## Construct example distributions
generation_time <- get_generation_time(disease = "SARS-CoV-2", source = "ganyani")
incubation_period <- get_incubation_period(disease = "SARS-CoV-2", source = "lauer")
reporting_delay <- EpiNow2::bootstrapped_dist_fit(rlnorm(100, log(6), 1), max_value = 30)

## Example case data
reported_cases <- EpiNow2::example_confirmed[1:40] 

## Report Rt along with forecasts
out <- epinow(reported_cases = reported_cases, generation_time = generation_time,
              delays = list(incubation_period, reporting_delay),
              stan_args = list(warmup = 200, cores = ifelse(interactive(), 4, 1)))

out
#> $estimates
#> $estimates$samples
#>                 variable        parameter time       date sample       value
#>      1:       infections       infections    1 2020-02-09      1    2.104578
#>      2:       infections       infections    2 2020-02-10      1    9.918045
#>      3:       infections       infections    3 2020-02-11      1   20.148745
#>      4:       infections       infections    4 2020-02-12      1   28.190642
#>      5:       infections       infections    5 2020-02-13      1   49.058274
#>     ---                                                                     
#> 214056: prior_infections prior_infections   56 2020-04-04      1 4607.088011
#> 214057: prior_infections prior_infections   57 2020-04-05      1 4562.019565
#> 214058: prior_infections prior_infections   58 2020-04-06      1 4517.391996
#> 214059: prior_infections prior_infections   59 2020-04-07      1 4473.200993
#> 214060: prior_infections prior_infections   60 2020-04-08      1 4429.442284
#>         strat     type
#>      1:  <NA> estimate
#>      2:  <NA> estimate
#>      3:  <NA> estimate
#>      4:  <NA> estimate
#>      5:  <NA> estimate
#>     ---               
#> 214056:  <NA> forecast
#> 214057:  <NA> forecast
#> 214058:  <NA> forecast
#> 214059:  <NA> forecast
#> 214060:  <NA> forecast
#> 
#> $estimates$summarised
#>            date       variable strat     type      bottom          top
#>   1: 2020-02-22              R  <NA> estimate    1.190113     1.865546
#>   2: 2020-02-23              R  <NA> estimate    1.331434     1.871549
#>   3: 2020-02-24              R  <NA> estimate    1.399129     1.816464
#>   4: 2020-02-25              R  <NA> estimate    1.464489     1.810915
#>   5: 2020-02-26              R  <NA> estimate    1.531088     1.828472
#>  ---                                                                  
#> 270: 2020-04-04 reported_cases  <NA> forecast  825.000000  9279.000000
#> 271: 2020-04-05 reported_cases  <NA> forecast 1077.000000 10940.000000
#> 272: 2020-04-06 reported_cases  <NA> forecast  618.000000 10933.000000
#> 273: 2020-04-07 reported_cases  <NA> forecast  579.000000 10074.000000
#> 274: 2020-04-08 reported_cases  <NA> forecast  476.000000  9825.000000
#>            lower       upper central_lower central_upper      median
#>   1:    1.399773    1.656284      1.439637      1.531139    1.529394
#>   2:    1.466085    1.674018      1.503920      1.580704    1.564417
#>   3:    1.490796    1.656797      1.566606      1.624793    1.600622
#>   4:    1.557524    1.687208      1.611329      1.657265    1.633066
#>   5:    1.582655    1.700511      1.651673      1.694200    1.665935
#>  ---                                                                
#> 270: 1680.000000 4571.000000   3126.000000   4131.000000 4255.500000
#> 271: 2020.000000 5418.000000   3666.000000   4809.000000 4761.000000
#> 272: 2231.000000 5715.000000   2299.000000   3451.000000 4577.500000
#> 273: 1332.000000 4231.000000   1866.000000   2940.000000 3898.000000
#> 274: 1259.000000 3781.000000   1381.000000   2242.000000 3332.500000
#>             mean           sd
#>   1:    1.539955 2.069887e-01
#>   2:    1.569134 1.671036e-01
#>   3:    1.601308 1.324607e-01
#>   4:    1.635035 1.056211e-01
#>   5:    1.668641 8.992658e-02
#>  ---                         
#> 270: 5047.191000 3.357203e+03
#> 271: 6003.977000 5.045123e+03
#> 272: 5920.518000 5.028147e+03
#> 273: 5537.903000 7.923904e+03
#> 274: 5150.376000 7.604268e+03
#> 
#> 
#> $estimated_reported_cases
#> $estimated_reported_cases$samples
#>              date sample cases  type
#>     1: 2020-02-22      1    53 gp_rt
#>     2: 2020-02-23      1   116 gp_rt
#>     3: 2020-02-24      1    36 gp_rt
#>     4: 2020-02-25      1   103 gp_rt
#>     5: 2020-02-26      1   156 gp_rt
#>    ---                              
#> 46996: 2020-04-04   1000  2364 gp_rt
#> 46997: 2020-04-05   1000  2248 gp_rt
#> 46998: 2020-04-06   1000  3128 gp_rt
#> 46999: 2020-04-07   1000  8209 gp_rt
#> 47000: 2020-04-08   1000  3898 gp_rt
#> 
#> $estimated_reported_cases$summarised
#>           date  type bottom   top lower upper central_lower central_upper
#>  1: 2020-02-22 gp_rt     18   100    35    67            46            57
#>  2: 2020-02-23 gp_rt     24   146    58   105            63            80
#>  3: 2020-02-24 gp_rt     37   198    62   126            89           111
#>  4: 2020-02-25 gp_rt     46   227    74   148           103           128
#>  5: 2020-02-26 gp_rt     41   230    83   151           119           144
#>  6: 2020-02-27 gp_rt     48   307    98   197           114           151
#>  7: 2020-02-28 gp_rt     89   477   166   317           168           223
#>  8: 2020-02-29 gp_rt     65   414   148   280           196           243
#>  9: 2020-03-01 gp_rt    131   645   178   378           293           363
#> 10: 2020-03-02 gp_rt    131   714   216   430           308           384
#> 11: 2020-03-03 gp_rt    132   744   233   477           264           351
#> 12: 2020-03-04 gp_rt    128   709   258   477           345           420
#> 13: 2020-03-05 gp_rt    157  1027   346   681           392           513
#> 14: 2020-03-06 gp_rt    278  1632   486  1002           675           860
#> 15: 2020-03-07 gp_rt    264  1459   511   994           652           833
#> 16: 2020-03-08 gp_rt    324  2031   701  1352           770           999
#> 17: 2020-03-09 gp_rt    512  2378   825  1563           936          1173
#> 18: 2020-03-10 gp_rt    536  2585   972  1741          1177          1443
#> 19: 2020-03-11 gp_rt    535  2470   796  1609          1070          1367
#> 20: 2020-03-12 gp_rt    613  3389  1225  2297          1155          1558
#> 21: 2020-03-13 gp_rt   1071  5014  1471  3003          2156          2682
#> 22: 2020-03-14 gp_rt    918  4513  1491  2837          1703          2177
#> 23: 2020-03-15 gp_rt    973  5645  1937  3681          2491          3115
#> 24: 2020-03-16 gp_rt   1151  6368  2063  4178          2374          3099
#> 25: 2020-03-17 gp_rt   1077  6311  2071  4094          2610          3273
#> 26: 2020-03-18 gp_rt   1365  6058  1882  3609          2563          3199
#> 27: 2020-03-19 gp_rt   1286  7301  2406  4682          2794          3580
#> 28: 2020-03-20 gp_rt   1880  9689  3691  6744          3828          4871
#> 29: 2020-03-21 gp_rt   1581  8268  2909  5400          3197          4059
#> 30: 2020-03-22 gp_rt   1709  9900  3481  6668          4086          5263
#> 31: 2020-03-23 gp_rt   1846 10288  3139  6367          5129          6367
#> 32: 2020-03-24 gp_rt   1511  8982  3188  6150          3167          4213
#> 33: 2020-03-25 gp_rt   1618  7836  2717  5125          3494          4354
#> 34: 2020-03-26 gp_rt   1601  9405  3103  5989          3959          4981
#> 35: 2020-03-27 gp_rt   2274 11857  4283  8057          4283          5588
#> 36: 2020-03-28 gp_rt   1795  9229  2708  5622          4380          5407
#> 37: 2020-03-29 gp_rt   1805 10743  3512  6752          3973          5085
#> 38: 2020-03-30 gp_rt   1607 10716  3254  6607          3923          5106
#> 39: 2020-03-31 gp_rt   1717  9463  3093  5900          3638          4627
#> 40: 2020-04-01 gp_rt   1233  8029  2356  4711          2583          3472
#> 41: 2020-04-02 gp_rt   1080  8654  2525  5217          2765          3732
#> 42: 2020-04-03 gp_rt   1284 11614  2695  6498          3179          4469
#> 43: 2020-04-04 gp_rt    825  9279  1680  4571          3126          4131
#> 44: 2020-04-05 gp_rt   1077 10940  2020  5418          3666          4809
#> 45: 2020-04-06 gp_rt    618 10933  2231  5715          2299          3451
#> 46: 2020-04-07 gp_rt    579 10074  1332  4231          1866          2940
#> 47: 2020-04-08 gp_rt    476  9825  1259  3781          1381          2242
#>           date  type bottom   top lower upper central_lower central_upper
#>     median     mean         sd
#>  1:   57.0   60.678   27.01218
#>  2:   85.0   91.670   41.47760
#>  3:  110.0  117.341   53.29979
#>  4:  126.0  133.783   61.07899
#>  5:  129.5  139.000   61.12140
#>  6:  175.0  187.931   86.16558
#>  7:  266.0  286.075  136.20108
#>  8:  238.0  258.636  121.47415
#>  9:  342.0  366.673  170.36690
#> 10:  381.5  418.444  195.69751
#> 11:  407.0  437.361  208.73650
#> 12:  409.0  443.529  203.49790
#> 13:  572.5  619.511  296.12467
#> 14:  826.0  912.488  462.79971
#> 15:  815.0  875.829  405.08167
#> 16: 1129.0 1243.271  603.95272
#> 17: 1256.5 1391.160  647.01949
#> 18: 1386.0 1482.637  667.17832
#> 19: 1405.0 1491.812  662.56437
#> 20: 1879.0 2023.209  926.20736
#> 21: 2673.0 2908.756 1325.18650
#> 22: 2450.0 2656.060 1182.45070
#> 23: 3126.5 3426.928 1567.43897
#> 24: 3565.5 3836.884 1749.45888
#> 25: 3551.5 3843.235 1748.11077
#> 26: 3190.5 3450.562 1560.67069
#> 27: 3966.5 4307.359 2080.74707
#> 28: 5374.0 5894.875 2760.88948
#> 29: 4449.0 4876.111 2295.83404
#> 30: 5536.5 5938.661 2823.16618
#> 31: 5646.0 6099.224 2929.93523
#> 32: 5062.5 5472.011 2576.84186
#> 33: 4369.5 4722.114 2087.23088
#> 34: 5164.5 5633.350 2640.13851
#> 35: 6378.0 7071.846 3409.32873
#> 36: 5038.0 5392.096 2535.73015
#> 37: 5671.0 6207.233 2934.55373
#> 38: 5559.5 6184.075 3218.14037
#> 39: 4873.0 5435.624 2607.27641
#> 40: 4071.5 4543.282 2385.49224
#> 41: 4505.5 5123.818 2956.89803
#> 42: 5743.0 6598.400 4043.86820
#> 43: 4255.5 5047.191 3357.20277
#> 44: 4761.0 6003.977 5045.12278
#> 45: 4577.5 5920.518 5028.14724
#> 46: 3898.0 5537.903 7923.90362
#> 47: 3332.5 5150.376 7604.26762
#>     median     mean         sd
#> 
#> 
#> $summary
#>                                  measure              estimate
#> 1: New confirmed cases by infection date    3746 (47 -- 12677)
#> 2:        Expected change in daily cases                Unsure
#> 3:            Effective reproduction no.      0.9 (0.3 -- 1.5)
#> 4:                        Rate of growth -0.04 (-0.22 -- 0.12)
#> 5:          Doubling/halving time (days)     -17 (5.6 -- -3.1)
#>     numeric_estimate
#> 1: <data.table[1x7]>
#> 2:              0.66
#> 3: <data.table[1x7]>
#> 4: <data.table[1x7]>
#> 5: <data.table[1x3]>
#> 
#> $plots
#> $plots$infections
#> 
#> $plots$reports
#> 
#> $plots$reff
#> 
#> $plots$growth_rate
#> 
#> $plots$summary
#> 
#> 

## For optional forecasting
if(requireNamespace("EpiSoon")){
   if(requireNamespace("forecastHybrid")){

   ## Report Rt along with forecasts
   out <- epinow(reported_cases = cases, samples = 1000, 
                 generation_time = generation_time,
                 delays = list(incubation_period, reporting_delay),
                 forecast_model = function(y, ...){
                   EpiSoon::forecastHybrid_model(
                     y = y[max(1, length(y) - 21):length(y)],
                     model_params = list(models = "aefz", weights = "equal"),
                     forecast_params = list(PI.combination = "mean"), ...)},
                  stan_args = list(warmup = 200, cores = ifelse(interactive(), 4, 1)))

out
   }
}
#> Loading required namespace: EpiSoon
#> Loading required namespace: forecastHybrid
#> Registered S3 method overwritten by 'quantmod':
#>   method            from
#>   as.zoo.data.frame zoo 
#> Error in data.table::setDT(reported_cases): object 'cases' not found

# }