don't estimate initial growth

given the very strong correlation with the number of initial infections, this reduces the estimation to one parameter and works out initial growth from that
epiforecasts · Jan 8, 2025 · 6225aa5 · 6225aa5
1 parent 746c3ba
commit 6225aa5
Show file tree

Hide file tree

Showing 2 changed files with 9 additions and 9 deletions.
diff --git a/inst/stan/estimate_infections.stan b/inst/stan/estimate_infections.stan
@@ -47,7 +47,6 @@ parameters {
   vector[fixed ? 0 : gp_type == 1 ? 2*M : M] eta;               // unconstrained noise
   // Rt
   array[estimate_r] real initial_infections;    // seed infections
-  array[estimate_r && seeding_time > 1 ? 1 : 0] real initial_growth; // seed growth rate
   array[bp_n > 0 ? 1 : 0] real<lower = 0> bp_sd; // standard deviation of breakpoint effect
   vector[bp_n] bp_effects;                   // Rt breakpoint effects
   // observation model
@@ -62,6 +61,12 @@ transformed parameters {
   vector[ot_h] reports;                                     // estimated reported cases
   vector[ot] obs_reports;                                   // observed estimated reported cases
   vector[estimate_r * (delay_type_max[gt_id] + 1)] gt_rev_pmf;
+  array[estimate_r && seeding_time > 1 ? 1 : 0] real initial_growth; // seed growth rate
+
+  if (num_elements(initial_growth) > 0) {
+    initial_growth[1] = prior_growth +
+      (prior_infections - initial_infections[1]) / seeding_time;
+  }
 
   // GP in noise - spectral densities
   profile("update gp") {
@@ -205,8 +210,7 @@ model {
     // priors on Rt
     profile("rt lp") {
       rt_lp(
-        initial_infections, initial_growth, bp_effects, bp_sd, bp_n,
-        seeding_time, prior_infections, prior_growth
+        initial_infections, bp_effects, bp_sd, bp_n, prior_infections
       );
     }
   }

diff --git a/inst/stan/functions/rt.stan b/inst/stan/functions/rt.stan
@@ -57,9 +57,8 @@ vector update_Rt(int t, real R0, vector noise, array[] int bps,
  * @param prior_infections Prior mean for initial infections
  * @param prior_growth Prior mean for initial growth rates
  */
-void rt_lp(array[] real initial_infections, array[] real initial_growth,
-           vector bp_effects, array[] real bp_sd, int bp_n, int seeding_time,
-           real prior_infections, real prior_growth) {
+void rt_lp(array[] real initial_infections, vector bp_effects,
+           array[] real bp_sd, int bp_n, real prior_infections) {
   //breakpoint effects on Rt
   if (bp_n > 0) {
     bp_sd[1] ~ normal(0, 0.1) T[0,];
@@ -68,7 +67,4 @@ void rt_lp(array[] real initial_infections, array[] real initial_growth,
   // initial infections
   initial_infections ~ normal(prior_infections, sqrt(prior_infections));
 
-  if (seeding_time > 1) {
-    initial_growth ~ normal(prior_growth, 0.2);
-  }
 }