Add intervention & restructure simulate_tree_from_pop()

R/intervention.R

@@ -0,0 +1,51 @@
+#' Reduce the mean of the offspring distribution
+#'
+#' @description
+#' `intvn_reduce_mean()` is a helper for the \code{simulate_*()} functions. It
+#' reduces/scales the mean of the offspring distribution in order to
+#' mimic the impact of a population-level intervention. Currently, it can only
+#' handle the poisson and negative binomial distributions and errors when other
+#' offspring distributions are specified alongside the `intvn_mean_reduction`
+#' argument.
+#'
+#' @inheritParams simulate_tree
+#' @param intvn_mean_reduction A number between 0
+#' and 1 for scaling/reducing the mean of `offspring_dist`. Serves as
+#' population-level intervention. `intvn_mean_reduction` = 0
+#' implies no intervention impact and `intvn_mean_reduction` = 1 implies full
+#' impact.
+#' @param pars_list Parameter(s) for poisson or negative binomial offspring
+#' distribution.
+#' @return List of the offspring distribution parameter(s) with the mean
+#' scaled by \code{1 - intvn_mean_reduction}.
+#' @details
+#' `intvn_reduce_mean()` scales the mean of the offspring distribution
+#' by \eqn{1 - {\sf intvn\_mean\_reduction}} so that the new mean is given as:
+#' \deqn{(1 - {\sf intvn\_mean\_reduction}) \times {\sf mean,}} This
+#' scaling when applied to the poisson and negative binomial offspring
+#' distributions corresponds to the population-level reduction of R0 as
+#' described in Lloyd-Smith et al, (2005). `intvn_reduce_mean()` is therefore
+#' only implemented for the aforementioned distributions and errors when other
+#' offspring distributions are specified along with the `intvn_mean_reduction`
+#' argument in the \code{simulate_*()} functions.
+#'
+#' @author James M. Azam
+#' @keywords internal
+#' @references Lloyd-Smith, J., Schreiber, S., Kopp, P. et al. Superspreading
+#' and the effect of individual variation on disease emergence. Nature 438,
+#' 355–359 (2005). \doi{10.1038/nature04153}
+intvn_reduce_mean <- function(intvn_mean_reduction, offspring_dist, pars_list) {
+  # Intervention only works for pois and nbinom
+  if (!offspring_dist %in% c("pois", "nbinom")) {
+    stop(
+      "`offspring_dist` must be one of c(\"pois\", \"nbinom\"), ",
+      "if intvn_mean_reduction is specified."
+    )
+  }
+  if (offspring_dist == "pois") {
+    pars_list$lambda <- (1 - intvn_mean_reduction) * pars_list$lambda
+  } else {
+    pars_list$mu <- (1 - intvn_mean_reduction) * pars_list$mu
+  }
+  return(pars_list)
+}

R/simulate.r

@@ -1,5 +1,6 @@
 #' Simulate transmission trees from an initial number of infections
 #'
+#' @inheritParams intvn_reduce_mean
 #' @param nchains Number of chains to simulate.
 #' @param offspring_dist Offspring distribution: a character string
 #' corresponding to the R distribution function (e.g., "pois" for Poisson,
@@ -97,6 +98,19 @@
 #'   serials_dist = function(x) 3,
 #'   lambda = 2
 #' )
+#'
+#' # Run model with intervention a 50% reduction in R0.
+#' chains_with_intvn <- simulate_tree(
+#'   nchains = 10,
+#'   statistic = "size",
+#'   offspring_dist = "pois",
+#'   intvn_mean_reduction = 0.5,
+#'   stat_max = 10,
+#'   serials_dist = function(x) 3,
+#'   lambda = 2
+#' )
+#'
+#' chains_with_intvn
 #' @references
 #' Lehtinen S, Ashcroft P, Bonhoeffer S. On the relationship
 #' between serial interval, infectiousness profile and generation time.
@@ -113,6 +127,7 @@
 #' 1186–1204. \doi{https://doi.org/10.3390/ijerph7031204}
 simulate_tree <- function(nchains, statistic = c("size", "length"),
                           offspring_dist, stat_max = Inf,
+                          intvn_mean_reduction = 0,
                           serials_dist, t0 = 0,
                           tf = Inf, ...) {
   statistic <- match.arg(statistic)
@@ -122,10 +137,29 @@ simulate_tree <- function(nchains, statistic = c("size", "length"),
   # check that offspring is properly specified
   check_offspring_valid(offspring_dist)
 
+  # Check that the intvn_mean_reduction is well specified
+  checkmate::assert_number(
+    intvn_mean_reduction,
+    lower = 0,
+    upper = 1
+  )
+
   # check that offspring function exists in base R
   roffspring_name <- paste0("r", offspring_dist)
   check_offspring_func_valid(roffspring_name)
 
+  # Gather offspring distribution parameters
+  pars <- list(...)
+
+  # Prepare interventions if specified
+  if (intvn_mean_reduction > 0) {
+    pars <- intvn_reduce_mean(
+      intvn_mean_reduction = intvn_mean_reduction,
+      offspring_dist = offspring_dist,
+      pars_list = pars
+    )
+  }
+
   if (!missing(serials_dist)) {
     check_serial_valid(serials_dist)
   } else if (!missing(tf)) {
@@ -155,7 +189,13 @@ simulate_tree <- function(nchains, statistic = c("size", "length"),
   # next, simulate n chains
   while (length(sim) > 0) {
     # simulate next generation
-    next_gen <- get(roffspring_name)(n = sum(n_offspring[sim]), ...)
+    next_gen <- do.call(
+      get(roffspring_name),
+      c(
+        list(n = sum(n_offspring[sim])),
+        pars
+      )
+    )
     if (any(next_gen %% 1 > 0)) {
       stop("Offspring distribution must return integers")
     }
@@ -244,6 +284,7 @@ simulate_tree <- function(nchains, statistic = c("size", "length"),
 #' Simulate transmission chains sizes/lengths
 #'
 #' @inheritParams simulate_tree
+#' @inheritParams intvn_reduce_mean
 #' @param stat_max A cut off for the chain statistic (size/length) being
 #' computed. Results above the specified value, are set to `Inf`.
 #' @inheritSection simulate_tree Calculating chain sizes and lengths
@@ -262,9 +303,22 @@ simulate_tree <- function(nchains, statistic = c("size", "length"),
 #'   stat_max = 10,
 #'   lambda = 2
 #' )
+#'
+#' # Run model with intervention a 50% reduction in R0.
+#' chain_summary_with_intvn <- simulate_summary(
+#'   nchains = 10,
+#'   statistic = "size",
+#'   offspring_dist = "pois",
+#'   intvn_mean_reduction = 0.5,
+#'   stat_max = 10,
+#'   lambda = 2
+#' )
+#'
+#' chain_summary_with_intvn
 #' @export
 simulate_summary <- function(nchains, statistic = c("size", "length"),
                              offspring_dist,
+                             intvn_mean_reduction = 0,
                              stat_max = Inf, ...) {
   statistic <- match.arg(statistic)
 
@@ -273,10 +327,29 @@ simulate_summary <- function(nchains, statistic = c("size", "length"),
   # check that offspring is properly specified
   check_offspring_valid(offspring_dist)
 
+  # Check that the intvn_mean_reduction is well specified
+  checkmate::assert_number(
+    intvn_mean_reduction,
+    lower = 0,
+    upper = 1
+  )
+
   # check that offspring function exists in base R
   roffspring_name <- paste0("r", offspring_dist)
   check_offspring_func_valid(roffspring_name)
 
+  # Gather offspring distribution parameters
+  pars <- list(...)
+
+  # Prepare interventions if specified
+  if (intvn_mean_reduction > 0) {
+    pars <- intvn_reduce_mean(
+      intvn_mean_reduction = intvn_mean_reduction,
+      offspring_dist = offspring_dist,
+      pars_list = pars
+    )
+  }
+
   # Initialisations
   stat_track <- rep(1, nchains) ## track length or size (depending on `stat`)
   n_offspring <- rep(1, nchains) ## current number of offspring
@@ -285,7 +358,13 @@ simulate_summary <- function(nchains, statistic = c("size", "length"),
   ## next, simulate nchains chains
   while (length(sim) > 0) {
     ## simulate next generation
-    next_gen <- get(roffspring_name)(n = sum(n_offspring[sim]), ...)
+    next_gen <- do.call(
+      get(roffspring_name),
+      c(
+        list(n = sum(n_offspring[sim])),
+        pars
+      )
+    )
     if (any(next_gen %% 1 > 0)) {
       stop("Offspring distribution must return integers")
     }
@@ -325,17 +404,12 @@ simulate_summary <- function(nchains, statistic = c("size", "length"),
 #' population
 #'
 #' @inheritParams simulate_tree
+#' @inheritParams intvn_mean_reduction
 #' @param pop The susceptible population size.
 #' @param offspring_dist Offspring distribution: a character string
 #' corresponding to the R distribution function (e.g., "pois" for Poisson,
 #' where \code{\link{rpois}} is the R function to generate Poisson random
 #' numbers). Only supports "pois" and "nbinom".
-#' @param offspring_mean The average number of secondary cases for each case.
-#' Same as \eqn{R_0}.
-#' @param offspring_disp The dispersion parameter of the number of
-#' secondary cases. Ignored if \code{offspring == "pois"}. Must be > 1 to
-#' avoid division by 0 when calculating the size. See details and
-#'  \code{?rnbinom} for details on the parameterisation in Ecology.
 #' @param initial_immune The number of initial immunes in the population.
 #' Must be less than `pop` - 1.
 #' @param t0 Start time; Defaults to 0.
@@ -352,15 +426,15 @@ simulate_summary <- function(nchains, statistic = c("size", "length"),
 #' than susceptibles at any point.
 #'
 #' The poisson model has mean, lambda, parametrised as:
-#' \deqn{{\sf lambda} = \dfrac{{\sf offspring\_mean} \times ({\sf pop} -
+#' \deqn{{\sf lambda} = \dfrac{{\sf lambda} \times ({\sf pop} -
 #' {\sf initial\_immune} - 1)}{{\sf pop}}}
 #'
 #' The negative binomial model, has mean, mu, parametrised as:
-#' \deqn{{\sf mu} = \dfrac{{\sf offspring\_mean} \times ({\sf pop} -
+#' \deqn{{\sf mu} = \dfrac{{\sf mu} \times ({\sf pop} -
 #' {\sf initial\_immune} - 1)}{{\sf pop}},}
 #' and dispersion, size, parametrised as:
-#' \deqn{{\sf size} = \dfrac{{\sf mu}}{{\sf offspring\_disp} - 1}.}
-#' This is why `offspring_disp` must be greater than 1.
+#' \deqn{{\sf size} = \dfrac{{\sf mu}}{{\sf size} - 1}.}
+#' This is why `size` must be greater than 1.
 #'
 #' # Differences with `simulate_tree()`
 #' `simulate_tree_from_pop()` has a couple of key differences from
@@ -380,66 +454,86 @@ simulate_summary <- function(nchains, statistic = c("size", "length"),
 #' simulate_tree_from_pop(
 #'   pop = 100,
 #'   offspring_dist = "pois",
-#'   offspring_mean = 0.5,
+#'   lambda = 0.5,
 #'   serials_dist = function(x) 3
 #' )
 #'
 #' # Simulate with negative binomial offspring
 #' simulate_tree_from_pop(
+#' pop = 100, offspring_dist = "nbinom",
+#' mu = 0.5,
+#' size = 1.1,
+#' serials_dist = function(x) 3
+#' )
+#'
+#' # Simulate with negative binomial offspring with intervention (50%
+#' # reduction in R0)
+#' simulate_tree_from_pop(
 #'   pop = 100,
 #'   offspring_dist = "nbinom",
-#'   offspring_mean = 0.5,
-#'   offspring_disp = 1.1,
+#'   intvn_mean_reduction = 0.5,
+#'   mu = 0.5,
+#'   size = 1.1,
 #'   serials_dist = function(x) 3
 #' )
 #' @export
 simulate_tree_from_pop <- function(pop,
                                    offspring_dist = c("pois", "nbinom"),
-                                   offspring_mean,
-                                   offspring_disp,
+                                   intvn_mean_reduction = 0,
                                    serials_dist,
                                    initial_immune = 0,
                                    t0 = 0,
-                                   tf = Inf) {
+                                   tf = Inf,
+                                   ...) {
   offspring_dist <- match.arg(offspring_dist)
 
-  if (offspring_dist == "pois") {
-    if (!missing(offspring_disp)) {
-      warning(sprintf(
-        "%s %s %s",
-        "'offspring_disp' is not used for",
-        "poisson offspring distribution.",
-        "Will be ignored."
-      ))
-    }
+  # Check that the intvn_mean_reduction is well specified
+  checkmate::assert_number(
+    intvn_mean_reduction,
+    lower = 0,
+    upper = 1
+  )
 
-    ## using a right truncated poisson distribution
+  # Gather offspring distribution parameters
+  pars <- list(...)
+
+  # Prepare interventions if specified
+  if (intvn_mean_reduction > 0) {
+    pars <- intvn_reduce_mean(
+      intvn_mean_reduction = intvn_mean_reduction,
+      offspring_dist = offspring_dist,
+      pars_list = pars
+    )
+  }
+
+  if (offspring_dist == "pois") {
+    ## Use a right truncated poisson distribution
     ## to avoid more cases than susceptibles
     offspring_func <- function(n, susc) {
       truncdist::rtrunc(
         n,
         spec = "pois",
-        lambda = offspring_mean * susc / pop,
+        lambda = pars$lambda * susc / pop,
         b = susc
       )
     }
   } else if (offspring_dist == "nbinom") {
-    if (missing(offspring_disp)) {
-      stop(sprintf("%s", "'offspring_disp' must be specified."))
-    } else if (offspring_disp <= 1) { ## dispersion coefficient
+    if (is.null(pars$size)) {
+      stop(sprintf("%s", "'size' must be specified."))
+    } else if (pars$size <= 1) { ## dispersion coefficient
       stop(sprintf(
         "%s %s %s",
         "Offspring distribution 'nbinom' requires",
-        "argument 'offspring_disp' > 1.",
+        "argument 'size' > 1.",
         "Use 'pois' if there is no overdispersion."
       ))
     }
     ## get distribution params from mean and dispersion
     offspring_func <- function(n, susc) {
       ## get distribution params from mean and dispersion
       ## see ?rnbinom for parameter definition
-      new_mn <- offspring_mean * susc / pop ## apply susceptibility
-      size <- new_mn / (offspring_disp - 1)
+      new_mn <- pars$mu * susc / pop ## apply susceptibility
+      size <- new_mn / (pars$size - 1)
 
       ## using a right truncated nbinom distribution
       ## to avoid more cases than susceptibles