Update growth calculation #213

seabbs · 2020-12-30T20:10:29Z

This updates the output growth rate to be directly calculated from infections (r = log(inf_(t-1)) - log(inf_t)) rather than via the growth rate approximation using Rt.

Adds calculate_growth and tests
Reduces growth rates output to be number of observations minus 1 (as can't be calculated for the first time step).

Everything looks as I would expect when unit testing and running the examples expect for this slightly odd artefact in the forecast portion that doesn't show up in the Rt or infections plot.

sbfnk · 2020-12-31T00:38:20Z

Makes sense. Is this with forecasts keeping R constant? I could possibly see a non-smooth change in R translating into a discontinuity in growth.

seabbs · 2020-12-31T15:05:10Z

Yes it is. Hmm.

hsbadr · 2020-12-31T22:39:03Z

R/extract.R

@@ -119,7 +119,7 @@ extract_parameter_samples <- function(stan_fit, data, reported_dates, reported_i
  out$growth_rate <- extract_parameter(
    "r",
    samples,
-    reported_dates
+    reported_dates[-1]


Instead of dropping the first time step, update calculate_growth to either get an approximate value from Rt or set it to NA .

hsbadr · 2020-12-31T22:39:51Z

inst/stan/estimate_infections.stan

@@ -118,21 +118,18 @@ model {
 generated quantities {
  int imputed_reports[ot_h];
  vector[estimate_r > 0 ? 0: ot_h] gen_R;
-  real r[ot_h];
+  real r[ot_h - 1];


This is related to approximating the growth rate for the 1st time step.

hsbadr · 2020-12-31T22:42:34Z

inst/stan/estimate_infections.stan

    // sample generation time
    real gt_mean_sample = normal_rng(gt_mean_mean, gt_mean_sd);
    real gt_sd_sample = normal_rng(gt_sd_mean, gt_sd_sd);
    // calculate Rt using infections and generation time
    gen_R = calculate_Rt(infections, seeding_time, gt_mean_sample, gt_mean_sample,
                         max_gt, rt_half_window);
-    // estimate growth from calculated Rt
-    r = R_to_growth(gen_R, gt_mean_sample, gt_sd_sample);


Now that calculate_growth and R_to_growth would have the same return value and vector length, you might keep both with a switch, if there's a value in keeping the Rt growth rate approximation (e.g., for diagnostics, not necessarily as a user option).

hsbadr · 2020-12-31T22:44:00Z

inst/stan/functions/generated_quantities.stan

+real[] calculate_growth(vector infections, int seeding_time) {
+  int t = num_elements(infections);
+  int ot = t - seeding_time;
+  vector[t] log_inf = log(infections); 
+  vector[ot] growth = log_inf[(seeding_time + 1):t] - log_inf[seeding_time:(t - 1)];
+  return(to_array_1d(growth));
+}


Update, if needed, according to the above comments.

hsbadr · 2020-12-31T22:44:23Z

inst/stan/simulate_infections.stan

@@ -27,7 +27,7 @@ generated quantities {
  matrix[n, t] infections; //latent infections
  matrix[n, t - seeding_time] reports; // observed cases
  int imputed_reports[n, t - seeding_time];
-  real r[n, t - seeding_time];
+  real r[n, t - seeding_time - 1];


This is related to approximating the growth rate for the 1st time step.

hsbadr · 2020-12-31T22:44:47Z

inst/stan/simulate_infections.stan

@@ -48,6 +48,6 @@ generated quantities {
    }
   // simulate reported cases
   imputed_reports[i] = report_rng(to_vector(reports[i]), rep_phi[i], model_type);
-   r[i] = R_to_growth(to_vector(R[i]), gt_mean[i, 1], gt_sd[i, 1]);
+   r[i] = calculate_growth(to_vector(infections[i]), seeding_time + 1);


Could use a switch as discussed above.

hsbadr · 2020-12-31T22:47:50Z

tests/testthat/test-stan-generated_quantities.R

+test_that("calculate_growth works as expected", {
+  skip_on_cran()
+  expect_equal(calculate_growth(rep(1, 5), 1), rep(0, 4))
+  expect_equal(round(calculate_growth(1:5, 2), 2), c(0.41, 0.29, 0.22))
+  expect_equal(round(calculate_growth(exp(0.4*1:5), 2), 2), rep(0.4, 3))
+  expect_error(calculate_growth(1:5, 6))
+  expect_error(calculate_growth(1:5, 0))
+})


Could also add tests to compare the outputs from calculate_growth and R_to_growth, and to check if there's discontinuities in the forecasts (maybe by comparing the difference between the first forecase time step and the previous day within a reasonable threshold).

See #213 (comment) #213 (comment) #213 (comment) For now not implementing comparison to the approximate growth rate as this seems quite a specific use case that could also be done outside the stan model. Also not implementing any approximate growth rate from seeding time - instead minimum seeding time is now set to 1, so the last seeding time is used to calculate the first growth rate.

sbfnk · 2022-07-27T16:10:30Z

Closing in favour of #307

See #213 (comment) #213 (comment) #213 (comment) For now not implementing comparison to the approximate growth rate as this seems quite a specific use case that could also be done outside the stan model. Also not implementing any approximate growth rate from seeding time - instead minimum seeding time is now set to 1, so the last seeding time is used to calculate the first growth rate.

@hsbadr

implements suggestions by @hsbadr See #213 (comment) #213 (comment) #213 (comment) For now not implementing comparison to the approximate growth rate as this seems quite a specific use case that could also be done outside the stan model. Also not implementing any approximate growth rate from seeding time - instead minimum seeding time is now set to 1, so the last seeding time is used to calculate the first growth rate.

@hsbadr

implements suggestions by @hsbadr See #213 (comment) #213 (comment) #213 (comment) For now not implementing comparison to the approximate growth rate as this seems quite a specific use case that could also be done outside the stan model. Also not implementing any approximate growth rate from seeding time - instead minimum seeding time is now set to 1, so the last seeding time is used to calculate the first growth rate.

seabbs added 4 commits December 30, 2020 18:43

switch to directly calculated growth rate

5038f42

tests for built package

3b45b9c

update to just estimate using within observation infections

f7fa03b

update stan

5aef61c

seabbs requested a review from sbfnk December 30, 2020 20:10

hsbadr reviewed Dec 31, 2020

View reviewed changes

sbfnk mentioned this pull request Jul 27, 2022

switch to directly calculated growth rate and extend model options #307

Closed

3 tasks

sbfnk closed this Jul 27, 2022

sbfnk mentioned this pull request Jan 18, 2023

switch to directly calculated growth rate and extend model options, take 2 #345

Closed

11 tasks

sbfnk deleted the update_growth branch March 12, 2024 16:29

sbfnk mentioned this pull request Mar 12, 2024

Directly calculated growth #610

Merged

7 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Update growth calculation #213

Update growth calculation #213

seabbs commented Dec 30, 2020 •

edited

Loading

sbfnk commented Dec 31, 2020

seabbs commented Dec 31, 2020

hsbadr Dec 31, 2020

hsbadr Dec 31, 2020

hsbadr Dec 31, 2020

hsbadr Dec 31, 2020

hsbadr Dec 31, 2020

hsbadr Dec 31, 2020

hsbadr Dec 31, 2020

sbfnk commented Jul 27, 2022

Update growth calculation #213

Update growth calculation #213

Conversation

seabbs commented Dec 30, 2020 • edited Loading

sbfnk commented Dec 31, 2020

seabbs commented Dec 31, 2020

hsbadr Dec 31, 2020

Choose a reason for hiding this comment

hsbadr Dec 31, 2020

Choose a reason for hiding this comment

hsbadr Dec 31, 2020

Choose a reason for hiding this comment

hsbadr Dec 31, 2020

Choose a reason for hiding this comment

hsbadr Dec 31, 2020

Choose a reason for hiding this comment

hsbadr Dec 31, 2020

Choose a reason for hiding this comment

hsbadr Dec 31, 2020

Choose a reason for hiding this comment

sbfnk commented Jul 27, 2022

seabbs commented Dec 30, 2020 •

edited

Loading