Merge pull request #97 from mrc-ide/particle

0.4.8 - calibrate reporting fraction passed through to likelihood

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: squire
 Type: Package
 Title: SEIR transmission model of COVID-19
-Version: 0.4.7
+Version: 0.4.8
 Authors@R: c(
   person("Patrick", "Walker", email = "[email protected]", role = c("aut")),
   person("OJ", "Watson", email = "[email protected]", role = c("aut", "cre")),

NEWS.md

@@ -1,4 +1,8 @@
 
+# squire 0.4.8
+
+* `calibrate` argument `reporting_fraction` implemented now in likelihood 
+
 # squire 0.4.7
 
 * `calibrate` can do 3D grid scans investigating Meff

R/calibration.R

@@ -213,9 +213,6 @@ calibrate <- function(data,
   # make the date definitely a date
   data$date <- as.Date(as.character(data$date))
 
-  # adjust for reporting fraction
-  data$deaths <- (data$deaths/reporting_fraction)
-
   # build model parameters
   model_params <- squire_model$parameter_func(country = country,
                                               population = population,
@@ -237,6 +234,13 @@ calibrate <- function(data,
                         date_hosp_bed_capacity_change = date_hosp_bed_capacity_change,
                         hosp_bed_capacity = hosp_bed_capacity)
 
+  # construct pars_obs for the user
+  pars_obs <-  list(phi_cases = reporting_fraction,
+                    k_cases = 2,
+                    phi_death = reporting_fraction,
+                    k_death = 2,
+                    exp_noise = 1e6)
+
   # construct scan
   if (Meff_min == Meff_max) {
     scan_results <- scan_R0_date(R0_min = R0_min,
@@ -253,6 +257,7 @@ calibrate <- function(data,
                                  date_ICU_bed_capacity_change = date_ICU_bed_capacity_change,
                                  date_hosp_bed_capacity_change = date_hosp_bed_capacity_change,
                                  squire_model = squire_model,
+                                 pars_obs = pars_obs,
                                  n_particles = n_particles)
   } else {
     scan_results <- scan_R0_date_Meff(R0_min = R0_min,
@@ -272,6 +277,7 @@ calibrate <- function(data,
                                       date_ICU_bed_capacity_change = date_ICU_bed_capacity_change,
                                       date_hosp_bed_capacity_change = date_hosp_bed_capacity_change,
                                       squire_model = squire_model,
+                                      pars_obs = pars_obs,
                                       n_particles = n_particles)
   }
 

R/estimation.R

@@ -85,7 +85,7 @@ scan_R0_date <- function(
 
   # Set up observation parameters that translate our model outputs to observations
   if (is.null(pars_obs)) {
-    pars_obs <-  list(phi_cases = 0.1,
+    pars_obs <-  list(phi_cases = 1,
                       k_cases = 2,
                       phi_death = 1,
                       k_death = 2,
@@ -282,7 +282,7 @@ scan_R0_date_Meff <- function(
 
   # Set up observation parameters that translate our model outputs to observations
   if (is.null(pars_obs)) {
-    pars_obs <-  list(phi_cases = 0.1,
+    pars_obs <-  list(phi_cases = 1,
                       k_cases = 2,
                       phi_death = 1,
                       k_death = 2,
@@ -905,8 +905,7 @@ plot_sample_grid_search <- function(x, what = "deaths") {
     base_plot <- base_plot +
       ggplot2::geom_line(ggplot2::aes(y=.data$ymin, x=as.Date(.data$date)), quants, linetype="dashed") +
       ggplot2::geom_line(ggplot2::aes(y=.data$ymax, x=as.Date(.data$date)), quants, linetype="dashed") +
-      ggplot2::geom_point(ggplot2::aes(y=.data$cases/x$scan_results$inputs$pars_obs$phi_cases,
-                                       x=as.Date(.data$date)), x$scan_results$inputs$data)
+      ggplot2::geom_point(ggplot2::aes(y=.data$cases, x=as.Date(.data$date)), x$scan_results$inputs$data)
 
   }
 
@@ -929,8 +928,7 @@ plot_sample_grid_search <- function(x, what = "deaths") {
     base_plot <- base_plot +
       ggplot2::geom_line(ggplot2::aes(y=.data$ymin, x=as.Date(.data$date)), quants, linetype="dashed") +
       ggplot2::geom_line(ggplot2::aes(y=.data$ymax, x=as.Date(.data$date)), quants, linetype="dashed") +
-      ggplot2::geom_point(ggplot2::aes(y=.data$deaths/x$scan_results$inputs$pars_obs$phi_death,
-                                       x=as.Date(.data$date)), x$scan_results$inputs$data)
+      ggplot2::geom_point(ggplot2::aes(y=.data$deaths,x=as.Date(.data$date)), x$scan_results$inputs$data)
 
   } else {
 

tests/testthat/test-calibration.R

@@ -737,3 +737,69 @@ test_that("R0 and Meff arge checking", {
 
 
 })
+
+
+#------------------------------------------------
+test_that("reporting fraction into pars_obs", {
+
+  data <- read.csv(squire_file("extdata/example.csv"),stringsAsFactors = FALSE)
+  interventions <- read.csv(squire_file("extdata/example_intervention.csv"))
+  int_unique <- interventions_unique(interventions)
+  country = "Algeria"
+  replicates = 2
+  R0_min = 2.6
+  R0_max = 2.6
+  R0_step = 0.1
+  first_start_date = "2020-02-01"
+  last_start_date = "2020-02-02"
+  day_step = 1
+  R0_change = int_unique$change
+  date_R0_change = as.Date(int_unique$dates_change)
+  date_contact_matrix_set_change = NULL
+  squire_model = explicit_model()
+  pars_obs = NULL
+  n_particles = 5
+
+  set.seed(93L)
+  out <- calibrate(
+    data = data,
+    R0_min = R0_min,
+    R0_max = R0_max,
+    R0_step = R0_step,
+    first_start_date = first_start_date,
+    last_start_date = last_start_date,
+    day_step = day_step,
+    squire_model = squire_model,
+    pars_obs = pars_obs,
+    n_particles = n_particles,
+    reporting_fraction = 0.1,
+    R0_change = R0_change,
+    date_R0_change = date_R0_change,
+    replicates = replicates,
+    country = country,
+    forecast = 0
+  )
+
+  out2 <- calibrate(
+    data = data,
+    R0_min = R0_min,
+    R0_max = R0_max,
+    R0_step = R0_step,
+    first_start_date = first_start_date,
+    last_start_date = last_start_date,
+    day_step = day_step,
+    squire_model = squire_model,
+    pars_obs = pars_obs,
+    n_particles = n_particles,
+    reporting_fraction = 1,
+    R0_change = R0_change,
+    date_R0_change = date_R0_change,
+    replicates = replicates,
+    country = country,
+    forecast = 0
+  )
+
+  index <- odin_index(out$model)
+  expect_true(sum(rowSums(out$output[,index$D,1]))  > sum(rowSums(out2$output[,index$D,1])))
+
+})