Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update growth calculation #213

Closed
wants to merge 4 commits into from
Closed

Update growth calculation #213

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5038f42
switch to directly calculated growth rate
seabbs Dec 30, 2020
3b45b9c
tests for built package
seabbs Dec 30, 2020
f7fa03b
update to just estimate using within observation infections
seabbs Dec 30, 2020
5aef61c
update stan
seabbs Dec 30, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion R/extract.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

)
if (data$week_effect == 1) {
out$day_of_week <- extract_parameter(
Expand Down
25 changes: 10 additions & 15 deletions R/stanmodels.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,22 +17,17 @@ Rcpp::loadModule("stan_fit4tune_inv_gamma_mod", what = TRUE)
# instantiate each stanmodel object
stanmodels <- sapply(stanmodels, function(model_name) {
# create C++ code for stan model
stan_file <- if (dir.exists("stan")) "stan" else file.path("inst", "stan")
stan_file <- if(dir.exists("stan")) "stan" else file.path("inst", "stan")
stan_file <- file.path(stan_file, paste0(model_name, ".stan"))
stanfit <- rstan::stanc_builder(stan_file,
allow_undefined = TRUE,
obfuscate_model_name = FALSE
)
stanfit$model_cpp <- list(
model_cppname = stanfit$model_name,
model_cppcode = stanfit$cppcode
)
allow_undefined = TRUE,
obfuscate_model_name = FALSE)
stanfit$model_cpp <- list(model_cppname = stanfit$model_name,
model_cppcode = stanfit$cppcode)
# create stanmodel object
methods::new(
Class = "stanmodel",
model_name = stanfit$model_name,
model_code = stanfit$model_code,
model_cpp = stanfit$model_cpp,
mk_cppmodule = function(x) get(paste0("model_", model_name))
)
methods::new(Class = "stanmodel",
model_name = stanfit$model_name,
model_code = stanfit$model_code,
model_cpp = stanfit$model_cpp,
mk_cppmodule = function(x) get(paste0("model_", model_name)))
})
11 changes: 4 additions & 7 deletions inst/stan/estimate_infections.stan
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

vector[ot] log_lik;
if (estimate_r){
// estimate growth from estimated Rt
r = R_to_growth(R, gt_mean[1], gt_sd[1]);
}else{
if (estimate_r == 0){
// 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);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

}
// estimate growth from infections
r = calculate_growth(infections, seeding_time + 1);
// simulate reported cases
imputed_reports = report_rng(reports, rep_phi, model_type);
// log likelihood of model
Expand Down
8 changes: 8 additions & 0 deletions inst/stan/functions/generated_quantities.stan
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,3 +44,11 @@ real[] R_to_growth(vector R, real gt_mean, real gt_sd) {
}
return(r);
}
// Calculate growth rate
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));
}
Comment on lines +48 to +54
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Update, if needed, according to the above comments.

4 changes: 2 additions & 2 deletions inst/stan/simulate_infections.stan
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

for (i in 1:n) {
// generate infections from Rt trace
infections[i] = to_row_vector(generate_infections(to_vector(R[i]), seeding_time,
Expand All @@ -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);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could use a switch as discussed above.

}
}
727 changes: 388 additions & 339 deletions src/stanExports_estimate_infections.h
View file
Open in desktop

Large diffs are not rendered by default.

144 changes: 72 additions & 72 deletions src/stanExports_estimate_secondary.h
View file
Open in desktop

Large diffs are not rendered by default.

102 changes: 51 additions & 51 deletions src/stanExports_estimate_truncation.h
View file
Open in desktop

Large diffs are not rendered by default.

454 changes: 253 additions & 201 deletions src/stanExports_simulate_infections.h
View file
Open in desktop

Large diffs are not rendered by default.

150 changes: 75 additions & 75 deletions src/stanExports_simulate_secondary.h
View file
Open in desktop

Large diffs are not rendered by default.

12 changes: 6 additions & 6 deletions src/stanExports_tune_inv_gamma.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,14 +77,14 @@ tail_delta(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& y,
current_statement_begin__ = 12;
stan::math::assign(beta, (logical_gte(beta, 1e6) ? stan::math::promote_scalar<local_scalar_t__>(1e6) : stan::math::promote_scalar<local_scalar_t__>(beta) ));
current_statement_begin__ = 13;
stan::model::assign(deltas,
stan::model::cons_list(stan::model::index_uni(1), stan::model::nil_index_list()),
(inv_gamma_cdf(get_base1(theta, 1, "theta", 1), alpha, beta) - 0.01),
stan::model::assign(deltas,
stan::model::cons_list(stan::model::index_uni(1), stan::model::nil_index_list()),
(inv_gamma_cdf(get_base1(theta, 1, "theta", 1), alpha, beta) - 0.01),
"assigning variable deltas");
current_statement_begin__ = 14;
stan::model::assign(deltas,
stan::model::cons_list(stan::model::index_uni(2), stan::model::nil_index_list()),
((1 - inv_gamma_cdf(get_base1(theta, 2, "theta", 1), alpha, beta)) - 0.01),
stan::model::assign(deltas,
stan::model::cons_list(stan::model::index_uni(2), stan::model::nil_index_list()),
((1 - inv_gamma_cdf(get_base1(theta, 2, "theta", 1), alpha, beta)) - 0.01),
"assigning variable deltas");
current_statement_begin__ = 15;
return stan::math::promote_scalar<fun_return_scalar_t__>(deltas);
Expand Down
16 changes: 16 additions & 0 deletions tests/testthat/test-stan-generated_quantities.R
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
context("estimate_infections")
if (!testthat:::on_cran()) {
files <- c("pmfs.stan", "infections.stan", "generated_quantities.stan")
suppressMessages(expose_stan_fns(files, target_dir = system.file("stan/functions", package = "EpiNow2")))
}


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))
})
Comment on lines +8 to +15
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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