diff --git a/R/borel.r b/R/borel.r
index 1a5ac8e3..f58dabd1 100644
--- a/R/borel.r
+++ b/R/borel.r
@@ -27,20 +27,19 @@ dborel <- function(x, mu, log = FALSE) {
 ##' Random numbers are generated by simulating from a Poisson branching process
 ##' @param n Number of random variates to generate.
 ##' @inheritParams dborel
-##' @param infinite A threshold above which the random number is set to
-##'  `Inf` as a form of truncation; A number coercible to integer. This
-##'  argument is important because `rborel()` uses `simulate_summary()`
-##'  internally with `offspring_dist = rpois`, so when `mu >= 1`, the
-##'  simulation could grow out of proportion and chain summaries (the random
-##'  numbers generated here) need to be truncated using
-##'  `infinite`.
+##' @param censor_at A stopping criterion; `<numeric>`. Defaults to `Inf`. A
+##' value above which the simulation ends and the random number is set to
+##' `Inf` (as a form of censoring). `rborel()` simulates chain sizes using
+##' [simulate_chain_stats()] with a Poisson offspring distribution, so if
+##' `mu >= 1`, the simulation could proceed unendingly. This parameter is used
+##' to prevent this.
 ##' @return A numeric vector of random numbers.
 ##' @author Sebastian Funk
 ##' @export
 ##' @examples
 ##' set.seed(32)
 ##' rborel(5, 1)
-rborel <- function(n, mu, infinite = Inf) {
+rborel <- function(n, mu, censor_at = Inf) {
   checkmate::assert_number(
     n, lower = 1, finite = TRUE, na.ok = FALSE
   )
@@ -52,7 +51,7 @@ rborel <- function(n, mu, infinite = Inf) {
     n_chains = n,
     offspring_dist = rpois,
     statistic = "size",
-    stat_max = infinite,
+    stat_threshold = censor_at,
     lambda = mu
   )
   out <- as.numeric(out)
@@ -72,7 +71,7 @@ rborel <- function(n, mu, infinite = Inf) {
 ##' @return Numeric vector of random numbers
 ##' @author Sebastian Funk
 ##' @export
-rgborel <- function(n, size, prob, mu, infinite = Inf) {
+rgborel <- function(n, size, prob, mu, censor_at = Inf) {
   ## This function was introduced to support estimating likelihoods using a
   ## Gamma-Borel mixture distribution. It is not actually called (it only needs)
   ## to exist and could be a dummy. However, the function is here included with
@@ -98,5 +97,7 @@ rgborel <- function(n, size, prob, mu, infinite = Inf) {
   ## first, sample from gamma
   x <- rgamma(n, shape = size, rate = size / mu)
   ## then, sample from borel
-  return(vapply(x, rborel, n = 1, infinite = infinite, FUN.VALUE = numeric(1)))
+  return(vapply(
+    x, rborel, n = 1, censor_at = censor_at, FUN.VALUE = numeric(1)
+  ))
 }
diff --git a/R/checks.R b/R/checks.R
index 7a2c08cb..539c4c10 100644
--- a/R/checks.R
+++ b/R/checks.R
@@ -28,7 +28,7 @@
     n_chains,
     statistic,
     offspring_dist,
-    stat_max,
+    stat_threshold,
     pop,
     percent_immune) {
   # Input checking
@@ -41,7 +41,7 @@
   # check that arguments related to the statistic are valid
   .check_statistic_args(
     statistic,
-    stat_max
+    stat_threshold
   )
   checkmate::assert(
     is.infinite(pop) ||
@@ -54,7 +54,7 @@
   invisible(NULL)
 }
 
-#' Check that the `statistic` and `stat_max` arguments are valid
+#' Check that the `statistic` and `stat_threshold` arguments are valid
 #'
 #' @inheritParams simulate_chains
 #' @description
@@ -64,15 +64,15 @@
 #' @return NULL; called for side effects
 #' @keywords internal
 .check_statistic_args <- function(statistic,
-                                  stat_max) {
+                                  stat_threshold) {
   checkmate::assert_choice(
     statistic,
     choices = c("size", "length")
   )
   checkmate::assert(
-    is.infinite(stat_max),
+    is.infinite(stat_threshold),
       checkmate::check_integerish(
-        stat_max,
+        stat_threshold,
         lower = 1,
         null.ok = FALSE
     ),
diff --git a/R/epichains.R b/R/epichains.R
index 58031a7e..f043b81b 100644
--- a/R/epichains.R
+++ b/R/epichains.R
@@ -19,7 +19,7 @@
                                n_chains,
                                statistic,
                                offspring_dist,
-                               stat_max,
+                               stat_threshold,
                                track_pop) {
   # Assemble the elements of the object
   obj <- sim_df
@@ -27,7 +27,7 @@
   attr(obj, "n_chains") <- n_chains
   attr(obj, "statistic") <- statistic
   attr(obj, "offspring_dist") <- offspring_dist
-  attr(obj, "stat_max") <- stat_max
+  attr(obj, "stat_threshold") <- stat_threshold
   attr(obj, "track_pop") <- track_pop
   return(obj)
 }
@@ -59,7 +59,7 @@
                       offspring_dist,
                       track_pop,
                       statistic = c("size", "length"),
-                      stat_max = Inf) {
+                      stat_threshold = Inf) {
   # Check that inputs are well specified
   checkmate::assert_data_frame(sim_df, min.cols = 3, min.rows = n_chains)
   checkmate::assert_integerish(
@@ -74,8 +74,8 @@
   .check_offspring_func_valid(offspring_dist)
   checkmate::assert_logical(track_pop, len = 1L)
   checkmate::assert(
-    is.infinite(stat_max),
-    checkmate::check_integerish(stat_max, lower = 1L)
+    is.infinite(stat_threshold),
+    checkmate::check_integerish(stat_threshold, lower = 1L)
   )
   # Create <epichains> object
   epichains <- .new_epichains(
@@ -83,7 +83,7 @@
     n_chains = n_chains,
     statistic = statistic,
     offspring_dist = offspring_dist,
-    stat_max = stat_max,
+    stat_threshold = stat_threshold,
     track_pop = track_pop
   )
 
@@ -107,22 +107,22 @@
 #' new `<epichains_summary>` object safely, use `epichains_summary()`.
 #'
 #' @param chains_summary A numeric `<vector>` of chain sizes and lengths.
+#' @inheritParams simulate_chain_stats
 #' @inheritParams .new_epichains
-#' @inheritParams simulate_chains
 #' @author James M. Azam
 #' @keywords internal
 .new_epichains_summary <- function(chains_summary,
                                   n_chains,
                                   statistic,
                                   offspring_dist,
-                                  stat_max) {
+                                  stat_threshold) {
   # Assemble the elements of the object
   obj <- chains_summary
   class(obj) <- c("epichains_summary", class(chains_summary))
   attr(obj, "n_chains") <- n_chains
   attr(obj, "statistic") <- statistic
   attr(obj, "offspring_dist") <- offspring_dist
-  attr(obj, "stat_max") <- stat_max
+  attr(obj, "stat_threshold") <- stat_threshold
   return(obj)
 }
 
@@ -144,7 +144,7 @@
                               n_chains,
                               offspring_dist,
                               statistic = c("size", "length"),
-                              stat_max = Inf) {
+                              stat_threshold = Inf) {
   # chain_summary can sometimes contain infinite values, so check
   # that finite elements are integerish.
   checkmate::check_integerish(
@@ -162,8 +162,8 @@
   statistic <- match.arg(statistic, c("size", "length"))
   .check_offspring_func_valid(offspring_dist)
   checkmate::assert(
-    is.infinite(stat_max),
-    checkmate::check_integerish(stat_max, lower = 1L)
+    is.infinite(stat_threshold),
+    checkmate::check_integerish(stat_threshold, lower = 1L)
   )
 
   # Create <epichains_summary> object
@@ -172,7 +172,7 @@
     n_chains = n_chains,
     statistic = statistic,
     offspring_dist = offspring_dist,
-    stat_max = stat_max
+    stat_threshold = stat_threshold
   )
 
   # Validate the created object
@@ -201,8 +201,8 @@ print.epichains <- function(x, ...) {
 #' prints the number of chains simulated, and the range of
 #' the statistic, represented as the maximum (`max_stat`) and minimum
 #' (`min_stat`). If the minimum or maximum is infinite, it is represented as
-#' `>= stat_max` where `stat_max` is the value of the censoring limit. See
-#' `?epichains_summary()` for the definition of `stat_max`.
+#' `>= stat_threshold` where `stat_threshold` is the value of the censoring
+#' limit. See `?epichains_summary()` for the definition of `stat_threshold`.
 #' @param ... Not used.
 #' @return Invisibly returns an `<epichains_summary>`. Called for
 #' side-effects.
@@ -296,7 +296,7 @@ format.epichains_summary <- function(x, ...) {
         ifelse(
           is.infinite(
             statistics[["max_stat"]]),
-            paste0(">=", attr(x, "stat_max")
+            paste0(">=", attr(x, "stat_threshold")
           ),
           statistics[["max_stat"]]
         )
@@ -306,7 +306,7 @@ format.epichains_summary <- function(x, ...) {
         ifelse(
           is.infinite(
             statistics[["min_stat"]]),
-          paste0(">=", attr(x, "stat_max")
+          paste0(">=", attr(x, "stat_threshold")
           ),
           statistics[["min_stat"]]
         )
@@ -344,7 +344,7 @@ format.epichains_summary <- function(x, ...) {
 #'   percent_immune = 0,
 #'   statistic = "size",
 #'   offspring_dist = rnbinom,
-#'   stat_max = 10,
+#'   stat_threshold = 10,
 #'   mu = 2,
 #'   size = 0.2
 #' )
@@ -360,7 +360,7 @@ format.epichains_summary <- function(x, ...) {
 #'   percent_immune = 0,
 #'   statistic = "size",
 #'   offspring_dist = rnbinom,
-#'   stat_max = 10,
+#'   stat_threshold = 10,
 #'   mu = 2,
 #'   size = 0.2
 #' )
@@ -392,11 +392,11 @@ summary.epichains <- function(object, ...) {
     }
   }
   # Get other required attributes from passed object
-  stat_max <- attr(object, "stat_max")
+  stat_threshold <- attr(object, "stat_threshold")
   offspring_dist <- attr(object, "offspring_dist")
 
   # Apply truncation
-  chain_summaries[chain_summaries >= stat_max] <- Inf
+  chain_summaries[chain_summaries >= stat_threshold] <- Inf
 
   # Return an <epichains_summary> object
   chain_summaries <- .epichains_summary(
@@ -404,7 +404,7 @@ summary.epichains <- function(object, ...) {
     n_chains = n_chains,
     statistic = statistic,
     offspring_dist = offspring_dist,
-    stat_max = stat_max
+    stat_threshold = stat_threshold
   )
   return(chain_summaries)
 }
@@ -536,7 +536,7 @@ summary.epichains_summary <- function(object, ...) {
 #'   n_chains = 10,
 #'   statistic = "size",
 #'   offspring_dist = rpois,
-#'   stat_max = 10,
+#'   stat_threshold = 10,
 #'   generation_time = function(n) rep(3, n),
 #'   lambda = 2
 #' )
@@ -557,7 +557,7 @@ head.epichains <- function(x, ...) {
 #'   n_chains = 10,
 #'   statistic = "size",
 #'   offspring_dist = rpois,
-#'   stat_max = 10,
+#'   stat_threshold = 10,
 #'   generation_time = function(n) rep(3, n),
 #'   lambda = 2
 #' )
@@ -590,7 +590,7 @@ tail.epichains <- function(x, ...) {
 #'   n_chains = 10,
 #'   statistic = "size",
 #'   offspring_dist = rpois,
-#'   stat_max = 10,
+#'   stat_threshold = 10,
 #'   generation_time = function(n) rep(3, n),
 #'   lambda = 2
 #' )
diff --git a/R/likelihood.R b/R/likelihood.R
index d4e1ea33..9b6fbd74 100644
--- a/R/likelihood.R
+++ b/R/likelihood.R
@@ -46,7 +46,7 @@
 #' )
 #' @export
 likelihood <- function(chains, statistic = c("size", "length"), offspring_dist,
-                       nsim_obs, obs_prob = 1, log = TRUE, stat_max = Inf,
+                       nsim_obs, obs_prob = 1, log = TRUE, stat_threshold = Inf,
                        exclude = NULL, individual = FALSE, ...) {
   statistic <- match.arg(statistic)
 
@@ -63,7 +63,7 @@ likelihood <- function(chains, statistic = c("size", "length"), offspring_dist,
   # check that arguments related to the statistic are valid
   .check_statistic_args(
     statistic,
-    stat_max
+    stat_threshold
   )
   .check_offspring_func_valid(offspring_dist)
   checkmate::assert_number(
@@ -90,22 +90,22 @@ likelihood <- function(chains, statistic = c("size", "length"), offspring_dist,
         length(chains),
         chains, obs_prob
       ),
-      stat_max
+      stat_threshold
     ), simplify = FALSE)
     stat_rep_vect <- unlist(stat_rep_list)
-    if (!is.finite(stat_max)) {
-      stat_max <- max(stat_rep_vect) + 1
+    if (!is.finite(stat_threshold)) {
+      stat_threshold <- max(stat_rep_vect) + 1
     }
   } else {
-    chains[chains >= stat_max] <- stat_max
+    chains[chains >= stat_threshold] <- stat_threshold
     stat_rep_vect <- chains
     stat_rep_list <- list(chains)
   }
 
   ## determine for which sizes to calculate the log-likelihood
   ## (for true chain size)
-  if (any(stat_rep_vect == stat_max)) {
-    calc_sizes <- seq_len(stat_max - 1)
+  if (any(stat_rep_vect == stat_threshold)) {
+    calc_sizes <- seq_len(stat_threshold - 1)
   } else {
     calc_sizes <- unique(c(stat_rep_vect, exclude))
   }
@@ -136,16 +136,16 @@ likelihood <- function(chains, statistic = c("size", "length"), offspring_dist,
             x = calc_sizes,
             offspring_dist = offspring_dist,
             statistic = statistic,
-            stat_max = stat_max
+            stat_threshold = stat_threshold
           ),
           pars
         )
       )
   }
 
-  ## assign probabilities to stat_max outbreak sizes
-  if (any(stat_rep_vect == stat_max)) {
-    likelihoods[stat_max] <- .complementary_logprob(likelihoods)
+  ## assign probabilities to stat_threshold outbreak sizes
+  if (any(stat_rep_vect == stat_threshold)) {
+    likelihoods[stat_threshold] <- .complementary_logprob(likelihoods)
   }
 
   if (!missing(exclude)) {
diff --git a/R/simulate.r b/R/simulate.r
index 512882c0..1a5ed1b7 100644
--- a/R/simulate.r
+++ b/R/simulate.r
@@ -22,7 +22,7 @@
 #' here are `rpois` for Poisson distributed offspring, `rnbinom` for negative
 #' binomial offspring, or custom functions.
 #' @param statistic The chain statistic to track as the
-#' stopping criteria for each chain being simulated when `stat_max` is not
+#' stopping criteria for each chain being simulated when `stat_threshold` is not
 #' `Inf`; A `<string>`. It can be one of:
 #' \itemize{
 #'   \item "size": the total number of cases produced by a chain before it goes
@@ -30,9 +30,13 @@
 #'   \item "length": the total number of generations reached by a chain before
 #'   it goes extinct.
 #' }
-#' @param stat_max The cut off for the chain statistic (size/length) being
-#' computed; A number coercible to integer. Results above `stat_max` are set to
-#' `stat_max`; An `<integer>`. Defaults to `Inf`.
+#' @param stat_threshold A stopping criterion for individual chain simulations;
+#' a positive number coercible to integer. When any chain's cumulative statistic
+#' reaches or surpasses `stat_threshold`, that chain ends. Defaults to `Inf`.
+#' For example, if `statistic = "size"` and `stat_threshold = 10`, then any
+#' chain that produces 10 or more cases will stop. Note that setting
+#' `stat_threshold` does not guarantee that all chains will stop at the same
+#' value.
 #' @param pop Population size; An `<Integer>`. Used alongside `percent_immune`
 #' to define the susceptible population. Defaults to `Inf`.
 #' @param percent_immune Percent of the population immune to
@@ -111,7 +115,7 @@
 #'   n_chains = 10,
 #'   statistic = "size",
 #'   offspring_dist = rpois,
-#'   stat_max = 10,
+#'   stat_threshold = 10,
 #'   generation_time = function(n) rep(3, n),
 #'   lambda = 2
 #' )
@@ -126,7 +130,7 @@
 #'   percent_immune = 0,
 #'   statistic = "size",
 #'   offspring_dist = rnbinom,
-#'   stat_max = 10,
+#'   stat_threshold = 10,
 #'   generation_time = function(n) rep(3, n),
 #'   mu = 2,
 #'   size = 0.2
@@ -150,7 +154,7 @@ simulate_chains <- function(n_chains,
                             statistic = c("size", "length"),
                             offspring_dist,
                             ...,
-                            stat_max = Inf,
+                            stat_threshold = Inf,
                             pop = Inf,
                             percent_immune = 0,
                             generation_time = NULL,
@@ -161,7 +165,7 @@ simulate_chains <- function(n_chains,
     n_chains = n_chains,
     statistic = statistic,
     offspring_dist = offspring_dist,
-    stat_max = stat_max,
+    stat_threshold = stat_threshold,
     pop = pop,
     percent_immune = percent_immune
   )
@@ -275,7 +279,7 @@ simulate_chains <- function(n_chains,
 
     # Find ongoing chains (those still infecting): those that have still
     # offspring and haven't reached the threshold statistic yet
-    chains_active <- which(n_offspring > 0 & stat_track < stat_max)
+    chains_active <- which(n_offspring > 0 & stat_track < stat_threshold)
     if (length(chains_active) > 0) {
       if (!missing(generation_time)) {
         ## only continue to simulate trees that don't go beyond tf
@@ -315,7 +319,7 @@ simulate_chains <- function(n_chains,
     n_chains = n_chains,
     statistic = statistic,
     offspring_dist = offspring_dist,
-    stat_max = stat_max,
+    stat_threshold = stat_threshold,
     track_pop = !missing(pop)
   )
   return(out)
@@ -336,9 +340,11 @@ simulate_chains <- function(n_chains,
 #' accepts population related inputs such as the population size (defaults
 #' to Inf) and percentage of the population initially immune (defaults to 0).
 #' @inheritParams simulate_chains
-#' @param stat_max A cut off for the chain statistic (size/length) being
-#' computed. A number coercible to integer. Results above the specified value,
-#' are set to `Inf`.
+#' @param stat_threshold A stopping criterion for individual chain simulations;
+#' a positive number coercible to integer. When any chain's cumulative statistic
+#' reaches or surpasses `stat_threshold`, that chain ends. It also serves as a
+#' censoring limit so that results above the specified value, are set to `Inf`.
+#' Defaults to `Inf`.
 #' @return An object of class `<epichains_summary>`, which is a numeric
 #' vector of chain sizes or lengths with extra attributes for storing the
 #' simulation parameters.
@@ -373,7 +379,7 @@ simulate_chains <- function(n_chains,
 #'   n_chains = 20,
 #'   statistic = "size",
 #'   offspring_dist = rpois,
-#'   stat_max = 10,
+#'   stat_threshold = 10,
 #'   lambda = 0.9
 #' )
 #' # Simulate chain sizes with a negative binomial distribution and assuming
@@ -385,7 +391,7 @@ simulate_chains <- function(n_chains,
 #'   n_chains = 20,
 #'   statistic = "size",
 #'   offspring_dist = rnbinom,
-#'   stat_max = 10,
+#'   stat_threshold = 10,
 #'   mu = 0.9,
 #'   size = 0.36
 #' )
@@ -394,7 +400,7 @@ simulate_chain_stats <- function(n_chains,
                              statistic = c("size", "length"),
                              offspring_dist,
                              ...,
-                             stat_max = Inf,
+                             stat_threshold = Inf,
                              pop = Inf,
                              percent_immune = 0) {
   # Check offspring and population-related arguments
@@ -402,7 +408,7 @@ simulate_chain_stats <- function(n_chains,
     n_chains = n_chains,
     statistic = statistic,
     offspring_dist = offspring_dist,
-    stat_max = stat_max,
+    stat_threshold = stat_threshold,
     pop = pop,
     percent_immune = percent_immune
   )
@@ -456,18 +462,19 @@ simulate_chain_stats <- function(n_chains,
     )
     # Update susceptible population
     susc_pop <- susc_pop - sum(n_offspring)
-    # only continue chains that have offspring and haven't reached stat_max
-    chains_active <- which(n_offspring > 0 & stat_track < stat_max)
+    # only continue chains that have offspring and haven't reached
+    # stat_threshold
+    chains_active <- which(n_offspring > 0 & stat_track < stat_threshold)
   }
 
-  stat_track[stat_track >= stat_max] <- Inf
+  stat_track[stat_track >= stat_threshold] <- Inf
 
   out <- .epichains_summary(
     chains_summary = stat_track,
     n_chains = n_chains,
     statistic = statistic,
     offspring_dist = offspring_dist,
-    stat_max = stat_max
+    stat_threshold = stat_threshold
   )
 
   return(out)
diff --git a/README.Rmd b/README.Rmd
index 3cd76bf4..daafcf46 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -104,7 +104,7 @@ sim_chains <- simulate_chains(
   n_chains = 20,
   statistic = "size",
   offspring_dist = rpois,
-  stat_max = 25,
+  stat_threshold = 25,
   generation_time = function(n) {rep(3, n)}, # constant generation time of 3
   lambda = 1 # mean of the Poisson distribution
 )
diff --git a/README.md b/README.md
index 55b32ba2..ac7ce70f 100644
--- a/README.md
+++ b/README.md
@@ -105,7 +105,7 @@ sim_chains <- simulate_chains(
   n_chains = 20,
   statistic = "size",
   offspring_dist = rpois,
-  stat_max = 25,
+  stat_threshold = 25,
   generation_time = function(n) {rep(3, n)}, # constant generation time of 3
   lambda = 1 # mean of the Poisson distribution
 )
diff --git a/man/aggregate.epichains.Rd b/man/aggregate.epichains.Rd
index 241b2900..d731f179 100644
--- a/man/aggregate.epichains.Rd
+++ b/man/aggregate.epichains.Rd
@@ -29,7 +29,7 @@ chains <- simulate_chains(
   n_chains = 10,
   statistic = "size",
   offspring_dist = rpois,
-  stat_max = 10,
+  stat_threshold = 10,
   generation_time = function(n) rep(3, n),
   lambda = 2
 )
diff --git a/man/dot-check_sim_args.Rd b/man/dot-check_sim_args.Rd
index 7d48ece3..de567a9b 100644
--- a/man/dot-check_sim_args.Rd
+++ b/man/dot-check_sim_args.Rd
@@ -8,7 +8,7 @@
   n_chains,
   statistic,
   offspring_dist,
-  stat_max,
+  stat_threshold,
   pop,
   percent_immune
 )
@@ -17,7 +17,7 @@
 \item{n_chains}{Number of chains to simulate.}
 
 \item{statistic}{The chain statistic to track as the
-stopping criteria for each chain being simulated when \code{stat_max} is not
+stopping criteria for each chain being simulated when \code{stat_threshold} is not
 \code{Inf}; A \verb{<string>}. It can be one of:
 \itemize{
 \item "size": the total number of cases produced by a chain before it goes
@@ -35,9 +35,13 @@ on to the random number generating functions. Examples that can be provided
 here are \code{rpois} for Poisson distributed offspring, \code{rnbinom} for negative
 binomial offspring, or custom functions.}
 
-\item{stat_max}{The cut off for the chain statistic (size/length) being
-computed; A number coercible to integer. Results above \code{stat_max} are set to
-\code{stat_max}; An \verb{<integer>}. Defaults to \code{Inf}.}
+\item{stat_threshold}{A stopping criterion for individual chain simulations;
+a positive number coercible to integer. When any chain's cumulative statistic
+reaches or surpasses \code{stat_threshold}, that chain ends. Defaults to \code{Inf}.
+For example, if \code{statistic = "size"} and \code{stat_threshold = 10}, then any
+chain that produces 10 or more cases will stop. Note that setting
+\code{stat_threshold} does not guarantee that all chains will stop at the same
+value.}
 
 \item{pop}{Population size; An \verb{<Integer>}. Used alongside \code{percent_immune}
 to define the susceptible population. Defaults to \code{Inf}.}
diff --git a/man/dot-check_statistic_args.Rd b/man/dot-check_statistic_args.Rd
index c54f8b08..2e3570ea 100644
--- a/man/dot-check_statistic_args.Rd
+++ b/man/dot-check_statistic_args.Rd
@@ -2,13 +2,13 @@
 % Please edit documentation in R/checks.R
 \name{.check_statistic_args}
 \alias{.check_statistic_args}
-\title{Check that the \code{statistic} and \code{stat_max} arguments are valid}
+\title{Check that the \code{statistic} and \code{stat_threshold} arguments are valid}
 \usage{
-.check_statistic_args(statistic, stat_max)
+.check_statistic_args(statistic, stat_threshold)
 }
 \arguments{
 \item{statistic}{The chain statistic to track as the
-stopping criteria for each chain being simulated when \code{stat_max} is not
+stopping criteria for each chain being simulated when \code{stat_threshold} is not
 \code{Inf}; A \verb{<string>}. It can be one of:
 \itemize{
 \item "size": the total number of cases produced by a chain before it goes
@@ -17,9 +17,13 @@ extinct.
 it goes extinct.
 }}
 
-\item{stat_max}{The cut off for the chain statistic (size/length) being
-computed; A number coercible to integer. Results above \code{stat_max} are set to
-\code{stat_max}; An \verb{<integer>}. Defaults to \code{Inf}.}
+\item{stat_threshold}{A stopping criterion for individual chain simulations;
+a positive number coercible to integer. When any chain's cumulative statistic
+reaches or surpasses \code{stat_threshold}, that chain ends. Defaults to \code{Inf}.
+For example, if \code{statistic = "size"} and \code{stat_threshold = 10}, then any
+chain that produces 10 or more cases will stop. Note that setting
+\code{stat_threshold} does not guarantee that all chains will stop at the same
+value.}
 }
 \value{
 NULL; called for side effects
diff --git a/man/dot-epichains.Rd b/man/dot-epichains.Rd
index 51b418f1..6f907dee 100644
--- a/man/dot-epichains.Rd
+++ b/man/dot-epichains.Rd
@@ -10,7 +10,7 @@
   offspring_dist,
   track_pop,
   statistic = c("size", "length"),
-  stat_max = Inf
+  stat_threshold = Inf
 )
 }
 \arguments{
@@ -32,7 +32,7 @@ binomial offspring, or custom functions.}
 \item{track_pop}{Was the susceptible population tracked? Logical.}
 
 \item{statistic}{The chain statistic to track as the
-stopping criteria for each chain being simulated when \code{stat_max} is not
+stopping criteria for each chain being simulated when \code{stat_threshold} is not
 \code{Inf}; A \verb{<string>}. It can be one of:
 \itemize{
 \item "size": the total number of cases produced by a chain before it goes
@@ -41,9 +41,13 @@ extinct.
 it goes extinct.
 }}
 
-\item{stat_max}{The cut off for the chain statistic (size/length) being
-computed; A number coercible to integer. Results above \code{stat_max} are set to
-\code{stat_max}; An \verb{<integer>}. Defaults to \code{Inf}.}
+\item{stat_threshold}{A stopping criterion for individual chain simulations;
+a positive number coercible to integer. When any chain's cumulative statistic
+reaches or surpasses \code{stat_threshold}, that chain ends. Defaults to \code{Inf}.
+For example, if \code{statistic = "size"} and \code{stat_threshold = 10}, then any
+chain that produces 10 or more cases will stop. Note that setting
+\code{stat_threshold} does not guarantee that all chains will stop at the same
+value.}
 }
 \value{
 An \verb{<epichains>} object.
diff --git a/man/dot-epichains_summary.Rd b/man/dot-epichains_summary.Rd
index f8ce0d0d..e103182d 100644
--- a/man/dot-epichains_summary.Rd
+++ b/man/dot-epichains_summary.Rd
@@ -9,7 +9,7 @@
   n_chains,
   offspring_dist,
   statistic = c("size", "length"),
-  stat_max = Inf
+  stat_threshold = Inf
 )
 }
 \arguments{
@@ -27,7 +27,7 @@ here are \code{rpois} for Poisson distributed offspring, \code{rnbinom} for nega
 binomial offspring, or custom functions.}
 
 \item{statistic}{The chain statistic to track as the
-stopping criteria for each chain being simulated when \code{stat_max} is not
+stopping criteria for each chain being simulated when \code{stat_threshold} is not
 \code{Inf}; A \verb{<string>}. It can be one of:
 \itemize{
 \item "size": the total number of cases produced by a chain before it goes
@@ -36,9 +36,11 @@ extinct.
 it goes extinct.
 }}
 
-\item{stat_max}{The cut off for the chain statistic (size/length) being
-computed; A number coercible to integer. Results above \code{stat_max} are set to
-\code{stat_max}; An \verb{<integer>}. Defaults to \code{Inf}.}
+\item{stat_threshold}{A stopping criterion for individual chain simulations;
+a positive number coercible to integer. When any chain's cumulative statistic
+reaches or surpasses \code{stat_threshold}, that chain ends. It also serves as a
+censoring limit so that results above the specified value, are set to \code{Inf}.
+Defaults to \code{Inf}.}
 }
 \value{
 An \verb{<epichains_summary>} object.
diff --git a/man/dot-new_epichains.Rd b/man/dot-new_epichains.Rd
index c92fbcae..13b5e054 100644
--- a/man/dot-new_epichains.Rd
+++ b/man/dot-new_epichains.Rd
@@ -9,7 +9,7 @@
   n_chains,
   statistic,
   offspring_dist,
-  stat_max,
+  stat_threshold,
   track_pop
 )
 }
@@ -21,7 +21,7 @@ for "time", and "chain_id".}
 \item{n_chains}{Number of chains to simulate.}
 
 \item{statistic}{The chain statistic to track as the
-stopping criteria for each chain being simulated when \code{stat_max} is not
+stopping criteria for each chain being simulated when \code{stat_threshold} is not
 \code{Inf}; A \verb{<string>}. It can be one of:
 \itemize{
 \item "size": the total number of cases produced by a chain before it goes
@@ -39,9 +39,13 @@ on to the random number generating functions. Examples that can be provided
 here are \code{rpois} for Poisson distributed offspring, \code{rnbinom} for negative
 binomial offspring, or custom functions.}
 
-\item{stat_max}{The cut off for the chain statistic (size/length) being
-computed; A number coercible to integer. Results above \code{stat_max} are set to
-\code{stat_max}; An \verb{<integer>}. Defaults to \code{Inf}.}
+\item{stat_threshold}{A stopping criterion for individual chain simulations;
+a positive number coercible to integer. When any chain's cumulative statistic
+reaches or surpasses \code{stat_threshold}, that chain ends. Defaults to \code{Inf}.
+For example, if \code{statistic = "size"} and \code{stat_threshold = 10}, then any
+chain that produces 10 or more cases will stop. Note that setting
+\code{stat_threshold} does not guarantee that all chains will stop at the same
+value.}
 
 \item{track_pop}{Was the susceptible population tracked? Logical.}
 }
diff --git a/man/dot-new_epichains_summary.Rd b/man/dot-new_epichains_summary.Rd
index 4294b6e3..5f5ab315 100644
--- a/man/dot-new_epichains_summary.Rd
+++ b/man/dot-new_epichains_summary.Rd
@@ -9,7 +9,7 @@
   n_chains,
   statistic,
   offspring_dist,
-  stat_max
+  stat_threshold
 )
 }
 \arguments{
@@ -18,7 +18,7 @@
 \item{n_chains}{Number of chains to simulate.}
 
 \item{statistic}{The chain statistic to track as the
-stopping criteria for each chain being simulated when \code{stat_max} is not
+stopping criteria for each chain being simulated when \code{stat_threshold} is not
 \code{Inf}; A \verb{<string>}. It can be one of:
 \itemize{
 \item "size": the total number of cases produced by a chain before it goes
@@ -36,9 +36,11 @@ on to the random number generating functions. Examples that can be provided
 here are \code{rpois} for Poisson distributed offspring, \code{rnbinom} for negative
 binomial offspring, or custom functions.}
 
-\item{stat_max}{The cut off for the chain statistic (size/length) being
-computed; A number coercible to integer. Results above \code{stat_max} are set to
-\code{stat_max}; An \verb{<integer>}. Defaults to \code{Inf}.}
+\item{stat_threshold}{A stopping criterion for individual chain simulations;
+a positive number coercible to integer. When any chain's cumulative statistic
+reaches or surpasses \code{stat_threshold}, that chain ends. It also serves as a
+censoring limit so that results above the specified value, are set to \code{Inf}.
+Defaults to \code{Inf}.}
 }
 \description{
 \code{new_epichains_summary()} constructs an \verb{<epichains_summary>} object from a
diff --git a/man/dot-offspring_ll.Rd b/man/dot-offspring_ll.Rd
index 19652d35..f4bf949b 100644
--- a/man/dot-offspring_ll.Rd
+++ b/man/dot-offspring_ll.Rd
@@ -20,7 +20,7 @@ here are \code{rpois} for Poisson distributed offspring, \code{rnbinom} for nega
 binomial offspring, or custom functions.}
 
 \item{statistic}{The chain statistic to track as the
-stopping criteria for each chain being simulated when \code{stat_max} is not
+stopping criteria for each chain being simulated when \code{stat_threshold} is not
 \code{Inf}; A \verb{<string>}. It can be one of:
 \itemize{
 \item "size": the total number of cases produced by a chain before it goes
diff --git a/man/head.epichains.Rd b/man/head.epichains.Rd
index 093eb188..c9df77d5 100644
--- a/man/head.epichains.Rd
+++ b/man/head.epichains.Rd
@@ -32,7 +32,7 @@ chains_pois_offspring <- simulate_chains(
   n_chains = 10,
   statistic = "size",
   offspring_dist = rpois,
-  stat_max = 10,
+  stat_threshold = 10,
   generation_time = function(n) rep(3, n),
   lambda = 2
 )
@@ -42,7 +42,7 @@ chains_pois_offspring <- simulate_chains(
   n_chains = 10,
   statistic = "size",
   offspring_dist = rpois,
-  stat_max = 10,
+  stat_threshold = 10,
   generation_time = function(n) rep(3, n),
   lambda = 2
 )
diff --git a/man/likelihood.Rd b/man/likelihood.Rd
index c24a784d..e58cb228 100644
--- a/man/likelihood.Rd
+++ b/man/likelihood.Rd
@@ -11,7 +11,7 @@ likelihood(
   nsim_obs,
   obs_prob = 1,
   log = TRUE,
-  stat_max = Inf,
+  stat_threshold = Inf,
   exclude = NULL,
   individual = FALSE,
   ...
@@ -21,7 +21,7 @@ likelihood(
 \item{chains}{A numeric vector of chain summaries (sizes/lengths).}
 
 \item{statistic}{The chain statistic to track as the
-stopping criteria for each chain being simulated when \code{stat_max} is not
+stopping criteria for each chain being simulated when \code{stat_threshold} is not
 \code{Inf}; A \verb{<string>}. It can be one of:
 \itemize{
 \item "size": the total number of cases produced by a chain before it goes
@@ -52,9 +52,11 @@ process is assumed to be constant.}
 \item{log}{Should the log-likelihoods be transformed to
 likelihoods? Logical. Defaults to \code{TRUE}.}
 
-\item{stat_max}{A cut off for the chain statistic (size/length) being
-computed. A number coercible to integer. Results above the specified value,
-are set to \code{Inf}.}
+\item{stat_threshold}{A stopping criterion for individual chain simulations;
+a positive number coercible to integer. When any chain's cumulative statistic
+reaches or surpasses \code{stat_threshold}, that chain ends. It also serves as a
+censoring limit so that results above the specified value, are set to \code{Inf}.
+Defaults to \code{Inf}.}
 
 \item{exclude}{A vector of indices of the sizes/lengths to exclude from the
 log-likelihood calculation.}
diff --git a/man/print.epichains_summary.Rd b/man/print.epichains_summary.Rd
index b993443c..6bdae92d 100644
--- a/man/print.epichains_summary.Rd
+++ b/man/print.epichains_summary.Rd
@@ -20,8 +20,8 @@ Prints a summary of the \verb{<epichains_summary>} object. In particular, it
 prints the number of chains simulated, and the range of
 the statistic, represented as the maximum (\code{max_stat}) and minimum
 (\code{min_stat}). If the minimum or maximum is infinite, it is represented as
-\verb{>= stat_max} where \code{stat_max} is the value of the censoring limit. See
-\code{?epichains_summary()} for the definition of \code{stat_max}.
+\verb{>= stat_threshold} where \code{stat_threshold} is the value of the censoring
+limit. See \code{?epichains_summary()} for the definition of \code{stat_threshold}.
 }
 \author{
 James M. Azam
diff --git a/man/rborel.Rd b/man/rborel.Rd
index 67ab8685..2c4138cf 100644
--- a/man/rborel.Rd
+++ b/man/rborel.Rd
@@ -4,20 +4,19 @@
 \alias{rborel}
 \title{Generate random numbers from the Borel distribution}
 \usage{
-rborel(n, mu, infinite = Inf)
+rborel(n, mu, censor_at = Inf)
 }
 \arguments{
 \item{n}{Number of random variates to generate.}
 
 \item{mu}{A non-negative number for the poisson mean.}
 
-\item{infinite}{A threshold above which the random number is set to
-\code{Inf} as a form of truncation; A number coercible to integer. This
-argument is important because \code{rborel()} uses \code{simulate_summary()}
-internally with \code{offspring_dist = rpois}, so when \code{mu >= 1}, the
-simulation could grow out of proportion and chain summaries (the random
-numbers generated here) need to be truncated using
-\code{infinite}.}
+\item{censor_at}{A stopping criterion; \verb{<numeric>}. Defaults to \code{Inf}. A
+value above which the simulation ends and the random number is set to
+\code{Inf} (as a form of censoring). \code{rborel()} simulates chain sizes using
+\code{\link[=simulate_chain_stats]{simulate_chain_stats()}} with a Poisson offspring distribution, so if
+\code{mu >= 1}, the simulation could proceed unendingly. This parameter is used
+to prevent this.}
 }
 \value{
 A numeric vector of random numbers.
diff --git a/man/rgborel.Rd b/man/rgborel.Rd
index d3258a25..fa5cc9b2 100644
--- a/man/rgborel.Rd
+++ b/man/rgborel.Rd
@@ -4,7 +4,7 @@
 \alias{rgborel}
 \title{Generate random numbers from a Gamma-Borel mixture distribution}
 \usage{
-rgborel(n, size, prob, mu, infinite = Inf)
+rgborel(n, size, prob, mu, censor_at = Inf)
 }
 \arguments{
 \item{n}{Number of random variates to generate.}
@@ -18,13 +18,12 @@ applications); A positive number.}
 
 \item{mu}{Mean; A positive number.}
 
-\item{infinite}{A threshold above which the random number is set to
-\code{Inf} as a form of truncation; A number coercible to integer. This
-argument is important because \code{rborel()} uses \code{simulate_summary()}
-internally with \code{offspring_dist = rpois}, so when \code{mu >= 1}, the
-simulation could grow out of proportion and chain summaries (the random
-numbers generated here) need to be truncated using
-\code{infinite}.}
+\item{censor_at}{A stopping criterion; \verb{<numeric>}. Defaults to \code{Inf}. A
+value above which the simulation ends and the random number is set to
+\code{Inf} (as a form of censoring). \code{rborel()} simulates chain sizes using
+\code{\link[=simulate_chain_stats]{simulate_chain_stats()}} with a Poisson offspring distribution, so if
+\code{mu >= 1}, the simulation could proceed unendingly. This parameter is used
+to prevent this.}
 }
 \value{
 Numeric vector of random numbers
diff --git a/man/simulate_chain_stats.Rd b/man/simulate_chain_stats.Rd
index 3e376ff7..2eacddd8 100644
--- a/man/simulate_chain_stats.Rd
+++ b/man/simulate_chain_stats.Rd
@@ -9,7 +9,7 @@ simulate_chain_stats(
   statistic = c("size", "length"),
   offspring_dist,
   ...,
-  stat_max = Inf,
+  stat_threshold = Inf,
   pop = Inf,
   percent_immune = 0
 )
@@ -18,7 +18,7 @@ simulate_chain_stats(
 \item{n_chains}{Number of chains to simulate.}
 
 \item{statistic}{The chain statistic to track as the
-stopping criteria for each chain being simulated when \code{stat_max} is not
+stopping criteria for each chain being simulated when \code{stat_threshold} is not
 \code{Inf}; A \verb{<string>}. It can be one of:
 \itemize{
 \item "size": the total number of cases produced by a chain before it goes
@@ -38,9 +38,11 @@ binomial offspring, or custom functions.}
 
 \item{...}{Parameters of the offspring distribution as required by R.}
 
-\item{stat_max}{A cut off for the chain statistic (size/length) being
-computed. A number coercible to integer. Results above the specified value,
-are set to \code{Inf}.}
+\item{stat_threshold}{A stopping criterion for individual chain simulations;
+a positive number coercible to integer. When any chain's cumulative statistic
+reaches or surpasses \code{stat_threshold}, that chain ends. It also serves as a
+censoring limit so that results above the specified value, are set to \code{Inf}.
+Defaults to \code{Inf}.}
 
 \item{pop}{Population size; An \verb{<Integer>}. Used alongside \code{percent_immune}
 to define the susceptible population. Defaults to \code{Inf}.}
@@ -116,7 +118,7 @@ simulate_chain_stats(
   n_chains = 20,
   statistic = "size",
   offspring_dist = rpois,
-  stat_max = 10,
+  stat_threshold = 10,
   lambda = 0.9
 )
 # Simulate chain sizes with a negative binomial distribution and assuming
@@ -128,7 +130,7 @@ simulate_chain_stats(
   n_chains = 20,
   statistic = "size",
   offspring_dist = rnbinom,
-  stat_max = 10,
+  stat_threshold = 10,
   mu = 0.9,
   size = 0.36
 )
diff --git a/man/simulate_chains.Rd b/man/simulate_chains.Rd
index 7ea31a8c..6962cf68 100644
--- a/man/simulate_chains.Rd
+++ b/man/simulate_chains.Rd
@@ -9,7 +9,7 @@ simulate_chains(
   statistic = c("size", "length"),
   offspring_dist,
   ...,
-  stat_max = Inf,
+  stat_threshold = Inf,
   pop = Inf,
   percent_immune = 0,
   generation_time = NULL,
@@ -21,7 +21,7 @@ simulate_chains(
 \item{n_chains}{Number of chains to simulate.}
 
 \item{statistic}{The chain statistic to track as the
-stopping criteria for each chain being simulated when \code{stat_max} is not
+stopping criteria for each chain being simulated when \code{stat_threshold} is not
 \code{Inf}; A \verb{<string>}. It can be one of:
 \itemize{
 \item "size": the total number of cases produced by a chain before it goes
@@ -41,9 +41,13 @@ binomial offspring, or custom functions.}
 
 \item{...}{Parameters of the offspring distribution as required by R.}
 
-\item{stat_max}{The cut off for the chain statistic (size/length) being
-computed; A number coercible to integer. Results above \code{stat_max} are set to
-\code{stat_max}; An \verb{<integer>}. Defaults to \code{Inf}.}
+\item{stat_threshold}{A stopping criterion for individual chain simulations;
+a positive number coercible to integer. When any chain's cumulative statistic
+reaches or surpasses \code{stat_threshold}, that chain ends. Defaults to \code{Inf}.
+For example, if \code{statistic = "size"} and \code{stat_threshold = 10}, then any
+chain that produces 10 or more cases will stop. Note that setting
+\code{stat_threshold} does not guarantee that all chains will stop at the same
+value.}
 
 \item{pop}{Population size; An \verb{<Integer>}. Used alongside \code{percent_immune}
 to define the susceptible population. Defaults to \code{Inf}.}
@@ -142,7 +146,7 @@ chains_pois_offspring <- simulate_chains(
   n_chains = 10,
   statistic = "size",
   offspring_dist = rpois,
-  stat_max = 10,
+  stat_threshold = 10,
   generation_time = function(n) rep(3, n),
   lambda = 2
 )
@@ -157,7 +161,7 @@ chains_nbinom_offspring <- simulate_chains(
   percent_immune = 0,
   statistic = "size",
   offspring_dist = rnbinom,
-  stat_max = 10,
+  stat_threshold = 10,
   generation_time = function(n) rep(3, n),
   mu = 2,
   size = 0.2
diff --git a/man/summary.epichains.Rd b/man/summary.epichains.Rd
index 327cd696..49be81cc 100644
--- a/man/summary.epichains.Rd
+++ b/man/summary.epichains.Rd
@@ -36,7 +36,7 @@ sim_chains_nbinom <- simulate_chains(
   percent_immune = 0,
   statistic = "size",
   offspring_dist = rnbinom,
-  stat_max = 10,
+  stat_threshold = 10,
   mu = 2,
   size = 0.2
 )
@@ -52,7 +52,7 @@ sim_summary_nbinom <- simulate_chain_stats(
   percent_immune = 0,
   statistic = "size",
   offspring_dist = rnbinom,
-  stat_max = 10,
+  stat_threshold = 10,
   mu = 2,
   size = 0.2
 )
diff --git a/tests/testthat/test-checks.R b/tests/testthat/test-checks.R
index 5b640cd1..38a9682a 100644
--- a/tests/testthat/test-checks.R
+++ b/tests/testthat/test-checks.R
@@ -5,7 +5,7 @@
         n_chains = 10,
         statistic = "size",
         offspring_dist = rpois,
-        stat_max = 10,
+        stat_threshold = 10,
         pop = 10,
         percent_immune = 0.1
     )
@@ -78,7 +78,7 @@ test_that("Smaller checker functions work", {
     expect_no_error(
         .check_statistic_args(
             statistic = "size",
-            stat_max = 10
+            stat_threshold = 10
         )
     )
 })
@@ -116,10 +116,10 @@ test_that(".check_sim_args() returns errors", {
         ),
         "object 'r' not found"
     )
-    # stat_max must be >= 1
+    # stat_threshold must be >= 1
     expect_error(
         .check_sim_args_default(
-            stat_max = 0
+            stat_threshold = 0
         ),
         "Assertion failed."
     )
diff --git a/tests/testthat/test-epichains.R b/tests/testthat/test-epichains.R
index eb1ae7e1..b79cf19c 100644
--- a/tests/testthat/test-epichains.R
+++ b/tests/testthat/test-epichains.R
@@ -80,7 +80,7 @@ test_that("Simulators return epichains objects", {
     n_chains = 10,
     statistic = "size",
     offspring_dist = rpois,
-    stat_max = 10,
+    stat_threshold = 10,
     generation_time = generation_time_fn,
     lambda = 2
   )
@@ -150,7 +150,7 @@ test_that("print.epichains works for simulation functions", {
     n_chains = 10,
     statistic = "size",
     offspring_dist = rpois,
-    stat_max = 10,
+    stat_threshold = 10,
     generation_time = generation_time_fn,
     lambda = 2
   )
@@ -162,23 +162,23 @@ test_that("print.epichains works for simulation functions", {
     statistic = "length",
     lambda = 0.9
   )
-  #' Simulate the case where Infs are produced and printed as ">= stat_max"
+  #' Simulate the case where Infs are produced and printed as ">= stat_threshold"
   set.seed(32)
   chain_lengths_with_Infs <- simulate_chain_stats(
     n_chains = 10,
     offspring_dist = rpois,
     statistic = "length",
     lambda = 1.1,
-    stat_max = 10
+    stat_threshold = 10
   )
-  #' Simulate the case where all are Infs printed as ">= stat_max"
+  #' Simulate the case where all are Infs printed as ">= stat_threshold"
   set.seed(32)
   chain_lengths_all_Infs <- simulate_chain_stats(
     n_chains = 2,
     offspring_dist = rpois,
     statistic = "length",
     lambda = 1.1,
-    stat_max = 10
+    stat_threshold = 10
   )
   #' Expectations
   expect_snapshot(susc_outbreak_raw)
@@ -276,7 +276,7 @@ test_that("validate_epichains works", {
     n_chains = 10,
     statistic = "size",
     offspring_dist = rpois,
-    stat_max = 10,
+    stat_threshold = 10,
     generation_time = generation_time_fn,
     lambda = 2
   )
@@ -338,7 +338,7 @@ test_that("is_chains_tree works", {
     n_chains = 10,
     statistic = "size",
     offspring_dist = rpois,
-    stat_max = 10,
+    stat_threshold = 10,
     generation_time = generation_time_fn,
     lambda = 2
   )
@@ -400,7 +400,7 @@ test_that("is_chains_summary works", {
     n_chains = 10,
     statistic = "size",
     offspring_dist = rpois,
-    stat_max = 10,
+    stat_threshold = 10,
     generation_time = generation_time_fn,
     lambda = 2
   )
@@ -436,7 +436,7 @@ test_that("aggregate.epichains method returns correct objects", {
     n_chains = 10,
     statistic = "size",
     offspring_dist = rpois,
-    stat_max = 10,
+    stat_threshold = 10,
     generation_time = generation_time_fn,
     lambda = 2
   )
@@ -467,7 +467,7 @@ test_that("aggregate.epichains method throws errors", {
         n_chains = 10,
         statistic = "size",
         offspring_dist = rpois,
-        stat_max = 10,
+        stat_threshold = 10,
         lambda = 2
       ),
       by = "time"
@@ -484,7 +484,7 @@ test_that("aggregate.epichains method is numerically correct", {
     n_chains = 10,
     statistic = "size",
     offspring_dist = rpois,
-    stat_max = 10,
+    stat_threshold = 10,
     lambda = 2
   )
   #' Simulate a tree of infections in an infinite population and with
@@ -493,7 +493,7 @@ test_that("aggregate.epichains method is numerically correct", {
     n_chains = 10,
     statistic = "size",
     offspring_dist = rpois,
-    stat_max = 10,
+    stat_threshold = 10,
     generation_time = generation_time_fn,
     lambda = 2
   )
@@ -532,7 +532,7 @@ test_that("head and tail print output as expected", {
     n_chains = 10,
     statistic = "size",
     offspring_dist = rpois,
-    stat_max = 10,
+    stat_threshold = 10,
     generation_time = generation_time_fn,
     lambda = 2
   )
diff --git a/tests/testthat/test-likelihood.R b/tests/testthat/test-likelihood.R
index abe5f92e..8bf3d5dd 100644
--- a/tests/testthat/test-likelihood.R
+++ b/tests/testthat/test-likelihood.R
@@ -26,7 +26,7 @@ test_that("Likelihoods can be calculated", {
         statistic = "size",
         offspring_dist = rpois,
         lambda = 0.5,
-        stat_max = 5
+        stat_threshold = 5
       ),
       0
     )
@@ -47,7 +47,7 @@ test_that("Likelihoods can be calculated", {
         statistic = "size",
         offspring_dist = rpois,
         lambda = 0.5,
-        stat_max = 5,
+        stat_threshold = 5,
         obs_prob = 0.5,
         nsim_obs = 1
       ),
diff --git a/tests/testthat/test-simulate.R b/tests/testthat/test-simulate.R
index 41984b25..72940556 100644
--- a/tests/testthat/test-simulate.R
+++ b/tests/testthat/test-simulate.R
@@ -27,7 +27,7 @@ test_that("simulate_chains has expected shape", {
     n_chains = 10,
     offspring_dist = rpois,
     statistic = "length",
-    stat_max = 10,
+    stat_threshold = 10,
     lambda = 2
   )
   # Simulate a tree of infections in an infinite population and with
@@ -60,7 +60,7 @@ test_that("simulate_chains has expected shape", {
     statistic = "size",
     offspring_dist = rpois,
     lambda = 2,
-    stat_max = 10,
+    stat_threshold = 10,
     generation_time = generation_time_fn,
     tf = 5
   )
@@ -328,7 +328,7 @@ test_that("simulate_chains produces expected snapshots", {
     n_chains = 10,
     offspring_dist = rpois,
     statistic = "size",
-    stat_max = 100,
+    stat_threshold = 100,
     lambda = 2,
     generation_time = generation_time_fn
   )
diff --git a/vignettes/epichains.Rmd b/vignettes/epichains.Rmd
index 4d8bb805..328ff784 100644
--- a/vignettes/epichains.Rmd
+++ b/vignettes/epichains.Rmd
@@ -303,7 +303,7 @@ sim_chains <- simulate_chains(
   n_chains = 10,
   statistic = "size",
   offspring_dist = rpois,
-  stat_max = 10,
+  stat_threshold = 10,
   generation_time = generation_time_fn,
   lambda = 0.9
 )
@@ -362,7 +362,7 @@ simulate_chain_stats_eg <- simulate_chain_stats(
   n_chains = 10,
   statistic = "size",
   offspring_dist = rpois,
-  stat_max = 10,
+  stat_threshold = 10,
   lambda = 0.9
 )
 
@@ -388,7 +388,7 @@ sim_chains <- simulate_chains(
   n_chains = 10,
   statistic = "size",
   offspring_dist = rpois,
-  stat_max = 10,
+  stat_threshold = 10,
   generation_time = generation_time_fn,
   lambda = 0.9
 )
@@ -402,7 +402,7 @@ simulate_chain_stats_eg <- simulate_chain_stats(
   n_chains = 10,
   statistic = "size",
   offspring_dist = rpois,
-  stat_max = 10,
+  stat_threshold = 10,
   lambda = 0.9
 )
 
@@ -451,7 +451,7 @@ sim_chains <- simulate_chains(
   n_chains = 10,
   statistic = "size",
   offspring_dist = rpois,
-  stat_max = 10,
+  stat_threshold = 10,
   generation_time = generation_time_fn,
   lambda = 0.9
 )
@@ -478,7 +478,7 @@ sim_chains <- simulate_chains(
   n_chains = 10,
   statistic = "size",
   offspring_dist = rpois,
-  stat_max = 1000,
+  stat_threshold = 1000,
   generation_time = generation_time_fn,
   lambda = 2
 )
diff --git a/vignettes/interventions.Rmd b/vignettes/interventions.Rmd
index a94c2db6..b18d857a 100644
--- a/vignettes/interventions.Rmd
+++ b/vignettes/interventions.Rmd
@@ -54,7 +54,7 @@ We simulate 200 chains tracking up to 99 infections:
 
 ```{r simulate_chains}
 sims <- simulate_chain_stats(
-  n_chains = 200, offspring_dist = rnbinom, stat_max = 99, mu = 1.2,
+  n_chains = 200, offspring_dist = rnbinom, stat_threshold = 99, mu = 1.2,
   size = 0.5, statistic = "size"
 )
 ```
@@ -80,7 +80,7 @@ For example, to reduce R by 25% at the population level we scale the `mu` parame
 
 ```{r simulate_chains_pop_control}
 sims <- simulate_chain_stats(
-  n_chains = 200, offspring_dist = rnbinom, stat_max = 99, mu = 0.9,
+  n_chains = 200, offspring_dist = rnbinom, stat_threshold = 99, mu = 0.9,
   size = 0.5, statistic = "size"
 )
 sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting.
@@ -118,7 +118,7 @@ Having defined this, we can generate simulations as before:
 
 ```{r simulate_chains_ind_control}
 sims <- simulate_chain_stats(
-  n_chains = 200, offspring_dist = rnbinom_ind, stat_max = 99, mu = 1.2,
+  n_chains = 200, offspring_dist = rnbinom_ind, stat_threshold = 99, mu = 1.2,
   size = 0.5, control = 0.25, statistic = "size"
 )
 sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting.
@@ -148,7 +148,7 @@ This can be likened to a disease control strategy where gatherings are limited t
 
 ```{r simulate_chains_truncated}
 sims <- simulate_chain_stats(
-  n_chains = 200, offspring_dist = rnbinom_truncated, stat_max = 99,
+  n_chains = 200, offspring_dist = rnbinom_truncated, stat_threshold = 99,
   mu = 1.2, size = 0.5, max =  10, statistic = "size"
 )
 sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting.
diff --git a/vignettes/projecting_incidence.Rmd b/vignettes/projecting_incidence.Rmd
index fa5aeded..1ab0568c 100644
--- a/vignettes/projecting_incidence.Rmd
+++ b/vignettes/projecting_incidence.Rmd
@@ -198,14 +198,14 @@ Lastly, since, we have specified that $R0 > 1$, it means the epidemic could
 potentially grow without end. Hence, we must specify an end point for the
 simulations. 
 
-`simulate_chains()` provides the `stat_max` argument for this purpose.
-Above `stat_max`, the simulation is cut off. If this value is 
+`simulate_chains()` provides the `stat_threshold` argument for this purpose.
+Above `stat_threshold`, the simulation is cut off. If this value is 
 not specified, it assumes a value of infinity. Here, we will
 assume a maximum chain size of $1000$.
 
-```{r stat_max_setup}
+```{r stat_threshold_setup}
 #' Maximum chain size allowed
-stat_max <- 1000
+stat_threshold <- 1000
 ```
 
 ## Modelling assumptions
@@ -222,7 +222,7 @@ This exercise makes the following simplifying assumptions:
 To summarise the whole set up so far, we are going to simulate 
 each chain `r sim_rep` times, projecting cases over
 `r projection_window` days after the first `r max(t0)` days, and 
-assuming that no outbreak size exceeds `r stat_max` cases.
+assuming that no outbreak size exceeds `r stat_threshold` cases.
 
 ## Running the simulations
 
@@ -239,7 +239,7 @@ sim_chain_sizes <- lapply(
       mu = mu,
       size = size,
       statistic = "size",
-      stat_max = stat_max,
+      stat_threshold = stat_threshold,
       generation_time = generation_time,
       t0 = t0,
       tf = tf