Issue #408: Fit the susceptible population #904

seabbs · 2024-12-19T18:20:57Z

Description

This PR closes #408 by allowing for fitting the susceptible population using the new distribution interface.

Maybe to do:

it does not add a inference unit test but we didn't have one before for susceptibility
This doesn't address if we should be reporting the susceptibility adjusted or unadjusted Rt.

Initial submission checklist

My PR is based on a package issue and I have explicitly linked it.
I have tested my changes locally (using devtools::test() and devtools::check()).
I have added or updated unit tests where necessary.
I have updated the documentation if required and rebuilt docs if yes (using devtools::document()).
I have followed the established coding standards (and checked using lintr::lint_package()).
I have added a news item linked to this PR.

After the initial Pull Request

I have reviewed Checks for this PR and addressed any issues as far as I am able.

seabbs · 2024-12-20T16:50:06Z

Maybe to do:

I am minded not to do any of these unless you feel strongly

sbfnk

This is really nice!

This doesn't address if we should be reporting the susceptibility adjusted or unadjusted Rt.

We already have a gen_R in generated quantities (for the nonmechanistic model) - so we could in principle model unadj_R and then always generate R in generated quantities?

R/opts.R

inst/stan/functions/infections.stan

sbfnk · 2024-12-20T16:57:20Z

Maybe to do:

I am minded not to do any of these unless you feel strongly

I don't feel strongly - not quite sure what's best to report anyway. Probably worth a new Issue though as we might be calling stuff R when it's not actually R.

R/opts.R

inst/stan/functions/infections.stan

…mple

github-actions · 2025-01-10T14:54:11Z

This is how benchmark results would change (along with a 95% confidence interval in relative change) if bafef0b is merged into main:

❗🐌default: 20.1s -> 22.5s [+1.1%, +22.44%]
✔️no_delays: 25.1s -> 23.8s [-15.01%, +4.71%]
✔️random_walk: 9.04s -> 13.3s [-32.12%, +126.71%]
✔️stationary: 11.6s -> 13.2s [-10.22%, +38.66%]
✔️uncertain: 34.9s -> 33.9s [-12.84%, +7.14%]
Further explanation regarding interpretation and methodology can be found in the documentation.

github-actions · 2025-01-10T16:57:07Z

This is how benchmark results would change (along with a 95% confidence interval in relative change) if 6184ab9 is merged into main:

✔️default: 22.1s -> 21s [-19.07%, +8.83%]
✔️no_delays: 25.2s -> 24.6s [-15.87%, +11.05%]
✔️random_walk: 9.17s -> 9.18s [-10.79%, +11.02%]
✔️stationary: 12s -> 13s [-11.45%, +27.35%]
✔️uncertain: 34.5s -> 35.4s [-7.75%, +13.15%]
Further explanation regarding interpretation and methodology can be found in the documentation.

sbfnk · 2025-01-13T11:43:43Z

R/opts.R

+#' Note that with "forecast", Rt estimates are unadjusted for susceptible
+#' depletion but posterior predictions are adjusted.


Suggested change

#' Note that with "forecast", Rt estimates are unadjusted for susceptible

#' depletion but posterior predictions are adjusted.

#' Note that with "all", Rt estimates are unadjusted for susceptible

#' depletion but posterior predictions of infections and reports are

#' adjusted.

sbfnk · 2025-01-13T11:46:25Z

R/simulate_infections.R

@@ -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))


can be a future issue but I think we want to be able to have susceptible depletion across full simulations, too.

sbfnk · 2025-01-13T11:46:50Z

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 use_pop;                       // use population size


Suggested change

int use_pop; // use population size

int use_pop; // use population size (0 = no; 1 = forecasts; 2 = all)

sbfnk · 2025-01-13T11:47:03Z

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


Suggested change

int use_pop; // use population size

int use_pop; // use population size (0 = no; 1 = forecasts)

sbfnk · 2025-01-13T11:51:25Z

inst/stan/functions/infections.stan

    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 == 1 || (use_pop == 2 && s <= nht)) {


Sorry I'm not still getting the logic. We have:
use_pop == 1: only adjust in forecast
use_pop == 2: adjust across all

So shouldn't this be

Suggested change

if (use_pop == 1 || (use_pop == 2 && s <= nht)) {

if ((use_pop == 1 && s > nht) || use_pop == 2)) {

sbfnk · 2025-01-13T12:09:59Z

tests/testthat/test-rt_opts.R

+  expect_warning(
+    rt_opts(pop = 1000),
+    "The `pop` argument of `rt_opts()` must be a `<dist_spec>` as of EpiNow2 1.7.0."
+  )


Suggested change

expect_warning(

rt_opts(pop = 1000),

"The `pop` argument of `rt_opts()` must be a `<dist_spec>` as of EpiNow2 1.7.0."

)

lifecycle::expect_deprecated(

rt_opts(pop = 1000),

"The `pop` argument of `rt_opts\\(\\)` must be a `<dist_spec>` as of EpiNow2 1.7.0."

)

sbfnk · 2025-01-13T12:14:42Z

inst/stan/functions/infections.stan

    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 == 1 || (use_pop == 2 && s <= nht)) {
      exp_adj_Rt = exp(-R[s] * infectiousness[s] / (pop - cum_infections[nht]));


Suggested change

exp_adj_Rt = exp(-R[s] * infectiousness[s] / (pop - cum_infections[nht]));

exp_adj_Rt = exp(-R[s] * infectiousness[s] / (pop - cum_infections[s]));

seabbs requested a review from sbfnk December 19, 2024 18:23

seabbs enabled auto-merge December 20, 2024 16:49

sbfnk requested changes Dec 20, 2024

View reviewed changes

R/opts.R Outdated Show resolved Hide resolved

inst/stan/functions/infections.stan Outdated Show resolved Hide resolved

This comment was marked as outdated.

Sign in to view

sbfnk reviewed Jan 7, 2025

View reviewed changes

R/opts.R Outdated Show resolved Hide resolved

inst/stan/functions/infections.stan Outdated Show resolved Hide resolved

seabbs and others added 9 commits January 10, 2025 10:41

change over opts and create interface

f372ff5

add support to stan

3e30b64

check works and update tests

288b424

update news

6b1cb51

update tests

c17fdcb

handle pop more carefully in a simulation setting

2967e05

update stan generate infection tests

86018b1

fix linting

07064bd

make it possible to either use pop out of sample or fit with it in sa…

3bc26d4

…mple

seabbs force-pushed the issue408 branch from a08f8f4 to 3bc26d4 Compare January 10, 2025 11:22

seabbs added 2 commits January 10, 2025 11:44

add support for in sample suscetpble

01f42c9

update docs

e21a5f8

This comment was marked as outdated.

Sign in to view

Document

101f1c5

seabbs mentioned this pull request Jan 10, 2025

Reporting adjusted or unadjusted Rt #922

Open

seabbs added 2 commits January 10, 2025 13:23

touch

27097cb

cleaner opt passing

bafef0b

seabbs and others added 3 commits January 10, 2025 15:46

update tests

841d238

Merge branch 'main' into issue408

1bedd6c

clean up whitespace

6184ab9

seabbs requested a review from sbfnk January 10, 2025 16:25

sbfnk requested changes Jan 13, 2025

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Issue #408: Fit the susceptible population #904

Issue #408: Fit the susceptible population #904

seabbs commented Dec 19, 2024 •

edited

Loading

seabbs commented Dec 20, 2024

sbfnk left a comment

sbfnk commented Dec 20, 2024

This comment was marked as outdated.

This comment was marked as outdated.

github-actions bot commented Jan 10, 2025

github-actions bot commented Jan 10, 2025

sbfnk Jan 13, 2025

sbfnk Jan 13, 2025

sbfnk Jan 13, 2025

sbfnk Jan 13, 2025

sbfnk Jan 13, 2025

sbfnk Jan 13, 2025

sbfnk Jan 13, 2025

		#' Note that with "forecast", Rt estimates are unadjusted for susceptible
		#' depletion but posterior predictions are adjusted.

-#' Note that with "forecast", Rt estimates are unadjusted for susceptible
-#' depletion but posterior predictions are adjusted.
+#' Note that with "all", Rt estimates are unadjusted for susceptible
+#' depletion but posterior predictions of infections and reports are
+#' adjusted.

	int use_pop; // use population size
	int use_pop; // use population size (0 = no; 1 = forecasts; 2 = all)

	if (use_pop == 1 \|\| (use_pop == 2 && s <= nht)) {
	if ((use_pop == 1 && s > nht) \|\| use_pop == 2)) {

	exp_adj_Rt = exp(-R[s] * infectiousness[s] / (pop - cum_infections[nht]));
	exp_adj_Rt = exp(-R[s] * infectiousness[s] / (pop - cum_infections[s]));

Issue #408: Fit the susceptible population #904

Are you sure you want to change the base?

Issue #408: Fit the susceptible population #904

Conversation

seabbs commented Dec 19, 2024 • edited Loading

Description

Initial submission checklist

After the initial Pull Request

seabbs commented Dec 20, 2024

sbfnk left a comment

Choose a reason for hiding this comment

sbfnk commented Dec 20, 2024

This comment was marked as outdated.

This comment was marked as outdated.

github-actions bot commented Jan 10, 2025

github-actions bot commented Jan 10, 2025

sbfnk Jan 13, 2025

Choose a reason for hiding this comment

sbfnk Jan 13, 2025

Choose a reason for hiding this comment

sbfnk Jan 13, 2025

Choose a reason for hiding this comment

sbfnk Jan 13, 2025

Choose a reason for hiding this comment

sbfnk Jan 13, 2025

Choose a reason for hiding this comment

sbfnk Jan 13, 2025

Choose a reason for hiding this comment

sbfnk Jan 13, 2025

Choose a reason for hiding this comment

seabbs commented Dec 19, 2024 •

edited

Loading