diff --git a/.lintr b/.lintr
index 7bedc8fa5..6588b4777 100644
--- a/.lintr
+++ b/.lintr
@@ -2,7 +2,16 @@ linters: linters_with_tags(
     tags = NULL, # include all linters
     implicit_integer_linter = NULL,
     extraction_operator_linter = NULL,
-    undesirable_function_linter = NULL,
+    keyword_quote_linter = NULL,
+    undesirable_function_linter(
+      fun = c(
+         # Base messaging
+         "message" = "use cli::cli_inform()",
+         "warning" = "use cli::cli_warn()",
+         "stop" = "use cli::cli_abort()",
+         "stopifnot" = "use cli::cli_abort()"
+      )
+    ),
     function_argument_linter = NULL,
     indentation_linter = NULL,
     object_name_linter = NULL,
diff --git a/DESCRIPTION b/DESCRIPTION
index 41fd5e0d8..0ad0fc252 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -100,6 +100,7 @@ Depends:
     R (>= 3.5.0)
 Imports:
     checkmate,
+    cli,
     data.table,
     futile.logger (>= 1.4),
     future,
@@ -148,7 +149,7 @@ Encoding: UTF-8
 Language: en-GB
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.2.9000
 NeedsCompilation: yes
 SystemRequirements: GNU make
        C++17
diff --git a/NAMESPACE b/NAMESPACE
index ff206f9aa..1da104c88 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -125,6 +125,11 @@ importFrom(checkmate,assert_string)
 importFrom(checkmate,assert_subset)
 importFrom(checkmate,test_data_frame)
 importFrom(checkmate,test_numeric)
+importFrom(cli,cli_abort)
+importFrom(cli,cli_inform)
+importFrom(cli,cli_warn)
+importFrom(cli,col_blue)
+importFrom(cli,col_red)
 importFrom(data.table,":=")
 importFrom(data.table,.N)
 importFrom(data.table,as.data.table)
diff --git a/NEWS.md b/NEWS.md
index 92c69cf07..7fa1e344f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -8,6 +8,7 @@
 
 ## Model changes
 
+- All functions now use the `{cli}` R package to signal errors, warnings, and messages. By @jamesmbaazam in #762 and reviewed by @seabbs.
 - `epinow()` now returns the "timing" output in a "time difference"" format that is easier to understand and work with. By @jamesmbaazam in #688 and reviewed by @sbfnk.
 - The interface for defining delay distributions has been generalised to also cater for continuous distributions
 - When defining probability distributions these can now be truncated using the `tolerance` argument
diff --git a/R/checks.R b/R/checks.R
index a43d8eb29..6a741a48f 100644
--- a/R/checks.R
+++ b/R/checks.R
@@ -61,6 +61,7 @@ check_reports_valid <- function(data,
 #' @param dist A `dist_spec` object.`
 #' @importFrom checkmate assert_class
 #' @importFrom rlang arg_match
+#' @importFrom cli cli_abort col_blue
 #' @return Called for its side effects.
 #' @keywords internal
 check_stan_delay <- function(dist) {
@@ -73,9 +74,12 @@ check_stan_delay <- function(dist) {
   if (
     !all(distributions %in% c("lognormal", "gamma", "fixed", "nonparametric"))
   ) {
-    stop(
-      "Distributions passed to the model need to be lognormal, gamma, fixed ",
-      "or nonparametric."
+    cli_abort(
+      c(
+       "!" = "Distributions passed to the model need to be
+        {col_blue(\"lognormal\")}, {col_blue(\"gamma\")},
+        {col_blue(\"fixed\")}, or {col_blue(\"nonparametric\")}."
+      )
     )
   }
   # Check that `dist` has parameters that are either numeric or normal
@@ -91,10 +95,13 @@ check_stan_delay <- function(dist) {
     }
   }))
   if (!all(numeric_or_normal)) {
-    stop(
-      "Delay distributions passed to the model need to have parameters that ",
-      "are either numeric or normally distributed with numeric parameters ",
-      "and infinite maximum."
+    cli_abort(
+      c(
+        "!" = "Delay distributions passed to the model need to have parameters
+        that are either {col_blue(\"numeric\")} or
+        {col_blue(\"normally distributed\")} with {col_blue(\"numeric\")}
+        parameters and {col_blue(\"infinite maximum\")}."
+      )
     )
   }
   if (is.null(attr(dist, "tolerance"))) {
@@ -103,9 +110,12 @@ check_stan_delay <- function(dist) {
   assert_numeric(attr(dist, "tolerance"), lower = 0, upper = 1)
   # Check that `dist` has a finite maximum
   if (any(is.infinite(max(dist))) && !(attr(dist, "tolerance") > 0)) {
-    stop(
-      "All distribution passed to the model need to have a finite maximum,",
-      "which can be achieved either by setting `max` or non-zero `tolerance`."
+    cli_abort(
+      c(
+        "i" = "All distribution passed to the model need to have a
+      {col_blue(\"finite maximum\")}, which can be achieved either by
+      setting {.var max} or non-zero {.var tolerance}."
+      )
     )
   }
 }
@@ -117,19 +127,21 @@ check_stan_delay <- function(dist) {
 #' @param pmf A probability mass function vector
 #' @param span The number of consecutive indices in the tail to check
 #' @param tol The value which to consider the tail as sparse
+#' @importFrom cli cli_warn col_blue
 #'
 #' @return Called for its side effects.
 #' @keywords internal
 check_sparse_pmf_tail <- function(pmf, span = 5, tol = 1e-6) {
   if (all(pmf[(length(pmf) - span + 1):length(pmf)] < tol)) {
-    warning(
-      sprintf(
-        "The PMF tail has %s consecutive values smaller than %s.",
-        span, tol
+    cli_warn(
+      c(
+        "!" = "The PMF tail has {col_blue(span)} consecutive value{?s} smaller
+        than {col_blue(tol)}.",
+        "i" = "This will drastically increase run times with very small
+        increases in accuracy. Consider increasing the tail values of the PMF."
       ),
-      " This will drastically increase run time with very small increases ",
-      "in accuracy. Consider increasing the tail values of the PMF.",
-      call. = FALSE
+      .frequency = "regularly",
+      .frequency_id = "sparse_pmf_tail"
     )
   }
 }
diff --git a/R/create.R b/R/create.R
index 55891230e..fb9677f50 100644
--- a/R/create.R
+++ b/R/create.R
@@ -246,6 +246,7 @@ create_future_rt <- function(future = c("latest", "project", "estimate"),
 #' breakpoints.
 #'
 #' @param horizon Numeric, forecast horizon.
+#' @importFrom cli cli_abort
 #'
 #' @seealso rt_settings
 #' @return A list of settings defining the time-varying reproduction number
@@ -285,7 +286,11 @@ create_rt_data <- function(rt = rt_opts(), breakpoints = NULL,
   # apply random walk
   if (rt$rw != 0) {
     if (is.null(breakpoints)) {
-      stop("breakpoints must be supplied when using random walk")
+      cli_abort(
+        c(
+          "!" = "breakpoints must be supplied when using random walk."
+        )
+      )
     }
 
     breakpoints <- seq_along(breakpoints)
diff --git a/R/deprecated.R b/R/deprecated.R
index 36234a48f..61a288e5b 100644
--- a/R/deprecated.R
+++ b/R/deprecated.R
@@ -595,13 +595,18 @@ dist_skel <- function(n, dist = FALSE, cum = TRUE, model,
 #' removed.
 #' @return A `<dist_spec>` where probability masses below the threshold level
 #' have been removed
+#' @importFrom cli cli_abort
 #' @keywords internal
 apply_tolerance <- function(x, tolerance) {
   lifecycle::deprecate_warn(
     "1.6.0", "apply_tolerance()", "bound_dist()"
   )
   if (!is(x, "dist_spec")) {
-    stop("Can only apply tolerance to distributions in a <dist_spec>.")
+    cli_abort(
+      c(
+        "!" = "Can only apply tolerance to distributions in a {.cls dist_spec}."
+      )
+    )
   }
   y <- lapply(x, function(x) {
     if (x$distribution == "nonparametric") {
diff --git a/R/dist_spec.R b/R/dist_spec.R
index ca3c52764..25a0a1aa4 100644
--- a/R/dist_spec.R
+++ b/R/dist_spec.R
@@ -135,6 +135,7 @@ discrete_pmf <- function(distribution =
 #' cannot be combined with other combinations of distributions.
 #'
 #' @param ... The delay distributions to combine
+#' @importFrom cli cli_abort
 #' @return Combined delay distributions (with class `<dist_spec>`)
 #' @method c dist_spec
 #' @export
@@ -155,8 +156,10 @@ c.dist_spec <- function(...) {
   dist_specs <- list(...)
   if (length(dist_specs) == 1) return(dist_specs[[1]])
   if (!(all(vapply(dist_specs, is, "dist_spec", FUN.VALUE = logical(1))))) {
-    stop(
-      "All distributions must be of class `<dist_spec>`."
+    cli_abort(
+      c(
+        "!" = "All distributions must be of class {.cls dist_spec}."
+      )
     )
   }
   convolutions <- vapply(
@@ -165,7 +168,11 @@ c.dist_spec <- function(...) {
   ## can only have one `multi_dist_spec`
   if (sum(convolutions) > 0) {
     if (sum(convolutions) > 1) {
-      stop("Can't convolve convolutions with other convolutions")
+      cli_abort(
+        c(
+          "!" = "Can't convolve convolutions with other convolutions"
+        )
+      )
     }
     ## preserve convolution attribute
     convolution_attributes <- attributes(dist_specs[[which(convolutions)]])
@@ -190,6 +197,7 @@ c.dist_spec <- function(...) {
 #' @param ignore_uncertainty Logical; whether to ignore any uncertainty in
 #'   parameters. If set to FALSE (the default) then the mean of any uncertain
 #'   parameters will be returned as NA.
+#' @importFrom cli cli_abort
 #' @method mean dist_spec
 #' @importFrom utils head
 #' @export
@@ -230,7 +238,11 @@ mean.dist_spec <- function(x, ..., ignore_uncertainty = FALSE) {
     } else if (dist == "fixed") {
       return(params$value)
     } else {
-      stop("Don't know how to calculate mean of ", dist, " distribution.")
+      cli_abort(
+        c(
+          "!" = "Don't know how to calculate mean of {dist} distribution."
+        )
+      )
     }
   }
 }
@@ -257,6 +269,7 @@ sd <- function(x, ...) {
 #' @param x The <dist_spec> to use
 #' @return A vector of standard deviations.
 #' @importFrom utils head
+#' @importFrom cli cli_abort
 #' @keywords internal
 #' @export
 #' @examples
@@ -292,9 +305,11 @@ sd.dist_spec <- function(x, ...) {
     } else if (x$distribution == "fixed") {
       0
     } else {
-      stop(
-        "Don't know how to calculate standard deviation of ",
-        x$distribution, " distribution."
+      cli_abort(
+        c(
+          "!" = "Don't know how to calculate standard deviation of
+        {x$distribution} distribution."
+        )
       )
     }
   }
@@ -370,6 +385,7 @@ discretise <- function(x, ...) {
 #' specified or parameters are uncertain). If `FALSE` then any distribution
 #' that cannot be discretised will be returned as is.
 #' @param ... ignored
+#' @importFrom cli cli_abort
 #' @return A `<dist_spec>` where all distributions with constant parameters are
 #'   nonparametric.
 #' @export
@@ -391,7 +407,12 @@ discretise.dist_spec <- function(x, strict = TRUE, ...) {
   }
   max_x <- max(x)
   if (is.infinite(max_x) && !(tolerance > 0) && strict) {
-    stop("Cannot discretise a distribution with infinite support.")
+    cli_abort(
+      c(
+        "!" = "Cannot discretise a distribution with infinite support.",
+        "i" = "Either set a finite maximum or a tolerance greater than 0."
+      )
+    )
   }
   if (get_distribution(x) == "nonparametric") {
     return(x)
@@ -411,8 +432,10 @@ discretise.dist_spec <- function(x, strict = TRUE, ...) {
       }
       return(y)
     } else if (strict) {
-      stop(
-        "Cannot discretise a distribution with uncertain parameters."
+      cli_abort(
+        c(
+          "!" = "Cannot discretise a distribution with uncertain parameters."
+        )
       )
     } else {
       return(x)
@@ -439,6 +462,7 @@ discretize <- discretise
 #' @return A `<dist_spec>` where consecutive nonparametric distributions
 #' have been convolved
 #' @importFrom stats convolve
+#' @importFrom cli cli_abort
 #' @export
 #' @examples
 #' # A fixed gamma distribution with mean 5 and sd 1.
@@ -451,7 +475,11 @@ discretize <- discretise
 #' collapse(discretise(dist1 + dist2))
 collapse <- function(x) {
   if (!is(x, "dist_spec")) {
-    stop("Can only convolve distributions in a <dist_spec>.")
+    cli_abort(
+      c(
+        "!" = "Can only convolve distributions in a {.cls dist_spec}."
+      )
+    )
   }
   ## get nonparametric distributions
   nonparametric <- vapply(
@@ -489,6 +517,7 @@ collapse <- function(x) {
 #' functions of fixed delay distributions combined in the passed <dist_spec>.
 #' @param x The `<dist_spec>` to use
 #' @param ... Not used
+#' @importFrom cli cli_abort
 #' @return invisible
 #' @method print dist_spec
 #' @export
@@ -578,6 +607,7 @@ print.dist_spec <- function(x, ...) {
 #' @param ... ignored
 #' @importFrom ggplot2 aes geom_col geom_step facet_wrap vars theme_bw
 #' @importFrom data.table data.table rbindlist
+#' @importFrom cli cli_abort
 #' @export
 #' @examples
 #' # A fixed lognormal distribution with mean 5 and sd 1.
@@ -626,9 +656,13 @@ plot.dist_spec <- function(x, samples = 50L, res = 1, cumulative = TRUE, ...) {
       }
       pmf_dt <- lapply(dists, function(y) {
         if (is.infinite(attr(y, "max"))) {
-          stop(
-            "Need to set a finite maximum, either as an argument to `plot()`",
-            "or when defining the distribution."
+          cli_abort(
+            c(
+              "!" = "All distributions in {.var x} must have a finite
+              maximum value.",
+              "i" = "You can set a finite maximum either as an
+              argument to {.fn plot} or when defining the distribution."
+            )
           )
         }
         x <- discrete_pmf(
@@ -676,6 +710,7 @@ plot.dist_spec <- function(x, samples = 50L, res = 1, cumulative = TRUE, ...) {
 #' @description `r lifecycle::badge("experimental")`
 #' @param x A composite `dist_spec` object
 #' @param i The index to extract
+#' @importFrom cli cli_abort
 #' @return A single `dist_spec` object
 #' @keywords internal
 #' @examples
@@ -693,7 +728,12 @@ plot.dist_spec <- function(x, samples = 50L, res = 1, cumulative = TRUE, ...) {
 #' }
 extract_single_dist <- function(x, i) {
   if (i > ndist(x)) {
-    stop("i can't be greater than the number of distributions.")
+    cli_abort(
+      c(
+        "!" = "i must be less than the number of distributions.",
+        "i" = "The number of distributions is {ndist(x)} whiles i is {i}."
+      )
+    )
   }
   if (ndist(x) == 1) {
     return(x)
@@ -968,11 +1008,17 @@ lower_bounds <- function(distribution) {
 #' @param tolerance Numeric; the desired tolerance level. Any part of the
 #' cumulative distribution function beyond 1 minus this tolerance level is
 #' removed. Default: `0`, i.e. use the full distribution.
+#' @importFrom cli cli_abort
 #' @return a `<dist_spec>` with relevant attributes set that define its bounds
 #' @export
 bound_dist <- function(x, max = Inf, tolerance = 0) {
   if (!is(x, "dist_spec")) {
-    stop("Can only get limit a <dist_spec>.")
+    cli_abort(
+      c(
+        "!" = "{.var x} must be of class {.cls dist_spec}.",
+        "i" = "It is currently of class {.cls class(x)}."
+      )
+    )
   }
   ## if it is a single nonparametric distribution we apply the bounds directly
   if (ndist(x) == 1 && get_distribution(x) == "nonparametric") {
@@ -1000,14 +1046,19 @@ bound_dist <- function(x, max = Inf, tolerance = 0) {
 #' @param params Given parameters (obtained using `as.list(environment())`)
 #' @return A character vector of parameters and their values.
 #' @inheritParams natural_params
+#' @importFrom cli cli_abort
 #' @keywords internal
 extract_params <- function(params, distribution) {
   params <- params[!vapply(params, inherits, "name", FUN.VALUE = TRUE)]
   n_params <- length(natural_params(distribution))
   if (length(params) != n_params) {
-    stop(
-      "Exactly ", n_params, " parameters of the ", distribution,
-      " distribution must be specified."
+    cli_abort(
+      c(
+        "!" = "Exactly {n_params} parameters of the {distribution}
+        distribution must be specified.",
+        "i" = "You have specified {length(params)} parameters, which is not
+        equal to {n_params}."
+      )
     )
   }
   return(params)
@@ -1023,6 +1074,7 @@ extract_params <- function(params, distribution) {
 #' @inheritParams extract_params
 #' @inheritParams bound_dist
 #' @importFrom purrr walk
+#' @importFrom cli cli_abort cli_warn
 #' @return A `dist_spec` of the given specification.
 #' @export
 #' @examples
@@ -1052,9 +1104,12 @@ new_dist_spec <- function(params, distribution, max = Inf, tolerance = 0) {
       for (param_name in names(params)) {
         lb <- lower_bounds(distribution)[param_name]
         if (is.numeric(params[[param_name]]) && params[[param_name]] < lb) {
-          stop(
-            "Parameter ", param_name, " is less than its lower bound ", lb,
-            "."
+          cli_abort(
+            c(
+              "!" = "Parameter {param_name} must be greater than its
+              lower bound {lb}.",
+              "i" = "It is currently set to less than the lower bound."
+            )
           )
         }
       }
@@ -1069,15 +1124,19 @@ new_dist_spec <- function(params, distribution, max = Inf, tolerance = 0) {
           return(is.na(sd_dist) || sd_dist > 0)
         }, logical(1))
         if (any(uncertain)) {
-          warning(
-            "Uncertain ", distribution, " distribution specified in terms of ",
-            "parameters that are not the \"natural\" parameters of the ",
-            "distribution (", toString(natural_params(distribution)),
-            "). Converting using a crude and very approximate method ",
-            "that is likely to produce biased results. If possible, ",
-            "it is preferable to specify the distribution directly ",
-            "in terms of the natural parameters."
+          #nolint start: duplicate_argument_linter
+          cli_warn(
+            c(
+              "!" = "Uncertain {distribution} distribution specified in
+              terms of parameters that are not the \"natural\" parameters of
+              the distribution {natural_params(distribution)}.",
+              "i" = "Converting using a crude and very approximate method
+            that is likely to produce biased results.",
+              "i" = "If possible it is preferable to specify the
+            distribution directly in terms of the natural parameters."
+            )
           )
+          #nolint end
         }
         ## generate natural parameters
         params <- convert_to_natural(params, distribution)
@@ -1109,6 +1168,7 @@ new_dist_spec <- function(params, distribution, max = Inf, tolerance = 0) {
 #' from a given set of parameters and distribution
 #' @param params A numerical named parameter vector
 #' @inheritParams natural_params
+#' @importFrom cli cli_abort
 #' @return A list with two elements, `params_mean` and `params_sd`, containing
 #' mean and sd of natural parameters.
 #' @keywords internal
@@ -1123,9 +1183,13 @@ convert_to_natural <- function(params, distribution) {
   ## unnatural parameter means
   ux <- lapply(params, mean)
   if (anyNA(ux)) {
-    stop(
-      "Cannot nest uncertainty in a distributions that is not specified with ",
-      "its natural parameters."
+    cli_abort(
+      c(
+        "!" = "Cannot nest uncertainty in a distributions that is not
+      specified with its natural parameters.",
+        "i" = "Specify the distribution in terms of its natural
+      parameters if you want to nest uncertainty."
+      )
     )
   }
   ## estimate relative uncertainty of parameters
@@ -1160,10 +1224,11 @@ convert_to_natural <- function(params, distribution) {
   ## sort
   x <- x[natural_params(distribution)]
   if (anyNA(names(x))) {
-    stop(
-      "Incompatible combination of parameters of a ", distribution,
-      " distribution specified:\n    ", toString(names(params)),
-      "."
+    cli_abort(
+      c(
+        "!" = "Incompatible combination of parameters of a {distribution}
+      distribution specified: {names(params)}."
+      )
     )
   }
   if (rel_unc > 0) {
@@ -1184,18 +1249,27 @@ convert_to_natural <- function(params, distribution) {
 ##' distribution). If `x` is a single distribution this is ignored and can be
 ##' left at its default value of `NULL`.
 ##' @param element The element, i.e. "parameters", "pmf" or "distribution".
+##' @importFrom cli cli_abort
 ##' @return The id to use.
 ##' @keywords internal
 get_element <- function(x, id = NULL, element) {
   if (!is.null(id) && id > ndist(x)) {
-    stop(
-      "`id` can't be greater than the number of distributions (", length(x),
-      ")."
+    cli_abort(
+      c(
+        "!" = "{.var id} cannot be greater than the number of distributions
+      ({length(x)}).",
+        "i" = "{.var id} currently has length {length(id)}."
+      )
     )
   }
   if (ndist(x) > 1) {
     if (is.null(id)) {
-      stop("`id` must be specified when `x` is a composite distribution.")
+      cli_abort(
+        c(
+         "!" = "{.var id} must be specified when {.var x} is a composite
+          distribution."
+        )
+      )
     }
     return(x[[id]][[element]])
   } else {
@@ -1207,6 +1281,7 @@ get_element <- function(x, id = NULL, element) {
 ##'
 ##' @inheritParams get_element
 ##' @description `r lifecycle::badge("experimental")`
+##' @importFrom cli cli_abort
 ##' @return A list of parameters of the distribution.
 ##' @export
 ##' @examples
@@ -1214,10 +1289,23 @@ get_element <- function(x, id = NULL, element) {
 ##' get_parameters(dist)
 get_parameters <- function(x, id = NULL) {
   if (!is(x, "dist_spec")) {
-    stop("Can only get parameters of a <dist_spec>.")
+    cli_abort(
+      c(
+        "!" = "Object must be of class {.cls dist_spec}",
+        "i" = "You have supplied an object of class {.cls {class(x)}}."
+      )
+    )
   }
   if (get_distribution(x, id) == "nonparametric") {
-    stop("Cannot get parameters of a nonparametric distribution.")
+    cli_abort(
+      c(
+        "!" = "To get parameters, distribution cannot not be
+        \"nonparametric\".",
+        "i" = "Distribution must be one of
+        {col_blue(\"gamma\")}, {col_blue(\"lognormal\")},
+        {col_blue(\"normal\")} or {col_blue(\"fixed\")}."
+      )
+    )
   }
   return(get_element(x, id, "parameters"))
 }
@@ -1227,16 +1315,26 @@ get_parameters <- function(x, id = NULL) {
 ##' @inheritParams get_element
 ##' @description `r lifecycle::badge("experimental")`
 ##' @return The pmf of the distribution
+##' @importFrom cli cli_abort
 ##' @export
 ##' @examples
 ##' dist <- discretise(Gamma(shape = 3, rate = 2, max = 10))
 ##' get_pmf(dist)
 get_pmf <- function(x, id = NULL) {
   if (!is(x, "dist_spec")) {
-    stop("Can only get pmf of a <dist_spec>.")
+    cli_abort(
+      c(
+        "!" = "Can only get pmf of a {.cls dist_spec}.",
+        "i" = "You have supplied an object of class {.cls {class(x)}}."
+      )
+    )
   }
   if (get_distribution(x, id) != "nonparametric") {
-    stop("Cannot get pmf of a parametric distribution.")
+    cli_abort(
+      c(
+        "!" = "To get PMF, distribution must be \"nonparametric\"."
+      )
+    )
   }
   return(get_element(x, id, "pmf"))
 }
@@ -1245,6 +1343,7 @@ get_pmf <- function(x, id = NULL) {
 ##'
 ##' @inheritParams get_element
 ##' @description `r lifecycle::badge("experimental")`
+##' @importFrom cli cli_abort
 ##' @return A character string naming the distribution (or "nonparametric")
 ##' @export
 ##' @examples
@@ -1252,7 +1351,12 @@ get_pmf <- function(x, id = NULL) {
 ##' get_distribution(dist)
 get_distribution <- function(x, id = NULL) {
   if (!is(x, "dist_spec")) {
-    stop("Can only get distribution of a <dist_spec>.")
+    cli_abort(
+      c(
+        "!" = "To get distribution of x, it must be a {.cls dist_spec}.",
+        "i" = "You have supplied an object of class {.cls {class(x)}}."
+      )
+    )
   }
   return(get_element(x, id, "distribution"))
 }
diff --git a/R/estimate_delay.R b/R/estimate_delay.R
index 5573e801d..11a5ad59f 100644
--- a/R/estimate_delay.R
+++ b/R/estimate_delay.R
@@ -24,6 +24,7 @@
 #' @return A stan fit of an interval censored distribution
 #' @export
 #' @inheritParams stan_opts
+#' @importFrom cli cli_warn col_blue
 #' @examples
 #' \donttest{
 #' # integer adjusted exponential model
@@ -50,10 +51,12 @@ dist_fit <- function(values = NULL, samples = 1000, cores = 1,
                      backend = "rstan") {
   if (samples < 1000) {
     samples <- 1000
-    warning(sprintf("%s %s", "`samples` must be at least 1000.",
-                    "Now setting it to 1000 internally."
-                    )
-            )
+    cli_warn(
+      c(
+        "!" = "{.var samples} must be at least {col_blue(\"1000\")}.",
+        "i" = "Now setting it to {col_blue(\"1000\")} internally."
+      )
+    )
   }
   # model parameters
   lows <- values - 1
@@ -149,6 +152,7 @@ dist_fit <- function(values = NULL, samples = 1000, cores = 1,
 #' @importFrom future.apply future_lapply
 #' @importFrom rstan extract
 #' @importFrom data.table data.table rbindlist
+#' @importFrom cli cli_abort col_blue
 #' @export
 #' @examples
 #' \donttest{
@@ -165,9 +169,14 @@ bootstrapped_dist_fit <- function(values, dist = "lognormal",
                                   bootstrap_samples = 250, max_value,
                                   verbose = FALSE) {
   if (!dist %in% c("gamma", "lognormal")) {
-    stop("Only lognormal and gamma distributions are supported")
+    cli_abort(
+      c(
+        "x" = "Unsupported distribution.",
+        "i" = "Only {col_blue(\"lognormal\")} and {col_blue(\"gamma\")}
+      distributions are supported"
+      )
+    )
   }
-
   if (samples < bootstraps) {
     samples <- bootstraps
   }
diff --git a/R/estimate_secondary.R b/R/estimate_secondary.R
index 87909779a..077825381 100644
--- a/R/estimate_secondary.R
+++ b/R/estimate_secondary.R
@@ -69,6 +69,7 @@
 #' @importFrom utils modifyList
 #' @importFrom checkmate assert_class assert_numeric assert_data_frame
 #' assert_logical
+#' @importFrom cli cli_abort col_red
 #' @examples
 #' \donttest{
 #' # set number of cores to use
@@ -201,8 +202,13 @@ estimate_secondary <- function(data,
   )
 
   if (burn_in >= nrow(reports)) {
-    stop("burn_in is greater or equal to the number of observations.
-         Some observations must be used in fitting")
+    cli_abort(
+      c(
+        "!" = "{.var burn_in} is {col_red(\"greater or equal to the number of
+        observations\")}.",
+        "i" = "Supply a {.var burn_in} less than the number of observations."
+      )
+    )
   }
   # observation and control data
   stan_data <- list(
@@ -281,6 +287,7 @@ estimate_secondary <- function(data,
 #' @export
 #' @inheritParams create_stan_args
 #' @importFrom data.table as.data.table
+#' @importFrom cli cli_inform cli_warn
 #' @examples
 #' priors <- data.frame(variable = "frac_obs", mean = 3, sd = 1)
 #' data <- list(obs_scale_mean = 4, obs_scale_sd = 3)
@@ -289,8 +296,9 @@ update_secondary_args <- function(data, priors, verbose = TRUE) {
   priors <- data.table::as.data.table(priors)
   if (!missing(priors) && !is.null(priors) && nrow(priors) > 0) {
     if (verbose) {
-      message(
-        "Replacing specified priors with those from the passed in prior dataframe" # nolint
+      cli_inform(
+        "Replacing specified priors with those passed through the
+        {.var prior} dataframe."
       )
     }
     # replace scaling if present in the prior
@@ -303,8 +311,11 @@ update_secondary_args <- function(data, priors, verbose = TRUE) {
     delay_params <- priors[grepl("delay_params", variable, fixed = TRUE)]
     if (nrow(delay_params) > 0) {
       if (is.null(data$delay_params_mean)) {
-        warning(
-          "Cannot replace delay distribution parameters as no default has been set" # nolint
+        cli_warn(
+          c(
+            "!" = "Cannot replace delay distribution parameters.",
+            "i" = "No default has been set."
+          )
         )
       }
       data$delay_params_mean <- as.array(signif(delay_params$mean, 3))
diff --git a/R/extract.R b/R/extract.R
index 331c100cf..3c6d04489 100644
--- a/R/extract.R
+++ b/R/extract.R
@@ -63,6 +63,7 @@ extract_static_parameter <- function(param, samples) {
 #' @param pars Any selection of parameters to extract
 #' @param include whether the parameters specified in `pars` should be included
 #' (`TRUE`, the default) or excluded (`FALSE`)
+#' @importFrom cli cli_abort
 #' @return List of data.tables with samples
 #' @export
 #'
@@ -76,7 +77,10 @@ extract_samples <- function(stan_fit, pars = NULL, include = TRUE) {
   }
   if (!inherits(stan_fit, "CmdStanMCMC") &&
       !inherits(stan_fit, "CmdStanFit")) {
-    stop("stan_fit must be a <stanfit>, <CmdStanMCMC> or <CmdStanFit> object")
+    cli_abort(
+      "{.var stan_fit} must be a {.cls stanfit}, {.cls CmdStanMCMC} or
+      {.cls CmdStanFit} object."
+    )
   }
 
   # extract sample from stan object
diff --git a/R/fit.R b/R/fit.R
index 6f2e53257..953f78996 100644
--- a/R/fit.R
+++ b/R/fit.R
@@ -163,6 +163,7 @@ fit_model_with_nuts <- function(args, future = FALSE, max_execution_time = Inf,
 #' @importFrom purrr safely
 #' @importFrom rstan vb
 #' @importFrom rlang abort
+#' @importFrom cli cli_abort
 #' @return A stan model object
 #' @keywords internal
 fit_model_approximate <- function(args, future = FALSE, id = "stan") {
@@ -192,7 +193,13 @@ fit_model_approximate <- function(args, future = FALSE, id = "stan") {
       if (method == "vb") {
         sample_func <- rstan::vb
       } else {
-        stop("Laplace approximation only available in the cmdstanr backend")
+        cli_abort(
+          c(
+            "!" = "Laplace approximation only available in the cmdstanr
+            backend.",
+            "i" = "You've supplied {.strong {method}}."
+          )
+        )
       }
     } else if (inherits(stan_args$object, "CmdStanModel")) {
       if (method == "vb") {
diff --git a/R/opts.R b/R/opts.R
index fe352100e..0c3bff03e 100644
--- a/R/opts.R
+++ b/R/opts.R
@@ -19,6 +19,7 @@
 #'   will be applied, i.e. any parameters in `dist` will be treated as a single
 #'   parameters.
 #' @inheritParams apply_default_tolerance
+#' @importFrom cli cli_warn cli_abort col_blue
 #' @return A `<generation_time_opts>` object summarising the input delay
 #' distributions.
 #' @seealso [convert_to_logmean()] [convert_to_logsd()]
@@ -52,22 +53,29 @@ gt_opts <- function(dist = Fixed(1), ...,
       !missing(disease) || !missing(source) || !missing(fixed) ||
       !missing(max) ||
       (!is(dist, "dist_spec"))) {
-    stop(
-      "The generation time distribution should be given to ",
-      "`generation_time_opts` using a `dist_spec`. ",
-      "This behaviour has changed from previous versions of `EpiNow2` and ",
-      "any code using it may need to be updated as any other ways of ",
-      "specifying the generation time are deprecated.",
-      "For examples and more ",
-      "information, see the relevant documentation pages using ",
-      "`?generation_time_opts`")
+    cli_abort(
+      c(
+        "!" = "The generation time distribution must be passed through
+      {.fn gt_opts} or {.fn generation_time_opts}. ",
+        "i" = "This behaviour has changed from previous versions of `EpiNow2`
+      and any code using it must be updated as any other ways of
+      specifying the generation time are deprecated. {col_blue(\"For
+      examples and more information, see the relevant documentation
+      pages using ?gt_opts\")}"
+      )
+    )
   }
   if (missing(dist)) {
-    warning(
-      "No generation time distribution given. Assuming a fixed generation ",
-      "time of 1 day, i.e. the reproduction number is the same as the daily ",
-      "growth rate. If this was intended then this warning can be silenced by ",
-      "setting `dist` explicitly to `Fixed(1)`."
+    #nolint start: duplicate_argument_linter
+    cli_warn(
+      c(
+        "!" = "No generation time distribution given.",
+        "i" = "Now using a fixed generation time of 1 day, i.e. the
+        reproduction number is the same as the daily growth rate.",
+        "i" = "If this was intended then this warning can be
+        silenced by setting {.var dist = Fixed(1)}'."
+      )
+      #nolint end
     )
   }
   ## apply default tolerance if `dist` is unconstrained
@@ -162,6 +170,7 @@ secondary_opts <- function(type = c("incidence", "prevalence"), ...) {
 #' @param ... deprecated; use `dist` instead
 #' @param fixed deprecated; use `dist` instead
 #' @inheritParams generation_time_opts
+#' @importFrom cli cli_abort
 #' @return A `<delay_opts>` object summarising the input delay distributions.
 #' @seealso [convert_to_logmean()] [convert_to_logsd()]
 #' [bootstrapped_dist_fit()] \code{\link{Distributions}}
@@ -184,20 +193,26 @@ delay_opts <- function(dist = Fixed(0), ..., fixed = FALSE,
                        default_tolerance = 0.001, weight_prior = TRUE) {
   dot_options <- list(...)
   if (!is(dist, "dist_spec") || !missing(fixed)) { ## could be old syntax
-    stop(
-      "Delay distributions must be of given either using a call to ",
-      "`dist_spec` or one of the `get_...` functions such as ",
-      "`get_incubation_period`. ",
-      "This behaviour has changed from previous versions of `EpiNow2` and ",
-      "any code using it may need to be updated as any other ways of ",
-      "specifying delays are deprecated. ",
-      "For examples and more ",
-      "information, see the relevant documentation pages using ",
-      "`?delay_opts`."
+    #nolint start: duplicate_argument_linter
+    cli_abort(
+      c(
+        "!" = "Delay distributions must be given using a call to
+        {.fn delay_opts}",
+        "i" = "This behaviour has changed from previous versions of `EpiNow2`
+      and any code using it must be updated as any other ways of specifying
+      delays are deprecated.",
+        "i" = "For examples and more information, see the relevant
+        documentation pages using {.code ?delay_opts}."
+      )
     )
+    #nolint end
   } else if (length(dot_options) > 0) {
     ## can be removed once dot options are hard deprecated
-    stop("Unknown named arguments passed to `delay_opts`")
+    cli_abort(
+      c(
+        "!" = "Unknown named arguments passed to {.fn delay_opts}."
+      )
+    )
   }
   ## apply default tolerance if `dist` is unconstrained
   dist <- apply_default_tolerance(
@@ -229,6 +244,7 @@ delay_opts <- function(dist = Fixed(0), ..., fixed = FALSE,
 #'   i.e. the truncation distribution will be treated as a single parameter.
 #'
 #' @inheritParams gt_opts
+#' @importFrom cli cli_abort
 #' @return A `<trunc_opts>` object summarising the input truncation
 #' distribution.
 #'
@@ -244,15 +260,18 @@ delay_opts <- function(dist = Fixed(0), ..., fixed = FALSE,
 trunc_opts <- function(dist = Fixed(0), default_tolerance = 0.001,
                        weight_prior = FALSE) {
   if (!is(dist, "dist_spec")) {
-    stop(
-      "Truncation distributions must be of given either using a call to ",
-      "`dist_spec` or one of the `get_...` functions. ",
-      "This behaviour has changed from previous versions of `EpiNow2` and ",
-      "any code using it may need to be updated as any other ways of ",
-      "specifying delays are deprecated and will be removed in ",
-      "the next version. For examples and more ",
-      "information, see the relevant documentation pages using ",
-      "`?trunc_opts`"
+    #nolint start: duplicate_argument_linter
+    cli_abort(
+      c(
+        "!" = "Truncation distributions must be given using a call to
+        {.fn trunc_opts}.",
+        "i" = "This behaviour has changed from previous versions of `EpiNow2`
+      and any code using it must be updated as any other ways of specifying
+      delays are deprecated.",
+        "i" = "For examples and more information, see the relevant
+        documentation pages using {.code ?trunc_opts}."
+      )
+    #nolint end
     )
   }
   ## apply default tolerance if `dist` is unconstrained
@@ -307,6 +326,7 @@ trunc_opts <- function(dist = Fixed(0), default_tolerance = 0.001,
 #' reproduction number.
 #' @inheritParams create_future_rt
 #' @importFrom rlang arg_match
+#' @importFrom cli cli_abort
 #' @export
 #' @examples
 #' # default settings
@@ -340,7 +360,13 @@ rt_opts <- function(prior = list(mean = 1, sd = 1),
   }
 
   if (!("mean" %in% names(rt$prior) && "sd" %in% names(rt$prior))) {
-    stop("prior must have both a mean and sd specified")
+    cli_abort(
+      c(
+        "!" = "{.var prior} must have both {.var mean} and {.var sd}
+        specified.",
+        "i" = "Did you forget to specify {.var mean} and/or {.var sd}?"
+      )
+    )
   }
   attr(rt, "class") <- c("rt_opts", class(rt))
   return(rt)
@@ -374,6 +400,7 @@ rt_opts <- function(prior = list(mean = 1, sd = 1),
 #' average to use when estimating Rt. This must be odd so that the central
 #' estimate is included.
 #' @importFrom rlang arg_match
+#' @importFrom cli cli_abort
 #'
 #' @return A `<backcalc_opts>` object of back calculation settings.
 #' @export
@@ -388,9 +415,12 @@ backcalc_opts <- function(prior = c("reports", "none", "infections"),
     rt_window = as.integer(rt_window)
   )
   if (backcalc$rt_window %% 2 == 0) {
-    stop(
-      "Rt rolling average window must be odd in order to include the current
-       estimate"
+    cli_abort(
+      c(
+        "!" = "{.var rt_window} must be odd in order to
+        include the current estimate.",
+        "i" = "You have supplied an even number."
+      )
     )
   }
   attr(backcalc, "class") <- c("backcalc_opts", class(backcalc))
@@ -453,6 +483,7 @@ backcalc_opts <- function(prior = c("reports", "none", "infections"),
 #' kernel. They are only used if `kernel` is set to "periodic".
 #'
 #' @importFrom rlang arg_match
+#' @importFrom cli cli_abort cli_warn
 #' @return A `<gp_opts>` object of settings defining the Gaussian process
 #' @export
 #' @examples
@@ -485,9 +516,12 @@ gp_opts <- function(basis_prop = 0.2,
 
   if (!missing(matern_type)) {
     if (!missing(matern_order) && matern_type != matern_order) {
-      stop(
-        "Incompatible `matern_order` and `matern_type`. ",
-        "Use `matern_order` only."
+      cli_abort(
+        c(
+          "!" = "{.var matern_order} and {.var matern_type} must be the same, if
+          both are supplied.",
+          "i" = "Rather only use {.var matern_order} only."
+        )
       )
     }
     matern_order <- matern_type
@@ -501,9 +535,11 @@ gp_opts <- function(basis_prop = 0.2,
   } else if (
       !(is.infinite(matern_order) || matern_order %in% c(1 / 2, 3 / 2, 5 / 2))
     ) {
-    warning(
-      "Uncommon Matern kernel order. Common orders are `1 / 2`, `3 / 2`,", # nolint
-      " and `5 / 2`" # nolint
+    cli_warn(
+      c(
+        "!" = "Uncommon Matern kernel order supplied.",
+        "i" = "Use one of `1 / 2`, `3 / 2`, or `5 / 2`" # nolint
+      )
     )
   }
 
@@ -565,6 +601,7 @@ gp_opts <- function(basis_prop = 0.2,
 #' @param return_likelihood Logical, defaults to `FALSE`. Should the likelihood
 #'   be returned by the model.
 #' @importFrom rlang arg_match
+#' @importFrom cli cli_inform cli_abort
 #' @return An `<obs_opts>` object of observation model settings.
 #' @export
 #' @examples
@@ -587,20 +624,29 @@ obs_opts <- function(family = c("negbin", "poisson"),
                      return_likelihood = FALSE) {
   na <- arg_match(na)
   if (na == "accumulate") {
-    message(
-      "Accumulating modelled values that correspond to NA values in the data ",
-      "by adding them to the next non-NA data point. This means that the ",
-      "first data point is not included in the likelihood but used only to ",
-      "reset modelled observations to zero. If the first data point should be ",
-      "included in the likelihood this can be achieved by adding a data point ",
-      "of arbitrary value before the first data point."
+    #nolint start: duplicate_argument_linter
+    cli_inform(
+      c(
+        "i" = "Accumulating modelled values that correspond to NA values in the
+      data by adding them to the next non-NA data point.",
+        "i" = "This means that the first data point is not included in the
+      likelihood but used only to reset modelled observations to zero.",
+        "i" = "{col_red('If the first data point should be included in the
+        likelihood this can be achieved by adding a data point of arbitrary
+        value before the first data point.')}"
+      ),
+      .frequency = "regularly",
+      .frequency_id = "obs_opts"
     )
+  #nolint end
   }
 
   if (length(phi) == 2 && is.numeric(phi)) {
-    stop(
-      "Specifying `phi` as a length 2 vector is deprecated. Mean and SD ",
-      "should be given as list elements."
+    cli_abort(
+      c(
+        "!" = "Specifying {.var phi} as a vector of length 2 is deprecated.",
+        "i" = "Mean and SD should be given as list elements."
+      )
     )
   }
   obs <- list(
@@ -620,7 +666,13 @@ obs_opts <- function(family = c("negbin", "poisson"),
       obs[[param]] <- list(mean = obs[[param]], sd = 0)
     }
     if (!(all(c("mean", "sd") %in% names(obs[[param]])))) {
-      stop("If specifying a ", param, " as list both a mean and sd are needed")
+      cli_abort(
+        c(
+          "!" = "Both a {.var mean} and {.var sd} are needed if specifying
+        {.strong {param}} as list.",
+          "i" = "Did you forget to specify {.var mean} and/or {.var sd}?"
+        )
+      )
     }
   }
 
@@ -673,6 +725,7 @@ obs_opts <- function(family = c("negbin", "poisson"),
 #' @param ... Additional parameters to pass to [rstan::sampling()] or
 #' [cmdstanr::sample()].
 #' @importFrom utils modifyList
+#' @importFrom cli cli_warn
 #' @return A list of arguments to pass to [rstan::sampling()] or
 #' [cmdstanr::sample()].
 #' @export
@@ -701,9 +754,11 @@ stan_sampling_opts <- function(cores = getOption("mc.cores", 1L),
   control_def <- list(adapt_delta = 0.9, max_treedepth = 12)
   control_def <- modifyList(control_def, control)
   if (any(c("iter", "iter_sampling") %in% names(dot_args))) {
-    warning(
-      "Number of samples should be specified using the `samples` and `warmup`",
-      "arguments rather than `iter` or `iter_sampliing` which will be ignored."
+    cli_warn(
+      "!" = "Number of samples must be specified using the {.var samples}
+      and {.var warmup} arguments rather than {.var iter} or
+      {.var iter_sampliing}.",
+      "i" = "Supplied {.var iter} or {.var iter_sampliing} will be ignored."
     )
   }
   dot_args$iter <- NULL
@@ -770,6 +825,7 @@ stan_vb_opts <- function(samples = 2000,
 #' @inheritParams stan_opts
 #' @inheritParams stan_vb_opts
 #' @param ... Additional parameters to pass to [cmdstanr::laplace()].
+#' @importFrom cli cli_abort col_blue
 #' @return A list of arguments to pass to [cmdstanr::laplace()].
 #' @export
 #' @examples
@@ -778,8 +834,12 @@ stan_laplace_opts <- function(backend = "cmdstanr",
                               trials = 10,
                               ...) {
   if (backend != "cmdstanr") {
-    stop(
-      "The Laplace algorithm is only available with the \"cmdstanr\" backend."
+    cli_abort(
+      c(
+        "!" = "Backend must be set to {col_blue(\"cmdstanr\")} to use
+        the Laplace algorithm.",
+        "i" = "Change {.var backend} to col_blue(\"cmdstanr\")}."
+      )
     )
   }
   opts <- list(trials = trials)
@@ -795,6 +855,7 @@ stan_laplace_opts <- function(backend = "cmdstanr",
 #' @inheritParams stan_opts
 #' @inheritParams stan_vb_opts
 #' @param ... Additional parameters to pass to [cmdstanr::laplace()].
+#' @importFrom cli cli_abort col_blue
 #' @return A list of arguments to pass to [cmdstanr::laplace()].
 #' @export
 #' @examples
@@ -804,9 +865,12 @@ stan_pathfinder_opts <- function(backend = "cmdstanr",
                                  trials = 10,
                                  ...) {
   if (backend != "cmdstanr") {
-    stop(
-      "The pathfinder algorithm is only available with the \"cmdstanr\" ",
-      "backend."
+    cli_abort(
+      c(
+        "!" = "Backend must be set to {col_blue(\"cmdstanr\")} to use
+        the pathfinder algorithm.",
+        "i" = "Change {.var backend} to col_blue(\"cmdstanr\")}."
+      )
     )
   }
   opts <- list(
@@ -847,6 +911,7 @@ stan_pathfinder_opts <- function(backend = "cmdstanr",
 #'   [stan_sampling_opts()] or [stan_vb_opts()], depending on the method
 #'
 #' @importFrom rlang arg_match
+#' @importFrom cli cli_abort cli_warn col_blue
 #' @return A `<stan_opts>` object of arguments to pass to the appropriate
 #' rstan functions.
 #' @export
@@ -869,17 +934,22 @@ stan_opts <- function(object = NULL,
   backend_passed <- !missing(backend)
   backend <- arg_match(backend)
   if (backend == "cmdstanr" && !requireNamespace("cmdstanr", quietly = TRUE)) {
-    stop(
-      "The `cmdstanr` package needs to be installed for using the ",
-      "\"cmdstanr\" backend."
+    cli_abort(
+      c(
+        "x" = "The {col_blue('cmdstanr')} R package is not installed.",
+        "i" = "Install it from {.url https://github.com/stan-dev/cmdstanr}
+        to use the {col_blue('cmdstanr')} backend."
+      )
     )
   }
   opts <- list()
   if (!is.null(object)) {
     if (backend_passed) {
-      warning(
-        "`backend` option will be ignored as a stan model object has been ",
-        "passed."
+      cli_warn(
+        c(
+          "!" = "{.var backend} option will be ignored as a stan model
+        object has been passed."
+        )
       )
     }
     if (inherits(object, "stanmodel")) {
@@ -887,7 +957,11 @@ stan_opts <- function(object = NULL,
     } else if (inherits(object, "CmdStanModel")) {
       backend <- "cmdstanr"
     } else {
-      stop("`object` must be a stan model object")
+      cli_abort(
+        c(
+          "!" = "{.var object} must be a stan model object."
+        )
+      )
     }
   } else {
     backend <- arg_match(backend, values = c("rstan", "cmdstanr"))
@@ -1002,23 +1076,31 @@ filter_opts <- function(opts, region) {
 #'   constrained by having a maximum or tolerance this is ignored.
 #' @param tolerance_set Logical; whether the default tolerance has been set by
 #'   the user; if yes and `dist` is constrained a warning is issued
+#' @importFrom cli cli_inform cli_warn
 #'
 #' @return A <dist_spec> with the default tolerance set if previously not
 #'   constrained
 #' @keywords internal
 apply_default_tolerance <- function(dist, default_tolerance, tolerance_set) {
   if (!is_constrained(dist)) {
-    message(
-      "Unconstrained distributon passed as a delay. ",
-      "Constraining with default tolerance ", default_tolerance,  ".\n",
-      "To silence this message, specify delay distributions with `max` or ",
-      "`tolerance`."
+    #nolint start: duplicate_argument_linter
+    cli_inform(
+      c(
+        "i" = "Unconstrained distributon passed as a delay. ",
+        "i" = "Constraining with default tolerance {default_tolerance}",
+        "i" = "To silence this message, specify delay distributions
+      with {.var max} or {.var tolerance}."
+      )
     )
+    #nolint end
     attr(dist, "tolerance") <- default_tolerance
   } else if (tolerance_set) {
-    warning(
-      "Ignoring given default tolerance as distribution is already ",
-      "constrained.")
+    cli_warn(
+      c(
+        "!" = "Ignoring given default tolerance.",
+        "i" = "Distribution is already constrained."
+      )
+    )
   }
   return(dist)
 }
diff --git a/R/regional_epinow.R b/R/regional_epinow.R
index da06b12a9..cc0296e89 100644
--- a/R/regional_epinow.R
+++ b/R/regional_epinow.R
@@ -275,6 +275,7 @@ clean_regions <- function(data, non_zero_points) {
       "Producing estimates for: %s regions",
       length(eval_regions)
     )
+    #nolint start: undesirable_function_linter
     message <- ifelse(length(orig_regions) == 0, 0,
       length(orig_regions)
     )
@@ -282,11 +283,13 @@ clean_regions <- function(data, non_zero_points) {
       "Regions excluded: %s regions",
       message
     )
+    #nolint end
   } else {
     futile.logger::flog.info(
       "Producing estimates for: %s",
       toString(eval_regions)
     )
+    #nolint start: undesirable_function_linter
     message <- ifelse(length(orig_regions) == 0, "none",
       toString(orig_regions)
     )
@@ -294,6 +297,7 @@ clean_regions <- function(data, non_zero_points) {
       "Regions excluded: %s",
       message
     )
+    #nolint end
   }
   # exclude zero regions
   reported_cases <- reported_cases[!is.na(region)][region %in% eval_regions]
diff --git a/R/setup.R b/R/setup.R
index 107a1ece2..2dfd8c96d 100644
--- a/R/setup.R
+++ b/R/setup.R
@@ -25,6 +25,7 @@
 #'
 #' @importFrom futile.logger flog.threshold flog.appender appender.tee
 #' @importFrom futile.logger appender.file flog.info
+#' @importFrom cli cli_inform col_blue
 #' @return Nothing
 #' @export
 setup_logging <- function(threshold = "INFO", file = NULL,
@@ -32,26 +33,31 @@ setup_logging <- function(threshold = "INFO", file = NULL,
   if (is.null(name)) {
     name <- "ROOT"
   }
-  message(sprintf(
-    "Logging threshold set at %s for the %s logger",
-    threshold, name
-  ))
+  cli_inform(
+    "Logging threshold set at {threshold} for the name logger"
+  )
   futile.logger::flog.threshold(threshold, name = name)
 
   if (!is.null(file)) {
     if (mirror_to_console) {
-      message(sprintf("Writing %s logs to the console and: %s", name, file))
+      cli_inform(
+        "Writing {col_blue(name)} logs to the console and: {.file {file}}."
+      )
       futile.logger::flog.appender(
         futile.logger::appender.tee(file), name = name
       )
     } else {
-      message(sprintf("Writing %s logs to: %s", name, file))
+      cli_inform(
+        "Writing {col_blue(name)} logs to: {.file {file}}."
+      )
       futile.logger::flog.appender(
         futile.logger::appender.file(file), name = name
       )
     }
   } else {
-    message(sprintf("Writing %s logs to the console", name))
+    cli_inform(
+      "Writing {col_blue(name)} logs to the console."
+    )
     futile.logger::flog.appender(futile.logger::appender.console(), name = name)
   }
   return(invisible(NULL))
@@ -131,6 +137,7 @@ setup_default_logging <- function(logs = tempdir(check = TRUE),
 #' @inheritParams regional_epinow
 #' @importFrom futile.logger flog.error flog.info flog.debug
 #' @importFrom future availableCores plan tweak
+#' @importFrom cli cli_abort
 #' @export
 #' @return Numeric number of cores to use per worker. If greater than 1 pass to
 #' `stan_args = list(cores = "output from setup future")` or use
@@ -140,11 +147,19 @@ setup_future <- function(data,
                          min_cores_per_worker = 4) {
   if (length(strategies) > 2 || length(strategies) == 0) {
     futile.logger::flog.error("1 or 2 strategies should be used")
-    stop("1 or 2 strategies should be used")
+    cli_abort(
+      c(
+        "!" = "{.var strategies} must either be of length 1 or 2."
+      )
+    )
   }
   if (is.null(data$region)) {
     futile.logger::flog.error("Reported cases must contain a region")
-    stop("Exactly 2 strategies should be used")
+    cli_abort(
+      c(
+        "!" = "Exactly 2 strategies should be used."
+      )
+    )
   }
   if (length(strategies) == 1) {
     workers <- future::availableCores()
diff --git a/R/simulate_infections.R b/R/simulate_infections.R
index 4b0974ed8..f7ade498b 100644
--- a/R/simulate_infections.R
+++ b/R/simulate_infections.R
@@ -40,6 +40,7 @@
 #' @importFrom checkmate assert_data_frame assert_date assert_numeric
 #'   assert_subset assert_integer
 #' @importFrom data.table data.table merge.data.table nafill rbindlist
+#' @importFrom cli cli_abort
 #' @return A data.table of simulated infections (variable `infections`) and
 #'   reported cases (variable `reported_cases`) by date.
 #' @export
@@ -134,10 +135,12 @@ simulate_infections <- function(estimates, R, initial_infections,
   ))
 
   if (length(data$delay_params_sd) > 0 && any(data$delay_params_sd > 0)) {
-    stop(
-      "Cannot simulate from uncertain parameters. Use the [fix_dist()] ",
-      "function to set the parameters of uncertain distributions either the ",
-      "mean or a randomly sampled value"
+    cli_abort(
+      c(
+        "!" = "Cannot simulate from uncertain parameters.",
+        "i" = "Use {.fn fix_dist} to set the parameters of uncertain
+        distributions using either the mean or a randomly sampled value."
+      )
     )
   }
   data$delay_params <- array(
@@ -150,9 +153,11 @@ simulate_infections <- function(estimates, R, initial_infections,
   ))
 
   if (data$obs_scale_sd > 0) {
-    stop(
-      "Cannot simulate from uncertain observation scaling; use fixed scaling ",
-      "instead."
+    cli_abort(
+      c(
+        "!" = "Cannot simulate from uncertain observation scaling.",
+        "i" = "Use fixed scaling instead."
+      )
     )
   }
   if (data$obs_scale) {
@@ -165,9 +170,11 @@ simulate_infections <- function(estimates, R, initial_infections,
 
   if (obs$family == "negbin") {
     if (data$phi_sd > 0) {
-      stop(
-        "Cannot simulate from uncertain overdispersion; use fixed ",
-        "overdispersion instead."
+      cli_abort(
+        c(
+          "!" = "Cannot simulate from uncertain overdispersion.",
+          "i" = "Use fixed overdispersion instead."
+        )
       )
     }
     data$rep_phi <- array(data$phi_mean, dim = c(1, 1))
@@ -253,6 +260,7 @@ simulate_infections <- function(estimates, R, initial_infections,
 #' @importFrom lubridate days
 #' @importFrom checkmate assert_class assert_names test_numeric test_data_frame
 #' assert_numeric assert_integerish assert_logical
+#' @importFrom cli cli_abort
 #' @return A list of output as returned by [estimate_infections()] but based on
 #' results from the specified scenario rather than fitting.
 #' @seealso [dist_spec()] [generation_time_opts()] [delay_opts()] [rt_opts()]
@@ -314,11 +322,15 @@ forecast_infections <- function(estimates,
   ## check inputs
   assert_class(estimates, "estimate_infections")
   assert_names(names(estimates), must.include = "fit")
-  stopifnot(
-    "R must either be a numeric vector or a data.frame" =
-    test_numeric(R, lower = 0, null.ok = TRUE) ||
-    test_data_frame(R, null.ok = TRUE)
-  )
+  if (!(test_numeric(R, lower = 0, null.ok = TRUE) ||
+    test_data_frame(R, null.ok = TRUE))) {
+    cli_abort(
+      c(
+        "!" = "R must either be a {.cls numeric} vector or
+        a {.cls data.frame}."
+      )
+    )
+  }
   if (test_data_frame(R)) {
     assert_names(names(R), must.include = c("date", "value"))
     assert_numeric(R$value, lower = 0)
diff --git a/R/simulate_secondary.R b/R/simulate_secondary.R
index 956f8470b..c639c3bd8 100644
--- a/R/simulate_secondary.R
+++ b/R/simulate_secondary.R
@@ -18,6 +18,7 @@
 #' @inheritParams estimate_secondary
 #' @importFrom checkmate assert_data_frame assert_date assert_numeric
 #'   assert_subset
+#' @importFrom cli cli_abort
 #' @return A data.table of simulated secondary observations (column `secondary`)
 #'   by date.
 #' @export
@@ -76,10 +77,12 @@ simulate_secondary <- function(primary,
   ))
 
   if (length(data$delay_params_sd) > 0 && any(data$delay_params_sd > 0)) {
-    stop(
-      "Cannot simulate from uncertain parameters. Use the [fix_dist()] ",
-      "function to set the parameters of uncertain distributions either the ",
-      "mean or a randomly sampled value"
+    cli_abort(
+      c(
+        "!" = "Cannot simulate from uncertain parameters.",
+        "i" = "Use {.fn fix_dist} to set the parameters of uncertain
+        distributions either using the mean or a randomly sampled value."
+      )
     )
   }
   data$delay_params <- array(
@@ -92,9 +95,11 @@ simulate_secondary <- function(primary,
   ))
 
   if (data$obs_scale_sd > 0) {
-    stop(
-      "Cannot simulate from uncertain observation scaling; use fixed scaling ",
-      "instead."
+    cli_abort(
+      c(
+        "!" = "Cannot simulate from uncertain observation scaling.",
+        "i" = "Use fixed scaling instead."
+      )
     )
   }
   if (data$obs_scale) {
@@ -107,9 +112,11 @@ simulate_secondary <- function(primary,
 
   if (obs$family == "negbin") {
     if (data$phi_sd > 0) {
-      stop(
-        "Cannot simulate from uncertain overdispersion; use fixed ",
-        "overdispersion instead."
+      cli_abort(
+        c(
+          "!" = "Cannot simulate from uncertain overdispersion.",
+          "i" = "Use fixed overdispersion instead."
+        )
       )
     }
     data$rep_phi <- array(data$phi_mean, dim = c(1, 1))
diff --git a/R/stan.R b/R/stan.R
index 204134793..07b49bf8d 100644
--- a/R/stan.R
+++ b/R/stan.R
@@ -14,6 +14,7 @@
 #'
 #' @param ... Additional arguments passed to [cmdstanr::cmdstan_model()].
 #'
+#' @importFrom cli cli_inform col_blue
 #' @return A `cmdstanr` model.
 #' @export
 epinow2_cmdstan_model <- function(model = "estimate_infections",
@@ -26,8 +27,8 @@ epinow2_cmdstan_model <- function(model = "estimate_infections",
     dir, paste0(model, ".stan")
   )
   if (verbose) {
-    message(sprintf("Using model %s.", model))
-    message(sprintf("dir is %s.", toString(dir)))
+    cli_inform("Using model {col_blue(model)}.")
+    cli_inform("{.var dir} is {.file {dir}}.")
   }
 
   monitor <- suppressMessages
@@ -91,6 +92,7 @@ epinow2_stan_model <- function(backend = c("rstan", "cmdstanr"),
 #'
 #' Internal function for dispatch to fitting with NUTS or VB.
 #' @inheritParams fit_model_with_nuts
+#' @importFrom cli cli_abort
 #' @keywords internal
 fit_model <- function(args, id = "stan") {
   if (args$method == "sampling") {
@@ -102,7 +104,13 @@ fit_model <- function(args, id = "stan") {
   } else if (args$method %in% c("vb", "laplace", "pathfinder")) {
     fit <- fit_model_approximate(args, id = id)
   } else {
-    stop("method ", args$method, " unknown")
+    cli_abort(
+      c(
+        "!" = "You supplied method {args$method}, which is unknown.",
+        "i" = "Use one of {col_blue(\"sampling\")}, {col_blue(\"vb\")},
+      {col_blue(\"laplace\")}, or {col_blue(\"pathfinder\")}."
+      )
+    )
   }
   return(fit)
 }
diff --git a/R/summarise.R b/R/summarise.R
index b340d4ba4..24a1c20fe 100644
--- a/R/summarise.R
+++ b/R/summarise.R
@@ -21,6 +21,7 @@
 #'
 #' @importFrom purrr safely map_chr map_dbl map_chr
 #' @importFrom data.table setorderv melt merge.data.table dcast
+#' @importFrom cli cli_abort
 #' @return A list of summary data
 #' @keywords internal
 summarise_results <- function(regions,
@@ -30,16 +31,21 @@ summarise_results <- function(regions,
                               region_scale = "Region") {
   if (is.null(results_dir)) {
     if (is.null(summaries)) {
-      stop(
-        "Either a results directory or a list of summary data frames must be",
-        " supplied"
+      cli_abort(
+        c(
+          "!" = "One of a {.var results_dir} or {.var summary}
+          must be supplied."
+        )
       )
     }
   } else {
     if (!is.null(summaries)) {
-      stop(
-        "Both a results directory and a list of summary data frames have been",
-        " supplied."
+      cli_abort(
+        c(
+          "!" = "Cannot supply both {.var results_dir} and {.var summary}.",
+          "i" = "Only one of {.var results_dir} or {.var summary} should be
+        supplied."
+        )
       )
     }
   }
@@ -160,6 +166,7 @@ summarise_results <- function(regions,
 #' @importFrom ggplot2 coord_cartesian guides guide_legend ggsave ggplot_build
 #' @importFrom data.table setDT fcase
 #' @importFrom futile.logger flog.info
+#' @importFrom cli cli_abort
 #' @examples
 #' # get example output from regional_epinow model
 #' regional_out <- readRDS(system.file(
@@ -191,7 +198,14 @@ regional_summary <- function(regional_output = NULL,
   }
 
   if (!is.null(results_dir) && !is.null(regional_output)) {
-    stop("Only one of results_dir and regional_output should be specified")
+    cli_abort(
+      c(
+        "!" = "Both {.var results_dir} and {.var regional_output} cannot be
+        specified.",
+        "i" = "Only supply one of {.var results_dir} or
+        {.var regional_output}."
+      )
+    )
   }
 
   if (is.null(regional_output)) {
@@ -421,6 +435,7 @@ regional_summary <- function(regional_output = NULL,
 #' region).
 #'
 #' @inheritParams get_regional_results
+#' @importFrom cli cli_abort
 #' @seealso regional_summary
 #' @return A list of summarised Rt, cases by date of infection and cases by
 #' date of report
@@ -431,7 +446,11 @@ summarise_key_measures <- function(regional_results = NULL,
                                    type = "region", date = "latest") {
   if (is.null(regional_results)) {
     if (is.null(results_dir)) {
-      stop("Missing results directory")
+      cli_abort(
+        c(
+          "!" = "{.var results_dir} must be specified."
+        )
+      )
     }
     timeseries <- EpiNow2::get_regional_results(
       results_dir = results_dir,
@@ -499,6 +518,7 @@ summarise_key_measures <- function(regional_results = NULL,
 #' @export
 #' @importFrom data.table data.table fwrite
 #' @importFrom purrr map safely map_vec
+#' @importFrom cli cli_abort
 #' @keywords internal
 #' @examples
 #' regional_out <- readRDS(system.file(
@@ -510,7 +530,12 @@ regional_runtimes <- function(regional_output = NULL,
                               target_date = NULL,
                               return_output = FALSE) {
   if (is.null(target_folder) && is.null(regional_output)) {
-    stop("Either an output should be passed in or a target folder specified")
+    cli_abort(
+      c(
+        "i" = "Either an output should be passed in or a target folder
+        specified."
+      )
+    )
   }
   if (is.null(target_folder)) {
     futile.logger::flog.info(
diff --git a/tests/testthat/test-checks.R b/tests/testthat/test-checks.R
index 0fe68715a..9a3fc6525 100644
--- a/tests/testthat/test-checks.R
+++ b/tests/testthat/test-checks.R
@@ -144,5 +144,5 @@ test_that("check_reports_valid errors for bad 'secondary' specifications", {
 
 test_that("check_sparse_pmf_tail throws a warning as expected", {
   pmf <- c(0.4, 0.30, 0.20, 0.05,  0.049995, 4.5e-06, rep(1e-7, 5))
-  expect_warning(check_sparse_pmf_tail(pmf), "consecutive values smaller than")
+  expect_warning(check_sparse_pmf_tail(pmf), "PMF tail has")
 })
diff --git a/tests/testthat/test-dist_spec.R b/tests/testthat/test-dist_spec.R
index 5a0f7308e..593bad907 100644
--- a/tests/testthat/test-dist_spec.R
+++ b/tests/testthat/test-dist_spec.R
@@ -309,8 +309,11 @@ test_that("delay distributions can be specified in different ways", {
   )})
 
 test_that("get functions report errors", {
-  expect_error(get_parameters("test"), "only get parameters")
-  expect_error(get_distribution(Gamma(mean = 4, sd = 1), 2), "can't be greater")
+  expect_error(get_parameters("test"), "Object must be of class")
+  expect_error(
+    get_distribution(Gamma(mean = 4, sd = 1), 2),
+    "cannot be greater than the number of distributions"
+  )
   expect_error(get_pmf(Gamma(mean = 4, sd = 1)), "parametric")
   expect_error(
     get_parameters(NonParametric(c(0.1, 0.3, 0.2, 0.1, 0.1))),
diff --git a/tests/testthat/test-gp_opts.R b/tests/testthat/test-gp_opts.R
index d91c18a87..cd848b75c 100644
--- a/tests/testthat/test-gp_opts.R
+++ b/tests/testthat/test-gp_opts.R
@@ -28,14 +28,14 @@ test_that("gp_opts warns for uncommon Matern kernel orders", {
 
 test_that("gp_opts handles deprecated matern_type parameter", {
   lifecycle::expect_deprecated(gp_opts(matern_type = 5 / 2))
-  gp <- gp_opts(matern_type = 5 / 2)
+  gp <- suppressWarnings(gp_opts(matern_type = 5 / 2))
   expect_equal(gp$matern_order, 5 / 2)
 })
 
 test_that("gp_opts stops for incompatible matern_order and matern_type", {
   expect_error(
-    gp_opts(matern_order = 3 / 2, matern_type = 5 / 2),
-    "Incompatible `matern_order` and `matern_type`"
+    suppressWarnings(gp_opts(matern_order = 3 / 2, matern_type = 5 / 2)),
+    "must be the same, if both are supplied."
   )
 })
 
diff --git a/tests/testthat/test-rt_opts.R b/tests/testthat/test-rt_opts.R
index 2d415ff8e..0a39be027 100644
--- a/tests/testthat/test-rt_opts.R
+++ b/tests/testthat/test-rt_opts.R
@@ -38,9 +38,9 @@ test_that("rt_opts sets use_breakpoints to TRUE when rw > 0", {
 
 test_that("rt_opts throws error for invalid prior", {
   expect_error(rt_opts(prior = list(mean = 1)),
-               "prior must have both a mean and sd specified")
+               "must have both")
   expect_error(rt_opts(prior = list(sd = 1)),
-               "prior must have both a mean and sd specified")
+               "must have both")
 })
 
 test_that("rt_opts validates gp_on argument", {