From f372ff51b7ca5f66bb34a55c2f45085ec5a2f3ce Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 19 Dec 2024 17:04:28 +0000 Subject: [PATCH 01/16] change over opts and create interface --- R/create.R | 2 +- R/opts.R | 22 ++++++++++++++++------ 2 files changed, 17 insertions(+), 7 deletions(-) diff --git a/R/create.R b/R/create.R index d65132f09..83045dfa6 100644 --- a/R/create.R +++ b/R/create.R @@ -263,7 +263,7 @@ create_rt_data <- function(rt = rt_opts(), breakpoints = NULL, breakpoints = breakpoints, future_fixed = as.numeric(future_rt$fixed), fixed_from = future_rt$from, - pop = rt$pop, + pop = as.integer(rt$pop != Fixed(0)), stationary = as.numeric(rt$gp_on == "R0"), future_time = horizon - future_rt$from ) diff --git a/R/opts.R b/R/opts.R index ca5c11dda..c3d87cff8 100644 --- a/R/opts.R +++ b/R/opts.R @@ -320,9 +320,9 @@ trunc_opts <- function(dist = Fixed(0), default_cdf_cutoff = 0.001, #' conservative estimate of break point changes (alter this by setting #' `gp = NULL`). #' -#' @param pop Integer, defaults to 0. Susceptible population initially present. -#' Used to adjust Rt estimates when otherwise fixed based on the proportion of -#' the population that is susceptible. When set to 0 no population adjustment +#' @param pop A `` giving the initial susceptible population size. +#' Used to adjust Rt estimates based on the proportion of the population that +#' is susceptible. Defaults to `Fixed(0)` which means no population adjustment #' is done. #' #' @param gp_on Character string, defaulting to "R_t-1". Indicates how the @@ -354,13 +354,12 @@ rt_opts <- function(prior = LogNormal(mean = 1, sd = 1), use_breakpoints = TRUE, future = "latest", gp_on = c("R_t-1", "R0"), - pop = 0) { + pop = Fixed(0)) { rt <- list( use_rt = use_rt, rw = rw, use_breakpoints = use_breakpoints, future = future, - pop = pop, gp_on = arg_match(gp_on) ) @@ -388,6 +387,17 @@ rt_opts <- function(prior = LogNormal(mean = 1, sd = 1), prior <- LogNormal(mean = prior$mean, sd = prior$sd) } + if (is.numeric(pop)) { + cli_warn( + c( + "!" = "Specifying {.var pop} as a numeric value is deprecated.", + "i" = "Use a {.cls dist_spec} instead, e.g. Fixed({pop})." + ) + ) + pop <- Fixed(pop) + } + rt$pop <- pop + if (rt$use_rt) { rt$prior <- prior } else { @@ -724,7 +734,7 @@ obs_opts <- function(family = c("negbin", "poisson"), cli_abort( c( "!" = "Specifying {.var phi} as a vector of length 2 is deprecated.", - "i" = "Mean and SD should be given as list elements." + "i" = "Use a {.cls dist_spec} instead, e.g. Normal(mean = {phi[1]}, sd = {phi[2]})." ) ) } From 3e30b646f385fb6df992906ecf9e54cfa2dbfeeb Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 19 Dec 2024 17:43:43 +0000 Subject: [PATCH 02/16] add support to stan --- R/create.R | 6 ++++-- inst/stan/data/estimate_infections_params.stan | 1 + inst/stan/data/rt.stan | 4 ++-- inst/stan/estimate_infections.stan | 6 +++++- inst/stan/functions/infections.stan | 8 ++++---- inst/stan/simulate_infections.stan | 8 +++++++- 6 files changed, 23 insertions(+), 10 deletions(-) diff --git a/R/create.R b/R/create.R index 83045dfa6..ef7c02237 100644 --- a/R/create.R +++ b/R/create.R @@ -263,7 +263,7 @@ create_rt_data <- function(rt = rt_opts(), breakpoints = NULL, breakpoints = breakpoints, future_fixed = as.numeric(future_rt$fixed), fixed_from = future_rt$from, - pop = as.integer(rt$pop != Fixed(0)), + use_pop = as.integer(rt$pop != Fixed(0)), stationary = as.numeric(rt$gp_on == "R0"), future_time = horizon - future_rt$from ) @@ -584,12 +584,14 @@ create_stan_data <- function(data, seeding_time, rt, gp, obs, backcalc, R0 = rt$prior, frac_obs = obs$scale, rep_phi = obs$phi, + pop = rt$pop, lower_bounds = c( alpha = 0, rho = 0, R0 = 0, frac_obs = 0, - rep_phi = 0 + rep_phi = 0, + pop = 0 ) ) ) diff --git a/inst/stan/data/estimate_infections_params.stan b/inst/stan/data/estimate_infections_params.stan index 85be5c1c9..a54ec9d32 100644 --- a/inst/stan/data/estimate_infections_params.stan +++ b/inst/stan/data/estimate_infections_params.stan @@ -3,3 +3,4 @@ int rho_id; // parameter id of rho (GP lengthscale) int R0_id; // parameter id of R0 int frac_obs_id; // parameter id of frac_obs int rep_phi_id; // parameter id of rep_phi_id +int pop_id; // parameter id of pop diff --git a/inst/stan/data/rt.stan b/inst/stan/data/rt.stan index 357330c19..709cb2ae9 100644 --- a/inst/stan/data/rt.stan +++ b/inst/stan/data/rt.stan @@ -5,5 +5,5 @@ array[t - seeding_time] int breakpoints; // when do breakpoints occur int future_fixed; // is underlying future Rt assumed to be fixed int fixed_from; // Reference date for when Rt estimation should be fixed - int pop; // Initial susceptible population - int gt_id; // id of generation time + int use_pop; // use population size + int gt_id; // id of generation time \ No newline at end of file diff --git a/inst/stan/estimate_infections.stan b/inst/stan/estimate_infections.stan index fe845b8fe..e079d345c 100644 --- a/inst/stan/estimate_infections.stan +++ b/inst/stan/estimate_infections.stan @@ -107,9 +107,13 @@ transformed parameters { frac_obs_id, params_fixed_lookup, params_variable_lookup, params_value, params ); + real pop = get_param( + pop_id, params_fixed_lookup, params_variable_lookup, params_value, + params + ); infections = generate_infections( R, seeding_time, gt_rev_pmf, initial_infections, initial_growth, pop, - future_time, obs_scale, frac_obs + use_pop, future_time, obs_scale, frac_obs ); } } else { diff --git a/inst/stan/functions/infections.stan b/inst/stan/functions/infections.stan index a599be8c2..f5e9f78c4 100644 --- a/inst/stan/functions/infections.stan +++ b/inst/stan/functions/infections.stan @@ -20,7 +20,7 @@ real update_infectiousness(vector infections, vector gt_rev_pmf, // generate infections by using Rt = Rt-1 * sum(reversed generation time pmf * infections) vector generate_infections(vector oR, int uot, vector gt_rev_pmf, array[] real initial_infections, array[] real initial_growth, - int pop, int ht, int obs_scale, real frac_obs) { + real pop, int use_pop, int ht, int obs_scale, real frac_obs) { // time indices and storage int ot = num_elements(oR); int nht = ot - ht; @@ -42,20 +42,20 @@ vector generate_infections(vector oR, int uot, vector gt_rev_pmf, } } // calculate cumulative infections - if (pop) { + if (use_pop) { cum_infections[1] = sum(infections[1:uot]); } // iteratively update infections for (s in 1:ot) { infectiousness[s] = update_infectiousness(infections, gt_rev_pmf, uot, s); - if (pop && s > nht) { + if (use_pop && s > nht) { exp_adj_Rt = exp(-R[s] * infectiousness[s] / (pop - cum_infections[nht])); exp_adj_Rt = exp_adj_Rt > 1 ? 1 : exp_adj_Rt; infections[s + uot] = (pop - cum_infections[s]) * (1 - exp_adj_Rt); }else{ infections[s + uot] = R[s] * infectiousness[s]; } - if (pop && s < ot) { + if (use_pop && s < ot) { cum_infections[s + 1] = cum_infections[s] + infections[s + uot]; } } diff --git a/inst/stan/simulate_infections.stan b/inst/stan/simulate_infections.stan index 43ea6cba9..f06006f6b 100644 --- a/inst/stan/simulate_infections.stan +++ b/inst/stan/simulate_infections.stan @@ -50,6 +50,12 @@ generated quantities { frac_obs_id, params_fixed_lookup, params_variable_lookup, params_value, params ); + + real pop = get_param( + pop_id, params_fixed_lookup, params_variable_lookup, + params_value, params + ); + for (i in 1:n) { // generate infections from Rt trace vector[delay_type_max[gt_id] + 1] gt_rev_pmf; @@ -62,7 +68,7 @@ generated quantities { infections[i] = to_row_vector(generate_infections( to_vector(R[i]), seeding_time, gt_rev_pmf, initial_infections[i], - initial_growth[i], pop, future_time, obs_scale, frac_obs[i] + initial_growth[i], pop, use_pop, future_time, obs_scale, frac_obs[i] )); if (delay_id) { From 288b424b90c34d6a893dca652af1b2cb1d509afa Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 19 Dec 2024 18:17:44 +0000 Subject: [PATCH 03/16] check works and update tests --- R/simulate_infections.R | 4 ++-- inst/stan/data/rt.stan | 2 +- inst/stan/data/simulation_rt.stan | 2 +- inst/stan/simulate_infections.stan | 4 ++-- man/EpiNow2-package.Rd | 2 +- man/rt_opts.Rd | 8 ++++---- man/simulate_infections.Rd | 8 ++++---- tests/testthat/test-create_rt_date.R | 4 ++-- tests/testthat/test-rt_opts.R | 14 ++++++++++---- vignettes/estimate_infections_options.Rmd.orig | 2 +- 10 files changed, 28 insertions(+), 22 deletions(-) diff --git a/R/simulate_infections.R b/R/simulate_infections.R index 6998a8171..25919d46a 100644 --- a/R/simulate_infections.R +++ b/R/simulate_infections.R @@ -69,7 +69,7 @@ simulate_infections <- function(estimates, R, initial_infections, CrIs = c(0.2, 0.5, 0.9), backend = "rstan", seeding_time = NULL, - pop = 0, ...) { + pop = Fixed(0), ...) { ## deprecated usage if (!missing(estimates)) { deprecate_stop( @@ -125,7 +125,7 @@ simulate_infections <- function(estimates, R, initial_infections, initial_infections = array(log_initial_infections, dim = c(1, 1)), initial_growth = array(initial_growth, dim = c(1, length(initial_growth))), R = array(R$R, dim = c(1, nrow(R))), - pop = pop + use_pop = as.integer(pop != Fixed(0)) ) data <- c(data, create_stan_delays( diff --git a/inst/stan/data/rt.stan b/inst/stan/data/rt.stan index 709cb2ae9..7ffff8a25 100644 --- a/inst/stan/data/rt.stan +++ b/inst/stan/data/rt.stan @@ -6,4 +6,4 @@ int future_fixed; // is underlying future Rt assumed to be fixed int fixed_from; // Reference date for when Rt estimation should be fixed int use_pop; // use population size - int gt_id; // id of generation time \ No newline at end of file + int gt_id; // id of generation time diff --git a/inst/stan/data/simulation_rt.stan b/inst/stan/data/simulation_rt.stan index 3beae3161..fdae11e2a 100644 --- a/inst/stan/data/simulation_rt.stan +++ b/inst/stan/data/simulation_rt.stan @@ -2,6 +2,6 @@ array[n, seeding_time > 1 ? 1 : 0] real initial_growth; //initial growth matrix[n, t - seeding_time] R; // reproduction number - int pop; // susceptible population + int use_pop; // use population size int gt_id; // id of generation time diff --git a/inst/stan/simulate_infections.stan b/inst/stan/simulate_infections.stan index f06006f6b..83fb14123 100644 --- a/inst/stan/simulate_infections.stan +++ b/inst/stan/simulate_infections.stan @@ -51,7 +51,7 @@ generated quantities { params_value, params ); - real pop = get_param( + vector[n] pop = get_param( pop_id, params_fixed_lookup, params_variable_lookup, params_value, params ); @@ -68,7 +68,7 @@ generated quantities { infections[i] = to_row_vector(generate_infections( to_vector(R[i]), seeding_time, gt_rev_pmf, initial_infections[i], - initial_growth[i], pop, use_pop, future_time, obs_scale, frac_obs[i] + initial_growth[i], pop[i], use_pop, future_time, obs_scale, frac_obs[i] )); if (delay_id) { diff --git a/man/EpiNow2-package.Rd b/man/EpiNow2-package.Rd index b0c820c1b..5e27df68a 100644 --- a/man/EpiNow2-package.Rd +++ b/man/EpiNow2-package.Rd @@ -43,7 +43,7 @@ Other contributors: \item Paul Mee \email{paul.mee@lshtm.ac.uk} [contributor] \item Peter Ellis \email{peter.ellis2013nz@gmail.com} [contributor] \item Pietro Monticone \email{pietro.monticone@edu.unito.it} [contributor] - \item Lloyd Chapman \email{lloyd.chapman1@lshtm.ac.uk} [contributor] + \item Lloyd Chapman \email{lloyd.chapman1@lshtm.ac.uk } [contributor] \item Andrew Johnson \email{andrew.johnson@arjohnsonau.com} [contributor] \item Kaitlyn Johnson \email{johnsonkaitlyne9@gmail.com} (\href{https://orcid.org/0000-0001-8011-0012}{ORCID}) [contributor] } diff --git a/man/rt_opts.Rd b/man/rt_opts.Rd index c48cfe21e..3155f9101 100644 --- a/man/rt_opts.Rd +++ b/man/rt_opts.Rd @@ -11,7 +11,7 @@ rt_opts( use_breakpoints = TRUE, future = "latest", gp_on = c("R_t-1", "R0"), - pop = 0 + pop = Fixed(0) ) } \arguments{ @@ -50,9 +50,9 @@ Rt = Rt-1 * GP), and applying the Gaussian process to a global mean (i.e Rt but the method relying on a global mean will revert to this for real time estimates, which may not be desirable.} -\item{pop}{Integer, defaults to 0. Susceptible population initially present. -Used to adjust Rt estimates when otherwise fixed based on the proportion of -the population that is susceptible. When set to 0 no population adjustment +\item{pop}{A \verb{} giving the initial susceptible population size. +Used to adjust Rt estimates based on the proportion of the population that +is susceptible. Defaults to \code{Fixed(0)} which means no population adjustment is done.} } \value{ diff --git a/man/simulate_infections.Rd b/man/simulate_infections.Rd index bc012013f..e0e34ef23 100644 --- a/man/simulate_infections.Rd +++ b/man/simulate_infections.Rd @@ -16,7 +16,7 @@ simulate_infections( CrIs = c(0.2, 0.5, 0.9), backend = "rstan", seeding_time = NULL, - pop = 0, + pop = Fixed(0), ... ) } @@ -74,9 +74,9 @@ R given) of \code{seeding_time} days is assumed to have followed exponential growth roughly in line with the growth rate implied by the first value of R.} -\item{pop}{Integer, defaults to 0. Susceptible population initially present. -Used to adjust Rt estimates when otherwise fixed based on the proportion of -the population that is susceptible. When set to 0 no population adjustment +\item{pop}{A \verb{} giving the initial susceptible population size. +Used to adjust Rt estimates based on the proportion of the population that +is susceptible. Defaults to \code{Fixed(0)} which means no population adjustment is done.} \item{...}{deprecated; only included for backward compatibility} diff --git a/tests/testthat/test-create_rt_date.R b/tests/testthat/test-create_rt_date.R index 3fc47b1cc..c63020f97 100644 --- a/tests/testthat/test-create_rt_date.R +++ b/tests/testthat/test-create_rt_date.R @@ -27,13 +27,13 @@ test_that("create_rt_data handles custom rt_opts correctly", { use_breakpoints = FALSE, future = "project", gp_on = "R0", - pop = 1000000 + pop = Normal(mean = 1000000, sd = 100) ) result <- create_rt_data(rt = custom_rt, horizon = 7) expect_equal(result$estimate_r, 0) - expect_equal(result$pop, 1000000) + expect_equal(result$use_pop, 1) expect_equal(result$stationary, 1) expect_equal(result$future_time, 7) }) diff --git a/tests/testthat/test-rt_opts.R b/tests/testthat/test-rt_opts.R index 69fb1183c..e35c5824c 100644 --- a/tests/testthat/test-rt_opts.R +++ b/tests/testthat/test-rt_opts.R @@ -19,7 +19,7 @@ test_that("rt_opts handles custom inputs correctly", { use_breakpoints = FALSE, future = "project", gp_on = "R0", - pop = 1000000 + pop = Normal(mean = 1000000, sd = 100) )) expect_null(result$prior) @@ -27,10 +27,17 @@ test_that("rt_opts handles custom inputs correctly", { expect_equal(result$rw, 7) expect_true(result$use_breakpoints) # Should be TRUE when rw > 0 expect_equal(result$future, "project") - expect_equal(result$pop, 1000000) + expect_equal(result$pop, Normal(mean = 1000000, sd = 100)) expect_equal(result$gp_on, "R0") }) +test_that("rt_opts warns when pop is passed as numeric", { + expect_warning( + rt_opts(pop = 1000), + "Specifying `pop` as a numeric value is deprecated" + ) +}) + test_that("rt_opts sets use_breakpoints to TRUE when rw > 0", { result <- rt_opts(rw = 3, use_breakpoints = FALSE) expect_true(result$use_breakpoints) @@ -59,8 +66,7 @@ test_that("rt_opts returns object of correct class", { }) test_that("rt_opts handles edge cases correctly", { - result <- rt_opts(rw = 0.1, pop = -1) + result <- rt_opts(rw = 0.1) expect_equal(result$rw, 0.1) - expect_equal(result$pop, -1) expect_true(result$use_breakpoints) }) diff --git a/vignettes/estimate_infections_options.Rmd.orig b/vignettes/estimate_infections_options.Rmd.orig index 0cd532013..622e47aa2 100644 --- a/vignettes/estimate_infections_options.Rmd.orig +++ b/vignettes/estimate_infections_options.Rmd.orig @@ -154,7 +154,7 @@ dep <- estimate_infections(reported_cases, delays = delay_opts(delay), rt = rt_opts( prior = rt_prior, - pop = 1000000, future = "latest" + pop = Normal(mean = 1000000, sd = 1000), future = "latest" ) ) # summarise results From 6b1cb519a228b9a1162ae9c02455957e8f4ec532 Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 19 Dec 2024 18:23:44 +0000 Subject: [PATCH 04/16] update news --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index a37bfc384..32724e466 100644 --- a/NEWS.md +++ b/NEWS.md @@ -15,6 +15,7 @@ estimate_infections() ``` +- Added support for fitting the susceptible population size. By @seabbs in #904 and reviewed by @sbfnk. - A bug was fixed where the initial growth was never estimated (i.e. the prior mean was always zero). By @sbfnk in #853 and reviewed by @seabbs. - A bug was fixed where an internal function for applying a default cdf cutoff failed due to a difference a vector length issue. By @jamesmbaazam in #858 and reviewed by @sbfnk. - All parameters have been changed to the new parameter interface. By @sbfnk in #871 and #890 and reviewed by @seabbs. From c17fdcb84eb4c41b6a42b89212944dc7c41638eb Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 19 Dec 2024 18:47:00 +0000 Subject: [PATCH 05/16] update tests --- tests/testthat/test-create_rt_date.R | 2 +- tests/testthat/test-rt_opts.R | 2 +- vignettes/EpiNow2.Rmd.orig | 3 ++- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/testthat/test-create_rt_date.R b/tests/testthat/test-create_rt_date.R index c63020f97..463c26261 100644 --- a/tests/testthat/test-create_rt_date.R +++ b/tests/testthat/test-create_rt_date.R @@ -7,7 +7,7 @@ test_that("create_rt_data returns expected default values", { expect_equal(result$breakpoints, numeric(0)) expect_equal(result$future_fixed, 1) expect_equal(result$fixed_from, 0) - expect_equal(result$pop, 0) + expect_equal(result$use_pop, 0) expect_equal(result$stationary, 0) expect_equal(result$future_time, 0) }) diff --git a/tests/testthat/test-rt_opts.R b/tests/testthat/test-rt_opts.R index e35c5824c..c8700d9df 100644 --- a/tests/testthat/test-rt_opts.R +++ b/tests/testthat/test-rt_opts.R @@ -7,7 +7,7 @@ test_that("rt_opts returns expected default values", { expect_equal(result$rw, 0) expect_true(result$use_breakpoints) expect_equal(result$future, "latest") - expect_equal(result$pop, 0) + expect_equal(result$pop, Fixed(0)) expect_equal(result$gp_on, "R_t-1") }) diff --git a/vignettes/EpiNow2.Rmd.orig b/vignettes/EpiNow2.Rmd.orig index 69d10bcdd..e46f53aa5 100644 --- a/vignettes/EpiNow2.Rmd.orig +++ b/vignettes/EpiNow2.Rmd.orig @@ -95,7 +95,8 @@ estimates <- epinow( 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)), - stan = stan_opts(cores = 4, control = list(adapt_delta = 0.99)), + stan = stan_opts( + cores = 4, control = list(adapt_delta = 0.99), backend = "cmdstanr"), verbose = interactive() ) names(estimates) From 2967e056d46cfece5d6887110e1ed31f192c60bb Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 19 Dec 2024 21:03:59 +0000 Subject: [PATCH 06/16] handle pop more carefully in a simulation setting --- R/simulate_infections.R | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/R/simulate_infections.R b/R/simulate_infections.R index 25919d46a..8f8531ef4 100644 --- a/R/simulate_infections.R +++ b/R/simulate_infections.R @@ -86,7 +86,6 @@ simulate_infections <- function(estimates, R, initial_infections, assert_numeric(R$R, lower = 0) assert_numeric(initial_infections, lower = 0) assert_numeric(day_of_week_effect, lower = 0, null.ok = TRUE) - assert_numeric(pop, lower = 0) if (!is.null(seeding_time)) { assert_integerish(seeding_time, lower = 1) } @@ -94,6 +93,7 @@ simulate_infections <- function(estimates, R, initial_infections, assert_class(truncation, "trunc_opts") assert_class(obs, "obs_opts") assert_class(generation_time, "generation_time_opts") + assert_class(pop, "dist_spec") ## create R for all dates modelled all_dates <- data.table(date = seq.Date(min(R$date), max(R$date), by = "day")) @@ -179,7 +179,8 @@ simulate_infections <- function(estimates, R, initial_infections, rho = NULL, R0 = NULL, frac_obs = obs$scale, - rep_phi = obs$phi + rep_phi = obs$phi, + pop = pop )) ## set empty params matrix - variable parameters not supported here data$params <- array(dim = c(1, 0)) From 86018b12ce1aaca22486e0c097268177835df902 Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Fri, 20 Dec 2024 14:25:21 +0000 Subject: [PATCH 07/16] update stan generate infection tests --- tests/testthat/test-stan-infections.R | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/tests/testthat/test-stan-infections.R b/tests/testthat/test-stan-infections.R index d08e06c01..0f8b631c0 100644 --- a/tests/testthat/test-stan-infections.R +++ b/tests/testthat/test-stan-infections.R @@ -25,35 +25,35 @@ gt_rev_pmf <- get_delay_rev_pmf( # test generate infections test_that("generate_infections works as expected", { expect_equal( - round(generate_infections(c(1, rep(1, 9)), 10, gt_rev_pmf, log(1000), 0, 0, 0, 0, 0), 0), + round(generate_infections(c(1, rep(1, 9)), 10, gt_rev_pmf, log(1000), 0, 0, 0, 0, 0, 0), 0), c(rep(1000, 10), 995, 996, rep(997, 8)) ) expect_equal( - round(generate_infections(c(1, rep(1.1, 9)), 10, gt_rev_pmf, log(20), 0.03, 0, 0, 0, 0), 0), + round(generate_infections(c(1, rep(1.1, 9)), 10, gt_rev_pmf, log(20), 0.03, 0, 0, 0, 0, 0), 0), c(20, 21, 21, 22, 23, 23, 24, 25, 25, 26, 24, 27, 28, 29, 30, 30, 31, 32, 33, 34) ) expect_equal( - round(generate_infections(c(1, rep(1.1, 9)), 10, gt_rev_pmf, log(100), 0, 0, 0, 0, 0), 0), + round(generate_infections(c(1, rep(1.1, 9)), 10, gt_rev_pmf, log(100), 0, 0, 0, 0, 0, 0), 0), c(rep(100, 10), 99, 110, 112, 115, 119, 122, 126, 130, 134, 138) ) expect_equal( - round(generate_infections(c(1, rep(1, 9)), 4, gt_rev_pmf, log(500), -0.02, 0, 0, 0, 0), 0), + round(generate_infections(c(1, rep(1, 9)), 4, gt_rev_pmf, log(500), -0.02, 0, 0, 0, 0, 0), 0), c(500, 490, 480, 471, 382, 403, 408, rep(409, 7)) ) expect_equal( - round(generate_infections(c(1, rep(1.1, 9)), 4, gt_rev_pmf, log(500), 0, 0, 0, 0, 0), 0), + round(generate_infections(c(1, rep(1.1, 9)), 4, gt_rev_pmf, log(500), 0, 0, 0, 0, 0, 0), 0), c(rep(500, 4), 394, 460, 475, 489, 505, 520, 536, 553, 570, 588) ) expect_equal( - round(generate_infections(c(1, rep(1, 9)), 1, gt_rev_pmf, log(40), numeric(0), 0, 0, 0, 0), 0), + round(generate_infections(c(1, rep(1, 9)), 1, gt_rev_pmf, log(40), numeric(0), 0, 0, 0, 0, 0), 0), c(40, 8, 11, 12, 12, rep(13, 6)) ) expect_equal( - round(generate_infections(c(1, rep(1.1, 9)), 1, gt_rev_pmf, log(100), 0.01, 0, 0, 0, 0), 0), + round(generate_infections(c(1, rep(1.1, 9)), 1, gt_rev_pmf, log(100), 0.01, 0, 0, 0, 0, 0), 0), c(100, 20, 31, 35, 36, 37, 38, 39, 41, 42, 43) ) expect_equal( - round(generate_infections(c(1, rep(1, 9)), 10, gt_rev_pmf, log(1000), 0, 100000, 4, 0, 0), 0), + round(generate_infections(c(1, rep(1, 9)), 10, gt_rev_pmf, log(1000), 0, 100000, 4, 0, 0, 0), 0), c(rep(1000, 10), 995, 996, rep(997, 4), 980, 965, 947, 926) ) }) From 07064bdac188f975de08d4c8518ffcbdde8a50bc Mon Sep 17 00:00:00 2001 From: Sam Date: Fri, 20 Dec 2024 16:44:50 +0000 Subject: [PATCH 08/16] fix linting --- R/opts.R | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/R/opts.R b/R/opts.R index c3d87cff8..e0a478d97 100644 --- a/R/opts.R +++ b/R/opts.R @@ -734,7 +734,10 @@ obs_opts <- function(family = c("negbin", "poisson"), cli_abort( c( "!" = "Specifying {.var phi} as a vector of length 2 is deprecated.", - "i" = "Use a {.cls dist_spec} instead, e.g. Normal(mean = {phi[1]}, sd = {phi[2]})." + "i" = paste0( + "Use a {.cls dist_spec} instead, e.g. Normal(mean = {phi[1]}, ", + "sd = {phi[2]})." + ) ) ) } From 3bc26d4c3efc83863d7bb8d6a9a87d8a1897b77e Mon Sep 17 00:00:00 2001 From: Sam Date: Fri, 20 Dec 2024 17:25:17 +0000 Subject: [PATCH 09/16] make it possible to either use pop out of sample or fit with it in sample --- R/create.R | 2 +- R/opts.R | 27 +++++++++++++++++++++------ inst/stan/functions/infections.stan | 2 +- man/rt_opts.Rd | 10 +++++++++- 4 files changed, 32 insertions(+), 9 deletions(-) diff --git a/R/create.R b/R/create.R index ef7c02237..97ae12fda 100644 --- a/R/create.R +++ b/R/create.R @@ -263,7 +263,7 @@ create_rt_data <- function(rt = rt_opts(), breakpoints = NULL, breakpoints = breakpoints, future_fixed = as.numeric(future_rt$fixed), fixed_from = future_rt$from, - use_pop = as.integer(rt$pop != Fixed(0)), + use_pop = as.integer(rt$pop != Fixed(0)) + as.integer(rt$estimate_pop), stationary = as.numeric(rt$gp_on == "R0"), future_time = horizon - future_rt$from ) diff --git a/R/opts.R b/R/opts.R index e0a478d97..6dc1f545f 100644 --- a/R/opts.R +++ b/R/opts.R @@ -324,6 +324,13 @@ trunc_opts <- function(dist = Fixed(0), default_cdf_cutoff = 0.001, #' Used to adjust Rt estimates based on the proportion of the population that #' is susceptible. Defaults to `Fixed(0)` which means no population adjustment #' is done. +#' +#' @param pop_within_horizon Logical, defaults to `FALSE`. Should the +#' susceptible population adjustment be applied within the time horizon of data +#' or just beyond it. Note that if `pop_within_horizon = TRUE` the Rt estimate +#' will be unadjusted for susceptible depletion but the resulting posterior +#' predictions for infections and cases will be adjusted for susceptible +#' depletion. #' #' @param gp_on Character string, defaulting to "R_t-1". Indicates how the #' Gaussian process, if in use, should be applied to Rt. Currently supported @@ -354,7 +361,8 @@ rt_opts <- function(prior = LogNormal(mean = 1, sd = 1), use_breakpoints = TRUE, future = "latest", gp_on = c("R_t-1", "R0"), - pop = Fixed(0)) { + pop = Fixed(0), + pop_within_horizon = FALSE) { rt <- list( use_rt = use_rt, rw = rw, @@ -388,15 +396,22 @@ rt_opts <- function(prior = LogNormal(mean = 1, sd = 1), } if (is.numeric(pop)) { - cli_warn( - c( - "!" = "Specifying {.var pop} as a numeric value is deprecated.", - "i" = "Use a {.cls dist_spec} instead, e.g. Fixed({pop})." - ) + lifecycle::deprecate_warn( + "1.7.0", + "rt_opts(pop = 'must be a ``')", + details = "For specifying a fixed population size, use `Fixed(pop)`" ) pop <- Fixed(pop) } rt$pop <- pop + if (isTRUE(pop_within_horizon) && pop == Fixed(0)) { + cli_abort( + c( + "!" = "pop_within_horizon = TRUE but pop is fixed at 0." + ) + ) + } + rt$estimate_pop <- TRUE if (rt$use_rt) { rt$prior <- prior diff --git a/inst/stan/functions/infections.stan b/inst/stan/functions/infections.stan index f5e9f78c4..2b626f4ad 100644 --- a/inst/stan/functions/infections.stan +++ b/inst/stan/functions/infections.stan @@ -48,7 +48,7 @@ vector generate_infections(vector oR, int uot, vector gt_rev_pmf, // iteratively update infections for (s in 1:ot) { infectiousness[s] = update_infectiousness(infections, gt_rev_pmf, uot, s); - if (use_pop && s > nht) { + if (use_pop == 2 && s > nht) { exp_adj_Rt = exp(-R[s] * infectiousness[s] / (pop - cum_infections[nht])); exp_adj_Rt = exp_adj_Rt > 1 ? 1 : exp_adj_Rt; infections[s + uot] = (pop - cum_infections[s]) * (1 - exp_adj_Rt); diff --git a/man/rt_opts.Rd b/man/rt_opts.Rd index 3155f9101..358117438 100644 --- a/man/rt_opts.Rd +++ b/man/rt_opts.Rd @@ -11,7 +11,8 @@ rt_opts( use_breakpoints = TRUE, future = "latest", gp_on = c("R_t-1", "R0"), - pop = Fixed(0) + pop = Fixed(0), + pop_within_horizon = FALSE ) } \arguments{ @@ -54,6 +55,13 @@ estimates, which may not be desirable.} Used to adjust Rt estimates based on the proportion of the population that is susceptible. Defaults to \code{Fixed(0)} which means no population adjustment is done.} + +\item{pop_within_horizon}{Logical, defaults to \code{FALSE}. Should the +susceptible population adjustment be applied within the time horizon of data +or just beyond it. Note that if \code{pop_within_horizon = TRUE} the Rt estimate +will be unadjusted for susceptible depletion but the resulting posterior +predictions for infections and cases will be adjusted for susceptible +depletion.} } \value{ An \verb{} object with settings defining the time-varying From 01f42c96ace216989491500dbdb78f0770f09941 Mon Sep 17 00:00:00 2001 From: Sam Abbott Date: Fri, 10 Jan 2025 11:44:05 +0000 Subject: [PATCH 10/16] add support for in sample suscetpble --- inst/stan/functions/infections.stan | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inst/stan/functions/infections.stan b/inst/stan/functions/infections.stan index 2b626f4ad..ee379e749 100644 --- a/inst/stan/functions/infections.stan +++ b/inst/stan/functions/infections.stan @@ -48,7 +48,7 @@ vector generate_infections(vector oR, int uot, vector gt_rev_pmf, // iteratively update infections for (s in 1:ot) { infectiousness[s] = update_infectiousness(infections, gt_rev_pmf, uot, s); - if (use_pop == 2 && s > nht) { + if (use_pop == 1 || (use_pop == 2 && s <= nht)) { exp_adj_Rt = exp(-R[s] * infectiousness[s] / (pop - cum_infections[nht])); exp_adj_Rt = exp_adj_Rt > 1 ? 1 : exp_adj_Rt; infections[s + uot] = (pop - cum_infections[s]) * (1 - exp_adj_Rt); From e21a5f8ce95f539a040fdb25781098769f0536db Mon Sep 17 00:00:00 2001 From: Sam Abbott Date: Fri, 10 Jan 2025 13:05:01 +0000 Subject: [PATCH 11/16] update docs --- R/opts.R | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/R/opts.R b/R/opts.R index 6dc1f545f..a37b3a63e 100644 --- a/R/opts.R +++ b/R/opts.R @@ -325,12 +325,11 @@ trunc_opts <- function(dist = Fixed(0), default_cdf_cutoff = 0.001, #' is susceptible. Defaults to `Fixed(0)` which means no population adjustment #' is done. #' -#' @param pop_within_horizon Logical, defaults to `FALSE`. Should the -#' susceptible population adjustment be applied within the time horizon of data -#' or just beyond it. Note that if `pop_within_horizon = TRUE` the Rt estimate -#' will be unadjusted for susceptible depletion but the resulting posterior -#' predictions for infections and cases will be adjusted for susceptible -#' depletion. +#' @param pop_period Character string, defaulting to "forecast". Controls when +#' susceptible population adjustment is applied. "forecast" only applies the +#' adjustment to forecasts while "all" applies it to both data and forecasts. +#' Note that with "forecast", Rt estimates are unadjusted for susceptible +#' depletion but posterior predictions are adjusted. #' #' @param gp_on Character string, defaulting to "R_t-1". Indicates how the #' Gaussian process, if in use, should be applied to Rt. Currently supported @@ -362,13 +361,14 @@ rt_opts <- function(prior = LogNormal(mean = 1, sd = 1), future = "latest", gp_on = c("R_t-1", "R0"), pop = Fixed(0), - pop_within_horizon = FALSE) { + pop_period = c("forecast", "all")) { rt <- list( use_rt = use_rt, rw = rw, use_breakpoints = use_breakpoints, future = future, - gp_on = arg_match(gp_on) + gp_on = arg_match(gp_on), + pop_period = arg_match(pop_period) ) # replace default settings with those specified by user @@ -396,18 +396,18 @@ rt_opts <- function(prior = LogNormal(mean = 1, sd = 1), } if (is.numeric(pop)) { - lifecycle::deprecate_warn( - "1.7.0", - "rt_opts(pop = 'must be a ``')", + lifecycle::deprecate_warn( + "1.7.0", + "rt_opts(pop = 'must be a ``')", details = "For specifying a fixed population size, use `Fixed(pop)`" ) pop <- Fixed(pop) } rt$pop <- pop - if (isTRUE(pop_within_horizon) && pop == Fixed(0)) { + if (rt$pop_period == "all" && pop == Fixed(0)) { cli_abort( c( - "!" = "pop_within_horizon = TRUE but pop is fixed at 0." + "!" = "pop_period = \"all\" but pop is fixed at 0." ) ) } From 101f1c56f9401787c0d41f38c6dc01e431606047 Mon Sep 17 00:00:00 2001 From: GitHub Actions Date: Fri, 10 Jan 2025 13:09:25 +0000 Subject: [PATCH 12/16] Document --- man/EpiNow2-package.Rd | 2 +- man/create_forecast_data.Rd | 2 +- man/create_stan_data.Rd | 2 +- man/epinow.Rd | 2 +- man/estimate_infections.Rd | 2 +- man/forecast_opts.Rd | 1 - man/regional_epinow.Rd | 2 +- man/rt_opts.Rd | 13 ++++++------- 8 files changed, 12 insertions(+), 14 deletions(-) diff --git a/man/EpiNow2-package.Rd b/man/EpiNow2-package.Rd index 5e27df68a..b0c820c1b 100644 --- a/man/EpiNow2-package.Rd +++ b/man/EpiNow2-package.Rd @@ -43,7 +43,7 @@ Other contributors: \item Paul Mee \email{paul.mee@lshtm.ac.uk} [contributor] \item Peter Ellis \email{peter.ellis2013nz@gmail.com} [contributor] \item Pietro Monticone \email{pietro.monticone@edu.unito.it} [contributor] - \item Lloyd Chapman \email{lloyd.chapman1@lshtm.ac.uk } [contributor] + \item Lloyd Chapman \email{lloyd.chapman1@lshtm.ac.uk} [contributor] \item Andrew Johnson \email{andrew.johnson@arjohnsonau.com} [contributor] \item Kaitlyn Johnson \email{johnsonkaitlyne9@gmail.com} (\href{https://orcid.org/0000-0001-8011-0012}{ORCID}) [contributor] } diff --git a/man/create_forecast_data.Rd b/man/create_forecast_data.Rd index 3ca7b0ccb..38a98a37e 100644 --- a/man/create_forecast_data.Rd +++ b/man/create_forecast_data.Rd @@ -9,7 +9,7 @@ create_forecast_data(forecast = forecast_opts(), data) \arguments{ \item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no -forecasting will be one.} +forecasting will be done.} \item{data}{A \verb{} of confirmed cases (confirm) by date (date). \code{confirm} must be numeric and \code{date} must be in date format. Optionally diff --git a/man/create_stan_data.Rd b/man/create_stan_data.Rd index e0f05ff26..e5680e01a 100644 --- a/man/create_stan_data.Rd +++ b/man/create_stan_data.Rd @@ -48,7 +48,7 @@ define the back calculation. Defaults to \code{\link[=backcalc_opts]{backcalc_op \item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no -forecasting will be one.} +forecasting will be done.} } \value{ A list of stan data diff --git a/man/epinow.Rd b/man/epinow.Rd index 8a62820a5..008759969 100644 --- a/man/epinow.Rd +++ b/man/epinow.Rd @@ -73,7 +73,7 @@ observation model. Defaults to \code{\link[=obs_opts]{obs_opts()}}.} \item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no -forecasting will be one.} +forecasting will be done.} \item{stan}{A list of stan options as generated by \code{\link[=stan_opts]{stan_opts()}}. Defaults to \code{\link[=stan_opts]{stan_opts()}}. Can be used to override \code{data}, \code{init}, and \code{verbose} diff --git a/man/estimate_infections.Rd b/man/estimate_infections.Rd index 992b65d46..4f5a97dde 100644 --- a/man/estimate_infections.Rd +++ b/man/estimate_infections.Rd @@ -69,7 +69,7 @@ observation model. Defaults to \code{\link[=obs_opts]{obs_opts()}}.} \item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no -forecasting will be one.} +forecasting will be done.} \item{stan}{A list of stan options as generated by \code{\link[=stan_opts]{stan_opts()}}. Defaults to \code{\link[=stan_opts]{stan_opts()}}. Can be used to override \code{data}, \code{init}, and \code{verbose} diff --git a/man/forecast_opts.Rd b/man/forecast_opts.Rd index 2bd02766b..cadd823aa 100644 --- a/man/forecast_opts.Rd +++ b/man/forecast_opts.Rd @@ -17,7 +17,6 @@ forecasts unless set explicitly here.} } \value{ A \verb{} object of forecast setting. -rstan functions. } \description{ \ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} diff --git a/man/regional_epinow.Rd b/man/regional_epinow.Rd index efaa0ee6b..1dfd319df 100644 --- a/man/regional_epinow.Rd +++ b/man/regional_epinow.Rd @@ -66,7 +66,7 @@ observation model. Defaults to \code{\link[=obs_opts]{obs_opts()}}.} \item{forecast}{A list of options as generated by \code{\link[=forecast_opts]{forecast_opts()}} defining the forecast opitions. Defaults to \code{\link[=forecast_opts]{forecast_opts()}}. If NULL then no -forecasting will be one.} +forecasting will be done.} \item{stan}{A list of stan options as generated by \code{\link[=stan_opts]{stan_opts()}}. Defaults to \code{\link[=stan_opts]{stan_opts()}}. Can be used to override \code{data}, \code{init}, and \code{verbose} diff --git a/man/rt_opts.Rd b/man/rt_opts.Rd index 358117438..8636fdb82 100644 --- a/man/rt_opts.Rd +++ b/man/rt_opts.Rd @@ -12,7 +12,7 @@ rt_opts( future = "latest", gp_on = c("R_t-1", "R0"), pop = Fixed(0), - pop_within_horizon = FALSE + pop_period = c("forecast", "all") ) } \arguments{ @@ -56,12 +56,11 @@ Used to adjust Rt estimates based on the proportion of the population that is susceptible. Defaults to \code{Fixed(0)} which means no population adjustment is done.} -\item{pop_within_horizon}{Logical, defaults to \code{FALSE}. Should the -susceptible population adjustment be applied within the time horizon of data -or just beyond it. Note that if \code{pop_within_horizon = TRUE} the Rt estimate -will be unadjusted for susceptible depletion but the resulting posterior -predictions for infections and cases will be adjusted for susceptible -depletion.} +\item{pop_period}{Character string, defaulting to "forecast". Controls when +susceptible population adjustment is applied. "forecast" only applies the +adjustment to forecasts while "all" applies it to both data and forecasts. +Note that with "forecast", Rt estimates are unadjusted for susceptible +depletion but posterior predictions are adjusted.} } \value{ An \verb{} object with settings defining the time-varying From 27097cb0ccb718273e9a6cf75fe37a050f8e1896 Mon Sep 17 00:00:00 2001 From: Sam Abbott Date: Fri, 10 Jan 2025 13:23:17 +0000 Subject: [PATCH 13/16] touch From bafef0b9c6dde4c1a6ef982532ba3be52d2fb639 Mon Sep 17 00:00:00 2001 From: Sam Abbott Date: Fri, 10 Jan 2025 13:53:43 +0000 Subject: [PATCH 14/16] cleaner opt passing --- R/create.R | 3 ++- R/opts.R | 1 - 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/create.R b/R/create.R index 97ae12fda..446336c37 100644 --- a/R/create.R +++ b/R/create.R @@ -263,7 +263,8 @@ create_rt_data <- function(rt = rt_opts(), breakpoints = NULL, breakpoints = breakpoints, future_fixed = as.numeric(future_rt$fixed), fixed_from = future_rt$from, - use_pop = as.integer(rt$pop != Fixed(0)) + as.integer(rt$estimate_pop), + use_pop = + as.integer(rt$pop != Fixed(0)) + as.integer(rt$pop_period == "all"), stationary = as.numeric(rt$gp_on == "R0"), future_time = horizon - future_rt$from ) diff --git a/R/opts.R b/R/opts.R index a37b3a63e..8c63589ab 100644 --- a/R/opts.R +++ b/R/opts.R @@ -411,7 +411,6 @@ rt_opts <- function(prior = LogNormal(mean = 1, sd = 1), ) ) } - rt$estimate_pop <- TRUE if (rt$use_rt) { rt$prior <- prior From 841d23881db18554d7d80c89882156f5ddbff59d Mon Sep 17 00:00:00 2001 From: Sam Abbott Date: Fri, 10 Jan 2025 15:46:58 +0000 Subject: [PATCH 15/16] update tests --- tests/testthat/test-rt_opts.R | 2 +- tests/testthat/test-stan-infections.R | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-rt_opts.R b/tests/testthat/test-rt_opts.R index c8700d9df..d76f0571b 100644 --- a/tests/testthat/test-rt_opts.R +++ b/tests/testthat/test-rt_opts.R @@ -34,7 +34,7 @@ test_that("rt_opts handles custom inputs correctly", { test_that("rt_opts warns when pop is passed as numeric", { expect_warning( rt_opts(pop = 1000), - "Specifying `pop` as a numeric value is deprecated" + "The `pop` argument of `rt_opts()` must be a `` as of EpiNow2 1.7.0." ) }) diff --git a/tests/testthat/test-stan-infections.R b/tests/testthat/test-stan-infections.R index 0f8b631c0..9d88a858a 100644 --- a/tests/testthat/test-stan-infections.R +++ b/tests/testthat/test-stan-infections.R @@ -53,7 +53,7 @@ test_that("generate_infections works as expected", { c(100, 20, 31, 35, 36, 37, 38, 39, 41, 42, 43) ) expect_equal( - round(generate_infections(c(1, rep(1, 9)), 10, gt_rev_pmf, log(1000), 0, 100000, 4, 0, 0, 0), 0), + round(generate_infections(c(1, rep(1, 9)), 10, gt_rev_pmf, log(1000), 0, 100000, 2, 0, 0, 0), 0), c(rep(1000, 10), 995, 996, rep(997, 4), 980, 965, 947, 926) ) }) From 6184ab9ca4b884b13adf87ae7611d71e7e0e1ada Mon Sep 17 00:00:00 2001 From: Sam Abbott Date: Fri, 10 Jan 2025 15:57:40 +0000 Subject: [PATCH 16/16] clean up whitespace --- R/opts.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/opts.R b/R/opts.R index 8c63589ab..c6e846bf8 100644 --- a/R/opts.R +++ b/R/opts.R @@ -324,7 +324,7 @@ trunc_opts <- function(dist = Fixed(0), default_cdf_cutoff = 0.001, #' Used to adjust Rt estimates based on the proportion of the population that #' is susceptible. Defaults to `Fixed(0)` which means no population adjustment #' is done. -#' +#' #' @param pop_period Character string, defaulting to "forecast". Controls when #' susceptible population adjustment is applied. "forecast" only applies the #' adjustment to forecasts while "all" applies it to both data and forecasts.