Expose dt/steps_per_day control

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: sircovid
 Title: SIR Model for COVID-19
-Version: 0.9.6
+Version: 0.9.7
 Authors@R:
     c(person(given = "Marc",
            family = "Baguelin",

R/carehomes.R

@@ -141,6 +141,8 @@ NULL
 ##'   rate is used for all age groups; if a vector of values if used it should
 ##'   have one value per age group.
 ##'
+##' @param steps_per_day Number of steps to use per day (integer, at least 1)
+##'
 ##' @param model_pcr_and_serology_user A value of 1 or 0 so switch on or off the
 ##'   flows out of T_PCR_neg and T_sero_neg and the corresponding cap on the
 ##'   number of individuals leaving the R compartments
@@ -276,9 +278,10 @@ carehomes_parameters <- function(start_date, region,
                                  vaccine_daily_doses = 0,
                                  waning_rate = 0,
                                  model_pcr_and_serology_user = 1,
+                                 steps_per_day = 4,
                                  exp_noise = 1e6) {
   ret <- sircovid_parameters_shared(start_date, region,
-                                    beta_date, beta_value)
+                                    beta_date, beta_value, steps_per_day)
 
   ## These are only used here, and are fixed
   carehome_occupancy <- 0.742

R/parameters.R

@@ -1,8 +1,14 @@
 ## These could be moved to be defaults within the models
 sircovid_parameters_shared <- function(start_date, region,
-                                       beta_date, beta_value) {
-  dt <- 0.25
+                                       beta_date, beta_value,
+                                       steps_per_day = 4) {
+  assert_scalar(steps_per_day)
+  assert_integer(steps_per_day)
+  if (steps_per_day < 1) {
+    stop("'steps_per_day' must be at least 1")
+  }
   assert_sircovid_date(start_date)
+  dt <- 1 / steps_per_day
   beta_step <- sircovid_parameters_beta(beta_date, beta_value %||% 0.08, dt)
   list(hosp_transmission = 0.1,
        ICU_transmission = 0.05,

R/support.R

@@ -364,3 +364,32 @@ combine_rt1 <- function(what, rt, samples) {
 
   Reduce(`+`, ret)
 }
+
+
+## adjust the rate of a Gamma/Erlang distribution, for a given time step, 
+## so that the mean of the discretized Gamma/Erlang  matches the mean
+## of the continuous distribution
+## this comes from the fact that the mean of the discretised Erlang distribution
+## is k * dt / (1 - exp(-gamma * dt))
+adjusted_gamma <- function(gamma, k, dt) {
+  if (is.null(k)) k <- 1 # exponential by default
+  true_mean <- k / gamma
+  tmp <- 1 - k*dt / true_mean
+  if(tmp < 0 ) {# mean is < dt so we can't achieve this
+   adj_gamma <- Inf # so setting the rate as high as possible 
+  } else {
+    adj_gamma <- -log(tmp) / dt
+  }
+  adj_gamma
+}
+
+
+adjust_all_gammas <- function(par) {
+  gamma_pars <- names(par)[grep("gamma", names(par))]
+  idx_gamma_pars <- lapply(gamma_pars, function(e) gsub("gamma_", "", e))
+  k_pars <- paste0("k_", idx_gamma_pars)
+  for (i in seq_along(gamma_pars)) {
+    par[[gamma_pars[i]]] <- adjusted_gamma(par[[gamma_pars[i]]], par[[k_pars[i]]], par$dt)
+  }
+  return(par)
+}

tests/testthat/test-carehomes-check.R

@@ -823,3 +823,58 @@ test_that("tots all summed correctly ", {
   expect_true(all(y$sero_pos == apply(y$T_sero_pos[4:13, , , 1, ], 3, sum)))
   expect_true(all(y$react_pos == apply(y$T_PCR_pos[2:18, , , 1, ], 3, sum)))
 })
+
+
+test_that("Can run the models with larger timestep", {
+  n <- c(1, 4)
+  pars <- lapply(n, function(n)
+    carehomes_parameters(0, "london", steps_per_day = n))
+  np <- 50 # use number large enough to get stable mean behaviour 
+
+  mod <- lapply(pars, function(p) carehomes$new(p, 0, np, seed = 1L))
+  end <- sircovid_date("2020-05-31")
+  for (i in seq_along(n)) {
+    info <- mod[[i]]$info()
+    initial <- carehomes_initial(info, np, pars[[i]])
+    mod[[i]]$set_state(initial$state, initial$step)
+    mod[[i]]$set_index(carehomes_index(info)$run)
+  }
+
+  f <- function(m, p) {
+    steps <- seq(0, length.out = 100, by = 1 / p$dt)
+    ## TODO: probably want to set a sensible index here:
+    m$transform_variables(dust::dust_iterate(m, steps))
+  }
+  y <- Map(f, mod, pars)
+
+  ## Below is some comparison to check that the outputs "agree" in
+  ## the broadest sense as they're stochastic and *should not*
+  ## disagree as we've made the integration much coarser.
+
+  ## compare mean attack rates
+  mean_AR_1 <- mean(y[[1]]$cum_infections[1, , 100]) # with dt = 1
+  mean_AR_2 <- mean(y[[2]]$cum_infections[1, , 100]) # with dt = 1/4
+  # compute relative error using dt = 1/4 as reference as more precise 
+  rel_error <- (mean_AR_1 - mean_AR_2) / mean_AR_2 
+  expect_true(rel_error < 1e-2)
+
+  ## compare mean peak size
+  y[[1]]$infections <-
+    sapply(1:np, function(i) diff(y[[1]]$cum_infections[1, i, ]))
+  y[[2]]$infections <-
+    sapply(1:np, function(i) diff(y[[2]]$cum_infections[1, i, ]))
+  mean_peak_1 <- mean(apply(y[[1]]$infections, 2, max)) # with dt = 1
+  mean_peak_2 <- mean(apply(y[[2]]$infections, 2, max)) # with dt = 1/4
+  rel_error <- (mean_peak_1 - mean_peak_2) / mean_peak_2 
+  # peak size is more variable than AR so using higher tolerance
+  expect_true(rel_error < 1e-1) 
+
+  # ## visual check
+  # matplot(t(apply(y[[1]]$infections, 1, quantile, probs = c(0.025, 0.5, 0.975))),
+  #         type = "l", lty = c(2, 1, 2), col = "black",
+  #         xlab = "Time",
+  #         ylab = "Incidence")
+  # matplot(t(apply(y[[2]]$infections, 1, quantile, probs = c(0.025, 0.5, 0.975))),
+  #         type = "l", lty = c(2, 1, 2), col = "blue", add = TRUE)
+
+})

tests/testthat/test-carehomes-rt.R

@@ -142,3 +142,98 @@ test_that("Can vary beta over time", {
   expect_true(all(res$Rt_all >= res$eff_Rt_all))
   expect_true(all(res$Rt_general >= res$eff_Rt_general))
 })
+
+test_that("Can compute Rt with larger timestep", {
+  ### Currently this test is failing: 
+  ### if we do not do the adjustment of parameters (using adjust_all_gammas)
+  ### then the epidemic trajectories are similar but R values don't match
+  ### if we do the adjustment of parameters (using adjust_all_gammas)
+  ### the the R match but not the epidemic trajectories
+
+  n <- c(1, 4, 100)
+  pars <- lapply(n, function(n)
+    carehomes_parameters(0, "london", steps_per_day = n))
+  np <- 50 # use number large enough to get stable mean behaviour 
+
+  ## adjust the gamma parameters to account for the discretisation
+  pars <- lapply(pars, adjust_all_gammas)
+
+  ## run model with different time steps
+  mod <- lapply(pars, function(p) carehomes$new(p, 0, np, seed = 1L))
+  end <- sircovid_date("2020-05-31")
+  for (i in seq_along(n)) {
+    info <- mod[[i]]$info()
+    initial <- carehomes_initial(info, np, pars[[i]])
+    mod[[i]]$set_state(initial$state, initial$step)
+    mod[[i]]$set_index(carehomes_index(info)$run)
+  }
+
+  index <- mod[[1]]$info()$index$S # should be the same for both models
+  f <- function(m, p) {
+    steps <- seq(0, length.out = 100, by = 1 / p$dt)
+    ## TODO: probably want to set a sensible index here:
+    dust::dust_iterate(m, steps, index)
+  }
+  y <- Map(f, mod, pars)
+
+  steps <-
+    lapply(1:3, function(i) seq(0, length.out = 100, by = 1 / pars[[i]]$dt))
+
+  rt_all_1 <- carehomes_Rt_trajectories(steps[[1]], y[[1]], pars[[1]])
+  rt_all_2 <- carehomes_Rt_trajectories(steps[[2]], y[[2]], pars[[2]])
+  rt_all_3 <- carehomes_Rt_trajectories(steps[[3]], y[[3]], pars[[3]])
+
+  ## Below is some comparison to check that the outputs "agree" in
+  ## the broadest sense as they're stochastic and *should not*
+  ## disagree as we've made the integration much coarser.
+
+  ## compare AR or rather n still susceptible
+  mean_S_final_1 <- mean(apply(y[[1]][, , 100], 2, sum))
+  mean_S_final_2 <- mean(apply(y[[2]][, , 100], 2, sum))
+  mean_S_final_3 <- mean(apply(y[[3]][, , 100], 2, sum))
+  # compute relative error using dt = 1/4 as reference
+  rel_error_1 <- (mean_S_final_1 - mean_S_final_2) / mean_S_final_2 
+  expect_true(rel_error_1 < 5e-2)
+  rel_error_3 <- (mean_S_final_3 - mean_S_final_2) / mean_S_final_2 
+  expect_true(rel_error_3 < 5e-2)
+
+  ## compare mean reproduction numbers 
+  mean_Rt_general_1 <- apply(rt_all_1$Rt_general, 1, mean) # with dt = 1
+  mean_Rt_general_2 <- apply(rt_all_2$Rt_general, 1, mean) # with dt = 1/4
+  mean_Rt_general_3 <- apply(rt_all_3$Rt_general, 1, mean) # with dt = 1/100
+  # compute relative error using dt = 1/4 as reference as more precise 
+  rel_error_1 <- (mean_Rt_general_1 - mean_Rt_general_2) / mean_Rt_general_2 
+  expect_true(all(rel_error_1 < 1e-2))
+  rel_error_3 <- (mean_Rt_general_3 - mean_Rt_general_2) / mean_Rt_general_2 
+  expect_true(all(rel_error_3 < 1e-2))
+
+  ## compare mean reproduction numbers 
+  mean_Rt_all_1 <- apply(rt_all_1$Rt_all, 1, mean) # with dt = 1
+  mean_Rt_all_2 <- apply(rt_all_2$Rt_all, 1, mean) # with dt = 1/4
+  # compute relative error using dt = 1/4 as reference as more precise 
+  rel_error <- (mean_Rt_all_1 - mean_Rt_all_2) / mean_Rt_all_2 
+  expect_true(all(rel_error < 1e-2))
+
+  ## compare mean reproduction numbers
+  mean_eff_Rt_general_1 <- apply(rt_all_1$eff_Rt_general, 1, mean) # with dt = 1
+  mean_eff_Rt_general_2 <- apply(rt_all_2$eff_Rt_general, 1, mean) # with dt = 1/4
+  # compute relative error using dt = 1/4 as reference as more precise 
+  rel_error <- (mean_eff_Rt_general_1 - mean_eff_Rt_general_2) / mean_eff_Rt_general_2 
+  expect_true(all(rel_error < 1e-2))
+
+  ## compare mean reproduction numbers
+  mean_eff_Rt_all_1 <- apply(rt_all_1$eff_Rt_all, 1, mean) # with dt = 1
+  mean_eff_Rt_all_2 <- apply(rt_all_2$eff_Rt_all, 1, mean) # with dt = 1/4
+  # compute relative error using dt = 1/4 as reference as more precise 
+  rel_error <- (mean_eff_Rt_all_1 - mean_eff_Rt_all_2) / mean_eff_Rt_all_2 
+  expect_true(all(rel_error < 1e-2))
+
+  # ## visual check
+  # matplot(t(apply(rt_all_1$eff_Rt_all, 1, quantile, probs = c(0.025, 0.5, 0.975))),
+  #         type = "l", lty = c(2, 1, 2), col = "black",
+  #         xlab = "Time",
+  #         ylab = "Effective Rt")
+  # matplot(t(apply(rt_all_2$eff_Rt_all, 1, quantile, probs = c(0.025, 0.5, 0.975))),
+  #         type = "l", lty = c(2, 1, 2), col = "blue", add = TRUE)
+
+})

tests/testthat/test-parameters.R

@@ -141,3 +141,19 @@ test_that("can expand beta", {
   beta3 <- sircovid_parameters_beta_expand(t3, beta)
   expect_equal(beta3, beta[1:65])
 })
+
+
+test_that("can change dt", {
+  expect_equal(sircovid_parameters_shared(0, "uk", NULL, NULL)$dt, 0.25)
+  expect_equal(sircovid_parameters_shared(0, "uk", NULL, NULL, 10)$dt, 0.1)
+  expect_equal(sircovid_parameters_shared(0, "uk", NULL, NULL, 1)$dt, 1)
+  expect_error(
+    sircovid_parameters_shared(0, "uk", NULL, NULL, 0),
+    "'steps_per_day' must be at least 1")
+  expect_error(
+    sircovid_parameters_shared(0, "uk", NULL, NULL, 0:1),
+    "'steps_per_day' must be a scalar")
+  expect_error(
+    sircovid_parameters_shared(0, "uk", NULL, NULL, 1.5),
+    "'steps_per_day' must be an integer")
+})

tests/testthat/test-utils.R

@@ -133,3 +133,47 @@ test_that("spread_integer", {
   expect_equal(spread_integer(0, 5), rep(0, 5))
   expect_equal(spread_integer(2, 5), c(1, 1, 0, 0, 0))
 })
+
+
+test_that("Gamma correction works", {
+  ## tests that adjusted_gamma improves match to mean of continuous distribution
+
+  ## compute the mean of a discretized Gamma/Erlang distribution
+  discretized_mean <- function(gamma, k, dt) {
+    max_t <- qgamma(0.999, shape = k, rate = gamma)
+    tt <- seq(0, max_t + 1, dt)
+    cdf <- pgamma(tt, shape = k, rate = gamma)
+    sum(tt[-length(tt)] * diff(cdf))
+  }
+
+  test_gamma_correction <- function(gamma, k, dt)
+  {
+    error_no_correction <- abs(discretized_mean(gamma, k, dt) - k / gamma)
+    error_with_correction <- abs(discretized_mean(adjusted_gamma(gamma, k, dt), k, dt) - k / gamma)
+    expect_true(error_no_correction > error_with_correction)
+  }
+
+  test_gamma_correction(gamma = 1 / 2.5, k = 1, dt = 1)
+  test_gamma_correction(gamma = 1 / 2.5, k = 1, dt = 1/4)
+  test_gamma_correction(gamma = 1 / 2.5, k = 1, dt = 1/100)
+
+  test_gamma_correction(gamma = 1 / 2.5, k = 2, dt = 1)
+  test_gamma_correction(gamma = 1 / 2.5, k = 2, dt = 1/4)
+  test_gamma_correction(gamma = 1 / 2.5, k = 2, dt = 1/100)
+
+  test_gamma_correction(gamma = 1 / 1, k = 1, dt = 1)
+  test_gamma_correction(gamma = 1 / 1, k = 1, dt = 1/4)
+  test_gamma_correction(gamma = 1 / 1, k = 1, dt = 1/100)
+
+  test_gamma_correction(gamma = 1 / 1, k = 2, dt = 1)
+  test_gamma_correction(gamma = 1 / 1, k = 2, dt = 1/4)
+  test_gamma_correction(gamma = 1 / 1, k = 2, dt = 1/100)
+
+  test_gamma_correction(gamma = 1 / 10, k = 1, dt = 1)
+  test_gamma_correction(gamma = 1 / 10, k = 1, dt = 1/4)
+  test_gamma_correction(gamma = 1 / 10, k = 1, dt = 1/100)
+
+  test_gamma_correction(gamma = 1 / 10, k = 2, dt = 1)
+  test_gamma_correction(gamma = 1 / 10, k = 2, dt = 1/4)
+  test_gamma_correction(gamma = 1 / 10, k = 2, dt = 1/100)
+})