Add epidemic risk vignette and probability_contain()

DESCRIPTION

@@ -34,7 +34,8 @@ Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.2.3
 Imports:
     checkmate,
-    stats
+    stats,
+    bpmodels
 Suggests: 
     epiparameter,
     distributional,
@@ -43,10 +44,12 @@ Suggests:
     bookdown,
     rmarkdown,
     ggplot2,
-    testthat (>= 3.0.0),
-    spelling
+    spelling,
+    ggtext,
+    testthat (>= 3.0.0)
 Remotes:
-    epiverse-trace/epiparameter
+    epiverse-trace/epiparameter,
+    epiverse-trace/bpmodels
 Config/testthat/edition: 3
 Config/Needs/website:
     epiverse-trace/epiversetheme

NAMESPACE

@@ -4,6 +4,7 @@ export(dpoislnorm)
 export(dpoisweibull)
 export(ppoislnorm)
 export(ppoisweibull)
+export(probability_contain)
 export(probability_epidemic)
 export(probability_extinct)
 export(proportion_cluster_size)

R/probability_contain.R

@@ -0,0 +1,68 @@
+#' Probability that an outbreak will be contained
+#'
+#' @description Containment is defined as the size of the transmission chain
+#' not reaching the `case_threshold` (default = 100).
+#'
+#' @inheritParams probability_epidemic
+#' @param c Control strength, 0 is no control measures, 1 is complete control.
+#' @param control_type Either `"population"` or `"individual"` for
+#' population-level or individual-level control measures.
+#' @param stochastic Whether to use a stochastic branching process model or the
+#' probability of extinction.
+#' @param ... arguments to be passed to [bpmodels::chain_sim()].
+#' @param case_threshold A number for the threshold of the number of cases below
+#' which the epidemic is considered contained.
+#'
+#' @return A number for the probability of containment
+#' @export
+#'
+#' @references
+#'
+#' Lloyd-Smith, J. O., Schreiber, S. J., Kopp, P. E., & Getz, W. M. (2005)
+#' Superspreading and the effect of individual variation on disease emergence.
+#' Nature, 438(7066), 355-359. <https://doi.org/10.1038/nature04153>
+#'
+#' @examples
+#' probability_contain(R = 1.5, k = 0.5, c = 1)
+probability_contain <- function(R, k, a = 1, c, # nolint
+                                control_type = c("population", "individual"),
+                                stochastic = TRUE,
+                                ...,
+                                case_threshold = 100) {
+  # check inputs
+  checkmate::assertNumber(R)
+  checkmate::assertNumber(k)
+  checkmate::assertCount(a)
+  checkmate::assertNumber(c, lower = 0, upper = 1)
+  checkmate::assertLogical(stochastic)
+  checkmate::assertNumber(case_threshold)
+
+  control_type <- match.arg(control_type)
+  if (control_type == "population") {
+    R <- (1 - c) * R # nolint
+  } else {
+    stop("individual-level controls not yet implemented", call. = FALSE)
+  }
+
+  if (a != 1) {
+    stop(
+      "Multiple introductions is not yet implemented for probability_contain",
+      call. = FALSE
+    )
+  }
+
+  if (stochastic) {
+    chain_size <- bpmodels::chain_sim(
+      n = 1e5,
+      offspring = "nbinom",
+      size = k,
+      mu = R,
+      infinite = case_threshold,
+      ...
+    )
+    prob_contain <- sum(!is.infinite(chain_size)) / length(chain_size)
+  } else {
+    prob_contain <- probability_extinct(R = R, k = k, a = a)
+  }
+  return(prob_contain)
+}

_pkgdown.yml

@@ -17,8 +17,8 @@ reference:
     contents:
     - ends_with("poisweibull")
 
-  - title: Probability of epidemic or extinction
-    desc: Probability of a disease causing an epidemic or going extinct
+  - title: Probability of epidemic, extinction or containment
+    desc: Probability of a disease causing an epidemic, going extinct or being contained
     contents:
     - starts_with("probability")
 
@@ -32,3 +32,7 @@ articles:
   navbar: Parameter estimation
   contents:
   - estimate_individual_level_transmission
+- title: Superspreading for decision-making
+  navbar: Superspreading for decision-making
+  contents:
+  - epidemic_risk

inst/WORDLIST

@@ -1,6 +1,8 @@
 Althaus
+aes
 Barré
 Boëlle
+bw
 Camacho
 CMD
 COVID
@@ -18,16 +20,27 @@ Getz
 Heleze
 Kopp
 Lifecycle
+linetype
 Liu
 Loucoubar
+lw
 Magassouba
+MERS
 N’Faly
 org
 Oumar
 Ousmane
+params
+pointrange
+prob
 Schreiber
 Soropogui
+viridis
+vline
 Wellcome
+xintercept
+ymax
+ymin
 al
 doi
 epiparameter

tests/testthat/test-probability_contain.R

@@ -0,0 +1,34 @@
+test_that("probability_contain works as expected", {
+  prob_contain <- probability_contain(R = 1.5, k = 0.5, a = 1, c = 0)
+  # larger tolerance for stochastic variance
+  expect_equal(prob_contain, 0.7672, tolerance = 1e-2)
+})
+
+test_that("probability_contain works as expected for deterministic", {
+  prob_contain <- probability_contain(
+    R = 1.5, k = 0.5, a = 1, c = 0, stochastic = FALSE
+  )
+  expect_equal(prob_contain, 0.768)
+})
+
+test_that("probability_contain works as expected for different threshold", {
+  prob_contain <- probability_contain(
+    R = 1.5, k = 0.5, a = 1, c = 0, case_threshold = 50
+  )
+  # larger tolerance for stochastic variance
+  expect_equal(prob_contain, 0.76255, tolerance = 1e-2)
+})
+
+test_that("probability_contain fails as expected", {
+  expect_error(
+    probability_contain(R = 1, k = 1, a = 2, c = 1),
+    regexp = "(Multiple introductions is not yet implemented)"
+  )
+
+  expect_error(
+    probability_contain(
+      R = 1, k = 1, a = 2, c = 1, control_type = "individual"
+    ),
+    regexp = "individual-level controls not yet implemented"
+  )
+})