From 2bc906f1876588f8939968c3ac52d53958efec4f Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Wed, 8 Nov 2023 14:29:52 +0000
Subject: [PATCH 1/6] remove `intvn_reduce_mean` argument

---
 R/intervention.R               |  51 ----------
 R/simulate.r                   |  90 ------------------
 man/intvn_reduce_mean.Rd       |  55 -----------
 man/simulate_summary.Rd        |  19 ----
 man/simulate_tree.Rd           |  20 ----
 man/simulate_tree_from_pop.Rd  |  18 ----
 tests/testthat/test-simulate.R | 169 ---------------------------------
 vignettes/epichains.Rmd        |  66 -------------
 8 files changed, 488 deletions(-)
 delete mode 100644 R/intervention.R
 delete mode 100644 man/intvn_reduce_mean.Rd

diff --git a/R/intervention.R b/R/intervention.R
deleted file mode 100644
index 3ab7339c..00000000
--- a/R/intervention.R
+++ /dev/null
@@ -1,51 +0,0 @@
-#' 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)
-}
diff --git a/R/simulate.r b/R/simulate.r
index 11f34afe..c5b81d1b 100644
--- a/R/simulate.r
+++ b/R/simulate.r
@@ -1,6 +1,5 @@
 #' 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,
@@ -98,19 +97,6 @@
 #'   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.
@@ -127,7 +113,6 @@
 #' 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)
@@ -137,13 +122,6 @@ 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)
@@ -151,15 +129,6 @@ simulate_tree <- function(nchains, statistic = c("size", "length"),
   # 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)) {
@@ -284,7 +253,6 @@ 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
@@ -303,22 +271,9 @@ 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)
 
@@ -327,13 +282,6 @@ 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)
@@ -341,15 +289,6 @@ simulate_summary <- function(nchains, statistic = c("size", "length"),
   # 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
@@ -404,7 +343,6 @@ 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,
@@ -465,21 +403,9 @@ simulate_summary <- function(nchains, statistic = c("size", "length"),
 #' 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",
-#'   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"),
-                                   intvn_mean_reduction = 0,
                                    serials_dist,
                                    initial_immune = 0,
                                    t0 = 0,
@@ -487,25 +413,9 @@ simulate_tree_from_pop <- function(pop,
                                    ...) {
   offspring_dist <- match.arg(offspring_dist)
 
-  # Check that the intvn_mean_reduction is well specified
-  checkmate::assert_number(
-    intvn_mean_reduction,
-    lower = 0,
-    upper = 1
-  )
-
   # 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
diff --git a/man/intvn_reduce_mean.Rd b/man/intvn_reduce_mean.Rd
deleted file mode 100644
index 080ec2a9..00000000
--- a/man/intvn_reduce_mean.Rd
+++ /dev/null
@@ -1,55 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/intervention.R
-\name{intvn_reduce_mean}
-\alias{intvn_reduce_mean}
-\title{Reduce the mean of the offspring distribution}
-\usage{
-intvn_reduce_mean(intvn_mean_reduction, offspring_dist, pars_list)
-}
-\arguments{
-\item{intvn_mean_reduction}{A number between 0
-and 1 for scaling/reducing the mean of \code{offspring_dist}. Serves as
-population-level intervention. \code{intvn_mean_reduction} = 0
-implies no intervention impact and \code{intvn_mean_reduction} = 1 implies full
-impact.}
-
-\item{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).}
-
-\item{pars_list}{Parameter(s) for poisson or negative binomial offspring
-distribution.}
-}
-\value{
-List of the offspring distribution parameter(s) with the mean
-scaled by \code{1 - intvn_mean_reduction}.
-}
-\description{
-\code{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 \code{intvn_mean_reduction}
-argument.
-}
-\details{
-\code{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). \code{intvn_reduce_mean()} is therefore
-only implemented for the aforementioned distributions and errors when other
-offspring distributions are specified along with the \code{intvn_mean_reduction}
-argument in the \code{simulate_*()} functions.
-}
-\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}
-}
-\author{
-James M. Azam
-}
-\keyword{internal}
diff --git a/man/simulate_summary.Rd b/man/simulate_summary.Rd
index 9c1283af..a453d680 100644
--- a/man/simulate_summary.Rd
+++ b/man/simulate_summary.Rd
@@ -8,7 +8,6 @@ simulate_summary(
   nchains,
   statistic = c("size", "length"),
   offspring_dist,
-  intvn_mean_reduction = 0,
   stat_max = Inf,
   ...
 )
@@ -29,12 +28,6 @@ corresponding to the R distribution function (e.g., "pois" for Poisson,
 where \code{\link{rpois}} is the R function to generate Poisson random
 numbers).}
 
-\item{intvn_mean_reduction}{A number between 0
-and 1 for scaling/reducing the mean of \code{offspring_dist}. Serves as
-population-level intervention. \code{intvn_mean_reduction} = 0
-implies no intervention impact and \code{intvn_mean_reduction} = 1 implies full
-impact.}
-
 \item{stat_max}{A cut off for the chain statistic (size/length) being
 computed. Results above the specified value, are set to \code{Inf}.}
 
@@ -105,18 +98,6 @@ simulate_summary(
   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
 }
 \seealso{
 \itemize{
diff --git a/man/simulate_tree.Rd b/man/simulate_tree.Rd
index cb3d0629..1d6722dc 100644
--- a/man/simulate_tree.Rd
+++ b/man/simulate_tree.Rd
@@ -9,7 +9,6 @@ simulate_tree(
   statistic = c("size", "length"),
   offspring_dist,
   stat_max = Inf,
-  intvn_mean_reduction = 0,
   serials_dist,
   t0 = 0,
   tf = Inf,
@@ -36,12 +35,6 @@ numbers).}
 computed. Results above the specified value, are set to this value.
 Defaults to \code{Inf}.}
 
-\item{intvn_mean_reduction}{A number between 0
-and 1 for scaling/reducing the mean of \code{offspring_dist}. Serves as
-population-level intervention. \code{intvn_mean_reduction} = 0
-implies no intervention impact and \code{intvn_mean_reduction} = 1 implies full
-impact.}
-
 \item{serials_dist}{The serial interval distribution function; the name
 of a user-defined named or anonymous function with only one argument \code{n},
 representing the number of serial intervals to generate. See details.}
@@ -128,19 +121,6 @@ chains <- simulate_tree(
   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
diff --git a/man/simulate_tree_from_pop.Rd b/man/simulate_tree_from_pop.Rd
index 030f9fe8..4f42d429 100644
--- a/man/simulate_tree_from_pop.Rd
+++ b/man/simulate_tree_from_pop.Rd
@@ -8,7 +8,6 @@ population}
 simulate_tree_from_pop(
   pop,
   offspring_dist = c("pois", "nbinom"),
-  intvn_mean_reduction = 0,
   serials_dist,
   initial_immune = 0,
   t0 = 0,
@@ -24,12 +23,6 @@ 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".}
 
-\item{intvn_mean_reduction}{A number between 0
-and 1 for scaling/reducing the mean of \code{offspring_dist}. Serves as
-population-level intervention. \code{intvn_mean_reduction} = 0
-implies no intervention impact and \code{intvn_mean_reduction} = 1 implies full
-impact.}
-
 \item{serials_dist}{The serial interval distribution function; the name
 of a user-defined named or anonymous function with only one argument \code{n},
 representing the number of serial intervals to generate. See details.}
@@ -137,17 +130,6 @@ 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",
-  intvn_mean_reduction = 0.5,
-  mu = 0.5,
-  size = 1.1,
-  serials_dist = function(x) 3
-)
 }
 \seealso{
 \itemize{
diff --git a/tests/testthat/test-simulate.R b/tests/testthat/test-simulate.R
index b6fd787e..4852fd9c 100644
--- a/tests/testthat/test-simulate.R
+++ b/tests/testthat/test-simulate.R
@@ -20,25 +20,6 @@ test_that("Simulators work", {
     size = 1.1,
     serials_dist = serial_func
   )
-  #' Simulate an outbreak from a susceptible population (pois) with
-  #' 50% R0 reduction
-  susc_outbreak_raw_intvn <- simulate_tree_from_pop(
-    pop = 100,
-    offspring_dist = "pois",
-    lambda = 1.5,
-    serials_dist = serial_func,
-    intvn_mean_reduction = 0.5
-  )
-  #' Simulate an outbreak from a susceptible population (nbinom) with
-  #' 50% R0 reduction
-  susc_outbreak_raw_intvn2 <- simulate_tree_from_pop(
-    pop = 100,
-    offspring_dist = "nbinom",
-    mu = 1.5,
-    size = 1.1,
-    serials_dist = serial_func,
-    intvn_mean_reduction = 0.5
-  )
   #' Simulate a tree of infections without serials
   tree_sim_raw <- simulate_tree(
     nchains = 2,
@@ -55,25 +36,6 @@ test_that("Simulators work", {
     serials_dist = function(x) 3,
     lambda = 2
   )
-  #' Simulate a tree of infections without serials and with 50% reduction
-  #' in R0
-  tree_sim_raw_intvn <- simulate_tree(
-    nchains = 2,
-    offspring_dist = "pois",
-    statistic = "length",
-    lambda = 0.9,
-    intvn_mean_reduction = 0.5
-  )
-  #' Simulate a tree of infections with nbinom offspring and with 50% reduction
-  #' in R0
-  tree_sim_raw_intvn2 <- simulate_tree(
-    nchains = 2,
-    offspring_dist = "nbinom",
-    statistic = "length",
-    mu = 0.9,
-    size = 1.1,
-    intvn_mean_reduction = 0.5
-  )
   #' Simulate chain statistics
   chain_summary_raw <- simulate_summary(
     nchains = 2,
@@ -81,37 +43,11 @@ test_that("Simulators work", {
     statistic = "length",
     lambda = 0.9
   )
-  #' Simulate chain statistics and with a 50% reduction in R0
-  chain_summary_raw_intvn <- simulate_summary(
-    nchains = 2,
-    offspring_dist = "pois",
-    statistic = "length",
-    lambda = 0.9,
-    intvn_mean_reduction = 0.5
-  )
-  #' Simulate chain statistics with nbinom offspring and with a 50% reduction
-  #' in R0
-  chain_summary_raw_intvn2 <- simulate_summary(
-    nchains = 2,
-    offspring_dist = "nbinom",
-    statistic = "length",
-    mu = 1.9,
-    size = 1.1,
-    intvn_mean_reduction = 0.5
-  )
   #' Expectations
   expect_length(
     chain_summary_raw,
     2
   )
-  expect_length(
-    chain_summary_raw_intvn,
-    2
-  )
-  expect_length(
-    chain_summary_raw_intvn2,
-    2
-  )
   expect_gte(
     nrow(tree_sim_raw),
     2
@@ -120,14 +56,6 @@ test_that("Simulators work", {
     nrow(tree_sim_raw2),
     2
   )
-  expect_identical(
-    nrow(tree_sim_raw_intvn),
-    3L
-  )
-  expect_identical(
-    nrow(tree_sim_raw_intvn2),
-    2L
-  )
   expect_gte(
     nrow(susc_outbreak_raw),
     1
@@ -136,14 +64,6 @@ test_that("Simulators work", {
     nrow(susc_outbreak_raw2),
     1
   )
-  expect_identical(
-    nrow(susc_outbreak_raw_intvn),
-    2L
-  )
-  expect_identical(
-    nrow(susc_outbreak_raw_intvn2),
-    1L
-  )
   expect_true(
     all(
       simulate_tree(
@@ -218,17 +138,6 @@ test_that("simulate_tree throws errors", {
     ),
     "must be specified"
   )
-  expect_error(
-    simulate_tree(
-      nchains = 2,
-      offspring_dist = "binom",
-      statistic = "length",
-      size = 1,
-      prob = 0.5,
-      intvn_mean_reduction = 0.5
-    ),
-    "must be one of"
-  )
 })
 
 test_that("simulate_summary throws errors", {
@@ -270,17 +179,6 @@ test_that("simulate_summary throws errors", {
     ),
     "character string"
   )
-  expect_error(
-    simulate_summary(
-      nchains = 2,
-      offspring_dist = "binom",
-      statistic = "length",
-      size = 1,
-      prob = 0.5,
-      intvn_mean_reduction = 0.5
-    ),
-    "must be one of"
-  )
 })
 
 test_that("simulate_tree_from_pop throws errors", {
@@ -333,18 +231,8 @@ test_that("simulate_tree is numerically correct", {
     statistic = "length",
     lambda = 0.9
   )
-  #' Simulate a tree of infections without serials and with 50% reduction
-  #' in R0
-  tree_sim_raw_intvn <- simulate_tree(
-    nchains = 2,
-    offspring_dist = "pois",
-    statistic = "length",
-    lambda = 0.9,
-    intvn_mean_reduction = 0.5
-  )
   #' summarise the results
   tree_sim_summary <- summary(tree_sim_raw)
-  tree_sim_intvn_summary <- summary(tree_sim_raw_intvn)
   #' Expectations
   expect_identical(
     tree_sim_summary$chains_run,
@@ -414,17 +302,8 @@ test_that("simulate_summary is numerically correct", {
     statistic = "length",
     lambda = 0.9
   )
-  #' Simulate chain statistics and with a 50% reduction in R0
-  chain_summary_raw_intvn <- simulate_summary(
-    nchains = 2,
-    offspring_dist = "pois",
-    statistic = "length",
-    lambda = 0.9,
-    intvn_mean_reduction = 0.5
-  )
   #' Summarise the results
   chain_summary_summaries <- summary(chain_summary_raw)
-  chain_summary_intvn_summaries <- summary(chain_summary_raw_intvn)
   #' Expectations
   expect_identical(
     chain_summary_summaries$chains_run,
@@ -442,22 +321,6 @@ test_that("simulate_summary is numerically correct", {
     as.vector(chain_summary_raw),
     c(1.00, 3.00)
   )
-  expect_identical(
-    chain_summary_intvn_summaries$chains_run,
-    2.00
-  )
-  expect_identical(
-    chain_summary_intvn_summaries$max_chain_stat,
-    2.00
-  )
-  expect_identical(
-    chain_summary_intvn_summaries$min_chain_stat,
-    1.00
-  )
-  expect_identical(
-    as.vector(chain_summary_raw_intvn),
-    c(2.00, 1.00)
-  )
 })
 
 test_that("simulate_tree_from_pop is numerically correct", {
@@ -469,19 +332,8 @@ test_that("simulate_tree_from_pop is numerically correct", {
     lambda = 0.9,
     serials_dist = serial_func
   )
-  #' Simulate an outbreak from a susceptible population (pois) with
-  #' 50% R0 reduction
-  set.seed(7)
-  susc_outbreak_raw_intvn <- simulate_tree_from_pop(
-    pop = 100,
-    offspring_dist = "pois",
-    lambda = 1.5,
-    serials_dist = serial_func,
-    intvn_mean_reduction = 0.5
-  )
   #' Summarise the results
   susc_outbreak_summary <- summary(susc_outbreak_raw)
-  susc_outbreak_summary_intvn <- summary(susc_outbreak_raw_intvn)
   #' Expectations
   expect_identical(
     susc_outbreak_summary$unique_ancestors,
@@ -512,25 +364,4 @@ test_that("simulate_tree_from_pop is numerically correct", {
     susc_outbreak_raw$time,
     0.00
   )
-  #' Expectations for intervention simulation
-  expect_identical(
-    susc_outbreak_summary_intvn$unique_ancestors,
-    12L
-  )
-  expect_identical(
-    round(
-      susc_outbreak_summary_intvn$max_time,
-      1
-    ),
-    72.1
-  )
-  expect_identical(
-    susc_outbreak_summary_intvn$max_generation,
-    10L
-  )
-  expect_null(susc_outbreak_summary_intvn$chains_run)
-  expect_identical(
-    sum(aggregate(susc_outbreak_raw_intvn, "time")$cases),
-    20L
-  )
 })
diff --git a/vignettes/epichains.Rmd b/vignettes/epichains.Rmd
index 3f8c0b90..7d361795 100644
--- a/vignettes/epichains.Rmd
+++ b/vignettes/epichains.Rmd
@@ -255,72 +255,6 @@ sim_tree_from_pop_eg <- simulate_tree_from_pop(
 head(sim_tree_from_pop_eg)
 ```
 
-#### Simulating chains with interventions
-
-All the `simulate_*()` functions can model interventions that reduce the $R_0$,
-using the `intvn_mean_reduction` argument. In general, these can be
-interpreted as population-level interventions.
-
-To illustrate this, we will use the previous examples for each function and specify
-a population-level intervention that reduces $R_0$ by $50\%$.
-
-Using `simulate_tree()`, we can specify an initial number of cases
-and a population level intervention, `intvn_mean_reduction`, that reduces $R_0$ by $50\%$.
-
-```{r}
-set.seed(123)
-# Define serial distribution
-serial_func <- function(x) {
-  return(3)
-}
-
-sim_tree_intvn_eg <- simulate_tree(
-  nchains = 10,
-  statistic = "size",
-  offspring_dist = "pois",
-  intvn_mean_reduction = 0.5,
-  stat_max = 10,
-  serials_dist = serial_func,
-  lambda = 0.9
-)
-
-head(sim_tree_intvn_eg)
-```
-
-Here is an example with `simulate_summary()`, modelling an intervention that reduces $R_0$ by $50\%$.
-```{r}
-simulate_summary_intvn_eg <- simulate_summary(
-  nchains = 10,
-  statistic = "size",
-  offspring_dist = "pois",
-  intvn_mean_reduction = 0.5,
-  stat_max = 10,
-  lambda = 0.9
-)
-
-# Print the results
-simulate_summary_intvn_eg
-```
-
-Finally, let's use `simulate_tree_from_pop()`.
-```{r}
-set.seed(7)
-# Define serial distribution
-serial_func <- function(x) {
-  return(3)
-}
-
-sim_tree_from_pop_intvn_eg <- simulate_tree_from_pop(
-  pop = 1000,
-  offspring_dist = "pois",
-  intvn_mean_reduction = 0.5,
-  lambda = 1,
-  serials_dist = serial_func
-)
-
-head(sim_tree_from_pop_intvn_eg)
-```
-
 ## Other functionalities
 
 ### Summarising

From 63f86366f0361c958df7e114e1a1072427444921 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Wed, 8 Nov 2023 15:25:40 +0000
Subject: [PATCH 2/6] add interventions vignette

---
 vignettes/interventions.Rmd | 162 ++++++++++++++++++++++++++++++++++++
 1 file changed, 162 insertions(+)
 create mode 100644 vignettes/interventions.Rmd

diff --git a/vignettes/interventions.Rmd b/vignettes/interventions.Rmd
new file mode 100644
index 00000000..f642fc5d
--- /dev/null
+++ b/vignettes/interventions.Rmd
@@ -0,0 +1,162 @@
+---
+title: "Modelling disease control interventions"
+author: "Sebastian Funk and James M. Azam"
+output:
+  bookdown::html_vignette2:
+    fig_caption: yes
+    code_folding: show
+pkgdown:
+  as_is: true
+bibliography: references.json
+link-citations: true
+vignette: >
+  %\VignetteIndexEntry{Modelling disease control interventions}
+  %\VignetteEncoding{UTF-8}
+  %\VignetteEngine{knitr::rmarkdown}
+editor_options: 
+  chunk_output_type: console
+---
+
+```{r include=FALSE}
+knitr::opts_chunk$set(
+  echo = TRUE,
+  message = FALSE,
+  warning = FALSE,
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+_epichains_ does not provide any direct functionality for studying public health interventions.
+However, the flexible simulation functionality that it includes can be used to consider some specific changes to the parameters that can be interpreted as the result of control measures.
+Here we investigate the effect on outbreak sizes, but the same approaches could be used for investigating chain lengths (using the `statistic` argument to `simulate_summary`) or the time progression of outbreaks (using the `simulate_tree` function).
+
+```{r}
+library("epichains")
+```
+
+As a base case we consider the spread of an infection with a negative binomial offspring distribution with mean 1.2 and overdispersion parameter 0.5.
+We simulate 200 chains tracking up to 99 infections:
+
+```{r simulate_chains}
+sims <- simulate_summary(
+  nchains = 200, offspring_dist = "nbinom", stat_max = 99, mu = 1.2, size = 0.5
+)
+```
+
+We then plot the resulting distribution of chain sizes
+```{r uncontrolled_chains_plot}
+library("ggplot2")
+sims[is.infinite(sims)] <- 100
+ggplot(data.frame(x = sims), aes(x = x)) +
+  geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") +
+  scale_x_continuous(breaks = c(0, 25, 50, 75, 100),
+                     labels = c(0, 25, 50, 75, ">99")) +
+  theme_bw()
+```
+
+# Reductions in the strength of transmission
+
+Following [@lloyd-smith2005] we consider two ways in which disease control interventions can reduce the reproduction number: _population-wide_ and _individual-specific_ control.
+
+## Population-wide control
+
+By population-level control we mean an intervention that reduces the mean number of offspring (i.e. the reproduction number) by a fixed proportion.
+For example, to reduce R by 25% at the population level we scale the `mu` parameter from 1.2 to 0.9:
+
+```{r simulate_chains_pop_control}
+sims <- simulate_summary(
+  nchains = 200, offspring_dist = "nbinom", stat_max = 99, mu = 0.9, size = 0.5
+)
+sims[is.infinite(sims)] <- 100
+ggplot(data.frame(x = sims), aes(x = x)) +
+  geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") +
+  scale_x_continuous(breaks = c(0, 25, 50, 75, 100),
+                     labels = c(0, 25, 50, 75, ">99")) +
+  theme_bw()
+```
+
+## Individual-level control.
+
+In simulating population-level control we now apply the same reduction as before (25%) but instead of assuming that the mean is reduced we apply this such that 25% of individuals do not transmit further at all, whereas the remaining 75% generate offspring as in the uncontrolled case.
+
+To do this, we can no longer use the standard negative binomial distribution that comes with R.
+Instead, we define a random generator from a modified negative binomial distribution that includes our individual-level control as a `control` argument indicating the level of individual-level control (0: no control; 1: full control):
+
+```{r nbinom_ind_control}
+rnbinom_ind <- function(n, ..., control = 0) {
+  ## initialise number of offspring to 0
+  offspring <- rep(0L, n)
+  ## for each individual, decide whether they transmit further
+  transmits <- rbinom(n = n, prob = 1 - control, size = 1)
+  ## check if anyone transmits further
+  if (sum(transmits) > 0L) {
+    ## for those that transmit, sample from negative binomial with given
+    ## parameters
+    offspring[which(transmits == 1L)] <- rnbinom(n = n, ...)
+  }
+  return(offspring)
+}
+```
+
+Having defined this, we can generate simulations as before:
+
+```{r simulate_chains_ind_control}
+sims <- simulate_summary(
+  nchains = 200, offspring_dist = "nbinom_ind", stat_max = 99, mu = 1.2, 
+  size = 0.5, control = 0.25
+)
+sims[is.infinite(sims)] <- 100
+ggplot(data.frame(x = sims), aes(x = x)) +
+  geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") +
+  scale_x_continuous(breaks = c(0, 25, 50, 75, 100),
+                     labels = c(0, 25, 50, 75, ">99")) +
+  theme_bw()
+```
+
+# Preventing superspreading events
+
+Another way of controlling a disease would be to prevent individuals from spreading to a large number of others, for example by preventing mass gatherings or, more generally, settings where superspreading events can occur.
+
+We can model this by truncating the offspring distribution at a certain size.
+This can be done, for example, using the [truncdist](https://cran.r-project.org/package=truncdist) R package.
+We use this to define a truncated negative binomial offspring distribution:
+
+```{r negbin_truncated}
+library("truncdist")
+rnbinom_truncated <- function(n, ..., max = Inf) {
+  return(rtrunc(n = n, spec = "nbinom", b = max, ...))
+}
+```
+
+We use this to simulate chains in a situation where the maximum of secondary cases that each infected person can generate is 10.
+
+```{r simulate_chains_truncated}
+sims <- simulate_summary(
+  nchains = 200, offspring_dist = "nbinom_truncated", stat_max = 99, mu = 1.2, 
+  size = 0.5, max =  10
+)
+sims[is.infinite(sims)] <- 100
+ggplot(data.frame(x = sims), aes(x = x)) +
+  geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") +
+  scale_x_continuous(breaks = c(0, 25, 50, 75, 100),
+                     labels = c(0, 25, 50, 75, ">99")) +
+  theme_bw()
+```
+
+# Truncating the generation interval
+
+Lastly, we consider a situation where the generation interval is shortened.
+We do not model this explicitly but instead consider the effect on the offspring distribution.
+
+For example, if our generation interval is from a gamma distribution with shape = 25 and rate = 5 (corresponding to a mean of 5 and standard deviation of 1), and we stop all transmission that would normally occur more than 6 days after infection, we can calculate the proportion of transmissions that are prevented as
+
+```{r truncate_gen_int}
+control <- 1 - pgamma(6, shape = 25, rate = 5)
+signif(control, 2)
+```
+
+# References
+
+In other words, this would prevent `r round(100 * control)`% of infections in this example.
+The value of `control` can be used in the examples above to study the effect on outbreak sizes.

From 766829a53d453043ffa1b24299a1651b1662b343 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Wed, 8 Nov 2023 15:25:53 +0000
Subject: [PATCH 3/6] update pkgdown

---
 _pkgdown.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/_pkgdown.yml b/_pkgdown.yml
index cc8cc70d..a84a41a0 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -10,6 +10,7 @@ articles:
   navbar: Package vignettes
   contents:
   - projecting_incidence
+  - interventions
 - title: Modelling guides and background
   navbar: Modelling guides and background
   contents:

From 2820bb04695dbb9596be9fbdf763ec8bd635eb57 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Wed, 8 Nov 2023 15:26:34 +0000
Subject: [PATCH 4/6] add news item

---
 NEWS.md | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/NEWS.md b/NEWS.md
index 228aab65..0a4161fe 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,8 @@
+# epichains 0.1.9999
+
+## Documentation
+* A vignette outlining how to simulate interventions has been added
+
 # epichains 0.1.0
 
 ## Package name change

From ce942e069d4c54594c30283f57b49ebd67e0562a Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Thu, 9 Nov 2023 09:55:00 +0000
Subject: [PATCH 5/6] Apply suggestions from code review

Co-authored-by: James Azam <james.m.azam@gmail.com>
---
 vignettes/interventions.Rmd | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/vignettes/interventions.Rmd b/vignettes/interventions.Rmd
index f642fc5d..0696a4ad 100644
--- a/vignettes/interventions.Rmd
+++ b/vignettes/interventions.Rmd
@@ -47,7 +47,7 @@ sims <- simulate_summary(
 We then plot the resulting distribution of chain sizes
 ```{r uncontrolled_chains_plot}
 library("ggplot2")
-sims[is.infinite(sims)] <- 100
+sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting.
 ggplot(data.frame(x = sims), aes(x = x)) +
   geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") +
   scale_x_continuous(breaks = c(0, 25, 50, 75, 100),
@@ -55,7 +55,7 @@ ggplot(data.frame(x = sims), aes(x = x)) +
   theme_bw()
 ```
 
-# Reductions in the strength of transmission
+# Reducing the strength of transmission
 
 Following [@lloyd-smith2005] we consider two ways in which disease control interventions can reduce the reproduction number: _population-wide_ and _individual-specific_ control.
 
@@ -68,7 +68,7 @@ For example, to reduce R by 25% at the population level we scale the `mu` parame
 sims <- simulate_summary(
   nchains = 200, offspring_dist = "nbinom", stat_max = 99, mu = 0.9, size = 0.5
 )
-sims[is.infinite(sims)] <- 100
+sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting.
 ggplot(data.frame(x = sims), aes(x = x)) +
   geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") +
   scale_x_continuous(breaks = c(0, 25, 50, 75, 100),
@@ -90,7 +90,7 @@ rnbinom_ind <- function(n, ..., control = 0) {
   ## for each individual, decide whether they transmit further
   transmits <- rbinom(n = n, prob = 1 - control, size = 1)
   ## check if anyone transmits further
-  if (sum(transmits) > 0L) {
+  if (any(transmits == 1L)) {
     ## for those that transmit, sample from negative binomial with given
     ## parameters
     offspring[which(transmits == 1L)] <- rnbinom(n = n, ...)
@@ -106,7 +106,7 @@ sims <- simulate_summary(
   nchains = 200, offspring_dist = "nbinom_ind", stat_max = 99, mu = 1.2, 
   size = 0.5, control = 0.25
 )
-sims[is.infinite(sims)] <- 100
+sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting.
 ggplot(data.frame(x = sims), aes(x = x)) +
   geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") +
   scale_x_continuous(breaks = c(0, 25, 50, 75, 100),
@@ -130,13 +130,14 @@ rnbinom_truncated <- function(n, ..., max = Inf) {
 ```
 
 We use this to simulate chains in a situation where the maximum of secondary cases that each infected person can generate is 10.
+This can be likened to a disease control strategy where gatherings are limited to 10 people.
 
 ```{r simulate_chains_truncated}
 sims <- simulate_summary(
   nchains = 200, offspring_dist = "nbinom_truncated", stat_max = 99, mu = 1.2, 
   size = 0.5, max =  10
 )
-sims[is.infinite(sims)] <- 100
+sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting.
 ggplot(data.frame(x = sims), aes(x = x)) +
   geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") +
   scale_x_continuous(breaks = c(0, 25, 50, 75, 100),
@@ -156,7 +157,8 @@ control <- 1 - pgamma(6, shape = 25, rate = 5)
 signif(control, 2)
 ```
 
-# References
 
 In other words, this would prevent `r round(100 * control)`% of infections in this example.
 The value of `control` can be used in the examples above to study the effect on outbreak sizes.
+
+# References

From e614d9890ba8e8361c03387930294a0260279ff3 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Thu, 9 Nov 2023 12:46:10 +0000
Subject: [PATCH 6/6] make lintr happy

---
 vignettes/interventions.Rmd | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/vignettes/interventions.Rmd b/vignettes/interventions.Rmd
index 0696a4ad..2176283f 100644
--- a/vignettes/interventions.Rmd
+++ b/vignettes/interventions.Rmd
@@ -31,8 +31,13 @@ _epichains_ does not provide any direct functionality for studying public health
 However, the flexible simulation functionality that it includes can be used to consider some specific changes to the parameters that can be interpreted as the result of control measures.
 Here we investigate the effect on outbreak sizes, but the same approaches could be used for investigating chain lengths (using the `statistic` argument to `simulate_summary`) or the time progression of outbreaks (using the `simulate_tree` function).
 
-```{r}
+```{r load_libraries}
+## main package
 library("epichains")
+## for plotting
+library("ggplot2")
+## for truncating the offspring distribution later
+library("truncdist")
 ```
 
 As a base case we consider the spread of an infection with a negative binomial offspring distribution with mean 1.2 and overdispersion parameter 0.5.
@@ -46,7 +51,6 @@ sims <- simulate_summary(
 
 We then plot the resulting distribution of chain sizes
 ```{r uncontrolled_chains_plot}
-library("ggplot2")
 sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting.
 ggplot(data.frame(x = sims), aes(x = x)) +
   geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") +
@@ -103,7 +107,7 @@ Having defined this, we can generate simulations as before:
 
 ```{r simulate_chains_ind_control}
 sims <- simulate_summary(
-  nchains = 200, offspring_dist = "nbinom_ind", stat_max = 99, mu = 1.2, 
+  nchains = 200, offspring_dist = "nbinom_ind", stat_max = 99, mu = 1.2,
   size = 0.5, control = 0.25
 )
 sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting.
@@ -123,7 +127,6 @@ This can be done, for example, using the [truncdist](https://cran.r-project.org/
 We use this to define a truncated negative binomial offspring distribution:
 
 ```{r negbin_truncated}
-library("truncdist")
 rnbinom_truncated <- function(n, ..., max = Inf) {
   return(rtrunc(n = n, spec = "nbinom", b = max, ...))
 }
@@ -134,7 +137,7 @@ This can be likened to a disease control strategy where gatherings are limited t
 
 ```{r simulate_chains_truncated}
 sims <- simulate_summary(
-  nchains = 200, offspring_dist = "nbinom_truncated", stat_max = 99, mu = 1.2, 
+  nchains = 200, offspring_dist = "nbinom_truncated", stat_max = 99, mu = 1.2,
   size = 0.5, max =  10
 )
 sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting.