From 4afb056241a96f2a8b0ba645cf2b5a463443cab3 Mon Sep 17 00:00:00 2001
From: Adam Howes <xwg3@cdc.gov>
Date: Thu, 14 Nov 2024 15:06:25 +0000
Subject: [PATCH] Issue #362: Reduce dependence on `cmdstanr` (#424)

* Switch default to rstan rather than cmdstanr

* Import BH

* Use match.arg as suggested by @jamesmbaazam

* Add RcppEigen to Suggests

* Change backend to cmdstanr in approximate inference vignette

* Use cmdstanr for Laplace test

* Remove erroneous output_dir arg

* drop backend entirely as an arg

* remove cmdstanr using vignettes from cran check

* update cmdstan install instructions

* update epidist docs

* update approx-inference vignette to make it clear we are using cmdstanr and not rstan

* switch ebola vignette over to cmdstanr

* update test fit returns

* update more fit uses in test

* fix helper function

* bring silence to laplace

* fix epidist

* get tests passing

* fix CRAN check

* use fewer cores

* drop additional deps

* drop duplicate prepping

* try and use linking with instead of imports

* Update vignettes/epidist.Rmd

Co-authored-by: Adam Howes <xwg3@cdc.gov>

* Update tests/testthat/test-int-latent_individual.R

Co-authored-by: Adam Howes <xwg3@cdc.gov>

* Update inst/generate_examples.R

Co-authored-by: Adam Howes <xwg3@cdc.gov>

---------

Co-authored-by: Sam <s.e.abbott12@gmail.com>
---
 .Rbuildignore                               |  4 +-
 DESCRIPTION                                 |  6 +-
 NAMESPACE                                   |  1 +
 R/{fit.R => epidist.R}                      | 11 +--
 R/utils.R                                   |  1 +
 README.Rmd                                  | 14 ++-
 inst/generate_examples.R                    | 10 +--
 man/epidist.Rd                              | 10 +--
 man/epidist.default.Rd                      |  9 +-
 man/epidist_family.Rd                       |  2 +-
 man/epidist_family_model.Rd                 |  2 +-
 man/epidist_family_model.default.Rd         |  2 +-
 man/epidist_family_prior.Rd                 |  2 +-
 man/epidist_family_prior.default.Rd         |  2 +-
 man/epidist_family_prior.lognormal.Rd       |  2 +-
 man/epidist_family_reparam.Rd               |  2 +-
 man/epidist_family_reparam.default.Rd       |  2 +-
 man/epidist_family_reparam.gamma.Rd         |  2 +-
 man/epidist_formula.Rd                      |  2 +-
 man/epidist_formula_model.default.Rd        |  2 +-
 man/epidist_model_prior.Rd                  |  2 +-
 man/epidist_model_prior.default.Rd          |  2 +-
 man/epidist_stancode.Rd                     |  2 +-
 man/epidist_stancode.default.Rd             |  2 +-
 tests/testthat/helper-functions.R           | 28 ++++++
 tests/testthat/setup.R                      | 29 ++++++
 tests/testthat/test-diagnostics.R           | 25 +++++-
 tests/testthat/test-family.R                |  3 -
 tests/testthat/test-formula.R               |  7 --
 tests/testthat/test-int-direct_model.R      | 17 ++--
 tests/testthat/test-int-latent_individual.R | 99 ++++-----------------
 tests/testthat/test-latent_gamma.R          | 18 ----
 tests/testthat/test-latent_lognormal.R      | 18 ----
 tests/testthat/test-postprocess.R           | 15 ----
 vignettes/approx-inference.Rmd              | 22 +++--
 vignettes/ebola.Rmd                         | 30 +++++--
 vignettes/epidist.Rmd                       | 22 ++---
 vignettes/faq.Rmd                           |  9 ++
 38 files changed, 202 insertions(+), 236 deletions(-)
 rename R/{fit.R => epidist.R} (82%)
 create mode 100644 tests/testthat/helper-functions.R

diff --git a/.Rbuildignore b/.Rbuildignore
index 5340855b4..55191a620 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -14,4 +14,6 @@
 ^doc$
 ^Meta$
 README.Rmd
-^pkgdown$
\ No newline at end of file
+^pkgdown$
+^vignettes/approx-inference\.Rmd$
+^vignettes/ebola\.Rmd$
diff --git a/DESCRIPTION b/DESCRIPTION
index 496571f62..21ceb1f3a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -31,7 +31,7 @@ Imports:
     stats,
     cli,
     checkmate,
-    rstan,
+    rstan (>= 2.26.0),
     dplyr
 Suggests:
     bookdown,
@@ -52,6 +52,10 @@ Suggests:
     patchwork,
     cmdstanr,
     priorsense
+LinkingTo: 
+    BH (>= 1.66.0),
+    Rcpp (>= 0.12.0),
+    RcppEigen (>= 0.3.3.3.0)
 Additional_repositories:
     https://stan-dev.r-universe.dev
 Config/Needs/website:
diff --git a/NAMESPACE b/NAMESPACE
index ad784791b..dba9c403e 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -68,3 +68,4 @@ importFrom(dplyr,filter)
 importFrom(dplyr,mutate)
 importFrom(dplyr,select)
 importFrom(stats,as.formula)
+importFrom(stats,setNames)
diff --git a/R/fit.R b/R/epidist.R
similarity index 82%
rename from R/fit.R
rename to R/epidist.R
index 4a8da7a2f..2496d4feb 100644
--- a/R/fit.R
+++ b/R/epidist.R
@@ -14,18 +14,15 @@
 #' @param prior One or more `brmsprior` objects created by [brms::set_prior()]
 #' or related functions. These priors are passed to [epidist_prior()] in the
 #' `prior` argument.
-#' @param backend Character string naming the package to use as the backend for
-#' fitting the Stan model. Options are `"rstan"` and `"cmdstanr"` (the default).
-#' This option is passed directly through to `fn`.
 #' @param fn The internal function to be called. By default this is
 #' [brms::brm()] which performs inference for the specified model. Other options
 #' are [brms::make_stancode()] which returns the Stan code for the specified
 #' model, or [brms::make_standata()] which returns the data passed to Stan.
 #' These two later options may be useful for model debugging and extensions.
-#' @param ... Additional arguments passed to method.
+#' @param ... Additional arguments passed to `fn` method.
 #' @family fit
 #' @export
-epidist <- function(data, formula, family, prior, backend, fn, ...) {
+epidist <- function(data, formula, family, prior, fn, ...) {
   UseMethod("epidist")
 }
 
@@ -38,7 +35,7 @@ epidist <- function(data, formula, family, prior, backend, fn, ...) {
 #' @export
 epidist.default <- function(data, formula = mu ~ 1,
                             family = "lognormal", prior = NULL,
-                            backend = "cmdstanr", fn = brms::brm, ...) {
+                            fn = brms::brm, ...) {
   epidist_validate_model(data)
   epidist_family <- epidist_family(data, family)
   epidist_formula <- epidist_formula(
@@ -52,7 +49,7 @@ epidist.default <- function(data, formula = mu ~ 1,
   )
   fit <- fn(
     formula = epidist_formula, family = epidist_family, prior = epidist_prior,
-    stanvars = epidist_stancode, backend = backend, data = data, ...
+    stanvars = epidist_stancode, data = data, ...
   )
   class(fit) <- c(class(fit), "epidist_fit")
   return(fit)
diff --git a/R/utils.R b/R/utils.R
index 6a881b09c..7d237e4e3 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -131,6 +131,7 @@
 #' @param new_names ...
 #' @param old_names ...
 #' @keywords internal
+#' @importFrom stats setNames
 .rename_columns <- function(df, new_names, old_names) {
   are_char <- is.character(new_names) & is.character(old_names)
   valid_new_names <- new_names[are_char]
diff --git a/README.Rmd b/README.Rmd
index 9887dac44..bcbdb32e1 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -76,14 +76,16 @@ remotes::install_github(
 
 </details>
 
-<details><summary>Installing CmdStan</summary>
+<details><summary>Installing CmdStan (optional)</summary>
 
-If you wish to do model fitting, you will need to install [CmdStan](https://mc-stan.org/users/interfaces/cmdstan), which also entails having a suitable C++ toolchain setup. We recommend using the [`cmdstanr` package](https://mc-stan.org/cmdstanr/). The Stan team provides instructions in the [_Getting started with
+By default `epidist` uses the `rstan` package for fitting models. If you wish to use the `cmdstanr` package instead, you will need to install [CmdStan](https://mc-stan.org/users/interfaces/cmdstan), which also entails having a suitable C++ toolchain setup. We recommend using the [`cmdstanr` package](https://mc-stan.org/cmdstanr/) to manage CmdStan. The Stan team provides instructions in the [_Getting started with
 `cmdstanr`_](https://mc-stan.org/cmdstanr/articles/cmdstanr.html) vignette, with other details and support at the [package site](https://mc-stan.org/cmdstanr/), but the brief version is:
 
 ```{r, eval = FALSE}
 # if you not yet installed `epidist`, or you installed it without `Suggests` dependencies
-install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
+install.packages(
+  "cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos"))
+)
 # once `cmdstanr` is installed:
 cmdstanr::install_cmdstan()
 ```
@@ -92,12 +94,6 @@ cmdstanr::install_cmdstan()
 
 </details>
 
-<!-- <details><summary>Alternative: Docker</summary> -->
-
-<!-- We also provide a [Docker](https://www.docker.com/get-started/) image with [all dependencies installed](https://github.com/orgs/epinowcast/packages/container/package/epidist). You can use this image to run `epidist` without installing dependencies. -->
-
-<!-- </details> -->
-
 ## Resources
 
 <details><summary>Organisation Website</summary>
diff --git a/inst/generate_examples.R b/inst/generate_examples.R
index beaf76f94..dd1fdb3a9 100644
--- a/inst/generate_examples.R
+++ b/inst/generate_examples.R
@@ -1,12 +1,4 @@
 source("tests/testthat/setup.R")
-set.seed(1)
-prep_obs <- as_latent_individual(sim_obs)
-fit <- epidist(prep_obs, seed = 1)
+
 saveRDS(fit, "inst/extdata/fit.rds")
-prep_obs_gamma <- as_latent_individual(sim_obs_gamma)
-fit_gamma <- epidist(
-  prep_obs_gamma,
-  family = stats::Gamma(link = "log"),
-  seed = 1
-)
 saveRDS(fit_gamma, "inst/extdata/fit_gamma.rds")
diff --git a/man/epidist.Rd b/man/epidist.Rd
index 89a9d661a..f6663d66c 100644
--- a/man/epidist.Rd
+++ b/man/epidist.Rd
@@ -1,10 +1,10 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/fit.R
+% Please edit documentation in R/epidist.R
 \name{epidist}
 \alias{epidist}
 \title{Fit epidemiological delay distributions using a \code{brms} interface}
 \usage{
-epidist(data, formula, family, prior, backend, fn, ...)
+epidist(data, formula, family, prior, fn, ...)
 }
 \arguments{
 \item{data}{A \code{data.frame} containing line list data.}
@@ -25,17 +25,13 @@ see \code{\link[=brmsfamily]{brmsfamily()}}.}
 or related functions. These priors are passed to \code{\link[=epidist_prior]{epidist_prior()}} in the
 \code{prior} argument.}
 
-\item{backend}{Character string naming the package to use as the backend for
-fitting the Stan model. Options are \code{"rstan"} and \code{"cmdstanr"} (the default).
-This option is passed directly through to \code{fn}.}
-
 \item{fn}{The internal function to be called. By default this is
 \code{\link[brms:brm]{brms::brm()}} which performs inference for the specified model. Other options
 are \code{\link[brms:stancode]{brms::make_stancode()}} which returns the Stan code for the specified
 model, or \code{\link[brms:standata]{brms::make_standata()}} which returns the data passed to Stan.
 These two later options may be useful for model debugging and extensions.}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 Fit epidemiological delay distributions using a \code{brms} interface
diff --git a/man/epidist.default.Rd b/man/epidist.default.Rd
index dc9a5a2b3..18353ca08 100644
--- a/man/epidist.default.Rd
+++ b/man/epidist.default.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/fit.R
+% Please edit documentation in R/epidist.R
 \name{epidist.default}
 \alias{epidist.default}
 \title{Default method used for interface using \code{brms}}
@@ -9,7 +9,6 @@
   formula = mu ~ 1,
   family = "lognormal",
   prior = NULL,
-  backend = "cmdstanr",
   fn = brms::brm,
   ...
 )
@@ -33,17 +32,13 @@ see \code{\link[=brmsfamily]{brmsfamily()}}.}
 or related functions. These priors are passed to \code{\link[=epidist_prior]{epidist_prior()}} in the
 \code{prior} argument.}
 
-\item{backend}{Character string naming the package to use as the backend for
-fitting the Stan model. Options are \code{"rstan"} and \code{"cmdstanr"} (the default).
-This option is passed directly through to \code{fn}.}
-
 \item{fn}{The internal function to be called. By default this is
 \code{\link[brms:brm]{brms::brm()}} which performs inference for the specified model. Other options
 are \code{\link[brms:stancode]{brms::make_stancode()}} which returns the Stan code for the specified
 model, or \code{\link[brms:standata]{brms::make_standata()}} which returns the data passed to Stan.
 These two later options may be useful for model debugging and extensions.}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 Default method used for interface using \code{brms}
diff --git a/man/epidist_family.Rd b/man/epidist_family.Rd
index af0200a8c..bb4bdc615 100644
--- a/man/epidist_family.Rd
+++ b/man/epidist_family.Rd
@@ -15,7 +15,7 @@ users to specify the link function to be applied on the response variable.
 If not specified, default links are used. For details of supported families
 see \code{\link[=brmsfamily]{brmsfamily()}}.}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 This function is used within \code{\link[=epidist]{epidist()}} to create a model specific custom
diff --git a/man/epidist_family_model.Rd b/man/epidist_family_model.Rd
index 747057069..08608529e 100644
--- a/man/epidist_family_model.Rd
+++ b/man/epidist_family_model.Rd
@@ -15,7 +15,7 @@ epidist_formula_model(data, formula, ...)
 \item{family}{Output of a call to \code{brms::brmsfamily()} with additional
 information as provided by \code{.add_dpar_info()}}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 
 \item{formula}{An object of class \link[stats:formula]{stats::formula} or \link[brms:brmsformula]{brms::brmsformula}
 (or one that can be coerced to those classes). A symbolic description of the
diff --git a/man/epidist_family_model.default.Rd b/man/epidist_family_model.default.Rd
index 51fc646ea..c7cde1af9 100644
--- a/man/epidist_family_model.default.Rd
+++ b/man/epidist_family_model.default.Rd
@@ -12,7 +12,7 @@
 \item{family}{Output of a call to \code{brms::brmsfamily()} with additional
 information as provided by \code{.add_dpar_info()}}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 Default method for defining a model specific family
diff --git a/man/epidist_family_prior.Rd b/man/epidist_family_prior.Rd
index 9912f0312..82ce122dd 100644
--- a/man/epidist_family_prior.Rd
+++ b/man/epidist_family_prior.Rd
@@ -13,7 +13,7 @@ users to specify the link function to be applied on the response variable.
 If not specified, default links are used. For details of supported families
 see \code{\link[=brmsfamily]{brmsfamily()}}.}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 This function contains \code{brms} prior distributions which are specific to
diff --git a/man/epidist_family_prior.default.Rd b/man/epidist_family_prior.default.Rd
index 19e20be81..eb2677626 100644
--- a/man/epidist_family_prior.default.Rd
+++ b/man/epidist_family_prior.default.Rd
@@ -19,7 +19,7 @@ model to be fitted. A formula must be provided for the distributional
 parameter \code{mu}, and may optionally be provided for other distributional
 parameters.}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 By default, we do not return any family specific prior distributions.
diff --git a/man/epidist_family_prior.lognormal.Rd b/man/epidist_family_prior.lognormal.Rd
index 1eedada97..2708ba9ff 100644
--- a/man/epidist_family_prior.lognormal.Rd
+++ b/man/epidist_family_prior.lognormal.Rd
@@ -19,7 +19,7 @@ model to be fitted. A formula must be provided for the distributional
 parameter \code{mu}, and may optionally be provided for other distributional
 parameters.}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 We suggest priors to overwrite the \code{brms} defaults for the lognormal family.
diff --git a/man/epidist_family_reparam.Rd b/man/epidist_family_reparam.Rd
index 23b3b08d2..27ecba0a4 100644
--- a/man/epidist_family_reparam.Rd
+++ b/man/epidist_family_reparam.Rd
@@ -13,7 +13,7 @@ users to specify the link function to be applied on the response variable.
 If not specified, default links are used. For details of supported families
 see \code{\link[=brmsfamily]{brmsfamily()}}.}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 Reparameterise an \code{epidist} family to align \code{brms} and Stan
diff --git a/man/epidist_family_reparam.default.Rd b/man/epidist_family_reparam.default.Rd
index 0feccf703..4fc4531de 100644
--- a/man/epidist_family_reparam.default.Rd
+++ b/man/epidist_family_reparam.default.Rd
@@ -13,7 +13,7 @@ users to specify the link function to be applied on the response variable.
 If not specified, default links are used. For details of supported families
 see \code{\link[=brmsfamily]{brmsfamily()}}.}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 Default method for families which do not require a reparameterisation
diff --git a/man/epidist_family_reparam.gamma.Rd b/man/epidist_family_reparam.gamma.Rd
index 2ac5dba69..70acfa454 100644
--- a/man/epidist_family_reparam.gamma.Rd
+++ b/man/epidist_family_reparam.gamma.Rd
@@ -13,7 +13,7 @@ users to specify the link function to be applied on the response variable.
 If not specified, default links are used. For details of supported families
 see \code{\link[=brmsfamily]{brmsfamily()}}.}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 Reparameterisation for the gamma family
diff --git a/man/epidist_formula.Rd b/man/epidist_formula.Rd
index 7b8ec6587..a1a94cad4 100644
--- a/man/epidist_formula.Rd
+++ b/man/epidist_formula.Rd
@@ -18,7 +18,7 @@ model to be fitted. A formula must be provided for the distributional
 parameter \code{mu}, and may optionally be provided for other distributional
 parameters.}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 This function is used within \code{\link[=epidist]{epidist()}} to create the formula object passed
diff --git a/man/epidist_formula_model.default.Rd b/man/epidist_formula_model.default.Rd
index 41e8419ae..7c8e7bbb2 100644
--- a/man/epidist_formula_model.default.Rd
+++ b/man/epidist_formula_model.default.Rd
@@ -15,7 +15,7 @@ model to be fitted. A formula must be provided for the distributional
 parameter \code{mu}, and may optionally be provided for other distributional
 parameters.}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 Default method for defining a model specific formula
diff --git a/man/epidist_model_prior.Rd b/man/epidist_model_prior.Rd
index ef04af4a6..df28ea733 100644
--- a/man/epidist_model_prior.Rd
+++ b/man/epidist_model_prior.Rd
@@ -9,7 +9,7 @@ epidist_model_prior(data, ...)
 \arguments{
 \item{data}{A \code{data.frame} containing line list data.}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 This function contains \code{brms} prior distributions which are specific to
diff --git a/man/epidist_model_prior.default.Rd b/man/epidist_model_prior.default.Rd
index 841a66e9d..27035bc07 100644
--- a/man/epidist_model_prior.default.Rd
+++ b/man/epidist_model_prior.default.Rd
@@ -15,7 +15,7 @@ model to be fitted. A formula must be provided for the distributional
 parameter \code{mu}, and may optionally be provided for other distributional
 parameters.}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 By default, we do not return any model specific prior distributions.
diff --git a/man/epidist_stancode.Rd b/man/epidist_stancode.Rd
index d86c8010d..31a58e5b5 100644
--- a/man/epidist_stancode.Rd
+++ b/man/epidist_stancode.Rd
@@ -9,7 +9,7 @@ epidist_stancode(data, ...)
 \arguments{
 \item{data}{A \code{data.frame} containing line list data.}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 This function is used within \code{\link[=epidist]{epidist()}} to create any custom Stan code which
diff --git a/man/epidist_stancode.default.Rd b/man/epidist_stancode.default.Rd
index 1aa781741..ef2f21f18 100644
--- a/man/epidist_stancode.default.Rd
+++ b/man/epidist_stancode.default.Rd
@@ -9,7 +9,7 @@
 \arguments{
 \item{data}{A \code{data.frame} containing line list data.}
 
-\item{...}{Additional arguments passed to method.}
+\item{...}{Additional arguments passed to \code{fn} method.}
 }
 \description{
 Default method for defining model specific Stan code
diff --git a/tests/testthat/helper-functions.R b/tests/testthat/helper-functions.R
new file mode 100644
index 000000000..7a1a5fd59
--- /dev/null
+++ b/tests/testthat/helper-functions.R
@@ -0,0 +1,28 @@
+on_ci <- function() {
+  isTRUE(as.logical(Sys.getenv("CI")))
+}
+
+not_on_cran <- function() {
+  identical(Sys.getenv("NOT_CRAN"), "true")
+}
+
+skip_on_local <- function() {
+  if (on_ci()) {
+    return(invisible(TRUE))
+  }
+  testthat::skip("Not on CI")
+}
+
+as_string_formula <- function(formula) {
+  form <- paste(deparse(formula), collapse = " ")
+  form <- gsub("\\s+", " ", form, perl = FALSE)
+  return(form)
+}
+
+extract_normal_parameters_brms <- function(prior) {
+  pattern <- "normal\\(([^,]+), ([^\\)]+)\\)"
+  match <- regmatches(prior, regexec(pattern, prior))
+  mean <- as.numeric(match[[1]][2])
+  sd <- as.numeric(match[[1]][3])
+  return(list(mean = mean, sd = sd))
+}
diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R
index 2a73f46a1..da29780cf 100644
--- a/tests/testthat/setup.R
+++ b/tests/testthat/setup.R
@@ -82,3 +82,32 @@ sim_obs_sex <- dplyr::bind_rows(sim_obs_sex_m, sim_obs_sex_f) |>
 
 # Temporary solution for classing time data
 sim_obs_sex <- as_epidist_linelist_time(sim_obs_sex)
+
+prep_obs <- as_latent_individual(sim_obs)
+prep_direct_obs <- as_direct_model(sim_obs)
+prep_obs_gamma <- as_latent_individual(sim_obs_gamma)
+prep_obs_sex <- as_latent_individual(sim_obs_sex)
+
+if (not_on_cran()) {
+  set.seed(1)
+  fit <- epidist(
+    data = prep_obs, seed = 1, chains = 2, cores = 2, silent = 2, refresh = 0,
+    backend = "cmdstanr"
+  )
+  fit_rstan <- epidist(
+    data = prep_obs, seed = 1, chains = 2, cores = 2, silent = 2, refresh = 0
+  )
+
+  fit_gamma <- epidist(
+    data = prep_obs_gamma, family = stats::Gamma(link = "log"),
+    seed = 1, chains = 2, cores = 2, silent = 2, refresh = 0,
+    backend = "cmdstanr"
+  )
+
+  fit_sex <- epidist(
+    data = prep_obs_sex,
+    formula = brms::bf(mu ~ 1 + sex, sigma ~ 1 + sex),
+    seed = 1, silent = 2, refresh = 0,
+    cores = 2, chains = 2, backend = "cmdstanr"
+  )
+}
diff --git a/tests/testthat/test-diagnostics.R b/tests/testthat/test-diagnostics.R
index e8e57496a..c1c127f07 100644
--- a/tests/testthat/test-diagnostics.R
+++ b/tests/testthat/test-diagnostics.R
@@ -1,8 +1,6 @@
 test_that("epidist_diagnostics", { # nolint: line_length_linter.
   skip_on_cran()
   set.seed(1)
-  prep_obs <- as_latent_individual(sim_obs)
-  fit <- epidist(data = prep_obs, seed = 1)
   diag <- epidist_diagnostics(fit)
   expected_names <- c(
     "time", "samples", "max_rhat", "divergent_transitions",
@@ -22,10 +20,31 @@ test_that("epidist_diagnostics", { # nolint: line_length_linter.
   expect_gt(diag$per_at_max_treedepth, 0)
 })
 
+test_that("epidist_diagnostics gives the same results for cmdstanr and rstan", {
+  skip_on_cran()
+  set.seed(1)
+  diag_cmdstanr <- epidist_diagnostics(fit)
+  diag_rstan <- epidist_diagnostics(fit_rstan)
+  expect_equal(colnames(diag_cmdstanr), colnames(diag_rstan))
+  expect_gt(diag_rstan$time, 0)
+  expect_gt(diag_rstan$samples, 0)
+  expect_gt(diag_rstan$max_rhat, 0.9)
+  expect_lt(diag_rstan$max_rhat, 1.1)
+  expect_gte(diag_rstan$divergent_transitions, 0)
+  expect_lt(diag_rstan$divergent_transitions, diag_rstan$samples)
+  expect_lt(diag_rstan$max_treedepth, 12)
+  expect_lte(diag_rstan$no_at_max_treedepth, diag_rstan$samples)
+  expect_lte(diag_rstan$per_at_max_treedepth, 1)
+  expect_gt(diag_rstan$per_at_max_treedepth, 0)
+})
+
 test_that("epidist_diagnostics gives an error when passed model fit using the Laplace algorithm", { # nolint: line_length_linter.
   skip_on_cran()
   set.seed(1)
   prep_obs <- as_latent_individual(sim_obs)
-  fit_laplace <- epidist(data = prep_obs, seed = 1, algorithm = "laplace")
+  fit_laplace <- epidist(
+    data = prep_obs, seed = 1, algorithm = "laplace", backend = "cmdstanr",
+    refresh = 0, silent = 2, show_messages = FALSE
+  )
   expect_error(epidist_diagnostics(fit_laplace))
 })
diff --git a/tests/testthat/test-family.R b/tests/testthat/test-family.R
index 0da6ae838..64d68810c 100644
--- a/tests/testthat/test-family.R
+++ b/tests/testthat/test-family.R
@@ -1,6 +1,3 @@
-prep_obs <- as_latent_individual(sim_obs)
-prep_obs_gamma <- as_latent_individual(sim_obs_gamma)
-
 test_that("epidist_family with default settings produces an object of the right class", { # nolint: line_length_linter.
   family <- epidist_family(prep_obs)
   expect_s3_class(family, "customfamily")
diff --git a/tests/testthat/test-formula.R b/tests/testthat/test-formula.R
index 89bf04086..d6410ae00 100644
--- a/tests/testthat/test-formula.R
+++ b/tests/testthat/test-formula.R
@@ -1,14 +1,7 @@
-prep_obs <- as_latent_individual(sim_obs)
 prep_obs_gamma <- as_latent_individual(sim_obs_gamma)
 
 family_lognormal <- epidist_family(prep_obs, family = brms::lognormal())
 
-as_string_formula <- function(formula) {
-  form <- paste(deparse(formula), collapse = " ")
-  form <- gsub("\\s+", " ", form, perl = FALSE)
-  return(form)
-}
-
 test_that("epidist_formula with default settings produces a brmsformula with the correct intercept only formula", { # nolint: line_length_linter.
   form <- epidist_formula(
     prep_obs, family = family_lognormal, formula = mu ~ 1
diff --git a/tests/testthat/test-int-direct_model.R b/tests/testthat/test-int-direct_model.R
index b366b938f..f1f49aefd 100644
--- a/tests/testthat/test-int-direct_model.R
+++ b/tests/testthat/test-int-direct_model.R
@@ -4,20 +4,16 @@
 # varying the input seed. Test failure at an unusually high rate does suggest
 # a potential code issue.
 
-prep_obs <- as_direct_model(sim_obs)
-
-test_that("epidist.epidist_direct_model Stan code has no syntax errors and compiles in the default case", { # nolint: line_length_linter.
+test_that("epidist.epidist_direct_model Stan code has no syntax errors in the default case", { # nolint: line_length_linter.
   skip_on_cran()
   stancode <- epidist(
-    data = prep_obs,
-    fn = brms::make_stancode,
-    output_dir = fs::dir_create(tempfile())
+    data = prep_direct_obs,
+    fn = brms::make_stancode
   )
   mod <- cmdstanr::cmdstan_model(
     stan_file = cmdstanr::write_stan_file(stancode), compile = FALSE
   )
   expect_true(mod$check_syntax())
-  expect_no_error(mod$compile())
 })
 
 test_that("epidist.epidist_direct_model fits and the MCMC converges in the default case", { # nolint: line_length_linter.
@@ -25,10 +21,11 @@ test_that("epidist.epidist_direct_model fits and the MCMC converges in the defau
   skip_on_cran()
   set.seed(1)
   fit <- epidist(
-    data = prep_obs,
+    data = prep_direct_obs,
     seed = 1,
-    silent = 2,
-    output_dir = fs::dir_create(tempfile())
+    silent = 2, refresh = 0,
+    cores = 2,
+    chains = 2
   )
   expect_s3_class(fit, "brmsfit")
   expect_s3_class(fit, "epidist_fit")
diff --git a/tests/testthat/test-int-latent_individual.R b/tests/testthat/test-int-latent_individual.R
index 236981ba3..04c249583 100644
--- a/tests/testthat/test-int-latent_individual.R
+++ b/tests/testthat/test-int-latent_individual.R
@@ -4,31 +4,18 @@
 # varying the input seed. Test failure at an unusually high rate does suggest
 # a potential code issue.
 
-prep_obs <- as_latent_individual(sim_obs)
-prep_obs_gamma <- as_latent_individual(sim_obs_gamma)
-
-test_that("epidist.epidist_latent_individual Stan code has no syntax errors and compiles in the default case", { # nolint: line_length_linter.
+test_that("epidist.epidist_latent_individual Stan code has no syntax errors in the default case", { # nolint: line_length_linter.
   skip_on_cran()
   stancode <- epidist(
     data = prep_obs,
-    fn = brms::make_stancode,
-    output_dir = fs::dir_create(tempfile())
+    fn = brms::make_stancode
   )
   mod <- cmdstanr::cmdstan_model(
     stan_file = cmdstanr::write_stan_file(stancode), compile = FALSE
   )
   expect_true(mod$check_syntax())
-  expect_no_error(mod$compile())
 })
 
-extract_normal_parameters_brms <- function(prior) {
-  pattern <- "normal\\(([^,]+), ([^\\)]+)\\)"
-  match <- regmatches(prior, regexec(pattern, prior))
-  mean <- as.numeric(match[[1]][2])
-  sd <- as.numeric(match[[1]][3])
-  return(list(mean = mean, sd = sd))
-}
-
 test_that("epidist.epidist_latent_individual samples from the prior according to marginal Kolmogorov-Smirnov tests in the default case.", { # nolint: line_length_linter.
   # Note: this test is stochastic. See note at the top of this script
   skip_on_cran()
@@ -38,8 +25,8 @@ test_that("epidist.epidist_latent_individual samples from the prior according to
     fn = brms::brm,
     sample_prior = "only",
     seed = 1,
-    silent = 2,
-    output_dir = fs::dir_create(tempfile())
+    silent = 2, refresh = 0,
+    cores = 2
   )
   pred <- predict_delay_parameters(prior_samples)
   family <- brms::lognormal()
@@ -69,13 +56,6 @@ test_that("epidist.epidist_latent_individual samples from the prior according to
 test_that("epidist.epidist_latent_individual fits and the MCMC converges in the default case", { # nolint: line_length_linter.
   # Note: this test is stochastic. See note at the top of this script
   skip_on_cran()
-  set.seed(1)
-  fit <- epidist(
-    data = prep_obs,
-    seed = 1,
-    silent = 2,
-    output_dir = fs::dir_create(tempfile())
-  )
   expect_s3_class(fit, "brmsfit")
   expect_s3_class(fit, "epidist_fit")
   expect_convergence(fit)
@@ -89,8 +69,9 @@ test_that("epidist.epidist_latent_individual fits, the MCMC converges, and the d
     data = prep_obs,
     formula = brms::bf(mu ~ 1, sigma = 1),
     seed = 1,
-    silent = 2,
-    output_dir = fs::dir_create(tempfile())
+    silent = 2, refresh = 0,
+    cores = 2,
+    chains = 2
   )
   expect_s3_class(fit_constant, "brmsfit")
   expect_s3_class(fit_constant, "epidist_fit")
@@ -99,7 +80,7 @@ test_that("epidist.epidist_latent_individual fits, the MCMC converges, and the d
   expect_true(all(sigma == 1))
 })
 
-test_that("epidist.epidist_latent_individual Stan code has no syntax errors and compiles with lognormal family as a string", { # nolint: line_length_linter.
+test_that("epidist.epidist_latent_individual Stan code has no syntax errors", { # nolint: line_length_linter.
   # Note: this test is stochastic. See note at the top of this script
   skip_on_cran()
   set.seed(1)
@@ -107,27 +88,19 @@ test_that("epidist.epidist_latent_individual Stan code has no syntax errors and
     data = prep_obs,
     family = "lognormal",
     seed = 1,
-    silent = 2,
-    output_dir = fs::dir_create(tempfile()),
+    silent = 2, refresh = 0,
     fn = brms::make_stancode
   )
   mod_string <- cmdstanr::cmdstan_model(
     stan_file = cmdstanr::write_stan_file(stancode_string), compile = FALSE
   )
   expect_true(mod_string$check_syntax())
-  expect_no_error(mod_string$compile())
 })
 
 test_that("epidist.epidist_latent_individual recovers the simulation settings for the delay distribution in the default case", { # nolint: line_length_linter.
   # Note: this test is stochastic. See note at the top of this script
   skip_on_cran()
   set.seed(1)
-  fit <- epidist(
-    data = prep_obs,
-    seed = 1,
-    silent = 2,
-    output_dir = fs::dir_create(tempfile())
-  )
   pred <- predict_delay_parameters(fit)
   # Unclear the extent to which we should expect parameter recovery here
   expect_equal(mean(pred$mu), meanlog, tolerance = 0.1)
@@ -140,7 +113,7 @@ test_that("epidist.epidist_latent_individual Stan code has no syntax errors and
     data = prep_obs_gamma,
     family = stats::Gamma(link = "log"),
     formula = mu ~ 1,
-    output_dir = fs::dir_create(tempfile()),
+    cores = 2,
     fn = brms::make_stancode
   )
   mod_gamma <- cmdstanr::cmdstan_model(
@@ -154,14 +127,6 @@ test_that("epidist.epidist_latent_individual fits and the MCMC converges in the
   # Note: this test is stochastic. See note at the top of this script
   skip_on_cran()
   set.seed(1)
-  fit_gamma <- epidist(
-    data = prep_obs_gamma,
-    family = stats::Gamma(link = "log"),
-    formula = mu ~ 1,
-    seed = 1,
-    silent = 2,
-    output_dir = fs::dir_create(tempfile())
-  )
   expect_s3_class(fit_gamma, "brmsfit")
   expect_s3_class(fit_gamma, "epidist_fit")
   expect_convergence(fit_gamma)
@@ -171,14 +136,6 @@ test_that("epidist.epidist_latent_individual recovers the simulation settings fo
   # Note: this test is stochastic. See note at the top of this script
   skip_on_cran()
   set.seed(1)
-  fit_gamma <- epidist(
-    data = prep_obs_gamma,
-    family = stats::Gamma(link = "log"),
-    formula = mu ~ 1,
-    seed = 1,
-    silent = 2,
-    output_dir = fs::dir_create(tempfile())
-  )
   draws_gamma <- posterior::as_draws_df(fit_gamma$fit)
   draws_gamma_mu <- exp(draws_gamma$Intercept)
   draws_gamma_shape <- exp(draws_gamma$Intercept_shape)
@@ -192,51 +149,27 @@ test_that("epidist.epidist_latent_individual recovers the simulation settings fo
   expect_lte(quantile_shape, 0.975)
 })
 
-test_that("epidist.epidist_latent_individual Stan code has no syntax errors and compiles for an alternative formula", { # nolint: line_length_linter.
+test_that("epidist.epidist_latent_individual Stan code has no syntax errors for an alternative formula", { # nolint: line_length_linter.
   skip_on_cran()
-  prep_obs$sex <- rbinom(n = nrow(prep_obs), size = 1, prob = 0.5)
   stancode_sex <- epidist(
-    data = prep_obs,
+    data = prep_obs_sex,
     formula = brms::bf(mu ~ 1 + sex, sigma ~ 1 + sex),
     fn = brms::make_stancode,
-    output_dir = fs::dir_create(tempfile())
+    cores = 2
   )
   mod_sex <- cmdstanr::cmdstan_model(
     stan_file = cmdstanr::write_stan_file(stancode_sex), compile = FALSE
   )
   expect_true(mod_sex$check_syntax())
-  expect_no_error(mod_sex$compile())
 })
 
-test_that("epidist.epidist_latent_individual recovers no sex effect when none is simulated", { # nolint: line_length_linter.
+test_that("epidist.epidist_latent_individual recovers a sex effect", { # nolint: line_length_linter.
   # Note: this test is stochastic. See note at the top of this script
   skip_on_cran()
   set.seed(1)
-  prep_obs$sex <- rbinom(n = nrow(prep_obs), size = 1, prob = 0.5)
-  fit_sex <- epidist(
-    data = prep_obs,
-    formula = brms::bf(mu ~ 1 + sex, sigma ~ 1 + sex),
-    seed = 1,
-    silent = 2,
-    output_dir = fs::dir_create(tempfile())
-  )
   draws <- posterior::as_draws_df(fit_sex$fit)
-  expect_equal(mean(draws$b_sex), 0, tolerance = 0.2)
-  expect_equal(mean(draws$b_sigma_sex), 0, tolerance = 0.2)
-})
-
-test_that("epidist.epidist_latent_individual fits and the MCMC converges for an alternative formula", { # nolint: line_length_linter.
-  # Note: this test is stochastic. See note at the top of this script
-  skip_on_cran()
-  set.seed(1)
-  prep_obs$sex <- rbinom(n = nrow(prep_obs), size = 1, prob = 0.5)
-  fit_sex <- epidist(
-    data = prep_obs,
-    formula = brms::bf(mu ~ 1 + sex, sigma ~ 1 + sex),
-    seed = 1,
-    silent = 2,
-    output_dir = fs::dir_create(tempfile())
-  )
+  expect_equal(mean(draws$b_sex), -0.73, tolerance = 0.2)
+  expect_equal(mean(draws$b_sigma_sex), 0.43, tolerance = 0.2)
   expect_s3_class(fit_sex, "brmsfit")
   expect_s3_class(fit_sex, "epidist_fit")
   expect_convergence(fit_sex)
diff --git a/tests/testthat/test-latent_gamma.R b/tests/testthat/test-latent_gamma.R
index 0a51dd69d..8ca99174b 100644
--- a/tests/testthat/test-latent_gamma.R
+++ b/tests/testthat/test-latent_gamma.R
@@ -1,8 +1,5 @@
 test_that("posterior_predict_latent_gamma outputs positive integers with length equal to draws", { # nolint: line_length_linter.
   skip_on_cran()
-  fit_gamma <- readRDS(
-    system.file("extdata/fit_gamma.rds", package = "epidist")
-  )
   prep <- brms::prepare_predictions(fit_gamma)
   i <- 1
   pred_i <- posterior_predict_latent_gamma(i = i, prep)
@@ -13,9 +10,6 @@ test_that("posterior_predict_latent_gamma outputs positive integers with length
 
 test_that("posterior_predict_latent_gamma errors for i out of bounds", { # nolint: line_length_linter.
   skip_on_cran()
-  fit_gamma <- readRDS(
-    system.file("extdata/fit_gamma.rds", package = "epidist")
-  )
   prep <- brms::prepare_predictions(fit_gamma)
   i_out_of_bounds <- length(prep$data$Y) + 1
   expect_error(posterior_predict_latent_gamma(i = i_out_of_bounds, prep))
@@ -23,9 +17,6 @@ test_that("posterior_predict_latent_gamma errors for i out of bounds", { # nolin
 
 test_that("posterior_predict_latent_gamma can generate predictions with no censoring", { # nolint: line_length_linter.
   skip_on_cran()
-  fit_gamma <- readRDS(
-    system.file("extdata/fit_gamma.rds", package = "epidist")
-  )
   draws <- data.frame(relative_obs_time = 1000, pwindow = 0, swindow = 0) |>
     tidybayes::add_predicted_draws(fit_gamma, ndraws = 100)
   expect_equal(draws$.draw, 1:100)
@@ -36,9 +27,6 @@ test_that("posterior_predict_latent_gamma can generate predictions with no censo
 
 test_that("posterior_predict_latent_gamma predicts delays for which the data is in the 95% credible interval", { # nolint: line_length_linter.
   skip_on_cran()
-  fit_gamma <- readRDS(
-    system.file("extdata/fit_gamma.rds", package = "epidist")
-  )
   prep <- brms::prepare_predictions(fit_gamma)
   prep$ndraws <- 1000 # Down from the 4000 for time saving
   q <- purrr::map_vec(seq_along(prep$data$Y), function(i) {
@@ -57,9 +45,6 @@ test_that("posterior_predict_latent_gamma predicts delays for which the data is
 
 test_that("posterior_epred_latent_gamma creates a array of non-negative numbers with the correct dimensions", { # nolint: line_length_linter.
   skip_on_cran()
-  fit_gamma <- readRDS(
-    system.file("extdata/fit_gamma.rds", package = "epidist")
-  )
   prep <- brms::prepare_predictions(fit_gamma)
   epred <- posterior_epred_latent_gamma(prep)
   expect_setequal(class(epred), c("matrix", "array"))
@@ -70,9 +55,6 @@ test_that("posterior_epred_latent_gamma creates a array of non-negative numbers
 
 test_that("log_lik_latent_gamma produces a vector with length ndraws of finite non-NA numbers", { # nolint: line_length_linter.
   skip_on_cran()
-  fit_gamma <- readRDS(
-    system.file("extdata/fit_gamma.rds", package = "epidist")
-  )
   prep <- brms::prepare_predictions(fit_gamma)
   i <- 1
   log_lik <- log_lik_latent_gamma(i, prep)
diff --git a/tests/testthat/test-latent_lognormal.R b/tests/testthat/test-latent_lognormal.R
index 4b35f1a61..351617701 100644
--- a/tests/testthat/test-latent_lognormal.R
+++ b/tests/testthat/test-latent_lognormal.R
@@ -1,8 +1,5 @@
 test_that("posterior_predict_latent_lognormal outputs positive integers with length equal to draws", { # nolint: line_length_linter.
   skip_on_cran()
-  fit <- readRDS(
-    system.file("extdata/fit.rds", package = "epidist")
-  )
   prep <- brms::prepare_predictions(fit)
   i <- 1
   pred_i <- posterior_predict_latent_lognormal(i = i, prep)
@@ -13,9 +10,6 @@ test_that("posterior_predict_latent_lognormal outputs positive integers with len
 
 test_that("posterior_predict_latent_lognormal errors for i out of bounds", { # nolint: line_length_linter.
   skip_on_cran()
-  fit <- readRDS(
-    system.file("extdata/fit.rds", package = "epidist")
-  )
   prep <- brms::prepare_predictions(fit)
   i_out_of_bounds <- length(prep$data$Y) + 1
   expect_error(posterior_predict_latent_lognormal(i = i_out_of_bounds, prep))
@@ -23,9 +17,6 @@ test_that("posterior_predict_latent_lognormal errors for i out of bounds", { # n
 
 test_that("posterior_predict_latent_lognormal can generate predictions with no censoring", { # nolint: line_length_linter.
   skip_on_cran()
-  fit <- readRDS(
-    system.file("extdata/fit.rds", package = "epidist")
-  )
   draws <- data.frame(relative_obs_time = 1000, pwindow = 0, swindow = 0) |>
     tidybayes::add_predicted_draws(fit, ndraws = 100)
   expect_equal(draws$.draw, 1:100)
@@ -36,9 +27,6 @@ test_that("posterior_predict_latent_lognormal can generate predictions with no c
 
 test_that("posterior_predict_latent_lognormal predicts delays for which the data is in the 95% credible interval", { # nolint: line_length_linter.
   skip_on_cran()
-  fit <- readRDS(
-    system.file("extdata/fit.rds", package = "epidist")
-  )
   prep <- brms::prepare_predictions(fit)
   prep$ndraws <- 1000 # Down from the 4000 for time saving
   q <- purrr::map_vec(seq_along(prep$data$Y), function(i) {
@@ -57,9 +45,6 @@ test_that("posterior_predict_latent_lognormal predicts delays for which the data
 
 test_that("posterior_epred_latent_lognormal creates a array of non-negative numbers with the correct dimensions", { # nolint: line_length_linter.
   skip_on_cran()
-  fit <- readRDS(
-    system.file("extdata/fit.rds", package = "epidist")
-  )
   prep <- brms::prepare_predictions(fit)
   epred <- posterior_epred_latent_lognormal(prep)
   expect_setequal(class(epred), c("matrix", "array"))
@@ -70,9 +55,6 @@ test_that("posterior_epred_latent_lognormal creates a array of non-negative numb
 
 test_that("log_lik_latent_lognormal produces a vector with length ndraws of finite non-NA numbers", { # nolint: line_length_linter.
   skip_on_cran()
-  fit <- readRDS(
-    system.file("extdata/fit.rds", package = "epidist")
-  )
   prep <- brms::prepare_predictions(fit)
   i <- 1
   log_lik <- log_lik_latent_lognormal(i, prep)
diff --git a/tests/testthat/test-postprocess.R b/tests/testthat/test-postprocess.R
index 0b9741ba9..8c0d205cd 100644
--- a/tests/testthat/test-postprocess.R
+++ b/tests/testthat/test-postprocess.R
@@ -1,13 +1,6 @@
 test_that("predict_delay_parameters works with NULL newdata and the latent lognormal model", { # nolint: line_length_linter.
   skip_on_cran()
   set.seed(1)
-  prep_obs <- as_latent_individual(sim_obs)
-  fit <- epidist(
-    data = prep_obs,
-    seed = 1,
-    silent = 2,
-    output_dir = fs::dir_create(tempfile())
-  )
   pred <- predict_delay_parameters(fit)
   expect_s3_class(pred, "lognormal_samples")
   expect_s3_class(pred, "data.frame")
@@ -21,14 +14,6 @@ test_that("predict_delay_parameters works with NULL newdata and the latent logno
 test_that("predict_delay_parameters accepts newdata arguments and prediction by sex recovers underlying parameters", { # nolint: line_length_linter.
   skip_on_cran()
   set.seed(1)
-  prep_obs_sex <- as_latent_individual(sim_obs_sex)
-  fit_sex <- epidist(
-    data = prep_obs_sex,
-    formula = brms::bf(mu ~ 1 + sex, sigma ~ 1 + sex),
-    seed = 1,
-    silent = 2,
-    output_dir = fs::dir_create(tempfile())
-  )
   pred_sex <- predict_delay_parameters(fit_sex, prep_obs_sex)
   expect_s3_class(pred_sex, "lognormal_samples")
   expect_s3_class(pred_sex, "data.frame")
diff --git a/vignettes/approx-inference.Rmd b/vignettes/approx-inference.Rmd
index c06d6bc63..5892fce9e 100644
--- a/vignettes/approx-inference.Rmd
+++ b/vignettes/approx-inference.Rmd
@@ -99,9 +99,16 @@ library(purrr)
 library(tidyr)
 library(tibble)
 library(tidybayes)
+library(cmdstanr)
 ```
 
-First, we begin by simulating data.
+To access the approximate inference methods used in this vignette we will need to use the `cmdstanr` backend for `brms` (we generally recommend using this backend for fitting models). To do this, we first need to install CmdStan (see the README for more details). We can check we have everything we need as follows:
+
+```{r}
+cmdstanr::cmdstan_version()
+```
+
+We can simulate data to use for fitting models.
 The example data simulation process follows that used in the [Getting started with epidist](https://epidist.epinowcast.org/articles/epidist.html#data) vignette, so we will not detail exactly what is happening here, but please consult that vignette if interested:
 
 ```{r}
@@ -134,7 +141,7 @@ linelist <- as_epidist_linelist_time(obs_cens_trunc_samp)
 data <- as_latent_individual(linelist)
 
 t <- proc.time()
-fit_hmc <- epidist(data = data, algorithm = "sampling")
+fit_hmc <- epidist(data = data, algorithm = "sampling", backend = "cmdstanr")
 time_hmc <- proc.time() - t
 ```
 
@@ -145,11 +152,15 @@ To match the four Markov chains of length 1000 in HMC above, we then draw 4000 s
 
 ```{r results='hide'}
 t <- proc.time()
-fit_laplace <- epidist(data = data, algorithm = "laplace", draws = 4000)
+fit_laplace <- epidist(
+  data = data, algorithm = "laplace", draws = 4000, backend = "cmdstanr"
+)
 time_laplace <- proc.time() - t
 
 t <- proc.time()
-fit_advi <- epidist(data = data, algorithm = "meanfield", draws = 4000)
+fit_advi <- epidist(
+  data = data, algorithm = "meanfield", draws = 4000, backend = "cmdstanr"
+)
 time_advi <- proc.time() - t
 ```
 
@@ -160,7 +171,8 @@ Although a `save_single_paths` option is available, which may have allowed recov
 ```{r}
 t <- proc.time()
 fit_pathfinder <- epidist(
-  data = data, algorithm = "pathfinder", draws = 4000, num_paths = 1
+  data = data, algorithm = "pathfinder", draws = 4000, num_paths = 1,
+  backend = "cmdstanr"
 )
 time_pathfinder <- proc.time() - t
 ```
diff --git a/vignettes/ebola.Rmd b/vignettes/ebola.Rmd
index 7f506f294..0fd5ce75b 100644
--- a/vignettes/ebola.Rmd
+++ b/vignettes/ebola.Rmd
@@ -53,10 +53,21 @@ library(tidybayes)
 library(modelr)
 library(patchwork)
 library(lubridate)
+library(cmdstanr)
 ```
 
 For users new to `epidist`, before reading this article we recommend beginning with the "[Getting started with `epidist`](http://epidist.epinowcast.org/articles/epidist.html)" vignette.
 
+# Using the cmdstanr backend
+
+As the models explored in this vignette are relatively complex, we recommend using the `cmdstanr` backend for fitting models as it is typically more performant than the default `rstan` backend.
+To use the `cmdstanr` backend, we first need to install CmdStan (see the README for more details).
+We can check we have everything we need as follows:
+
+```{r}
+cmdstanr::cmdstan_version()
+```
+
 # Data preparation
 
 We begin by loading the Ebola line list data:
@@ -115,7 +126,9 @@ These shapefiles were obtained from the Database of Global Administrative Areas
 
 ```{r message=FALSE, warning=FALSE}
 sf <- st_read(
-  system.file("extdata/gadm41_SLE_shp/gadm41_SLE_2.shp", package = "epidist"),
+  system.file(
+    "extdata", "gadm41_SLE_shp", "gadm41_SLE_2.shp", package = "epidist"
+  ),
   quiet = TRUE
 )
 
@@ -231,8 +244,9 @@ fit <- epidist(
   family = lognormal(),
   algorithm = "sampling",
   refresh = 0,
-  silent = 2,
-  seed = 1
+  silent = 2, refresh = 0,
+  seed = 1,
+  backend = "cmdstanr"
 )
 ```
 
@@ -255,8 +269,9 @@ fit_sex <- epidist(
   family = lognormal(),
   algorithm = "sampling",
   refresh = 0,
-  silent = 2,
-  seed = 1
+  silent = 2, refresh = 0,
+  seed = 1,
+  backend = "cmdstanr"
 )
 ```
 
@@ -284,8 +299,9 @@ fit_sex_district <- epidist(
   family = lognormal(),
   algorithm = "sampling",
   refresh = 0,
-  silent = 2,
-  seed = 1
+  silent = 2, refresh = 0,
+  seed = 1,
+  backend = "cmdstanr"
 )
 ```
 
diff --git a/vignettes/epidist.Rmd b/vignettes/epidist.Rmd
index cbfd3972f..62082587a 100644
--- a/vignettes/epidist.Rmd
+++ b/vignettes/epidist.Rmd
@@ -28,9 +28,6 @@ knitr::opts_chunk$set(
   warning = FALSE,
   error = FALSE
 )
-
-# Required for R-CMD-as-cran-check to pass
-cmdstanr::set_cmdstan_path()
 ```
 
 Many quantities in epidemiology can be thought of as the time between two events or "delays".
@@ -57,7 +54,7 @@ To run this vignette yourself, along with the `epidist` package, you will need t
 
 ```{r load-requirements}
 library(epidist)
-library(purrr)
+library(brms)
 library(ggplot2)
 library(gt)
 library(dplyr)
@@ -157,7 +154,6 @@ This means that rather than exact event times, we observe event times within an
 Here we suppose that the interval is daily, meaning that only the date of the primary or secondary event, not the exact event time, is reported (Figure \@ref(fig:cens)):
 
 ```{r}
-# observe_process() should be renamed and include a "daily" i.e. 1 argument
 obs_cens <- obs |> observe_process()
 ```
 
@@ -179,7 +175,6 @@ This is called (right) truncation, and biases the observation process towards sh
 
 ```{r}
 obs_time <- 25
-# filter_obs_by_obs_time() should be renamed to refer to stime
 obs_cens_trunc <- filter_obs_by_obs_time(obs_cens, obs_time = obs_time)
 ```
 
@@ -269,20 +264,25 @@ The parameters of the model are inferred using Bayesian inference.
 In particular, we use the the No-U-Turn Sampler (NUTS) Markov chain Monte Carlo (MCMC) algorithm via the [`brms`](https://paul-buerkner.github.io/brms/) R package [@brms].
 
 ```{r}
-fit <- epidist(data = data, cores = 4, refresh = 0)
+fit <- epidist(data = data, chains = 2, cores = 2, refresh = 0)
 ```
 
 The `fit` object is a `brmsfit` object containing MCMC samples from each of the parameters (Table \@ref(tab:pars)) in the model.
 Users familiar with Stan and `brms`, can work with `fit` directly.
 Any tool that supports `brms` fitted model objects will be compatible with `fit`.
 
-(ref:pars) All of the parameters that are included in the model. Many of these parameters (e.g. `swindow` and `pwindow`) the so called latent variables in the model, and have lengths corresponding to the `sample_size`.
+**Note that here we use the default `rstan` backend but we generally recommend using the `cmdstanr` backend for faster sampling and additional features. This can be set using `backend = "cmdstanr"` after following the installing `cmdstan` instructions in the README.**
+
+(ref:pars) All of the parameters that are included in the model. Many of these parameters (e.g. `swindow` and `pwindow`) are the so called latent variables in the model, and have lengths corresponding to the `sample_size`. We extracted the model parameters using `brms::variables()` and removed the indices.
 
 ```{r pars}
-pars <- fit$fit@par_dims |>
-  map(.f = function(x) if (identical(x, integer(0))) return(1) else return(x))
+pars <- fit |>
+  variables(reserved = FALSE) |>
+  gsub(pattern = "\\[\\d+\\]", replacement = "")
 
-data.frame("Parameter" = names(pars), "Length" = unlist(pars)) |>
+data.frame(
+  "Parameter" = unique(pars), "Length" = table(pars)
+) |>
   gt() |>
   tab_caption("(ref:pars)")
 ```
diff --git a/vignettes/faq.Rmd b/vignettes/faq.Rmd
index d8fc2605b..f93a80a2a 100644
--- a/vignettes/faq.Rmd
+++ b/vignettes/faq.Rmd
@@ -223,4 +223,13 @@ For example, for the `latent_individual` model, currently methods are implemente
 If you are using another family, consider [submitting a pull request](https://github.com/epinowcast/epidist/pulls) to implement these methods!
 In doing so, you may find it useful to use the [`primarycensored`](https://primarycensored.epinowcast.org/) package.
 
+# How can I use the `cmdstanr` backend?
+
+The `cmdstanr` backend is typically more performant than the default `rstan` backend.
+To use the `cmdstanr` backend, we first need to install CmdStan (see the README for more details). We can check we have everything we need as follows:
+
+```{r}
+cmdstanr::cmdstan_version()
+```
+
 ## Bibliography