don't defer negative binomial to Poisson at low levels of overdispers…

…ion (#432)

DESCRIPTION

@@ -2,7 +2,7 @@ Type: Package
 Package: EpiNow2
 Title: Estimate Real-Time Case Counts and Time-Varying
     Epidemiological Parameters
-Version: 1.3.6.9000
+Version: 1.3.6.9001
 Authors@R: 
     c(person(given = "Sam",
              family = "Abbott",

NEWS.md

@@ -22,6 +22,7 @@ This release is in development. For a stable release install 1.3.5 from CRAN.
 * The internal distribution interface has been streamlined to reduce code duplication. By @sbfnk in #363 and reviewed by @seabbs.
 * A small bug has been fixed where the seeding time was too long. When a single delay is used this shortens the seeding time by one day and when more delays are used it shortens the seeding time by n days where n is the number of delays used e.g. for two parametric delays it's two days. By @sbfnk in #413 and reviewed by @seabbs.
 * Some tuning was done to speed up the renewal model. By @sbfnk in #416 and reviewed by @seabbs.
+* An approximation of the negative binomial by the Poisson at low levels of overdispersion was disabled as it led to parameter identification issues. By @sbfnk in #432 and reviewed by @seabbs.
 
 # EpiNow2 1.3.5
 

inst/stan/functions/observation_model.stan

@@ -54,50 +54,40 @@ void truncation_lp(real[] truncation_mean, real[] truncation_sd,
 void report_lp(int[] cases, vector reports,
                real[] rep_phi, real phi_mean, real phi_sd,
                int model_type, real weight) {
-  real sqrt_phi = 1e5;
   if (model_type) {
-    // the reciprocal overdispersion parameter (phi)
+    real sqrt_phi; // the reciprocal overdispersion parameter (phi)
     rep_phi[model_type] ~ normal(phi_mean, phi_sd) T[0,];
     sqrt_phi = 1 / sqrt(rep_phi[model_type]);
-    // defer to poisson if phi is large, to avoid overflow or
-    // if poisson specified
-  }
-  if (sqrt_phi > 1e4) {
     if (weight == 1) {
-      cases ~ poisson(reports);
-    }else{
-      target += poisson_lpmf(cases | reports) * weight;
+      cases ~ neg_binomial_2(reports, sqrt_phi);
+    } else {
+      target += neg_binomial_2_lpmf(cases | reports, sqrt_phi) * weight;
     }
   } else {
     if (weight == 1) {
-      cases ~ neg_binomial_2(reports, sqrt_phi);
-    }else{
-      target += neg_binomial_2_lpmf(cases | reports, sqrt_phi);
+      cases ~ poisson(reports);
+    } else {
+      target += poisson_lpmf(cases | reports) * weight;
     }
   }
-
 }
 // update log likelihood (as above but not vectorised and returning log likelihood)
 vector report_log_lik(int[] cases, vector reports,
                       real[] rep_phi, int model_type, real weight) {
   int t = num_elements(reports);
   vector[t] log_lik;
-  real sqrt_phi = 1e5;
-  if (model_type) {
-  // the reciprocal overdispersion parameter (phi)
-  sqrt_phi = 1 / sqrt(rep_phi[model_type]);
-  }
 
   // defer to poisson if phi is large, to avoid overflow
-  if (sqrt_phi > 1e4) {
+  if (model_type == 0) {
     for (i in 1:t) {
       log_lik[i] = poisson_lpmf(cases[i] | reports[i]) * weight;
     }
   } else {
+    real sqrt_phi = 1 / sqrt(rep_phi[model_type]);
     for (i in 1:t) {
       log_lik[i] = neg_binomial_2_lpmf(cases[i] | reports[i], sqrt_phi) * weight;
     }
-    }
+  }
   return(log_lik);
 }
 // sample reported cases from the observation model

tests/testthat/test-estimate_secondary.R

@@ -23,7 +23,6 @@ cases[
 # with a secondary case
 inc <- estimate_secondary(cases[1:60],
   obs = obs_opts(scale = list(mean = 0.2, sd = 0.2), week_effect = FALSE),
-  control = list(adapt_delta = 0.98),
   verbose = FALSE
 )
 
@@ -86,14 +85,12 @@ test_that("estimate_secondary can recover simulated parameters", {
     inc_posterior[, median], c(1.8, 0.5, 0.4),
     tolerance = 0.1
   )
-  # These tests currently fail indicating the model is not recovering the
-  # simulated parameters.
-  # expect_equal(
-  #   prev_posterior[, mean], c(1.6, 0.8, 0.3), tolerance = 0.1
-  # )
-  # expect_equal(
-  #   prev_posterior[, median], c(1.6, 0.8, 0.3), tolerance = 0.1
-  # )
+  expect_equal(
+    prev_posterior[, mean], c(1.6, 0.8, 0.3), tolerance = 0.2
+  )
+  expect_equal(
+    prev_posterior[, median], c(1.6, 0.8, 0.3), tolerance = 0.2
+  )
 })
 
 test_that("forecast_secondary can return values from simulated data and plot