From 457d4acb12a49c2e0bac7e1cde659aeb9afd7839 Mon Sep 17 00:00:00 2001 From: jamesaazam Date: Thu, 12 Dec 2024 13:55:37 +0000 Subject: [PATCH 01/18] Use natural parameters in all examples --- R/epinow.R | 6 +++++- R/estimate_infections.R | 6 +++++- R/regional_epinow.R | 15 ++++++++++----- R/simulate_infections.R | 7 ++++++- man/epinow.Rd | 6 +++++- man/estimate_infections.Rd | 6 +++++- man/forecast_infections.Rd | 7 ++++++- man/regional_epinow.Rd | 15 ++++++++++----- vignettes/EpiNow2.Rmd.orig | 19 ++++++++++++++++--- vignettes/epinow.Rmd.orig | 11 +++++++++-- .../estimate_infections_options.Rmd.orig | 15 +++++++++++---- .../estimate_infections_workflow.Rmd.orig | 6 +++--- 12 files changed, 91 insertions(+), 28 deletions(-) diff --git a/R/epinow.R b/R/epinow.R index 6ca1fd078..b9edc5538 100644 --- a/R/epinow.R +++ b/R/epinow.R @@ -56,7 +56,11 @@ #' ) #' # set an example reporting delay. In practice this should use an estimate #' # from the literature or be estimated from data -#' reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10) +#' reporting_delay <- LogNormal( +#' meanlog = Fixed(2), +#' sdlog = Fixed(1), +#' max = 10 +#' ) #' #' # example case data #' reported_cases <- example_confirmed[1:40] diff --git a/R/estimate_infections.R b/R/estimate_infections.R index 1d5b5f004..da59d7156 100644 --- a/R/estimate_infections.R +++ b/R/estimate_infections.R @@ -101,7 +101,11 @@ #' ) #' # set an example reporting delay. In practice this should use an estimate #' # from the literature or be estimated from data -#' reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10) +#' reporting_delay <- LogNormal( +#' meanlog = Fixed(2), +#' sdlog = Fixed(1), +#' max = 10 +#' ) #' #' # for more examples, see the "estimate_infections examples" vignette #' def <- estimate_infections(reported_cases, diff --git a/R/regional_epinow.R b/R/regional_epinow.R index a5f08d062..e88e8ebf7 100644 --- a/R/regional_epinow.R +++ b/R/regional_epinow.R @@ -79,11 +79,16 @@ #' data = cases, #' generation_time = gt_opts(example_generation_time), #' delays = delay_opts(example_incubation_period + example_reporting_delay), -#' rt = rt_opts(prior = LogNormal(mean = 2, sd = 0.2)), -#' stan = stan_opts( -#' samples = 100, warmup = 200 -#' ), -#' verbose = interactive() +#' rt = rt_opts( +#' prior = LogNormal( +#' meanlog = Fixed(2), +#' sdlog = Fixed(0.2) +#' ), +#' stan = stan_opts( +#' samples = 100, warmup = 200 +#' ), +#' verbose = interactive() +#' ) #' ) #' options(old_opts) #' } diff --git a/R/simulate_infections.R b/R/simulate_infections.R index f3b66cc3f..beebe3cc5 100644 --- a/R/simulate_infections.R +++ b/R/simulate_infections.R @@ -277,7 +277,12 @@ simulate_infections <- function(estimates, R, initial_infections, #' est <- estimate_infections(reported_cases, #' generation_time = generation_time_opts(example_generation_time), #' delays = delay_opts(example_incubation_period + example_reporting_delay), -#' rt = rt_opts(prior = LogNormal(mean = 2, sd = 0.1), rw = 7), +#' rt = rt_opts(prior = LogNormal( +#' meanlog = Fixed(2), +#' sdlog = Fixed(0.1) +#' ), +#' rw = 7 +#' ), #' obs = obs_opts(scale = Normal(mean = 0.1, sd = 0.01)), #' gp = NULL, horizon = 0 #' ) diff --git a/man/epinow.Rd b/man/epinow.Rd index bfe6dc102..c1973bdb7 100644 --- a/man/epinow.Rd +++ b/man/epinow.Rd @@ -157,7 +157,11 @@ incubation_period <- LogNormal( ) # set an example reporting delay. In practice this should use an estimate # from the literature or be estimated from data -reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10) +reporting_delay <- LogNormal( + meanlog = Fixed(2), + sdlog = Fixed(1), + max = 10 +) # example case data reported_cases <- example_confirmed[1:40] diff --git a/man/estimate_infections.Rd b/man/estimate_infections.Rd index 0372966a0..e65d0f6d3 100644 --- a/man/estimate_infections.Rd +++ b/man/estimate_infections.Rd @@ -146,7 +146,11 @@ incubation_period <- LogNormal( ) # set an example reporting delay. In practice this should use an estimate # from the literature or be estimated from data -reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10) +reporting_delay <- LogNormal( + meanlog = Fixed(2), + sdlog = Fixed(1), + max = 10 +) # for more examples, see the "estimate_infections examples" vignette def <- estimate_infections(reported_cases, diff --git a/man/forecast_infections.Rd b/man/forecast_infections.Rd index 24c5b2c5d..23ceac66d 100644 --- a/man/forecast_infections.Rd +++ b/man/forecast_infections.Rd @@ -65,7 +65,12 @@ reported_cases <- example_confirmed[1:50] est <- estimate_infections(reported_cases, generation_time = generation_time_opts(example_generation_time), delays = delay_opts(example_incubation_period + example_reporting_delay), - rt = rt_opts(prior = LogNormal(mean = 2, sd = 0.1), rw = 7), + rt = rt_opts(prior = LogNormal( + meanlog = Fixed(2), + sdlog = Fixed(0.1) + ), + rw = 7 + ), obs = obs_opts(scale = Normal(mean = 0.1, sd = 0.01)), gp = NULL, horizon = 0 ) diff --git a/man/regional_epinow.Rd b/man/regional_epinow.Rd index 6eca0c83a..20b9ee734 100644 --- a/man/regional_epinow.Rd +++ b/man/regional_epinow.Rd @@ -156,11 +156,16 @@ def <- regional_epinow( data = cases, generation_time = gt_opts(example_generation_time), delays = delay_opts(example_incubation_period + example_reporting_delay), - rt = rt_opts(prior = LogNormal(mean = 2, sd = 0.2)), - stan = stan_opts( - samples = 100, warmup = 200 - ), - verbose = interactive() + rt = rt_opts( + prior = LogNormal( + meanlog = Fixed(2), + sdlog = Fixed(0.2) + ), + stan = stan_opts( + samples = 100, warmup = 200 + ), + verbose = interactive() + ) ) options(old_opts) } diff --git a/vignettes/EpiNow2.Rmd.orig b/vignettes/EpiNow2.Rmd.orig index 69d10bcdd..f2df29bc8 100644 --- a/vignettes/EpiNow2.Rmd.orig +++ b/vignettes/EpiNow2.Rmd.orig @@ -50,7 +50,11 @@ If data was not available we could instead specify an informed estimate of the l 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*. ```{r} -reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10) +reporting_delay <- LogNormal( + meanlog = Fixed(2), + sdlog = Fixed(1), + max = 10 +) reporting_delay ``` @@ -94,7 +98,10 @@ estimates <- epinow( data = reported_cases, generation_time = gt_opts(example_generation_time), delays = delay_opts(example_incubation_period + reporting_delay), - rt = rt_opts(prior = LogNormal(mean = 2, sd = 0.2)), + rt = rt_opts(prior = LogNormal( + meanlog = Fixed(2), + sdlog = Fixed(0.2) + )), stan = stan_opts(cores = 4, control = list(adapt_delta = 0.99)), verbose = interactive() ) @@ -148,7 +155,13 @@ estimates <- regional_epinow( data = reported_cases, generation_time = gt_opts(example_generation_time), delays = delay_opts(example_incubation_period + reporting_delay), - rt = rt_opts(prior = LogNormal(mean = 2, sd = 0.2), rw = 7), + rt = rt_opts( + prior = LogNormal( + meanlog = Fixed(2), + sdlog = Fixed(0.2) + ), + rw = 7 + ), gp = NULL, stan = stan_opts(cores = 4, warmup = 250, samples = 1000) ) diff --git a/vignettes/epinow.Rmd.orig b/vignettes/epinow.Rmd.orig index d8960bfbd..ab7726db6 100644 --- a/vignettes/epinow.Rmd.orig +++ b/vignettes/epinow.Rmd.orig @@ -38,9 +38,16 @@ This should be replaced with parameters relevant to the system that is being stu library("EpiNow2") options(mc.cores = 4) reported_cases <- example_confirmed[1:60] -reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10) +reporting_delay <- LogNormal( + meanlog = Fixed(2), + sdlog = Fixed(1), + max = 10 +) delay <- example_incubation_period + reporting_delay -rt_prior <- LogNormal(mean = 2, sd = 0.1) +rt_prior <- LogNormal( + meanlog = Fixed(2), + sdlog = Fixed(0.1) +) ``` We can then run the `epinow()` function with the same arguments as `estimate_infections()`. diff --git a/vignettes/estimate_infections_options.Rmd.orig b/vignettes/estimate_infections_options.Rmd.orig index 25706703e..08067ca2d 100644 --- a/vignettes/estimate_infections_options.Rmd.orig +++ b/vignettes/estimate_infections_options.Rmd.orig @@ -73,7 +73,11 @@ For the reporting delay, we use a lognormal distribution with mean of 2 days and Note that the mean and standard deviation must be converted to the log scale, which can be done using the `convert_log_logmean()` function. ```{r reporting_delay} -reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10) +reporting_delay <- LogNormal( + meanlog = Fixed(2), + sdlog = Fixed(1), + max = 10 +) reporting_delay ``` @@ -97,7 +101,10 @@ example_generation_time Lastly we need to choose a prior for the initial value of the reproduction number. This is assumed by the model to be normally distributed and we can set the mean and the standard deviation. We decide to set the mean to 2 and the standard deviation to 1. ```{r initial_r} -rt_prior <- LogNormal(mean = 2, sd = 0.1) +rt_prior <- LogNormal( + meanlog = Fixed(2), + sdlog = Fixed(0.1) +) ``` # Running the model @@ -172,8 +179,8 @@ Here, instead of doing so we assume that we know about truncation with mean of 1 ```{r define_truncation} trunc_dist <- LogNormal( - mean = Normal(0.5, 0.1), - sd = Normal(0.5, 0.1), + meanlog = Normal(0.5, 0.1), + sdlog = Normal(0.5, 0.1), max = 3 ) trunc_dist diff --git a/vignettes/estimate_infections_workflow.Rmd.orig b/vignettes/estimate_infections_workflow.Rmd.orig index 2c8f48a19..59485840a 100644 --- a/vignettes/estimate_infections_workflow.Rmd.orig +++ b/vignettes/estimate_infections_workflow.Rmd.orig @@ -79,7 +79,7 @@ In all cases, the distributions given can be *fixed* (i.e. have no uncertainty) For example, to define a fixed gamma distribution with mean 3, standard deviation (sd) 1 and maximum value 10, you would write ```{r} -fixed_gamma <- Gamma(mean = 3, sd = 1, max = 10) +fixed_gamma <- Gamma(shape = 3, rate = 1, max = 10) fixed_gamma ``` @@ -93,7 +93,7 @@ If distributions are variable, the values with uncertainty are treated as [prior For example, to define a variable gamma distribution where uncertainty in the mean is given by a normal distribution with mean 3 and sd 2, and uncertainty in the standard deviation is given by a normal distribution with mean 1 and sd 0.1, with a maximum value 10, you would write ```{r} -uncertain_gamma <- Gamma(mean = Normal(3, 2), sd = Normal(1, 0.1), max = 10) +uncertain_gamma <- Gamma(shape = Normal(3, 2), rate = Normal(1, 0.1), max = 10) uncertain_gamma ``` @@ -209,7 +209,7 @@ It can be changed using the `rt_opts()` function. For example, if the user believes that at the very start of the data the reproduction number was 2, with uncertainty in this belief represented by a standard deviation of 1, they would use ```{r results = 'hide'} -rt_prior <- LogNormal(mean = 2, sd = 1) +rt_prior <- LogNormal(meanlog = Fixed(2), sdlog = Fixed(1)) rt_opts(prior = rt_prior) ``` From 9b9dbf176e4b5fe4690529327c79e4a1105599db Mon Sep 17 00:00:00 2001 From: jamesaazam Date: Thu, 12 Dec 2024 13:58:43 +0000 Subject: [PATCH 02/18] Add NEWS item --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index d8edd70b6..8b3c46aeb 100644 --- a/NEWS.md +++ b/NEWS.md @@ -29,6 +29,7 @@ - Brought the docs on `alpha_sd` up to date with the code change from prior PR #853. By @zsusswein in #862 and reviewed by @jamesmbaazam. - The `...` argument in `estimate_secondary()` has been removed because it was not used. By @jamesmbaazam in #894 and reviewed by @. +- All examples now use the natural parameters of distributions rather than the mean and standard deviation to remove warnings and encourage best practice. By @jamesmbaazam in #893 and reviewed by @. # EpiNow2 1.6.1 From ed39b8de40ea2396fbd87470eec121616bd2083f Mon Sep 17 00:00:00 2001 From: jamesaazam Date: Thu, 12 Dec 2024 14:48:47 +0000 Subject: [PATCH 03/18] Remove statement about warning --- vignettes/estimate_infections_workflow.Rmd.orig | 4 ---- 1 file changed, 4 deletions(-) diff --git a/vignettes/estimate_infections_workflow.Rmd.orig b/vignettes/estimate_infections_workflow.Rmd.orig index 59485840a..a250ae32d 100644 --- a/vignettes/estimate_infections_workflow.Rmd.orig +++ b/vignettes/estimate_infections_workflow.Rmd.orig @@ -102,10 +102,6 @@ which looks like this when plotted ```{r plot_uncertain_gamma} plot(uncertain_gamma) ``` - -Note the warning about parameters. -We used the mean and standard deviation to define this distribution with uncertain parameters, but it would be better to use the "natural" parameters of the gamma distribution, shape and rate, for example using the values estimate and reported back after calling the previous command. - There are various ways the specific delay distributions mentioned below might be obtained. Often, they will come directly from the existing literature reviewed by the user and studies conducted elsewhere. Sometimes it might be possible to obtain them from existing databases, e.g. using the [epiparameter](https://github.com/epiverse-trace/epiparameter) R package. From 19be97ebb89e99ccc261261971c593e7e008673c Mon Sep 17 00:00:00 2001 From: jamesaazam Date: Thu, 12 Dec 2024 15:13:39 +0000 Subject: [PATCH 04/18] Remove use of Fixed() for specifying fixed parameters --- R/epinow.R | 4 ++-- R/estimate_infections.R | 4 ++-- R/regional_epinow.R | 4 ++-- R/simulate_infections.R | 4 ++-- vignettes/EpiNow2.Rmd.orig | 12 ++++++------ vignettes/epinow.Rmd.orig | 8 ++++---- vignettes/estimate_infections_options.Rmd.orig | 8 ++++---- vignettes/estimate_infections_workflow.Rmd.orig | 2 +- 8 files changed, 23 insertions(+), 23 deletions(-) diff --git a/R/epinow.R b/R/epinow.R index b9edc5538..e90bb0fa8 100644 --- a/R/epinow.R +++ b/R/epinow.R @@ -57,8 +57,8 @@ #' # set an example reporting delay. In practice this should use an estimate #' # from the literature or be estimated from data #' reporting_delay <- LogNormal( -#' meanlog = Fixed(2), -#' sdlog = Fixed(1), +#' meanlog = 2, +#' sdlog = 1, #' max = 10 #' ) #' diff --git a/R/estimate_infections.R b/R/estimate_infections.R index da59d7156..d59e0b988 100644 --- a/R/estimate_infections.R +++ b/R/estimate_infections.R @@ -102,8 +102,8 @@ #' # set an example reporting delay. In practice this should use an estimate #' # from the literature or be estimated from data #' reporting_delay <- LogNormal( -#' meanlog = Fixed(2), -#' sdlog = Fixed(1), +#' meanlog = 2, +#' sdlog = 1, #' max = 10 #' ) #' diff --git a/R/regional_epinow.R b/R/regional_epinow.R index e88e8ebf7..38d001c98 100644 --- a/R/regional_epinow.R +++ b/R/regional_epinow.R @@ -81,8 +81,8 @@ #' delays = delay_opts(example_incubation_period + example_reporting_delay), #' rt = rt_opts( #' prior = LogNormal( -#' meanlog = Fixed(2), -#' sdlog = Fixed(0.2) +#' meanlog = 2, +#' sdlog = 0.2 #' ), #' stan = stan_opts( #' samples = 100, warmup = 200 diff --git a/R/simulate_infections.R b/R/simulate_infections.R index beebe3cc5..ac43b1400 100644 --- a/R/simulate_infections.R +++ b/R/simulate_infections.R @@ -278,8 +278,8 @@ simulate_infections <- function(estimates, R, initial_infections, #' generation_time = generation_time_opts(example_generation_time), #' delays = delay_opts(example_incubation_period + example_reporting_delay), #' rt = rt_opts(prior = LogNormal( -#' meanlog = Fixed(2), -#' sdlog = Fixed(0.1) +#' meanlog = 2, +#' sdlog = 0.1 #' ), #' rw = 7 #' ), diff --git a/vignettes/EpiNow2.Rmd.orig b/vignettes/EpiNow2.Rmd.orig index f2df29bc8..23cfe67a8 100644 --- a/vignettes/EpiNow2.Rmd.orig +++ b/vignettes/EpiNow2.Rmd.orig @@ -51,8 +51,8 @@ To demonstrate, we choose a lognormal distribution with mean 2, standard deviati ```{r} reporting_delay <- LogNormal( - meanlog = Fixed(2), - sdlog = Fixed(1), + meanlog = 2, + sdlog = 1, max = 10 ) reporting_delay @@ -99,8 +99,8 @@ estimates <- epinow( generation_time = gt_opts(example_generation_time), delays = delay_opts(example_incubation_period + reporting_delay), rt = rt_opts(prior = LogNormal( - meanlog = Fixed(2), - sdlog = Fixed(0.2) + meanlog = 2, + sdlog = 0.2 )), stan = stan_opts(cores = 4, control = list(adapt_delta = 0.99)), verbose = interactive() @@ -157,8 +157,8 @@ estimates <- regional_epinow( delays = delay_opts(example_incubation_period + reporting_delay), rt = rt_opts( prior = LogNormal( - meanlog = Fixed(2), - sdlog = Fixed(0.2) + meanlog = 2, + sdlog = 0.2 ), rw = 7 ), diff --git a/vignettes/epinow.Rmd.orig b/vignettes/epinow.Rmd.orig index ab7726db6..10d79f032 100644 --- a/vignettes/epinow.Rmd.orig +++ b/vignettes/epinow.Rmd.orig @@ -39,14 +39,14 @@ library("EpiNow2") options(mc.cores = 4) reported_cases <- example_confirmed[1:60] reporting_delay <- LogNormal( - meanlog = Fixed(2), - sdlog = Fixed(1), + meanlog = 2, + sdlog = 1, max = 10 ) delay <- example_incubation_period + reporting_delay rt_prior <- LogNormal( - meanlog = Fixed(2), - sdlog = Fixed(0.1) + meanlog = 2, + sdlog = 0.1 ) ``` diff --git a/vignettes/estimate_infections_options.Rmd.orig b/vignettes/estimate_infections_options.Rmd.orig index 08067ca2d..068ffad34 100644 --- a/vignettes/estimate_infections_options.Rmd.orig +++ b/vignettes/estimate_infections_options.Rmd.orig @@ -74,8 +74,8 @@ Note that the mean and standard deviation must be converted to the log scale, wh ```{r reporting_delay} reporting_delay <- LogNormal( - meanlog = Fixed(2), - sdlog = Fixed(1), + meanlog = 2, + sdlog = 1, max = 10 ) reporting_delay @@ -102,8 +102,8 @@ Lastly we need to choose a prior for the initial value of the reproduction numbe ```{r initial_r} rt_prior <- LogNormal( - meanlog = Fixed(2), - sdlog = Fixed(0.1) + meanlog = 2, + sdlog = 0.1 ) ``` diff --git a/vignettes/estimate_infections_workflow.Rmd.orig b/vignettes/estimate_infections_workflow.Rmd.orig index a250ae32d..e4e294bc2 100644 --- a/vignettes/estimate_infections_workflow.Rmd.orig +++ b/vignettes/estimate_infections_workflow.Rmd.orig @@ -205,7 +205,7 @@ It can be changed using the `rt_opts()` function. For example, if the user believes that at the very start of the data the reproduction number was 2, with uncertainty in this belief represented by a standard deviation of 1, they would use ```{r results = 'hide'} -rt_prior <- LogNormal(meanlog = Fixed(2), sdlog = Fixed(1)) +rt_prior <- LogNormal(meanlog = 2, sdlog = 1) rt_opts(prior = rt_prior) ``` From 36cff84d61eb7ed79b954262d5ef02feec665a23 Mon Sep 17 00:00:00 2001 From: jamesaazam Date: Thu, 12 Dec 2024 15:15:44 +0000 Subject: [PATCH 05/18] Redocument --- man/epinow.Rd | 4 ++-- man/estimate_infections.Rd | 4 ++-- man/forecast_infections.Rd | 4 ++-- man/regional_epinow.Rd | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/man/epinow.Rd b/man/epinow.Rd index c1973bdb7..40b295a02 100644 --- a/man/epinow.Rd +++ b/man/epinow.Rd @@ -158,8 +158,8 @@ incubation_period <- LogNormal( # set an example reporting delay. In practice this should use an estimate # from the literature or be estimated from data reporting_delay <- LogNormal( - meanlog = Fixed(2), - sdlog = Fixed(1), + meanlog = 2, + sdlog = 1, max = 10 ) diff --git a/man/estimate_infections.Rd b/man/estimate_infections.Rd index e65d0f6d3..0c1686397 100644 --- a/man/estimate_infections.Rd +++ b/man/estimate_infections.Rd @@ -147,8 +147,8 @@ incubation_period <- LogNormal( # set an example reporting delay. In practice this should use an estimate # from the literature or be estimated from data reporting_delay <- LogNormal( - meanlog = Fixed(2), - sdlog = Fixed(1), + meanlog = 2, + sdlog = 1, max = 10 ) diff --git a/man/forecast_infections.Rd b/man/forecast_infections.Rd index 23ceac66d..b836dc323 100644 --- a/man/forecast_infections.Rd +++ b/man/forecast_infections.Rd @@ -66,8 +66,8 @@ est <- estimate_infections(reported_cases, generation_time = generation_time_opts(example_generation_time), delays = delay_opts(example_incubation_period + example_reporting_delay), rt = rt_opts(prior = LogNormal( - meanlog = Fixed(2), - sdlog = Fixed(0.1) + meanlog = 2, + sdlog = 0.1 ), rw = 7 ), diff --git a/man/regional_epinow.Rd b/man/regional_epinow.Rd index 20b9ee734..ecf613441 100644 --- a/man/regional_epinow.Rd +++ b/man/regional_epinow.Rd @@ -158,8 +158,8 @@ def <- regional_epinow( delays = delay_opts(example_incubation_period + example_reporting_delay), rt = rt_opts( prior = LogNormal( - meanlog = Fixed(2), - sdlog = Fixed(0.2) + meanlog = 2, + sdlog = 0.2 ), stan = stan_opts( samples = 100, warmup = 200 From 172743a95bcdfde34a4e7865fe2631a9308a395b Mon Sep 17 00:00:00 2001 From: jamesaazam Date: Thu, 12 Dec 2024 17:27:23 +0000 Subject: [PATCH 06/18] Fix missing comma --- R/regional_epinow.R | 12 ++++++------ man/regional_epinow.Rd | 12 ++++++------ 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/R/regional_epinow.R b/R/regional_epinow.R index 38d001c98..5df23def5 100644 --- a/R/regional_epinow.R +++ b/R/regional_epinow.R @@ -83,12 +83,12 @@ #' prior = LogNormal( #' meanlog = 2, #' sdlog = 0.2 -#' ), -#' stan = stan_opts( -#' samples = 100, warmup = 200 -#' ), -#' verbose = interactive() -#' ) +#' ) +#' ), +#' stan = stan_opts( +#' samples = 100, warmup = 200 +#' ), +#' verbose = interactive() #' ) #' options(old_opts) #' } diff --git a/man/regional_epinow.Rd b/man/regional_epinow.Rd index ecf613441..ad7cd3bbf 100644 --- a/man/regional_epinow.Rd +++ b/man/regional_epinow.Rd @@ -160,12 +160,12 @@ def <- regional_epinow( prior = LogNormal( meanlog = 2, sdlog = 0.2 - ), - stan = stan_opts( - samples = 100, warmup = 200 - ), - verbose = interactive() - ) + ) + ), + stan = stan_opts( + samples = 100, warmup = 200 + ), + verbose = interactive() ) options(old_opts) } From e049f74a81dd94744bdef5ba3e6b20a43a8fd0b8 Mon Sep 17 00:00:00 2001 From: James Azam Date: Fri, 13 Dec 2024 11:22:08 +0000 Subject: [PATCH 07/18] Change uncertain parameters to numeric --- vignettes/estimate_infections_options.Rmd.orig | 4 ++-- vignettes/estimate_infections_workflow.Rmd.orig | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/vignettes/estimate_infections_options.Rmd.orig b/vignettes/estimate_infections_options.Rmd.orig index 068ffad34..f49776ada 100644 --- a/vignettes/estimate_infections_options.Rmd.orig +++ b/vignettes/estimate_infections_options.Rmd.orig @@ -179,8 +179,8 @@ Here, instead of doing so we assume that we know about truncation with mean of 1 ```{r define_truncation} trunc_dist <- LogNormal( - meanlog = Normal(0.5, 0.1), - sdlog = Normal(0.5, 0.1), + meanlog = 0.5, + sdlog = 0.1, max = 3 ) trunc_dist diff --git a/vignettes/estimate_infections_workflow.Rmd.orig b/vignettes/estimate_infections_workflow.Rmd.orig index e4e294bc2..d3762f53e 100644 --- a/vignettes/estimate_infections_workflow.Rmd.orig +++ b/vignettes/estimate_infections_workflow.Rmd.orig @@ -93,7 +93,7 @@ If distributions are variable, the values with uncertainty are treated as [prior For example, to define a variable gamma distribution where uncertainty in the mean is given by a normal distribution with mean 3 and sd 2, and uncertainty in the standard deviation is given by a normal distribution with mean 1 and sd 0.1, with a maximum value 10, you would write ```{r} -uncertain_gamma <- Gamma(shape = Normal(3, 2), rate = Normal(1, 0.1), max = 10) +uncertain_gamma <- Gamma(shape = 3, rate = 1, max = 10) uncertain_gamma ``` From bd759983930f705c38371b59737398ce814d5911 Mon Sep 17 00:00:00 2001 From: James Azam Date: Fri, 13 Dec 2024 11:35:51 +0000 Subject: [PATCH 08/18] Revert to using mean and sd --- R/epinow.R | 4 ++-- R/estimate_infections.R | 6 +----- R/regional_epinow.R | 7 +------ R/simulate_infections.R | 7 +------ vignettes/EpiNow2.Rmd.orig | 19 +++---------------- vignettes/epinow.Rmd.orig | 11 ++--------- .../estimate_infections_options.Rmd.orig | 11 ++--------- .../estimate_infections_workflow.Rmd.orig | 4 ++-- 8 files changed, 14 insertions(+), 55 deletions(-) diff --git a/R/epinow.R b/R/epinow.R index e90bb0fa8..ef6b03cea 100644 --- a/R/epinow.R +++ b/R/epinow.R @@ -57,8 +57,8 @@ #' # set an example reporting delay. In practice this should use an estimate #' # from the literature or be estimated from data #' reporting_delay <- LogNormal( -#' meanlog = 2, -#' sdlog = 1, +#' mean = 2, +#' sd = 1, #' max = 10 #' ) #' diff --git a/R/estimate_infections.R b/R/estimate_infections.R index d59e0b988..1d5b5f004 100644 --- a/R/estimate_infections.R +++ b/R/estimate_infections.R @@ -101,11 +101,7 @@ #' ) #' # set an example reporting delay. In practice this should use an estimate #' # from the literature or be estimated from data -#' reporting_delay <- LogNormal( -#' meanlog = 2, -#' sdlog = 1, -#' max = 10 -#' ) +#' reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10) #' #' # for more examples, see the "estimate_infections examples" vignette #' def <- estimate_infections(reported_cases, diff --git a/R/regional_epinow.R b/R/regional_epinow.R index 5df23def5..a5f08d062 100644 --- a/R/regional_epinow.R +++ b/R/regional_epinow.R @@ -79,12 +79,7 @@ #' data = cases, #' generation_time = gt_opts(example_generation_time), #' delays = delay_opts(example_incubation_period + example_reporting_delay), -#' rt = rt_opts( -#' prior = LogNormal( -#' meanlog = 2, -#' sdlog = 0.2 -#' ) -#' ), +#' rt = rt_opts(prior = LogNormal(mean = 2, sd = 0.2)), #' stan = stan_opts( #' samples = 100, warmup = 200 #' ), diff --git a/R/simulate_infections.R b/R/simulate_infections.R index ac43b1400..f3b66cc3f 100644 --- a/R/simulate_infections.R +++ b/R/simulate_infections.R @@ -277,12 +277,7 @@ simulate_infections <- function(estimates, R, initial_infections, #' est <- estimate_infections(reported_cases, #' generation_time = generation_time_opts(example_generation_time), #' delays = delay_opts(example_incubation_period + example_reporting_delay), -#' rt = rt_opts(prior = LogNormal( -#' meanlog = 2, -#' sdlog = 0.1 -#' ), -#' rw = 7 -#' ), +#' rt = rt_opts(prior = LogNormal(mean = 2, sd = 0.1), rw = 7), #' obs = obs_opts(scale = Normal(mean = 0.1, sd = 0.01)), #' gp = NULL, horizon = 0 #' ) diff --git a/vignettes/EpiNow2.Rmd.orig b/vignettes/EpiNow2.Rmd.orig index 23cfe67a8..69d10bcdd 100644 --- a/vignettes/EpiNow2.Rmd.orig +++ b/vignettes/EpiNow2.Rmd.orig @@ -50,11 +50,7 @@ If data was not available we could instead specify an informed estimate of the l 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*. ```{r} -reporting_delay <- LogNormal( - meanlog = 2, - sdlog = 1, - max = 10 -) +reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10) reporting_delay ``` @@ -98,10 +94,7 @@ estimates <- epinow( data = reported_cases, generation_time = gt_opts(example_generation_time), delays = delay_opts(example_incubation_period + reporting_delay), - rt = rt_opts(prior = LogNormal( - meanlog = 2, - sdlog = 0.2 - )), + rt = rt_opts(prior = LogNormal(mean = 2, sd = 0.2)), stan = stan_opts(cores = 4, control = list(adapt_delta = 0.99)), verbose = interactive() ) @@ -155,13 +148,7 @@ estimates <- regional_epinow( data = reported_cases, generation_time = gt_opts(example_generation_time), delays = delay_opts(example_incubation_period + reporting_delay), - rt = rt_opts( - prior = LogNormal( - meanlog = 2, - sdlog = 0.2 - ), - rw = 7 - ), + rt = rt_opts(prior = LogNormal(mean = 2, sd = 0.2), rw = 7), gp = NULL, stan = stan_opts(cores = 4, warmup = 250, samples = 1000) ) diff --git a/vignettes/epinow.Rmd.orig b/vignettes/epinow.Rmd.orig index 10d79f032..d8960bfbd 100644 --- a/vignettes/epinow.Rmd.orig +++ b/vignettes/epinow.Rmd.orig @@ -38,16 +38,9 @@ This should be replaced with parameters relevant to the system that is being stu library("EpiNow2") options(mc.cores = 4) reported_cases <- example_confirmed[1:60] -reporting_delay <- LogNormal( - meanlog = 2, - sdlog = 1, - max = 10 -) +reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10) delay <- example_incubation_period + reporting_delay -rt_prior <- LogNormal( - meanlog = 2, - sdlog = 0.1 -) +rt_prior <- LogNormal(mean = 2, sd = 0.1) ``` We can then run the `epinow()` function with the same arguments as `estimate_infections()`. diff --git a/vignettes/estimate_infections_options.Rmd.orig b/vignettes/estimate_infections_options.Rmd.orig index f49776ada..23eca5459 100644 --- a/vignettes/estimate_infections_options.Rmd.orig +++ b/vignettes/estimate_infections_options.Rmd.orig @@ -73,11 +73,7 @@ For the reporting delay, we use a lognormal distribution with mean of 2 days and Note that the mean and standard deviation must be converted to the log scale, which can be done using the `convert_log_logmean()` function. ```{r reporting_delay} -reporting_delay <- LogNormal( - meanlog = 2, - sdlog = 1, - max = 10 -) +reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10) reporting_delay ``` @@ -101,10 +97,7 @@ example_generation_time Lastly we need to choose a prior for the initial value of the reproduction number. This is assumed by the model to be normally distributed and we can set the mean and the standard deviation. We decide to set the mean to 2 and the standard deviation to 1. ```{r initial_r} -rt_prior <- LogNormal( - meanlog = 2, - sdlog = 0.1 -) +rt_prior <- LogNormal(mean = 2, sd = 0.1) ``` # Running the model diff --git a/vignettes/estimate_infections_workflow.Rmd.orig b/vignettes/estimate_infections_workflow.Rmd.orig index d3762f53e..02b71bc33 100644 --- a/vignettes/estimate_infections_workflow.Rmd.orig +++ b/vignettes/estimate_infections_workflow.Rmd.orig @@ -79,7 +79,7 @@ In all cases, the distributions given can be *fixed* (i.e. have no uncertainty) For example, to define a fixed gamma distribution with mean 3, standard deviation (sd) 1 and maximum value 10, you would write ```{r} -fixed_gamma <- Gamma(shape = 3, rate = 1, max = 10) +fixed_gamma <- Gamma(mean = 3, sd = 1, max = 10) fixed_gamma ``` @@ -205,7 +205,7 @@ It can be changed using the `rt_opts()` function. For example, if the user believes that at the very start of the data the reproduction number was 2, with uncertainty in this belief represented by a standard deviation of 1, they would use ```{r results = 'hide'} -rt_prior <- LogNormal(meanlog = 2, sdlog = 1) +rt_prior <- LogNormal(mean = 2, sd = 1) rt_opts(prior = rt_prior) ``` From 4960f7f81f8d58b6788e4f2f4bcb8c7b1b7b8120 Mon Sep 17 00:00:00 2001 From: James Azam Date: Fri, 13 Dec 2024 11:36:44 +0000 Subject: [PATCH 09/18] Use numeric instead of uncertain parameters --- vignettes/estimate_infections_options.Rmd.orig | 4 ++-- vignettes/estimate_infections_workflow.Rmd.orig | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/vignettes/estimate_infections_options.Rmd.orig b/vignettes/estimate_infections_options.Rmd.orig index 23eca5459..d23566fbb 100644 --- a/vignettes/estimate_infections_options.Rmd.orig +++ b/vignettes/estimate_infections_options.Rmd.orig @@ -172,8 +172,8 @@ Here, instead of doing so we assume that we know about truncation with mean of 1 ```{r define_truncation} trunc_dist <- LogNormal( - meanlog = 0.5, - sdlog = 0.1, + mean = 0.5, + sd = 0.1, max = 3 ) trunc_dist diff --git a/vignettes/estimate_infections_workflow.Rmd.orig b/vignettes/estimate_infections_workflow.Rmd.orig index 02b71bc33..3b67c9722 100644 --- a/vignettes/estimate_infections_workflow.Rmd.orig +++ b/vignettes/estimate_infections_workflow.Rmd.orig @@ -93,7 +93,7 @@ If distributions are variable, the values with uncertainty are treated as [prior For example, to define a variable gamma distribution where uncertainty in the mean is given by a normal distribution with mean 3 and sd 2, and uncertainty in the standard deviation is given by a normal distribution with mean 1 and sd 0.1, with a maximum value 10, you would write ```{r} -uncertain_gamma <- Gamma(shape = 3, rate = 1, max = 10) +uncertain_gamma <- Gamma(mean = 3, sd = 0.1, max = 10) uncertain_gamma ``` From a34e03d4787ed14be8dec8827cc82393e9f534b4 Mon Sep 17 00:00:00 2001 From: James Azam Date: Fri, 13 Dec 2024 11:37:47 +0000 Subject: [PATCH 10/18] Revert indentation --- R/epinow.R | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/R/epinow.R b/R/epinow.R index ef6b03cea..6ca1fd078 100644 --- a/R/epinow.R +++ b/R/epinow.R @@ -56,11 +56,7 @@ #' ) #' # set an example reporting delay. In practice this should use an estimate #' # from the literature or be estimated from data -#' reporting_delay <- LogNormal( -#' mean = 2, -#' sd = 1, -#' max = 10 -#' ) +#' reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10) #' #' # example case data #' reported_cases <- example_confirmed[1:40] From d5c76120456035c7dc0ded0384b59004722ad9c1 Mon Sep 17 00:00:00 2001 From: jamesaazam Date: Fri, 13 Dec 2024 12:09:56 +0000 Subject: [PATCH 11/18] Use natural parameters to specify uncertain distributions --- R/dist_spec.R | 88 +++++++++++++++++++++++++++---------- R/opts.R | 10 +++-- man/Distributions.Rd | 3 +- man/c.dist_spec.Rd | 7 ++- man/collapse.Rd | 9 +++- man/delay_opts.Rd | 6 ++- man/discretise.Rd | 9 +++- man/extract_single_dist.Rd | 7 ++- man/fix_parameters.Rd | 9 ++-- man/generation_time_opts.Rd | 4 +- man/is_constrained.Rd | 9 +++- man/max.dist_spec.Rd | 9 +++- man/mean.dist_spec.Rd | 7 ++- man/plot.dist_spec.Rd | 7 ++- man/plus-.dist_spec.Rd | 7 ++- man/print.dist_spec.Rd | 5 ++- 16 files changed, 142 insertions(+), 54 deletions(-) diff --git a/R/dist_spec.R b/R/dist_spec.R index 5b9a155ba..a28cb2141 100644 --- a/R/dist_spec.R +++ b/R/dist_spec.R @@ -116,9 +116,12 @@ discrete_pmf <- function(distribution = #' ) #' dist1 + dist1 #' -#' # An uncertain gamma distribution with mean 3 and sd 2 +#' # An uncertain gamma distribution with shape and rate normally distributed +#' # as Normal(3, 0.5) and Normal(2, 0.5) respectively #' dist2 <- Gamma( -#' mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20 +#' shape = Normal(3, 0.5), +#' rate = Normal(2, 0.5), +#' max = 20 #' ) #' dist1 + dist2 `+.dist_spec` <- function(e1, e2) { @@ -197,9 +200,12 @@ discrete_pmf <- function(distribution = #' ) #' dist1 + dist1 #' -#' # An uncertain gamma distribution with mean 3 and sd 2 +#' # An uncertain gamma distribution with shape and rate normally distributed +#' # as Normal(3, 0.5) and Normal(2, 0.5) respectively #' dist2 <- Gamma( -#' mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20 +#' shape = Normal(3, 0.5), +#' rate = Normal(2, 0.5), +#' max = 20 #' ) #' c(dist1, dist2) c.dist_spec <- function(...) { @@ -257,9 +263,12 @@ c.dist_spec <- function(...) { #' dist1 <- LogNormal(mean = 5, sd = 1, max = 20) #' mean(dist1) #' -#' # An uncertain gamma distribution with mean 3 and sd 2 +#' # An uncertain gamma distribution with shape and rate normally distributed +#' # as Normal(3, 0.5) and Normal(2, 0.5) respectively #' dist2 <- Gamma( -#' mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20 +#' shape = Normal(3, 0.5), +#' rate = Normal(2, 0.5), +#' max = 20 #' ) #' mean(dist2) #' @@ -393,8 +402,13 @@ sd.default <- function(x, ...) { #' dist1 <- Gamma(mean = 5, sd = 1, max = 20) #' max(dist1) #' -#' # An uncertain lognormal distribution with mean 3 and sd 2 -#' dist2 <- LogNormal(mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20) +#' # An uncertain lognormal distribution with meanlog and sdlog normally +#' # distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively +#' dist2 <- LogNormal( +#' meanlog = Normal(3, 0.5), +#' sdlog = Normal(2, 0.5), +#' max = 20 +#' ) #' max(dist2) #' #' # The max the sum of two distributions @@ -447,8 +461,13 @@ discretise <- function(x, ...) { #' # A fixed gamma distribution with mean 5 and sd 1. #' dist1 <- Gamma(mean = 5, sd = 1, max = 20) #' -#' # An uncertain lognormal distribution with mean 3 and sd 2 -#' dist2 <- LogNormal(mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20) +#' # An uncertain lognormal distribution with meanlog and sdlog normally +#' # distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively +#' dist2 <- LogNormal( +#' meanlog = Normal(3, 0.5), +#' sdlog = Normal(2, 0.5), +#' max = 20 +#' ) #' #' # The maxf the sum of two distributions #' discretise(dist1 + dist2, strict = FALSE) @@ -531,8 +550,13 @@ collapse <- function(x, ...) { #' # A fixed gamma distribution with mean 5 and sd 1. #' dist1 <- Gamma(mean = 5, sd = 1, max = 20) #' -#' # An uncertain lognormal distribution with mean 3 and sd 2 -#' dist2 <- LogNormal(mean = 3, sd = 2, max = 20) +#' # An uncertain lognormal distribution with meanlog and sdlog normally +#' # distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively +#' dist2 <- LogNormal( +#' meanlog = Normal(3, 0.5), +#' sdlog = Normal(2, 0.5), +#' max = 20 +#' ) #' #' # The maxf the sum of two distributions #' collapse(discretise(dist1 + dist2)) @@ -588,9 +612,10 @@ collapse.multi_dist_spec <- function(x, ...) { #' dist1 <- LogNormal(mean = 1.5, sd = 0.5, max = 20) #' print(dist1) #' -#' # An uncertain gamma distribution with mean 3 and sd 2 +#' # An uncertain gamma distribution with shape and rate normally distributed +#' # as Normal(3, 0.5) and Normal(2, 0.5) respectively #' dist2 <- Gamma( -#' mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20 +#' shape = Normal(3, 0.5), rate = Normal(2, 0.5), max = 20 #' ) #' print(dist2) print.dist_spec <- function(x, ...) { @@ -680,9 +705,12 @@ print.dist_spec <- function(x, ...) { #' # Plot discretised distribution with 0.01 day discretisation window #' plot(dist1, res = 0.01, cumulative = FALSE) #' -#' # An uncertain gamma distribution with mean 3 and sd 2 +#' # An uncertain gamma distribution with shape and rate normally distributed +#' # as Normal(3, 0.5) and Normal(2, 0.5) respectively #' dist2 <- Gamma( -#' mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20 +#' shape = Normal(3, 0.5), +#' rate = Normal(2, 0.5), +#' max = 20 #' ) #' plot(dist2) #' @@ -781,9 +809,12 @@ plot.dist_spec <- function(x, samples = 50L, res = 1, cumulative = TRUE, ...) { #' @examples #' dist1 <- LogNormal(mean = 1.6, sd = 0.5, max = 20) #' -#' # An uncertain gamma distribution with mean 3 and sd 2 +#' # An uncertain gamma distribution with shape and rate normally distributed +#' # as Normal(3, 0.5) and Normal(2, 0.5) respectively #' dist2 <- Gamma( -#' mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20 +#' shape = Normal(3, 0.5), +#' rate = Normal(2, 0.5), +#' max = 20 #' ) #' #' # Multiple distributions @@ -828,9 +859,12 @@ fix_parameters <- function(x, ...) { #' @importFrom rlang arg_match #' @method fix_parameters dist_spec #' @examples -#' # An uncertain gamma distribution with mean 3 and sd 2 -#' dist <- LogNormal( -#' meanlog = Normal(3, 0.5), sdlog = Normal(2, 0.5), max = 20 +#' # An uncertain gamma distribution with shape and rate normally distributed +#' # as Normal(3, 0.5) and Normal(2, 0.5) respectively +#' dist <- Gamma( +#' shape = Normal(3, 0.5), +#' rate = Normal(2, 0.5), +#' max = 20 #' ) #' #' fix_parameters(dist) @@ -891,8 +925,13 @@ is_constrained <- function(x, ...) { #' # A fixed gamma distribution with mean 5 and sd 1. #' dist1 <- Gamma(mean = 5, sd = 1, max = 20) #' -#' # An uncertain lognormal distribution with mean 3 and sd 2 -#' dist2 <- LogNormal(mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20) +#' # An uncertain lognormal distribution with meanlog and sdlog normally +#' # distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively +#' dist2 <- LogNormal( +#' meanlog = Normal(3, 0.5), +#' sdlog = Normal(2, 0.5), +#' max = 20 +#' ) #' #' # both distributions are constrained and therefore so is the sum #' is_constrained(dist1 + dist2) @@ -952,7 +991,8 @@ is_constrained.multi_dist_spec <- function(x, ...) { #' @examples #' LogNormal(mean = 4, sd = 1) #' LogNormal(mean = 4, sd = 1, max = 10) -#' LogNormal(mean = Normal(4, 1), sd = 1, max = 10) +#' # If specifying uncertain parameters, use the natural parameters +#' LogNormal(meanlog = Normal(4, 1), sdlog = 1, max = 10) LogNormal <- function(meanlog, sdlog, mean, sd, ...) { params <- as.list(environment()) return(new_dist_spec(params, "lognormal", ...)) diff --git a/R/opts.R b/R/opts.R index 04d9e5577..613a6452f 100644 --- a/R/opts.R +++ b/R/opts.R @@ -42,8 +42,8 @@ #' # An uncertain gamma distributed generation time #' generation_time_opts( #' Gamma( -#' mean = Normal(mean = 3, sd = 1), -#' sd = Normal(mean = 2, sd = 0.5), +#' shape = Normal(mean = 3, sd = 1), +#' rate = Normal(mean = 2, sd = 0.5), #' max = 14 #' ) #' ) @@ -186,7 +186,11 @@ secondary_opts <- function(type = c("incidence", "prevalence"), ...) { #' delay_opts() #' #' # A single delay that has uncertainty -#' delay <- LogNormal(mean = Normal(1, 0.2), sd = Normal(0.5, 0.1), max = 14) +#' delay <- LogNormal( +#' meanlog = Normal(1, 0.2), +#' sdlog = Normal(0.5, 0.1), +#' max = 14 +#' ) #' delay_opts(delay) #' #' # A single delay without uncertainty diff --git a/man/Distributions.Rd b/man/Distributions.Rd index 583b0193d..6db8e8911 100644 --- a/man/Distributions.Rd +++ b/man/Distributions.Rd @@ -77,7 +77,8 @@ probability mass function is given directly as a numeric vector. \examples{ LogNormal(mean = 4, sd = 1) LogNormal(mean = 4, sd = 1, max = 10) -LogNormal(mean = Normal(4, 1), sd = 1, max = 10) +# If specifying uncertain parameters, use the natural parameters +LogNormal(meanlog = Normal(4, 1), sdlog = 1, max = 10) Gamma(mean = 4, sd = 1) Gamma(shape = 16, rate = 4) Gamma(shape = Normal(16, 2), rate = Normal(4, 1)) diff --git a/man/c.dist_spec.Rd b/man/c.dist_spec.Rd index 45cffbf55..ae603e9e2 100644 --- a/man/c.dist_spec.Rd +++ b/man/c.dist_spec.Rd @@ -27,9 +27,12 @@ dist1 <- LogNormal( ) dist1 + dist1 -# An uncertain gamma distribution with mean 3 and sd 2 +# An uncertain gamma distribution with shape and rate normally distributed +# as Normal(3, 0.5) and Normal(2, 0.5) respectively dist2 <- Gamma( - mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20 + shape = Normal(3, 0.5), + rate = Normal(2, 0.5), + max = 20 ) c(dist1, dist2) } diff --git a/man/collapse.Rd b/man/collapse.Rd index 97ed0d1c2..e3d67be00 100644 --- a/man/collapse.Rd +++ b/man/collapse.Rd @@ -25,8 +25,13 @@ in the . # A fixed gamma distribution with mean 5 and sd 1. dist1 <- Gamma(mean = 5, sd = 1, max = 20) -# An uncertain lognormal distribution with mean 3 and sd 2 -dist2 <- LogNormal(mean = 3, sd = 2, max = 20) +# An uncertain lognormal distribution with meanlog and sdlog normally +# distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively +dist2 <- LogNormal( + meanlog = Normal(3, 0.5), + sdlog = Normal(2, 0.5), + max = 20 +) # The maxf the sum of two distributions collapse(discretise(dist1 + dist2)) diff --git a/man/delay_opts.Rd b/man/delay_opts.Rd index b28a47b72..dc383a541 100644 --- a/man/delay_opts.Rd +++ b/man/delay_opts.Rd @@ -44,7 +44,11 @@ functions. delay_opts() # A single delay that has uncertainty -delay <- LogNormal(mean = Normal(1, 0.2), sd = Normal(0.5, 0.1), max = 14) +delay <- LogNormal( + meanlog = Normal(1, 0.2), + sdlog = Normal(0.5, 0.1), + max = 14 +) delay_opts(delay) # A single delay without uncertainty diff --git a/man/discretise.Rd b/man/discretise.Rd index e6c2cda1a..a76646d4a 100644 --- a/man/discretise.Rd +++ b/man/discretise.Rd @@ -42,8 +42,13 @@ from the difference of two censored events. # A fixed gamma distribution with mean 5 and sd 1. dist1 <- Gamma(mean = 5, sd = 1, max = 20) -# An uncertain lognormal distribution with mean 3 and sd 2 -dist2 <- LogNormal(mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20) +# An uncertain lognormal distribution with meanlog and sdlog normally +# distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively +dist2 <- LogNormal( + meanlog = Normal(3, 0.5), + sdlog = Normal(2, 0.5), + max = 20 +) # The maxf the sum of two distributions discretise(dist1 + dist2, strict = FALSE) diff --git a/man/extract_single_dist.Rd b/man/extract_single_dist.Rd index e3407cf6b..e6761e96c 100644 --- a/man/extract_single_dist.Rd +++ b/man/extract_single_dist.Rd @@ -20,9 +20,12 @@ A single \code{dist_spec} object \examples{ dist1 <- LogNormal(mean = 1.6, sd = 0.5, max = 20) -# An uncertain gamma distribution with mean 3 and sd 2 +# An uncertain gamma distribution with shape and rate normally distributed +# as Normal(3, 0.5) and Normal(2, 0.5) respectively dist2 <- Gamma( - mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20 + shape = Normal(3, 0.5), + rate = Normal(2, 0.5), + max = 20 ) # Multiple distributions diff --git a/man/fix_parameters.Rd b/man/fix_parameters.Rd index 04a8a4136..a43df5bbd 100644 --- a/man/fix_parameters.Rd +++ b/man/fix_parameters.Rd @@ -25,9 +25,12 @@ If the given \verb{} has any uncertainty, it is removed and the corresponding distribution converted into a fixed one. } \examples{ -# An uncertain gamma distribution with mean 3 and sd 2 -dist <- LogNormal( - meanlog = Normal(3, 0.5), sdlog = Normal(2, 0.5), max = 20 +# An uncertain gamma distribution with shape and rate normally distributed +# as Normal(3, 0.5) and Normal(2, 0.5) respectively +dist <- Gamma( + shape = Normal(3, 0.5), + rate = Normal(2, 0.5), + max = 20 ) fix_parameters(dist) diff --git a/man/generation_time_opts.Rd b/man/generation_time_opts.Rd index 1200101c7..8306d7a61 100644 --- a/man/generation_time_opts.Rd +++ b/man/generation_time_opts.Rd @@ -79,8 +79,8 @@ generation_time_opts(Gamma(mean = 3, sd = 2, max = 14)) # An uncertain gamma distributed generation time generation_time_opts( Gamma( - mean = Normal(mean = 3, sd = 1), - sd = Normal(mean = 2, sd = 0.5), + shape = Normal(mean = 3, sd = 1), + rate = Normal(mean = 2, sd = 0.5), max = 14 ) ) diff --git a/man/is_constrained.Rd b/man/is_constrained.Rd index 7a097ab6d..39803ee62 100644 --- a/man/is_constrained.Rd +++ b/man/is_constrained.Rd @@ -23,8 +23,13 @@ Logical; TRUE if \code{x} is constrained # A fixed gamma distribution with mean 5 and sd 1. dist1 <- Gamma(mean = 5, sd = 1, max = 20) -# An uncertain lognormal distribution with mean 3 and sd 2 -dist2 <- LogNormal(mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20) +# An uncertain lognormal distribution with meanlog and sdlog normally +# distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively +dist2 <- LogNormal( + meanlog = Normal(3, 0.5), + sdlog = Normal(2, 0.5), + max = 20 +) # both distributions are constrained and therefore so is the sum is_constrained(dist1 + dist2) diff --git a/man/max.dist_spec.Rd b/man/max.dist_spec.Rd index c81f811ef..516a0cc3a 100644 --- a/man/max.dist_spec.Rd +++ b/man/max.dist_spec.Rd @@ -25,8 +25,13 @@ in parameters) dist1 <- Gamma(mean = 5, sd = 1, max = 20) max(dist1) -# An uncertain lognormal distribution with mean 3 and sd 2 -dist2 <- LogNormal(mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20) +# An uncertain lognormal distribution with meanlog and sdlog normally +# distributed as Normal(3, 0.5) and Normal(2, 0.5) respectively +dist2 <- LogNormal( + meanlog = Normal(3, 0.5), + sdlog = Normal(2, 0.5), + max = 20 +) max(dist2) # The max the sum of two distributions diff --git a/man/mean.dist_spec.Rd b/man/mean.dist_spec.Rd index 5448ca1b8..8be98f2d0 100644 --- a/man/mean.dist_spec.Rd +++ b/man/mean.dist_spec.Rd @@ -25,9 +25,12 @@ distributions combined in the passed . dist1 <- LogNormal(mean = 5, sd = 1, max = 20) mean(dist1) -# An uncertain gamma distribution with mean 3 and sd 2 +# An uncertain gamma distribution with shape and rate normally distributed +# as Normal(3, 0.5) and Normal(2, 0.5) respectively dist2 <- Gamma( - mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20 + shape = Normal(3, 0.5), + rate = Normal(2, 0.5), + max = 20 ) mean(dist2) diff --git a/man/plot.dist_spec.Rd b/man/plot.dist_spec.Rd index 27698649a..7e1a8d19e 100644 --- a/man/plot.dist_spec.Rd +++ b/man/plot.dist_spec.Rd @@ -33,9 +33,12 @@ plot(dist1) # Plot discretised distribution with 0.01 day discretisation window plot(dist1, res = 0.01, cumulative = FALSE) -# An uncertain gamma distribution with mean 3 and sd 2 +# An uncertain gamma distribution with shape and rate normally distributed +# as Normal(3, 0.5) and Normal(2, 0.5) respectively dist2 <- Gamma( - mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20 + shape = Normal(3, 0.5), + rate = Normal(2, 0.5), + max = 20 ) plot(dist2) diff --git a/man/plus-.dist_spec.Rd b/man/plus-.dist_spec.Rd index 16e967faa..091b7e536 100644 --- a/man/plus-.dist_spec.Rd +++ b/man/plus-.dist_spec.Rd @@ -26,9 +26,12 @@ dist1 <- LogNormal( ) dist1 + dist1 -# An uncertain gamma distribution with mean 3 and sd 2 +# An uncertain gamma distribution with shape and rate normally distributed +# as Normal(3, 0.5) and Normal(2, 0.5) respectively dist2 <- Gamma( - mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20 + shape = Normal(3, 0.5), + rate = Normal(2, 0.5), + max = 20 ) dist1 + dist2 } diff --git a/man/print.dist_spec.Rd b/man/print.dist_spec.Rd index c3ea43cfb..b6925af9f 100644 --- a/man/print.dist_spec.Rd +++ b/man/print.dist_spec.Rd @@ -24,9 +24,10 @@ functions of fixed delay distributions combined in the passed . dist1 <- LogNormal(mean = 1.5, sd = 0.5, max = 20) print(dist1) -# An uncertain gamma distribution with mean 3 and sd 2 +# An uncertain gamma distribution with shape and rate normally distributed +# as Normal(3, 0.5) and Normal(2, 0.5) respectively dist2 <- Gamma( - mean = Normal(3, 0.5), sd = Normal(2, 0.5), max = 20 + shape = Normal(3, 0.5), rate = Normal(2, 0.5), max = 20 ) print(dist2) } From 78f0800be185b468a3752981cc9ad222820c9705 Mon Sep 17 00:00:00 2001 From: jamesaazam Date: Fri, 13 Dec 2024 12:10:50 +0000 Subject: [PATCH 12/18] Redoc after revert to mean and sd for numeric parameters --- man/epinow.Rd | 6 +----- man/estimate_infections.Rd | 6 +----- man/forecast_infections.Rd | 7 +------ man/regional_epinow.Rd | 7 +------ 4 files changed, 4 insertions(+), 22 deletions(-) diff --git a/man/epinow.Rd b/man/epinow.Rd index 40b295a02..bfe6dc102 100644 --- a/man/epinow.Rd +++ b/man/epinow.Rd @@ -157,11 +157,7 @@ incubation_period <- LogNormal( ) # set an example reporting delay. In practice this should use an estimate # from the literature or be estimated from data -reporting_delay <- LogNormal( - meanlog = 2, - sdlog = 1, - max = 10 -) +reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10) # example case data reported_cases <- example_confirmed[1:40] diff --git a/man/estimate_infections.Rd b/man/estimate_infections.Rd index 0c1686397..0372966a0 100644 --- a/man/estimate_infections.Rd +++ b/man/estimate_infections.Rd @@ -146,11 +146,7 @@ incubation_period <- LogNormal( ) # set an example reporting delay. In practice this should use an estimate # from the literature or be estimated from data -reporting_delay <- LogNormal( - meanlog = 2, - sdlog = 1, - max = 10 -) +reporting_delay <- LogNormal(mean = 2, sd = 1, max = 10) # for more examples, see the "estimate_infections examples" vignette def <- estimate_infections(reported_cases, diff --git a/man/forecast_infections.Rd b/man/forecast_infections.Rd index b836dc323..24c5b2c5d 100644 --- a/man/forecast_infections.Rd +++ b/man/forecast_infections.Rd @@ -65,12 +65,7 @@ reported_cases <- example_confirmed[1:50] est <- estimate_infections(reported_cases, generation_time = generation_time_opts(example_generation_time), delays = delay_opts(example_incubation_period + example_reporting_delay), - rt = rt_opts(prior = LogNormal( - meanlog = 2, - sdlog = 0.1 - ), - rw = 7 - ), + rt = rt_opts(prior = LogNormal(mean = 2, sd = 0.1), rw = 7), obs = obs_opts(scale = Normal(mean = 0.1, sd = 0.01)), gp = NULL, horizon = 0 ) diff --git a/man/regional_epinow.Rd b/man/regional_epinow.Rd index ad7cd3bbf..6eca0c83a 100644 --- a/man/regional_epinow.Rd +++ b/man/regional_epinow.Rd @@ -156,12 +156,7 @@ def <- regional_epinow( data = cases, generation_time = gt_opts(example_generation_time), delays = delay_opts(example_incubation_period + example_reporting_delay), - rt = rt_opts( - prior = LogNormal( - meanlog = 2, - sdlog = 0.2 - ) - ), + rt = rt_opts(prior = LogNormal(mean = 2, sd = 0.2)), stan = stan_opts( samples = 100, warmup = 200 ), From 3a534bc07c37d21c5841524629e52389ef3b77d9 Mon Sep 17 00:00:00 2001 From: jamesaazam Date: Fri, 13 Dec 2024 12:16:59 +0000 Subject: [PATCH 13/18] Update NEWS --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 8b3c46aeb..594403165 100644 --- a/NEWS.md +++ b/NEWS.md @@ -29,7 +29,7 @@ - Brought the docs on `alpha_sd` up to date with the code change from prior PR #853. By @zsusswein in #862 and reviewed by @jamesmbaazam. - The `...` argument in `estimate_secondary()` has been removed because it was not used. By @jamesmbaazam in #894 and reviewed by @. -- All examples now use the natural parameters of distributions rather than the mean and standard deviation to remove warnings and encourage best practice. By @jamesmbaazam in #893 and reviewed by @. +- All examples now use the natural parameters of distributions rather than the mean and standard deviation when specifying uncertain distributions. This is to eliminate warnings and encourage best practice. By @jamesmbaazam in #893 and reviewed by @sbfnk. # EpiNow2 1.6.1 From 6b0140f56a8abf60fff7d312762263c71e6203c7 Mon Sep 17 00:00:00 2001 From: jamesaazam Date: Fri, 13 Dec 2024 13:49:28 +0000 Subject: [PATCH 14/18] Fix example in collapse.dist_spec --- R/dist_spec.R | 2 +- man/collapse.Rd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/dist_spec.R b/R/dist_spec.R index a28cb2141..9c82a4c30 100644 --- a/R/dist_spec.R +++ b/R/dist_spec.R @@ -559,7 +559,7 @@ collapse <- function(x, ...) { #' ) #' #' # The maxf the sum of two distributions -#' collapse(discretise(dist1 + dist2)) +#' collapse(discretise(dist1 + dist2, strict = FALSE)) collapse.dist_spec <- function(x, ...) { return(x) } diff --git a/man/collapse.Rd b/man/collapse.Rd index e3d67be00..9715d9c33 100644 --- a/man/collapse.Rd +++ b/man/collapse.Rd @@ -34,5 +34,5 @@ dist2 <- LogNormal( ) # The maxf the sum of two distributions -collapse(discretise(dist1 + dist2)) +collapse(discretise(dist1 + dist2, strict = FALSE)) } From 07e94bce64e3555b9a3da10a820079719a3287ab Mon Sep 17 00:00:00 2001 From: jamesaazam Date: Fri, 13 Dec 2024 13:53:35 +0000 Subject: [PATCH 15/18] Align code with text --- vignettes/estimate_infections_options.Rmd.orig | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/estimate_infections_options.Rmd.orig b/vignettes/estimate_infections_options.Rmd.orig index d23566fbb..0cd532013 100644 --- a/vignettes/estimate_infections_options.Rmd.orig +++ b/vignettes/estimate_infections_options.Rmd.orig @@ -173,7 +173,7 @@ Here, instead of doing so we assume that we know about truncation with mean of 1 ```{r define_truncation} trunc_dist <- LogNormal( mean = 0.5, - sd = 0.1, + sd = 0.5, max = 3 ) trunc_dist From d9265c3fe6cafebf0a009ddda3970ce84fe23f27 Mon Sep 17 00:00:00 2001 From: jamesaazam Date: Fri, 13 Dec 2024 13:54:02 +0000 Subject: [PATCH 16/18] Use natural params of Gamma in vignette example --- vignettes/estimate_infections_workflow.Rmd.orig | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/vignettes/estimate_infections_workflow.Rmd.orig b/vignettes/estimate_infections_workflow.Rmd.orig index 3b67c9722..6f28d6cec 100644 --- a/vignettes/estimate_infections_workflow.Rmd.orig +++ b/vignettes/estimate_infections_workflow.Rmd.orig @@ -90,10 +90,10 @@ plot(fixed_gamma) ``` If distributions are variable, the values with uncertainty are treated as [prior probability densities](https://en.wikipedia.org/wiki/Prior_probability) in the Bayesian inference framework used by _EpiNow2_, i.e. they are estimated as part of the inference. -For example, to define a variable gamma distribution where uncertainty in the mean is given by a normal distribution with mean 3 and sd 2, and uncertainty in the standard deviation is given by a normal distribution with mean 1 and sd 0.1, with a maximum value 10, you would write +For example, to define a variable gamma distribution where uncertainty in the shape is given by a normal distribution with mean 3 and sd 2, and uncertainty in the rate is given by a normal distribution with mean 1 and sd 0.1, with a maximum value 10, you would write ```{r} -uncertain_gamma <- Gamma(mean = 3, sd = 0.1, max = 10) +uncertain_gamma <- Gamma(shape = Normal(3, 2), rate = Normal(1, 0.1), max = 10) uncertain_gamma ``` From 2f2b3bb633c2fea66746f5696daa5ae666524b78 Mon Sep 17 00:00:00 2001 From: James Azam Date: Tue, 17 Dec 2024 11:17:36 +0000 Subject: [PATCH 17/18] Update LogNormal example Co-authored-by: Sebastian Funk --- R/dist_spec.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/dist_spec.R b/R/dist_spec.R index 9c82a4c30..cce218198 100644 --- a/R/dist_spec.R +++ b/R/dist_spec.R @@ -992,7 +992,7 @@ is_constrained.multi_dist_spec <- function(x, ...) { #' LogNormal(mean = 4, sd = 1) #' LogNormal(mean = 4, sd = 1, max = 10) #' # If specifying uncertain parameters, use the natural parameters -#' LogNormal(meanlog = Normal(4, 1), sdlog = 1, max = 10) +#' LogNormal(meanlog = Normal(1.5, 0.5), sdlog = 0.25, max = 10) LogNormal <- function(meanlog, sdlog, mean, sd, ...) { params <- as.list(environment()) return(new_dist_spec(params, "lognormal", ...)) From 89e89a7345fcd9b4bb4eb0827058655fe016e117 Mon Sep 17 00:00:00 2001 From: jamesaazam Date: Tue, 17 Dec 2024 11:20:57 +0000 Subject: [PATCH 18/18] Generate docs --- man/Distributions.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/man/Distributions.Rd b/man/Distributions.Rd index 6db8e8911..add503e54 100644 --- a/man/Distributions.Rd +++ b/man/Distributions.Rd @@ -78,7 +78,7 @@ probability mass function is given directly as a numeric vector. LogNormal(mean = 4, sd = 1) LogNormal(mean = 4, sd = 1, max = 10) # If specifying uncertain parameters, use the natural parameters -LogNormal(meanlog = Normal(4, 1), sdlog = 1, max = 10) +LogNormal(meanlog = Normal(1.5, 0.5), sdlog = 0.25, max = 10) Gamma(mean = 4, sd = 1) Gamma(shape = 16, rate = 4) Gamma(shape = Normal(16, 2), rate = Normal(4, 1))