Use of {epiparameter} in Rt pipeline

[ntncmch]

Here is the R script by @annecori and @ntncmch that attempts to use the {epiparameter} package in a typical Rt estimation pipeline, using both {EpiEstim} and {EpiNow2}.

The 10 problems discovered are marked with # TODO-1,2,etc. It remains to be discussed whether these issues should be posted on {epiparameter} or in other packages.

# EpiParameter workshop
library(epiparameter)
library(epireview)
library(simulist)
library(here)
library(tidyverse)
library(EpiEstim)
library(incidence2)
library(EpiNow2)

# downloaded from https://royalsocietypublishing.org/doi/suppl/10.1098/rstb.2016.0308
data <- read_csv(here("data/epiparameter/rstb20160308supp1.csv"))

# select Sierra Leone
data_sle <- data %>% filter(Country == "Sierra Leone")

# show data
ggplot(data_sle, aes(x = DateOnset)) + geom_histogram(binwidth = 1)

# transform LL into incidence TS
inc_sle <-
  incidence2::incidence(data_sle, date_index = "DateOnset") %>%
  complete(
    date_index = seq(min(date_index), max(date_index), 1),
    fill = list(count = 0, count_variable = "DateOnset")
  ) # add missing dates with 0 counts
# TODO-1: it would be nice to have an option in incidence 
# to add missing 0 when interval is specified
# e.g. incidence2::incidence(data_sle, date_index = "DateOnset", interval = "day", complete = TRUE)

# get Ebola serial interval from EpiReview via Epiparameter

ebola <- epireview::load_epidata("ebola")
params <- ebola[["params"]]

# we use Faye study for our estimates
# this is an arbitrary choice, we looked for study with both mean and SD estimate
# TODO-2: would be nice to extract some consensus estimates
param_evd_faye <- params %>% filter(covidence_id == 4966)

# get serial interval (2 lines: mean + sd)
# TODO-3: in epireview mean and sd should be combined in a single 
# entry as for Zika (parameter 1, paramater 2)
serial_interval_faye <-
  param_evd_faye %>% filter(parameter_type == "Human delay - serial interval")

serial_interval_faye %>% select(
  pathogen,
  article_label,
  parameter_type,
  parameter_value_type,
  parameter_value,
  parameter_unit
)

# create epidist object from epireview
# several manual steps needed
# TODO-4: ideally we should get prob_dist directly 
# e.g. epiparameter:::as_epidist(serial_interval_faye, prob_dist = "gamma")
evd_SI_ep <- epiparameter:::as_epidist(serial_interval_faye)
evd_SI_ep_param <-
  convert_summary_stats_to_params("gamma",
                                  mean = evd_SI_ep$summary_stats$mean,
                                  sd = evd_SI_ep$summary_stats$sd)
evd_SI_ep$prob_dist <-
  distributional::dist_gamma(shape = evd_SI_ep_param$shape,
                             rate = 1 / evd_SI_ep_param$scale)
# TODO-5: bug in dist_gamma: error when inputing shape and scale only

is_parameterised(evd_SI_ep)
plot(evd_SI_ep, day_range = 0:50)
evd_SI_pmf <- density(discretize(evd_SI_ep), at = 0:50)


# EpiEstim example --------------------------------------------------------

# epiestim v1
# TODO-6: epiestim should accept incidence2 objects
res1 <- estimate_R(incid = inc_sle$count,
                   method = "non_parametric_si",
                   config = make_config(list(si_distr = evd_SI_pmf)))
plot(res1)

# epiestim v2
res2 <- estimate_R(incid = inc_sle$count,
                   method = "parametric",
                   config = make_config(
                     list(
                       mean_si = evd_SI_ep$summary_stats$mean,
                       std_si = evd_SI_ep$summary_stats$sd
                     )
                   ))
plot(res2)

# epiestim v3
# TODO-7: epiestim should be able to automatically create the make_config
# from an epidist object depending on the method chosen
p_std_mean_si <- 1 - (1 - evd_SI_ep$summary_stats$mean_ci / 100) / 2
std_mean_si <-
  diff(evd_SI_ep$summary_stats$mean_ci_limits) / qnorm(p = p_std_mean_si) /
  2

p_std_std_si <- 1 - (1 - evd_SI_ep$summary_stats$sd_ci / 100) / 2
std_std_si <-
  diff(evd_SI_ep$summary_stats$sd_ci_limits) / qnorm(p = p_std_std_si) /
  2

res3 <- estimate_R(incid = inc_sle$count,
                   method = "uncertain_si",
                   config = make_config(
                     list(
                       mean_si = evd_SI_ep$summary_stats$mean,
                       std_mean_si = std_mean_si,
                       min_mean_si = min(evd_SI_ep$summary_stats$mean_ci_limits),
                       max_mean_si = max(evd_SI_ep$summary_stats$mean_ci_limits),
                       std_si = evd_SI_ep$summary_stats$sd,
                       std_std_si = std_std_si,
                       min_std_si = min(evd_SI_ep$summary_stats$sd_ci_limits),
                       max_std_si = max(evd_SI_ep$summary_stats$sd_ci_limits),
                       n1 = 10,
                       n2 = 100
                     )
                   ))
plot(res3)
# ggsave(here("output/epiparameter/epiestime_v3.png"), width = 5, height = 7)

# EpiNow2 example ---------------------------------------------------------

# in addition to the generation time we need also the delay in outcome
# we use the serial interval as an approximation for the generation time
# we use the incubation period as an approximation for the delay in outcome
param_evd_faye_incub <-
  param_evd_faye %>% filter(parameter_type == "Human delay - incubation period")

evd_IP_ep <- epiparameter:::as_epidist(param_evd_faye_incub)
evd_IP_ep_param <-
  convert_summary_stats_to_params("gamma",
                                  mean = evd_IP_ep$summary_stats$mean,
                                  sd = evd_IP_ep$summary_stats$sd)
evd_IP_ep$prob_dist <-
  distributional::dist_gamma(shape = evd_IP_ep_param$shape,
                             rate = 1 / evd_IP_ep_param$scale)

# need to convert into epinow2 format
evd_IP_epinow <- EpiNow2::Gamma(
  mean = evd_IP_ep$summary_stats$mean,
  sd = evd_IP_ep$summary_stats$sd,
  max = 20
)
# TODO-8: write function that compute a max based on a 
# given percentile (e.g. max such that 99% of the density)
plot(evd_IP_epinow)

# need to convert into epinow2 format
evd_SI_epinow <- EpiNow2::Gamma(
  mean = evd_SI_ep$summary_stats$mean,
  sd = evd_SI_ep$summary_stats$sd,
  max = 30
)
plot(evd_SI_epinow)

# epinow2 don't take incidence2 objects
# need to rename "date" and "count" variables
# TODO-9: epinow accept incidence2 objects
inc_sle_epinow <-
  inc_sle %>% rename(date = date_index, confirm = count) %>% select(-count_variable) %>%
  filter(date > "2014-05-01", date < "2014-08-01") # reduce epidemic length due to long running time of epinow

ggplot(inc_sle_epinow, aes(x = date, y = confirm)) + geom_line()

## NOTE this is slow
# TODO-10: generation_time_opts and delay_opts should be able to
# take "epidist" object without converting them with Gamma()
R_epinow <- epinow(
  # cases
  data = inc_sle_epinow,
  # delays
  generation_time = generation_time_opts(evd_SI_epinow),
  # delays
  delays = delay_opts(evd_IP_epinow),
  # computation
  stan = stan_opts(
    cores = 4,
    samples = 1000,
    chains = 3,
    control = list(adapt_delta = 0.99)
  )
)

plot(R_epinow)
# ggsave(here("output/epiparameter/epinow2.png"), width = 5, height = 7)

[sbfnk]

This is great, thanks for sharing! Just a couple of comments relating to EpiNow2:

# TODO-8: write function that compute a max based on a 
# given percentile (e.g. max such that 99% of the density)

See epiforecasts/EpiNow2#533

# TODO-9: epinow accept incidence2 objects

We're thinking of making the column names more flexible in EpiNow2, which would address this to some degree. We could also add an S3 method to accept incidence2 objects but I'm not sure it's worth the overhead if all it does is rename some columns (happy to convinced otherwise).
epiforecasts/EpiNow2#505

## NOTE this is slow

This is largely a consequence of the settings used. We're hoping to clarify this in the future:
epiforecasts/EpiNow2#629

[TimTaylor]

TODO-1

# transform LL into incidence TS
inc_sle <-
  incidence2::incidence(data_sle, date_index = "DateOnset") %>%
  complete(
    date_index = seq(min(date_index), max(date_index), 1),
    fill = list(count = 0, count_variable = "DateOnset")
  ) # add missing dates with 0 counts
# TODO-1: it would be nice to have an option in incidence 
# to add missing 0 when interval is specified
# e.g. incidence2::incidence(data_sle, date_index = "DateOnset", interval = "day", complete = TRUE)

currently you can do

inc_sle <-
  incidence2::incidence(data_sle, date_index = "DateOnset") %>%
  incidence2::complete_dates()

to achieve the same.

That said , this has been requested before though so I will add a complete_dates option to incidence itself and get it in to the next release. There are some complications due to people using POSIXct instead of dates that I need to account for but should hopefully get something out in the next couple of weeks.

TODO-9

# TODO-9: epinow accept incidence2 objects

I've various thoughts here but I currently lean towards one of:

• having it in incidence2 itself (i.e. useful EpiNow2 methods brought in via suggested packages)
• a package which pulls together incidence2 and other useful packages and combines them.

Either way, happy to help out here as needed.

[TimTaylor]

TODO-6

# epiestim v1
# TODO-6: epiestim should accept incidence2 objects

Similar thoughts to above. Happy to work this through with you either way.

[TimTaylor]

@annecori , @ntncmch - I've just release incidence 2.3.0 to CRAN (binaries waiting to be built but source there) and you can now do:

inc_sle <-
    incidence2::incidence(data_sle, date_index = "DateOnset", complete_dates = TRUE)

hopefully that crosses off TODO-1.

[ntncmch]

This is great, thanks! I was not aware of the function incidence2::complete_dates but I agree it's nice to have the option directly in incidence2::incidence.

[joshwlambert]

Thanks for contributing this @ntncmch and @annecori, and for pointing out the areas where improvements were needed. It's also great to see other package developers incorporating the suggestions.

TODO-4

TODO-4: ideally we should get prob_dist directly
e.g. epiparameter:::as_epidist(serial_interval_faye, prob_dist = "gamma")

This has been addressed in epiverse-trace/epiparameter#335 and is available by reinstalling the package.

Here is a run through with the updates:

library(epiparameter)
library(epireview)
#> Loading required package: epitrix
#> Loading required package: ggplot2
#> Loading required package: ggforce

library(tidyverse)

# get Ebola serial interval from EpiReview via Epiparameter

ebola <- epireview::load_epidata("ebola")
#> Warning: One or more parsing issues, call `problems()` on your data frame for details,
#> e.g.:
#>   dat <- vroom(...)
#>   problems(dat)
#> Warning in load_epidata_raw(pathogen, "outbreak"): No data found for ebola
#> Warning: One or more parsing issues, call `problems()` on your data frame for details,
#> e.g.:
#>   dat <- vroom(...)
#>   problems(dat)
#> Warning in epireview::load_epidata("ebola"): No outbreaks information found for
#> ebola
#> Data loaded for ebola

params <- ebola[["params"]]

# we use Faye study for our estimates
# this is an arbitrary choice, we looked for study with both mean and SD estimate
param_evd_faye <- params %>% filter(covidence_id == 4966)

# get serial interval (2 lines: mean + sd)
serial_interval_faye <-
  param_evd_faye %>% filter(parameter_type == "Human delay - serial interval")

serial_interval_faye %>% select(
  pathogen,
  article_label,
  parameter_type,
  parameter_value_type,
  parameter_value,
  parameter_unit
)
#> # A tibble: 2 × 6
#>   pathogen    article_label  parameter_type parameter_value_type parameter_value
#>   <chr>       <chr>          <chr>          <chr>                          <dbl>
#> 1 Ebola virus Faye 2015 (a)… Human delay -… Mean                            14.2
#> 2 Ebola virus Faye 2015 (a)… Human delay -… Standard Deviation               7.1
#> # ℹ 1 more variable: parameter_unit <chr>

# create epidist object from epireview
evd_SI_ep <- epiparameter:::as_epidist(serial_interval_faye, prob_dist = "gamma")
#> Using Faye (2015). "<title not available>." _<journal not available>_. 
#> To retrieve the citation use the 'get_citation' function
#> Warning: Cannot create full citation for epidemiological parameters without bibliographic information 
#>  see ?as_epidist for help.

is_parameterised(evd_SI_ep)
#> [1] TRUE

plot(evd_SI_ep)

evd_SI_pmf <- density(discretize(evd_SI_ep), at = 0:50)

^{Created on 2024-06-28 with reprex v2.1.0}