Merge pull request #168 from mrc-ide/sero_fitting

0.6.7 sero fitting in

Author: OJWatson

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: squire
 Type: Package
 Title: SEIR transmission model of COVID-19
-Version: 0.6.6
+Version: 0.6.7
 Authors@R: c(
   person("OJ", "Watson", email = "o.watson15@imperial.ac.uk", role = c("aut", "cre")),
   person("Patrick", "Walker", email = "patrick.walker06@imperial.ac.uk", role = c("aut")),

NAMESPACE

@@ -56,6 +56,7 @@ importFrom(odin,odin)
 importFrom(rlang,.data)
 importFrom(stats,cor)
 importFrom(stats,cov)
+importFrom(stats,dbinom)
 importFrom(stats,dnbinom)
 importFrom(stats,median)
 importFrom(stats,plogis)

NEWS.md

@@ -1,3 +1,8 @@
+# squire 0.6.7
+
+* `pmcmc` can now be used to fit to serology data (deterministic model only) 
+by passing `sero_df` and `sero_det` as names list elements of `pars_obs`
+
 # squire 0.6.6
 
 * `projections` can be used now for `nimue` models

R/particle.R

@@ -650,6 +650,7 @@ scale_log_weights <- function(log_weights) {
 #'
 #' @return Results from particle filter
 #'
+#' @importFrom stats dbinom
 run_deterministic_comparison <- function(data,
                                          squire_model,
                                          model_params,
@@ -700,7 +701,7 @@ run_deterministic_comparison <- function(data,
   Ds[Ds < 0] <- 0
   deaths <- data$deaths[-1]
 
-  # calculate ll
+  # calculate ll for deaths
   if (obs_params$treated_deaths_only) {
 
     Ds_heathcare <- diff(rowSums(out[,index$D_get]))
@@ -713,14 +714,57 @@ run_deterministic_comparison <- function(data,
 
   }
 
+  # calculate ll for the seroprevalence
+  lls <- 0
+  if("sero_df" %in% obs_params && "sero_det" %in% obs_params) {
+
+    sero_df <- obs_params$sero_df
+    sero_det <- obs_params$sero_det
+
+    # were there actually seroprevalence data points to compare against
+    if(nrow(sero_df) > 0) {
+
+      sero_at_date <- function(date, symptoms, det, dates, N) {
+
+        di <- which(dates == date)
+        if(length(di) > 0) {
+          to_sum <- tail(symptoms[seq_len(di)], length(det))
+          min(sum(rev(to_sum)*head(det, length(to_sum)), na.rm=TRUE)/N, 0.99)
+        } else {
+          0
+        }
+
+      }
+
+    # get symptom incidence
+    symptoms <- rowSums(out[,index$E2]) * model_params$gamma_E
+
+    # dates of incidence, pop size and dates of sero surveys
+    dates <- data$date[[1]] + seq_len(nrow(out)) - 1L
+    N <- sum(model_params$population)
+    sero_dates <- list(sero_df$date_end, sero_df$date_start, sero_df$date_start + as.integer((sero_df$date_end - sero_df$date_start)/2))
+    unq_sero_dates <- unique(c(sero_df$date_end, sero_df$date_start, sero_df$date_start + as.integer((sero_df$date_end - sero_df$date_start)/2)))
+    det <- obs_params$sero_det
+
+    # estimate model seroprev
+    sero_model <- vapply(unq_sero_dates, sero_at_date, numeric(1), symptoms, det, dates, N)
+    sero_model_mat <- do.call(cbind,lapply(sero_dates, function(x) {sero_model[match(x, unq_sero_dates)]}))
+
+    # likelihood of model obvs
+    lls <- rowMeans(dbinom(sero_df$sero_pos, sero_df$samples, sero_model_mat, log = TRUE))
+
+    }
+
+  }
+
   # format the out object
   date <- data$date[[1]] + seq_len(nrow(out)) - 1L
   rownames(out) <- as.character(date)
   attr(out, "date") <- date
 
   # format similar to particle_filter nomenclature
   pf_results <- list()
-  pf_results$log_likelihood <- sum(ll)
+  pf_results$log_likelihood <- sum(ll) + sum(lls)
 
   # single returns final state
   if (save_history) {

R/pmcmc.R

@@ -293,8 +293,8 @@ pmcmc <- function(data,
   }
   assert_logical(unlist(pars_discrete))
   assert_list(pars_obs)
-  assert_eq(names(pars_obs), c("phi_cases", "k_cases", "phi_death", "k_death", "exp_noise"))
-  assert_numeric(unlist(pars_obs))
+  assert_in(c("phi_cases", "k_cases", "phi_death", "k_death", "exp_noise"), names(pars_obs))
+  assert_numeric(unlist(pars_obs[c("phi_cases", "k_cases", "phi_death", "k_death", "exp_noise")]))
 
   # mcmc items
   assert_pos_int(n_mcmc)
@@ -1315,6 +1315,12 @@ calc_loglikelihood <- function(pars, data, squire_model, model_params,
   R0 <- pars[["R0"]]
   start_date <- pars[["start_date"]]
 
+  # reporting fraction par if in pars list
+  if("rf" %in% names(pars)) {
+    assert_numeric(pars[["rf"]])
+    pars_obs$phi_death <- pars[["rf"]]
+  }
+
   #----------------..
   # more assertions
   #----------------..

R/pmcmc_object.R

@@ -49,10 +49,17 @@ plot_pmcmc_sample  <- function(x, what = "deaths") {
 
     base_plot <- plot(x, "deaths", ci = FALSE, replicates = TRUE, x_var = "date",
                       date_0 = max(x$pmcmc_results$inputs$data$date))
+
+    if ("rf" %in% names(x$replicate_parameters)) {
+      rf <- mean(x$replicate_parameters$rf)
+    } else {
+      rf <- x$pmcmc_results$inputs$pars_obs$phi_death
+    }
+
     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$pmcmc_results$inputs$pars_obs$phi_death,
+      ggplot2::geom_point(ggplot2::aes(y=.data$deaths/rf,
                                        x=as.Date(.data$date)), x$pmcmc_results$inputs$data) +
       ggplot2::theme(legend.position = "top")
 

tests/testthat/test-pmcmc.R

@@ -1643,3 +1643,114 @@ test_that("pmcmc restarting covariance/scaling", {
 
 
 })
+
+
+
+#------------------------------------------------
+test_that("sero fitting works", {
+
+  set.seed(12)
+
+  Sys.setenv("SQUIRE_PARALLEL_DEBUG" = "TRUE")
+  data <- read.csv(squire_file("extdata/example.csv"),stringsAsFactors = FALSE)
+  interventions <- read.csv(squire_file("extdata/example_intervention.csv"))
+  int_unique <- interventions_unique(interventions)
+  reporting_fraction = 1
+  country = "Algeria"
+  pars_init = list('start_date'     = as.Date("2020-02-07"),
+                   'R0'             = 3,
+                   'Meff'           = 2,
+                   "rf"             = 0.25) # correct rf for the data
+  pars_min = list('start_date'      = as.Date("2020-02-01"),
+                  'R0'              = 1e-10,
+                  'Meff'            = 0.1,
+                  "rf"              = 0.1)
+  pars_max = list('start_date'      = as.Date("2020-02-20"),
+                  'R0'              = 5,
+                  'Meff'            = 5,
+                  "rf"              = 1)
+  pars_discrete = list('start_date' = TRUE,
+                       'R0'         = FALSE,
+                       'Meff'       = FALSE,
+                       'rf'       = FALSE)
+  pars_obs = list(phi_cases = 0.1,
+                  k_cases = 2,
+                  phi_death = 1,
+                  k_death = 2,
+                  exp_noise = 1e6)
+
+  steps_per_day = 1
+  R0_change = int_unique$change
+  date_R0_change = as.Date(int_unique$dates_change)
+  date_contact_matrix_set_change = NULL
+  squire_model = squire:::deterministic_model()
+  n_particles = 2
+  # proposal kernel covriance
+  proposal_kernel <- matrix(0.5, ncol=length(pars_init), nrow = length(pars_init))
+  diag(proposal_kernel) <- 1
+  rownames(proposal_kernel) <- colnames(proposal_kernel) <- names(pars_init)
+
+  sero_df <- data.frame("samples" = 1000, "sero_pos" = 10,
+                        "date_start" = as.Date("2020-04-15"),
+                        "date_end" = as.Date("2020-04-19"))
+  # seroconversion data from brazeay report 34
+  prob_conversion <-  cumsum(dgamma(0:300,shape = 5, rate = 1/2))/max(cumsum(dgamma(0:300,shape = 5, rate = 1/2)))
+  sero_det <- cumsum(dweibull(0:300, 3.669807, scale = 143.7046))
+  sero_det <- prob_conversion-sero_det
+  sero_det[sero_det < 0] <- 0
+  sero_det <- sero_det/max(sero_det)
+
+  pars_obs$sero_df <- sero_df
+  pars_obs$sero_det <- sero_det
+
+  Sys.setenv("SQUIRE_PARALLEL_DEBUG"=TRUE)
+  out <- pmcmc(data = data,
+               n_mcmc = 5,
+               log_likelihood = NULL,
+               log_prior = NULL,
+               n_particles = 2,
+               steps_per_day = steps_per_day,
+               output_proposals = FALSE,
+               n_chains = 1,
+               replicates = 20,
+               burnin = 5,
+               squire_model = squire_model,
+               pars_init = pars_init,
+               pars_min = pars_min,
+               pars_max = pars_max,
+               pars_discrete = pars_discrete,
+               pars_obs = pars_obs,
+               proposal_kernel = proposal_kernel,
+               R0_change = R0_change,
+               date_R0_change = date_R0_change,
+               country = country)
+
+  pars_init$rf <- 1
+
+  out2 <- pmcmc(data = data,
+               n_mcmc = 5,
+               log_likelihood = NULL,
+               log_prior = NULL,
+               n_particles = 2,
+               steps_per_day = steps_per_day,
+               output_proposals = FALSE,
+               n_chains = 1,
+               replicates = 20,
+               burnin = 5,
+               squire_model = squire_model,
+               pars_init = pars_init,
+               pars_min = pars_min,
+               pars_max = pars_max,
+               pars_discrete = pars_discrete,
+               pars_obs = pars_obs,
+               proposal_kernel = proposal_kernel,
+               R0_change = R0_change,
+               date_R0_change = date_R0_change,
+               country = country)
+
+  expect_gt(sum(out$pmcmc_results$results$log_likelihood),
+            sum(out2$pmcmc_results$results$log_likelihood))
+
+  expect_s3_class(plot(out, what = "deaths", particle_fit = TRUE), "gg")
+
+})