add: how to reuse distribution delays to quantify the time-varying reproduction number (Rt)?

[avallecam]

from epiverse-trace/tutorials#104

# howto epinow2 -----------------------------------------------------------

library(epiparameter)
library(EpiNow2)
library(tidyverse)

# cases -------------------------------------------------------------------

example_confirmed

# delay: generation time --------------------------------------------------------------

covid_serialint <- 
  epiparameter::epidist_db(
    disease = "covid",
    epi_dist = "serial",
    author = "Nishiura",
    single_epidist = T
  )

covid_serialint

covid_serialint_discrete <- 
  epiparameter::discretise(covid_serialint)

covid_serialint_discrete_max <- 
  covid_serialint_discrete$prob_dist$q(p = 0.999)

serial_interval_covid <- 
  dist_spec(
    mean = covid_serialint$summary_stats$mean,
    sd = covid_serialint$summary_stats$sd,
    max = covid_serialint_discrete_max,
    distribution = "lognormal"
  )

# delay: incubation period -------------------------------------------------------------------

covid_incubation <- epiparameter::epidist_db(
  disease = "covid",
  epi_dist = "incubation",
  author = "Natalie",
  single_epidist = T
)

covid_incubation

covid_incubation_discrete <- epiparameter::discretise(covid_incubation)

incubation_time_covid <- dist_spec(
  mean = covid_incubation$summary_stats$mean,
  sd = covid_incubation$summary_stats$sd,
  max = covid_incubation_discrete$prob_dist$q(p = 0.999),
  distribution = "lognormal"
)

# epinow -------------------------------------------------------------------

epinow_estimates <- epinow(
  # cases
  reported_cases = example_confirmed[1:60],
  # delays
  generation_time = generation_time_opts(serial_interval_covid),
  delays = delay_opts(incubation_time_covid),
  # computation
  stan = stan_opts(
    cores = 4, samples = 1000, chains = 3,
    control = list(adapt_delta = 0.99)
  )
)

base::plot(epinow_estimates)

updated after #34 (comment)

[avallecam]

Thank you @seabbs for catching that up! This came from the source, so I also fixed it there epiverse-trace/tutorials@f14a2e2

[avallecam]

updated with the version in EpiNow2@dist-interface

epiverse-trace/episoap#130 (comment)

[avallecam]

update:

• we'll show first EpiNow2::dist_spec(mean,sd,max,distribution) then EpiNow2::dist_spec(pmf = ...).
  • EpiNow2::dist_spec(mean,sd,max,distribution) (or EpiNow2::Gamma() or EpiNow2::LogNormal()) works for fixed and uncertain distributions.
  • EpiNow2::dist_spec(pmf = ...) (or EpiNow2::pmf()) cover other distributions (like weibull) but only for fixed distributions.
• about the EpiNow2@dist-interface feature branch:
  • after it is merged, we'll keep the order of these two scenarios.
  • EpiNow2::LogNormal() or EpiNow2::Gamma() provide a safer interface than EpiNow2::dist_spec(). Their arguments avoid the misspecification of inputs to arguments as reported in this issue comment for {episoap}.
• This will be reassessed after having any <epidist> to <dist_spec> conversion method.

in detail:

There are two ways to get valid inputs to EpiNow2::epinow():

• in tutorials/case studies:
  • get parameters with epiparameter::get_parameters() from continuous distribution
  • get discretised distribution from continuous distribution
  • get a maximum value
  • plug them in to EpiNow2::dist_spec() or (EpiNow2::LogNormal() or EpiNow2::Gamma(), after branch is merged)
  • example in reprex below
• in episoap:
  • get discretised distribution from continuous distribution
  • get a maximum value
  • get sequence of quantiles with seq()
  • get density values with $d() from sequence of quantiles
  • plug them into EpiNow2::dist_spec(pmf = ...) (or EpiNow2::pmf(), after branch is merged)
  • example in this {episoap} commit and with dist-interface in a reprex comment for {episoap}

episoap gets benefit from distribution-independent workflow, making that second alternative more flexible for different inputs. But in learning materials, we'll prioritize the display of the distribution functions because it reduces one step to get to a goal, which includes the possibility of adding uncertainty, as written in epiverse-trace/epiparameter#218 (comment).

reprex:

# howto epinow2 -----------------------------------------------------------

library(epiparameter)
library(EpiNow2)
library(tidyverse)

# cases -------------------------------------------------------------------

example_confirmed # assumption: covid data

# delay: generation time --------------------------------------------------------------

covid_serialint <- 
  epiparameter::epidist_db(
    disease = "covid",
    epi_dist = "serial",
    author = "Nishiura",
    single_epidist = T
  )

covid_serialint

covid_serialint_parameters <- epiparameter::get_parameters(covid_serialint)

covid_serialint_parameters

covid_serialint_discrete <- 
  epiparameter::discretise(covid_serialint)

covid_serialint_discrete_max <- 
  covid_serialint_discrete$prob_dist$q(p = 0.999)

serial_interval_covid <- 
  dist_spec(
    mean = covid_serialint_parameters["meanlog"],
    sd = covid_serialint_parameters["sdlog"],
    max = covid_serialint_discrete_max, #must for epinow()
    distribution = "lognormal"
  )

serial_interval_covid

# delay: incubation period -------------------------------------------------------------------

covid_incubation <- epiparameter::epidist_db(
  disease = "covid",
  epi_dist = "incubation",
  author = "Natalie",
  single_epidist = T
)

covid_incubation

covid_incubation_parameters <- epiparameter::get_parameters(covid_incubation)

covid_incubation_parameters

covid_incubation_discrete <- epiparameter::discretise(covid_incubation)

incubation_time_covid <- dist_spec(
  mean = covid_incubation_parameters["meanlog"],
  sd = covid_incubation_parameters["sdlog"],
  max = covid_incubation_discrete$prob_dist$q(p = 0.999), #must for epinow()
  distribution = "lognormal"
)

incubation_time_covid

# epinow -------------------------------------------------------------------

epinow_estimates <- epinow(
  # cases
  reported_cases = example_confirmed[1:60],
  # delays
  generation_time = generation_time_opts(serial_interval_covid),
  delays = delay_opts(incubation_time_covid),
  # computation
  stan = stan_opts(
    cores = 4, samples = 1000, chains = 3,
    control = list(adapt_delta = 0.99)
  )
)

summary(epinow_estimates)
plot(epinow_estimates)

[avallecam]

reference using incidence2 + epiparameter + epinow2 WHO-Collaboratory/collaboratory-epiparameter-community#18