diff --git a/DESCRIPTION b/DESCRIPTION
index a569f8fc..9b8a93e3 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -55,6 +55,7 @@ Suggests:
     loo (>= 2.0.0),
     qs2,
     rmarkdown,
-    testthat (>= 2.1.0),
+    testthat (>= 3.3.0),
     Rcpp
 VignetteBuilder: knitr
+Config/testthat/edition: 3
diff --git a/tests/testthat/_snaps/cpp_opts.md b/tests/testthat/_snaps/cpp_opts.md
new file mode 100644
index 00000000..7d166fa1
--- /dev/null
+++ b/tests/testthat/_snaps/cpp_opts.md
@@ -0,0 +1,8 @@
+# validate_cpp_options works
+
+    STAN_OPENCL set to FALSE Since this is a non-empty value, it will result in the corresponding ccp option being turned ON. To turn this option off, use cpp_options = list(stan_opencl = NULL).
+
+# exe_info cpp_options comparison works
+
+    Recompiling is recommended due to missing exe_info.
+
diff --git a/tests/testthat/_snaps/data/process-data-float-rounding.json b/tests/testthat/_snaps/data/process-data-float-rounding.json
new file mode 100644
index 00000000..5861f67f
--- /dev/null
+++ b/tests/testthat/_snaps/data/process-data-float-rounding.json
@@ -0,0 +1,4 @@
+{
+  "a": 3,
+  "b": 2.0
+}
diff --git a/tests/testthat/_snaps/data/process-data-int-matrix.json b/tests/testthat/_snaps/data/process-data-int-matrix.json
new file mode 100644
index 00000000..ad68d520
--- /dev/null
+++ b/tests/testthat/_snaps/data/process-data-int-matrix.json
@@ -0,0 +1,7 @@
+{
+  "k": [
+    [18, 18, 16],
+    [13, 9, 6],
+    [4, 4, 4]
+  ]
+}
diff --git a/tests/testthat/_snaps/data/process-data-int-real.json b/tests/testthat/_snaps/data/process-data-int-real.json
new file mode 100644
index 00000000..0d644740
--- /dev/null
+++ b/tests/testthat/_snaps/data/process-data-int-real.json
@@ -0,0 +1,4 @@
+{
+  "a": 1,
+  "b": 2.0
+}
diff --git a/tests/testthat/_snaps/data/process-data-large-real.json b/tests/testthat/_snaps/data/process-data-large-real.json
new file mode 100644
index 00000000..8e82b0a6
--- /dev/null
+++ b/tests/testthat/_snaps/data/process-data-large-real.json
@@ -0,0 +1,4 @@
+{
+  "a": 1,
+  "b": 1774000000.0
+}
diff --git a/tests/testthat/_snaps/example.md b/tests/testthat/_snaps/example.md
new file mode 100644
index 00000000..55bce5f3
--- /dev/null
+++ b/tests/testthat/_snaps/example.md
@@ -0,0 +1,47 @@
+# print_example_program outputs stay stable
+
+    Code
+      cat(print_example_program("schools"))
+    Output
+      data {
+        int<lower=1> J;
+        vector<lower=0>[J] sigma;
+        vector[J] y;
+      }
+      parameters {
+        real mu;
+        real<lower=0> tau;
+        vector[J] theta;
+      }
+      model {
+        target += normal_lpdf(tau | 0, 10);
+        target += normal_lpdf(mu | 0, 10);
+        target += normal_lpdf(theta | mu, tau);
+        target += normal_lpdf(y | theta, sigma);
+      }
+
+---
+
+    Code
+      cat(print_example_program("schools_ncp"))
+    Output
+      data {
+        int<lower=1> J;
+        vector<lower=0>[J] sigma;
+        vector[J] y;
+      }
+      parameters {
+        real mu;
+        real<lower=0> tau;
+        vector[J] theta_raw;
+      }
+      transformed parameters {
+        vector[J] theta = mu + tau * theta_raw;
+      }
+      model {
+        target += normal_lpdf(tau | 0, 10);
+        target += normal_lpdf(mu | 0, 10);
+        target += normal_lpdf(theta_raw | 0, 1);
+        target += normal_lpdf(y | theta, sigma);
+      }
+
diff --git a/tests/testthat/_snaps/fit-gq.md b/tests/testthat/_snaps/fit-gq.md
new file mode 100644
index 00000000..28199208
--- /dev/null
+++ b/tests/testthat/_snaps/fit-gq.md
@@ -0,0 +1,79 @@
+# print() method works after gq
+
+    variable mean median sd mad q5 q95
+    y_rep[1]
+    y_rep[2]
+    y_rep[3]
+    y_rep[4]
+    y_rep[5]
+    y_rep[6]
+    y_rep[7]
+    y_rep[8]
+    y_rep[9]
+    y_rep[10]
+    
+    # showing 10 of 11 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)
+    variable mean median sd mad q5 q95
+    y_rep[1]
+    
+    # showing 1 of 11 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)
+    variable mad
+    y_rep[1]
+    y_rep[2]
+    y_rep[3]
+    y_rep[4]
+    y_rep[5]
+    y_rep[6]
+    y_rep[7]
+    y_rep[8]
+    y_rep[9]
+    y_rep[10]
+    
+    # showing 10 of 11 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)
+
+---
+
+    variable mean median sd mad q5 q95
+    y_rep[1]
+    y_rep[2]
+    
+    # showing 2 of 11 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)
+    variable mean median sd mad q5 q95
+    y_rep[1]
+    y_rep[2]
+    y_rep[3]
+    y_rep[4]
+    y_rep[5]
+    y_rep[6]
+    y_rep[7]
+    y_rep[8]
+    y_rep[9]
+    y_rep[10]
+    sum_y
+
+---
+
+    variable mean median sd mad q5 q95
+    y_rep[1]
+    y_rep[2]
+    
+    # showing 2 of 10 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)
+    variable mean median sd mad q5 q95
+    y_rep[1]
+    y_rep[2]
+    y_rep[3]
+    y_rep[4]
+    y_rep[5]
+    y_rep[6]
+    y_rep[7]
+    y_rep[8]
+    y_rep[9]
+    y_rep[10]
+
+---
+
+    variable mean median sd mad q5 q95
+    y_rep[1]
+    sum_y
+    y_rep[3]
+
diff --git a/tests/testthat/_snaps/fit-laplace.md b/tests/testthat/_snaps/fit-laplace.md
new file mode 100644
index 00000000..e3a75ab5
--- /dev/null
+++ b/tests/testthat/_snaps/fit-laplace.md
@@ -0,0 +1,17 @@
+# summary() and print() methods works after laplace
+
+        variable   mean median   sd  mad     q5    q95
+     lp__        -65.92 -65.65 1.36 1.20 -68.50 -64.28
+     lp_approx__  -2.04  -1.75 1.40 1.21  -4.68  -0.35
+     alpha         0.39   0.38 0.22 0.23   0.03   0.74
+     beta[1]      -0.65  -0.65 0.25 0.26  -1.05  -0.25
+     beta[2]      -0.28  -0.28 0.22 0.23  -0.63   0.09
+     beta[3]       0.65   0.66 0.26 0.26   0.21   1.07
+
+---
+
+     variable   mean median   sd  mad     q5    q95
+         lp__ -65.92 -65.65 1.36 1.20 -68.50 -64.28
+    
+     # showing 1 of 6 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)
+
diff --git a/tests/testthat/_snaps/fit-mcmc.md b/tests/testthat/_snaps/fit-mcmc.md
new file mode 100644
index 00000000..19b6d705
--- /dev/null
+++ b/tests/testthat/_snaps/fit-mcmc.md
@@ -0,0 +1,85 @@
+# print() method works after mcmc
+
+    variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
+    lp__
+    alpha
+    beta[1]
+    beta[2]
+    beta[3]
+    variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
+    lp__
+    
+    # showing 1 of 5 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)
+    variable ess_sd
+    lp__
+    alpha
+    beta[1]
+    beta[2]
+    beta[3]
+
+---
+
+    variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
+    lp__
+    mu
+    tau
+    theta_raw[1]
+    theta_raw[2]
+    theta_raw[3]
+    theta_raw[4]
+    theta_raw[5]
+    theta_raw[6]
+    theta_raw[7]
+    
+    # showing 10 of 19 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)
+    variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
+    lp__
+    mu
+    
+    # showing 2 of 19 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)
+    variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
+    lp__
+    mu
+    tau
+    theta_raw[1]
+    theta_raw[2]
+    theta_raw[3]
+    theta_raw[4]
+    theta_raw[5]
+    theta_raw[6]
+    theta_raw[7]
+    theta_raw[8]
+    theta[1]
+    theta[2]
+    theta[3]
+    theta[4]
+    theta[5]
+    theta[6]
+    theta[7]
+    theta[8]
+
+---
+
+    variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
+    theta[1]
+    theta[2]
+    
+    # showing 2 of 8 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)
+    variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
+    theta[1]
+    theta[2]
+    theta[3]
+    theta[4]
+    theta[5]
+    theta[6]
+    theta[7]
+    theta[8]
+
+---
+
+    variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
+    theta[1]
+    tau
+    mu
+    theta_raw[3]
+
diff --git a/tests/testthat/_snaps/fit-mle.md b/tests/testthat/_snaps/fit-mle.md
new file mode 100644
index 00000000..9eb3a1b3
--- /dev/null
+++ b/tests/testthat/_snaps/fit-mle.md
@@ -0,0 +1,20 @@
+# print() method works after optimization
+
+     variable estimate
+      lp__      -63.92
+      alpha       0.36
+      beta[1]    -0.63
+      beta[2]    -0.26
+      beta[3]     0.65
+
+---
+
+     variable estimate
+         lp__   -63.92
+    
+     # showing 1 of 5 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)
+
+---
+
+    Can't find the following variable(s): unknown
+
diff --git a/tests/testthat/_snaps/fit-vb.md b/tests/testthat/_snaps/fit-vb.md
new file mode 100644
index 00000000..f66d66ca
--- /dev/null
+++ b/tests/testthat/_snaps/fit-vb.md
@@ -0,0 +1,71 @@
+# print() method works after vb
+
+        variable   mean median   sd  mad     q5    q95
+     lp__        -66.75 -66.18 2.16 1.79 -70.59 -64.40
+     lp_approx__  -2.04  -1.74 1.39 1.26  -4.91  -0.40
+     alpha         0.30   0.29 0.21 0.20  -0.03   0.64
+     beta[1]      -0.59  -0.58 0.37 0.37  -1.18   0.02
+     beta[2]      -0.25  -0.25 0.22 0.23  -0.64   0.09
+     beta[3]       0.60   0.61 0.27 0.28   0.13   1.02
+     variable   mean median   sd  mad     q5    q95
+         lp__ -66.75 -66.18 2.16 1.79 -70.59 -64.40
+    
+     # showing 1 of 6 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)
+
+---
+
+         variable   mean median   sd  mad     q5    q95
+     lp__         -47.27 -46.75 2.69 2.15 -51.93 -43.98
+     lp_approx__   -5.05  -4.68 2.27 2.08  -9.18  -1.98
+     mu             6.55   6.52 3.82 3.78   0.24  12.67
+     tau            3.87   3.05 2.81 2.00   1.01   9.44
+     theta_raw[1]   0.44   0.44 0.98 0.97  -1.15   1.97
+     theta_raw[2]   0.14   0.12 0.92 0.89  -1.39   1.63
+     theta_raw[3]  -0.21  -0.18 1.00 1.01  -1.86   1.38
+     theta_raw[4]   0.11   0.12 0.88 0.88  -1.35   1.55
+     theta_raw[5]  -0.11  -0.12 1.00 1.03  -1.66   1.57
+     theta_raw[6]   0.04   0.03 0.78 0.81  -1.26   1.31
+    
+     # showing 10 of 20 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)
+         variable   mean median   sd  mad     q5    q95
+     lp__         -47.27 -46.75 2.69 2.15 -51.93 -43.98
+     lp_approx__   -5.05  -4.68 2.27 2.08  -9.18  -1.98
+     mu             6.55   6.52 3.82 3.78   0.24  12.67
+     tau            3.87   3.05 2.81 2.00   1.01   9.44
+     theta_raw[1]   0.44   0.44 0.98 0.97  -1.15   1.97
+     theta_raw[2]   0.14   0.12 0.92 0.89  -1.39   1.63
+     theta_raw[3]  -0.21  -0.18 1.00 1.01  -1.86   1.38
+     theta_raw[4]   0.11   0.12 0.88 0.88  -1.35   1.55
+     theta_raw[5]  -0.11  -0.12 1.00 1.03  -1.66   1.57
+     theta_raw[6]   0.04   0.03 0.78 0.81  -1.26   1.31
+     theta_raw[7]   0.34   0.33 0.98 0.99  -1.24   1.96
+     theta_raw[8]   0.10   0.10 0.97 0.97  -1.47   1.68
+     theta[1]       8.21   7.82 6.46 5.22  -1.00  17.39
+     theta[2]       7.23   7.00 6.04 4.97  -2.01  16.92
+     theta[3]       5.75   5.92 6.03 5.26  -3.41  14.73
+     theta[4]       6.90   6.98 5.90 5.03  -2.08  16.23
+     theta[5]       6.03   5.78 6.50 5.13  -3.18  16.17
+     theta[6]       6.74   6.84 5.36 4.86  -1.99  15.91
+     theta[7]       7.80   7.57 6.21 5.41  -1.44  18.38
+     theta[8]       7.07   6.85 5.94 5.30  -2.19  17.13
+
+---
+
+     variable   mean median   sd  mad     q5    q95
+     theta[1]   8.21   7.82 6.46 5.22  -1.00  17.39
+     theta[2]   7.23   7.00 6.04 4.97  -2.01  16.92
+     theta[3]   5.75   5.92 6.03 5.26  -3.41  14.73
+     theta[4]   6.90   6.98 5.90 5.03  -2.08  16.23
+     theta[5]   6.03   5.78 6.50 5.13  -3.18  16.17
+     theta[6]   6.74   6.84 5.36 4.86  -1.99  15.91
+     theta[7]   7.80   7.57 6.21 5.41  -1.44  18.38
+     theta[8]   7.07   6.85 5.94 5.30  -2.19  17.13
+     tau        3.87   3.05 2.81 2.00   1.01   9.44
+     lp__     -47.27 -46.75 2.69 2.15 -51.93 -43.98
+    
+     # showing 10 of 11 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)
+
+---
+
+    Can't find the following variable(s): unknown
+
diff --git a/tests/testthat/_snaps/install.md b/tests/testthat/_snaps/install.md
new file mode 100644
index 00000000..a0fcc57a
--- /dev/null
+++ b/tests/testthat/_snaps/install.md
@@ -0,0 +1,52 @@
+# install_cmdstan() errors if invalid version or URL
+
+    Download of CmdStan failed with error: cannot open URL 'https://github.com/stan-dev/cmdstan/releases/download/v2.35.5/cmdstan-2.35.5.tar.gz'
+    Please check if the supplied version number is valid.
+
+---
+
+    Download of CmdStan failed with error: cannot open URL 'https://github.com/stan-dev/cmdstan/releases/download/v2.35.5/cmdstan-2.35.5.tar.gz'
+    Please check if the supplied release URL is valid.
+
+---
+
+    https://github.com/stan-dev/cmdstan/releases/tag/v2.24.0 is not a .tar.gz archive!cmdstanr supports installing from .tar.gz archives only.
+
+# toolchain checks on Unix work
+
+    Code
+      check_unix_cpp_compiler()
+    Condition
+      Error:
+      ! C++ compiler missing.
+
+---
+
+    Code
+      check_unix_make()
+    Condition
+      Error:
+      ! make missing.
+
+# check_rtools4x_windows_toolchain reports missing Rtools and make
+
+    Code
+      check_rtools4x_windows_toolchain()
+    Condition
+      Error:
+      ! 
+      Rtools<version> was not found but is required to run CmdStan with R version <version>.
+      Please install or reinstall the appropriate Rtools version for this R installation,
+      restart R, and then run cmdstanr::check_cmdstan_toolchain().
+
+# check_rtools4x_windows_toolchain validates install path and empty candidates
+
+    
+    Rtools44 is installed in a path with spaces or brackets, which is not supported.
+    Please reinstall the appropriate Rtools version for this R installation to a valid path,
+    restart R, and then run cmdstanr::check_cmdstan_toolchain().
+
+# check_cmdstan_toolchain(fix = TRUE) is deprecated
+
+    The 'fix' argument is deprecated and will be removed in a future release.
+
diff --git a/tests/testthat/_snaps/json/json-always-decimal-false.json b/tests/testthat/_snaps/json/json-always-decimal-false.json
new file mode 100644
index 00000000..756b0338
--- /dev/null
+++ b/tests/testthat/_snaps/json/json-always-decimal-false.json
@@ -0,0 +1,4 @@
+{
+  "a": 1,
+  "b": 2
+}
diff --git a/tests/testthat/_snaps/json/json-always-decimal-true.json b/tests/testthat/_snaps/json/json-always-decimal-true.json
new file mode 100644
index 00000000..0d644740
--- /dev/null
+++ b/tests/testthat/_snaps/json/json-always-decimal-true.json
@@ -0,0 +1,4 @@
+{
+  "a": 1,
+  "b": 2.0
+}
diff --git a/tests/testthat/_snaps/json/json-boolean.json b/tests/testthat/_snaps/json/json-boolean.json
new file mode 100644
index 00000000..0d8a3bf3
--- /dev/null
+++ b/tests/testthat/_snaps/json/json-boolean.json
@@ -0,0 +1,3 @@
+{
+  "N": [1, 0, 1]
+}
diff --git a/tests/testthat/_snaps/json/json-df-matrix.json b/tests/testthat/_snaps/json/json-df-matrix.json
new file mode 100644
index 00000000..a5132dcf
--- /dev/null
+++ b/tests/testthat/_snaps/json/json-df-matrix.json
@@ -0,0 +1,7 @@
+{
+  "X": [
+    [1, 0.2],
+    [2, 0.3],
+    [3, 0.4]
+  ]
+}
diff --git a/tests/testthat/_snaps/json/json-factor.json b/tests/testthat/_snaps/json/json-factor.json
new file mode 100644
index 00000000..cf77867f
--- /dev/null
+++ b/tests/testthat/_snaps/json/json-factor.json
@@ -0,0 +1,3 @@
+{
+  "N": [1, 2, 3, 3, 2, 1]
+}
diff --git a/tests/testthat/_snaps/json/json-integer.json b/tests/testthat/_snaps/json/json-integer.json
new file mode 100644
index 00000000..872c8b0a
--- /dev/null
+++ b/tests/testthat/_snaps/json/json-integer.json
@@ -0,0 +1,3 @@
+{
+  "N": [1, 2, 3, 4]
+}
diff --git a/tests/testthat/_snaps/json/json-matrix-lists.json b/tests/testthat/_snaps/json/json-matrix-lists.json
new file mode 100644
index 00000000..22456268
--- /dev/null
+++ b/tests/testthat/_snaps/json/json-matrix-lists.json
@@ -0,0 +1,12 @@
+{
+  "M": [
+    [
+      [1, 3],
+      [2, 4]
+    ],
+    [
+      [5, 6],
+      [7, 8]
+    ]
+  ]
+}
diff --git a/tests/testthat/_snaps/json/json-table-array.json b/tests/testthat/_snaps/json/json-table-array.json
new file mode 100644
index 00000000..e8e864e1
--- /dev/null
+++ b/tests/testthat/_snaps/json/json-table-array.json
@@ -0,0 +1,28 @@
+{
+  "x": [
+    [
+      [5, 0, 0, 0],
+      [0, 0, 0, 0],
+      [0, 0, 0, 0],
+      [0, 0, 0, 0]
+    ],
+    [
+      [0, 0, 0, 0],
+      [0, 5, 0, 0],
+      [0, 0, 0, 0],
+      [0, 0, 0, 0]
+    ],
+    [
+      [0, 0, 0, 0],
+      [0, 0, 0, 0],
+      [0, 0, 5, 0],
+      [0, 0, 0, 0]
+    ],
+    [
+      [0, 0, 0, 0],
+      [0, 0, 0, 0],
+      [0, 0, 0, 0],
+      [0, 0, 0, 5]
+    ]
+  ]
+}
diff --git a/tests/testthat/_snaps/json/json-table-matrix.json b/tests/testthat/_snaps/json/json-table-matrix.json
new file mode 100644
index 00000000..b10267e3
--- /dev/null
+++ b/tests/testthat/_snaps/json/json-table-matrix.json
@@ -0,0 +1,8 @@
+{
+  "x": [
+    [5, 0, 0, 0],
+    [0, 5, 0, 0],
+    [0, 0, 5, 0],
+    [0, 0, 0, 5]
+  ]
+}
diff --git a/tests/testthat/_snaps/json/json-table-vector.json b/tests/testthat/_snaps/json/json-table-vector.json
new file mode 100644
index 00000000..63e4c2b4
--- /dev/null
+++ b/tests/testthat/_snaps/json/json-table-vector.json
@@ -0,0 +1,3 @@
+{
+  "x": [5, 5, 5, 5]
+}
diff --git a/tests/testthat/_snaps/json/json-unboxing.json b/tests/testthat/_snaps/json/json-unboxing.json
new file mode 100644
index 00000000..e7dc5060
--- /dev/null
+++ b/tests/testthat/_snaps/json/json-unboxing.json
@@ -0,0 +1,3 @@
+{
+  "N": 10
+}
diff --git a/tests/testthat/_snaps/json/json-vector-lists.json b/tests/testthat/_snaps/json/json-vector-lists.json
new file mode 100644
index 00000000..9f9bb3e3
--- /dev/null
+++ b/tests/testthat/_snaps/json/json-vector-lists.json
@@ -0,0 +1,6 @@
+{
+  "N": [
+    [1, 2, 3],
+    [4, 5, 6]
+  ]
+}
diff --git a/tests/testthat/_snaps/model-code-print.md b/tests/testthat/_snaps/model-code-print.md
new file mode 100644
index 00000000..0bc92b13
--- /dev/null
+++ b/tests/testthat/_snaps/model-code-print.md
@@ -0,0 +1,40 @@
+# code() and print() methods work
+
+    data {
+      int<lower=0> N;
+      array[N] int<lower=0, upper=1> y;
+    }
+    parameters {
+      real<lower=0, upper=1> theta;
+    }
+    model {
+      theta ~ beta(1, 1); // uniform prior on interval 0,1
+      y ~ bernoulli(theta);
+    }
+
+---
+
+    data {
+      int<lower=0> N;
+      array[N] int<lower=0, upper=1> y;
+    }
+    parameters {
+      real<lower=0, upper=1> theta;
+    }
+    model {
+      theta ~ beta(1, 1); // uniform prior on interval 0,1
+      y ~ bernoulli(theta);
+    }
+
+# code() warns and print() errors if only exe and no Stan file
+
+    '$code()' will return NULL because the 'CmdStanModel' was not created with a Stan file.
+
+---
+
+    '$print()' cannot be used because the 'CmdStanModel' was not created with a Stan file.
+
+# check_syntax() errors if only exe and no Stan file
+
+    '$check_syntax()' cannot be used because the 'CmdStanModel' was not created with a Stan file.
+
diff --git a/tests/testthat/_snaps/model-compile.md b/tests/testthat/_snaps/model-compile.md
new file mode 100644
index 00000000..7f71e315
--- /dev/null
+++ b/tests/testthat/_snaps/model-compile.md
@@ -0,0 +1,21 @@
+# name in STANCFLAGS is set correctly
+
+    Code
+      cat(trim_stanc_invocations(out), sep = "\n")
+    Output
+      bin/stanc --name='bernoulli_model' --o
+
+---
+
+    Code
+      cat(trim_stanc_invocations(out), sep = "\n")
+    Output
+      bin/stanc --name='bernoulli2_model' --o
+
+# STANCFLAGS from get_cmdstan_flags() are included in compile output
+
+    Code
+      cat(trim_stanc_invocations(out), sep = "\n")
+    Output
+      bin/stanc --name='bernoulli_model'  --O1  --warn-pedantic --o
+
diff --git a/tests/testthat/_snaps/path.md b/tests/testthat/_snaps/path.md
new file mode 100644
index 00000000..9b05ce9e
--- /dev/null
+++ b/tests/testthat/_snaps/path.md
@@ -0,0 +1,28 @@
+# Setting bad path leads to warning (can't find directory)
+
+    Path not set. Can't find directory: BAD_PATH
+
+# Setting bad path from env leads to warning (can't find directory)
+
+    Can't find directory specified by environment variable 'CMDSTAN'. Path not set.
+
+# Existing CMDSTAN env path with no install resets cached state
+
+    No CmdStan installation found in the path specified by the environment variable 'CMDSTAN'.
+
+# Getting missing path leads to error (path not set)
+
+    CmdStan path has not been set yet. See ?set_cmdstan_path.
+
+# cmdstan_version() behaves correctly when version is not set
+
+    CmdStan path has not been set yet. See ?set_cmdstan_path.
+
+# Warning message is thrown if can't detect version number
+
+    Can't find CmdStan makefile to detect version number. Path may not point to valid installation.
+
+# Setting path rejects unsupported CmdStan versions
+
+    CmdStan path not set. CmdStan v2.34.0 is no longer supported. cmdstanr now requires CmdStan v2.35.0 or newer.
+
diff --git a/tests/testthat/_snaps/utils.md b/tests/testthat/_snaps/utils.md
new file mode 100644
index 00000000..862b6a06
--- /dev/null
+++ b/tests/testthat/_snaps/utils.md
@@ -0,0 +1,16 @@
+# check_ebfmi and computing ebfmi works
+
+    E-BFMI not computed because it is undefined for posterior chains of length less than 3.
+
+---
+
+    E-BFMI not computed because it is undefined for posterior chains of length less than 3.
+
+---
+
+    E-BFMI not computed because the 'energy__' diagnostic could not be located.
+
+---
+
+    E-BFMI not computed because the 'energy__' diagnostic could not be located.
+
diff --git a/tests/testthat/helper-custom-expectations.R b/tests/testthat/helper-custom-expectations.R
index 47a6e38b..731bd13d 100644
--- a/tests/testthat/helper-custom-expectations.R
+++ b/tests/testthat/helper-custom-expectations.R
@@ -1,17 +1,27 @@
 #' @param ... arguments passed to mod$compile()
 expect_compilation <- function(mod, ...) {
+  mtime_check_enabled <- FALSE
   if(length(mod$exe_file()) > 0 && file.exists(mod$exe_file())) {
+    original_mtime <- file.mtime(mod$exe_file())
+    suppressWarnings(Sys.setFileTime(mod$exe_file(), original_mtime - 10))
     before_mtime <- file.mtime(mod$exe_file())
+    mtime_check_enabled <- before_mtime < original_mtime
   } else {
     before_mtime <- NULL
   }
-  expect_interactive_message(mod$compile(...), "Compiling Stan program...")
+  rlang::with_interactive(value = TRUE, {
+    expect_message(mod$compile(...), "Compiling Stan program...")
+  })
   if(length(mod$exe_file()) == 0 || !file.exists(mod$exe_file())) {
     fail(sprint("Model executable '%s' does not exist after compilation.", mod$exe_file()))
   }
-  if(!is.null(before_mtime)) {
+  if(!is.null(before_mtime) && mtime_check_enabled) {
     after_mtime <- file.mtime(mod$exe_file())
-    expect(before_mtime != after_mtime, sprintf("Exe file '%s' has NOT changed, despite expecting (re)compilation", mod$exe_file()))
+    expect_gt(
+      after_mtime,
+      before_mtime,
+      sprintf("Exe file '%s' has NOT changed, despite expecting (re)compilation", mod$exe_file())
+    )
   }
   invisible(mod)
 }
@@ -20,13 +30,19 @@ expect_compilation <- function(mod, ...) {
 #' @param constructor_call a call returning a CmdStanModel object that should have been compiled
 #' @return the newly created model
 expect_call_compilation <- function(constructor_call) {
-  before_time <- Sys.time()
-  mod <- expect_interactive_message(constructor_call, "Compiling Stan program...")
+  constructor_call <- substitute(constructor_call)
+  before_time <- Sys.time() - 10
+  rlang::with_interactive(value = TRUE, {
+    expect_message(
+      mod <- rlang::eval_bare(constructor_call, parent.frame()),
+      "Compiling Stan program..."
+    )
+  })
   if(length(mod$exe_file()) == 0 || !file.exists(mod$exe_file())) {
     fail(sprint("Model executable '%s' does not exist after compilation.", mod$exe_file()))
   }
   after_mtime <- file.mtime(mod$exe_file())
-  expect(before_time <= after_mtime, sprintf("Exe file '%s' has old timestamp, despite expecting (re)compilation", mod$exe_file()))
+  expect_true(before_time <= after_mtime, sprintf("Exe file '%s' has old timestamp, despite expecting (re)compilation", mod$exe_file()))
   invisible(mod)
 }
 
@@ -38,9 +54,11 @@ expect_no_recompilation <- function(mod, ...) {
   }
 
   before_mtime <- file.mtime(mod$exe_file())
-  expect_interactive_message(mod$compile(...), "Model executable is up to date!")
+  rlang::with_interactive(value = TRUE, {
+    expect_message(mod$compile(...), "Model executable is up to date!")
+  })
   after_mtime <- file.mtime(mod$exe_file())
-  expect(before_mtime == after_mtime, sprintf("Model executable '%s' has changed, despite expecting no recompilation", mod$exe_file()))
+  expect_true(before_mtime == after_mtime, sprintf("Model executable '%s' has changed, despite expecting no recompilation", mod$exe_file()))
   invisible(mod)
 }
 
@@ -108,3 +126,43 @@ expect_equal_ignore_order <- function(object, expected, ...) {
 }
 
 expect_not_true <- function(...) expect_false(isTRUE(...))
+
+transform_print_snapshot <- function(x) {
+  vapply(x, function(line) {
+    line <- trimws(line)
+    if (!nzchar(line)) {
+      return(line)
+    }
+    if (grepl("^variable\\b", line)) {
+      return(gsub("\\s+", " ", line))
+    }
+    if (grepl("^# showing", line) ||
+        grepl("^Can't find the following variable\\(s\\):", line)) {
+      return(line)
+    }
+    sub("\\s+.*$", "", line)
+  }, character(1))
+}
+
+transform_unix_toolchain_snapshot <- function(x) {
+  x <- gsub("A suitable C\\+\\+ compiler was not found\\..*", "C++ compiler missing.", x)
+  x <- gsub("A C\\+\\+ compiler was not found\\..*", "C++ compiler missing.", x)
+  gsub("The 'make' tool was not found\\..*", "make missing.", x)
+}
+
+transform_r_version_snapshot <- function(x) {
+  x <- gsub("Rtools[0-9]+", "Rtools<version>", x)
+  gsub("R version [0-9]+\\.[0-9]+\\.[0-9]+", "R version <version>", x)
+}
+
+capture_print_snapshot <- function(code) {
+  code <- substitute(code)
+  value <- NULL
+  lines <- utils::capture.output(
+    value <- rlang::eval_bare(code, parent.frame())
+  )
+  list(
+    value = value,
+    lines = transform_print_snapshot(lines)
+  )
+}
diff --git a/tests/testthat/helper-mock-cli.R b/tests/testthat/helper-mock-cli.R
index 60a9e52d..799e8d1d 100644
--- a/tests/testthat/helper-mock-cli.R
+++ b/tests/testthat/helper-mock-cli.R
@@ -1,8 +1,9 @@
 real_wcr <- wsl_compatible_run
 
 with_mocked_cli <- function(code, compile_ret, info_ret) {
-  with_mocked_bindings(
-    code,
+  code <- substitute(code)
+  caller <- parent.frame()
+  local_mocked_bindings(
     wsl_compatible_run = function(command, args, ...) {
       if (
         !is.null(command)
@@ -17,8 +18,11 @@ with_mocked_cli <- function(code, compile_ret, info_ret) {
       } else {
         real_wcr(command = command, args = args, ...)
       }
-    }
+    },
+    .package = "cmdstanr",
+    .env = caller
   )
+  rlang::eval_bare(code, env = caller)
 }
 
 ######## Mock Compile Expectations #######
diff --git a/tests/testthat/helper-models.R b/tests/testthat/helper-models.R
index f1e248b8..759d72ef 100644
--- a/tests/testthat/helper-models.R
+++ b/tests/testthat/helper-models.R
@@ -1,3 +1,12 @@
+local({
+  stan_files <- dir(test_path("resources", "stan"), pattern = "\\.stan$", full.names = TRUE)
+  exe_files <- cmdstanr:::cmdstan_ext(cmdstanr:::strip_ext(stan_files))
+  existing_exe_files <- exe_files[file.exists(exe_files)]
+  if (length(existing_exe_files) > 0) {
+    unlink(existing_exe_files, force = TRUE)
+  }
+})
+
 testing_data <- function(name) {
   if (file.exists(test_path("resources", "data", paste0(name, ".data.rds")))) {
     readRDS(test_path("resources", "data", paste0(name, ".data.rds")))
diff --git a/tests/testthat/test-cpp_opts.R b/tests/testthat/test-cpp_opts.R
index e02d1c9d..aea550c4 100644
--- a/tests/testthat/test-cpp_opts.R
+++ b/tests/testthat/test-cpp_opts.R
@@ -34,7 +34,7 @@ test_that("validate_cpp_options works", {
       abc = FALSE
     )
   )
-  expect_warning(validate_cpp_options(list(STAN_OPENCL = FALSE)))
+  expect_snapshot_warning(validate_cpp_options(list(STAN_OPENCL = FALSE)))
 })
 
 
@@ -64,8 +64,5 @@ test_that("exe_info cpp_options comparison works", {
   ))
 
   # no exe_info -> no recompile based on cpp info
-  expect_warning(
-    expect_true(exe_info_reflects_cpp_options(list(), list())),
-    "Recompiling is recommended"
-  )
+  expect_snapshot_warning(expect_true(exe_info_reflects_cpp_options(list(), list())))
 })
diff --git a/tests/testthat/test-csv.R b/tests/testthat/test-csv.R
index 89c0faf3..0922c0f0 100644
--- a/tests/testthat/test-csv.R
+++ b/tests/testthat/test-csv.R
@@ -1,5 +1,3 @@
-context("read_cmdstan_csv")
-
 set_cmdstan_path()
 fit_bernoulli_optimize <- testing_fit("bernoulli", method = "optimize", seed = 1234)
 fit_bernoulli_variational <- testing_fit("bernoulli", method = "variational", seed = 123)
@@ -519,9 +517,10 @@ test_that("returning time works for read_cmdstan_csv", {
 
 test_that("time from read_cmdstan_csv matches time from fit$time()", {
   fit <- fit_bernoulli_thin_1
-  expect_equivalent(
+  expect_equal(
     read_cmdstan_csv(fit$output_files())$time$chains,
-    fit$time()$chains
+    fit$time()$chains,
+    ignore_attr = TRUE
   )
 })
 
diff --git a/tests/testthat/test-data.R b/tests/testthat/test-data.R
index b024f5e6..7b8f506f 100644
--- a/tests/testthat/test-data.R
+++ b/tests/testthat/test-data.R
@@ -1,5 +1,3 @@
-context("data-utils")
-
 set_cmdstan_path()
 fit <- testing_fit("bernoulli", method = "sample", seed = 123)
 fit_vb <- testing_fit("bernoulli", method = "variational", seed = 123)
@@ -328,27 +326,9 @@ test_that("process_data() corrrectly casts integers and floating point numbers",
   ")
   mod <- cmdstan_model(stan_file, compile = FALSE)
   test_file <- process_data(list(a = 1, b = 2), model_variables = mod$variables())
-  expect_match(
-    "  \"a\": 1,",
-    readLines(test_file)[2],
-    fixed = TRUE
-  )
-  expect_match(
-    "  \"b\": 2.0",
-    readLines(test_file)[3],
-    fixed = TRUE
-  )
+  expect_snapshot_file(test_file, "process-data-int-real.json")
   test_file <- process_data(list(a = 1L, b = 1774000000), model_variables = mod$variables())
-  expect_match(
-    "  \"a\": 1,",
-    readLines(test_file)[2],
-    fixed = TRUE
-  )
-  expect_match(
-    "  \"b\": 1774000000.0",
-    readLines(test_file)[3],
-    fixed = TRUE
-  )
+  expect_snapshot_file(test_file, "process-data-large-real.json")
 
   stan_file <- write_stan_file("
   data {
@@ -357,16 +337,7 @@ test_that("process_data() corrrectly casts integers and floating point numbers",
   ")
   mod <- cmdstan_model(stan_file, compile = FALSE)
   test_file <- process_data(list(k = matrix(c(18, 18, 16, 13, 9, 6, 4, 4, 4), nrow=3, ncol=3, byrow=T)), model_variables = mod$variables())
-  expect_match(
-    "  \"k\": [",
-    readLines(test_file)[2],
-    fixed = TRUE
-  )
-  expect_match(
-    "    [18, 18, 16],",
-    readLines(test_file)[3],
-    fixed = TRUE
-  )
+  expect_snapshot_file(test_file, "process-data-int-matrix.json")
 })
 
 test_that("process_data warns on int coercion", {
@@ -412,9 +383,5 @@ test_that("Floating-point differences do not cause truncation towards 0", {
   a <- 10*(3-2.7)
   expect_false(is.integer(a))
   test_file <- process_data(list(a = a, b = 2.0), model_variables = mod$variables())
-  expect_match(
-    "  \"a\": 3,",
-    readLines(test_file)[2],
-    fixed = TRUE
-  )
+  expect_snapshot_file(test_file, "process-data-float-rounding.json")
 })
diff --git a/tests/testthat/test-example.R b/tests/testthat/test-example.R
index d8c91910..5af93eb4 100644
--- a/tests/testthat/test-example.R
+++ b/tests/testthat/test-example.R
@@ -1,4 +1,15 @@
-context("cmdstanr_example")
+stan_program <- "
+  data {
+    int<lower=0> N;
+    array[N] int<lower=0,upper=1> y;
+  }
+  parameters {
+    real<lower=0,upper=1> theta;
+  }
+  model {
+    y ~ bernoulli(theta);
+  }
+  "
 
 test_that("cmdstanr_example works", {
   fit_mcmc <- cmdstanr_example("logistic", chains = 2, force_recompile = TRUE)
@@ -13,25 +24,12 @@ test_that("cmdstanr_example works", {
 
   fit_vb <- cmdstanr_example("logistic", method = "variational")
   checkmate::expect_r6(fit_vb, "CmdStanVB")
-
-  expect_output(print_example_program("schools"), "vector[J] theta", fixed=TRUE)
-  expect_output(print_example_program("schools_ncp"), "vector[J] theta_raw", fixed=TRUE)
 })
 
-
-# used in multiple tests below
-stan_program <- "
-  data {
-    int<lower=0> N;
-    array[N] int<lower=0,upper=1> y;
-  }
-  parameters {
-    real<lower=0,upper=1> theta;
-  }
-  model {
-    y ~ bernoulli(theta);
-  }
-  "
+test_that("print_example_program outputs stay stable", {
+  expect_snapshot(cat(print_example_program("schools")))
+  expect_snapshot(cat(print_example_program("schools_ncp")))
+})
 
 test_that("write_stan_file writes Stan file correctly", {
   skip_if_not_installed("rlang")
@@ -47,17 +45,16 @@ test_that("write_stan_file writes Stan file correctly", {
 })
 
 test_that("write_stan_file writes to specified directory and filename", {
-  dir <- file.path(test_path(), "answers")
+  dir <- withr::local_tempdir()
+  explicit_dir <- withr::local_tempdir()
   expect_equal(dirname(f1 <- write_stan_file(stan_program, dir = dir, basename = "pasta")),
                absolute_path(dir))
   expect_equal(f2 <- write_stan_file(stan_program, dir = dir, basename = "fruit.stan"),
                absolute_path(file.path(dir, "fruit.stan")))
   expect_equal(f3 <- write_stan_file(stan_program, dir = dir, basename = "vegetable"),
                absolute_path(file.path(dir, "vegetable.stan"))) # should add .stan extension if missing
-  expect_equal(f4 <- write_stan_file(stan_program, dir = tempdir(), basename = "test"),
-               absolute_path(file.path(tempdir(), "test.stan")))
-
-  try(file.remove(f1, f2, f3, f4), silent = TRUE)
+  expect_equal(f4 <- write_stan_file(stan_program, dir = explicit_dir, basename = "test"),
+               absolute_path(file.path(explicit_dir, "test.stan")))
 })
 
 test_that("write_stan_file creates dir if necessary", {
@@ -69,7 +66,7 @@ test_that("write_stan_file creates dir if necessary", {
 
 test_that("write_stan_file by default creates the same file for the same Stan model", {
   skip_if_not_installed("rlang")
-  dir <- file.path(test_path(), "answers")
+  dir <- withr::local_tempdir()
 
   f1 <- write_stan_file(stan_program, dir = dir)
   mtime1 <- file.info(f1)$mtime
@@ -95,24 +92,16 @@ test_that("write_stan_file by default creates the same file for the same Stan mo
 
   mtime5 <- file.info(f5)$mtime
   expect_true(mtime1 < mtime5)
-
-
-  try(file.remove(f1, f2, f4), silent = TRUE)
 })
 
 test_that("cmdstanr_write_stan_file_dir option works", {
   base_dir <- tempdir()
-  test_dir <- file.path(base_dir, "option_test")
-  if (!dir.exists(test_dir)) {
-    dir.create(test_dir)
-  }
-  options("cmdstanr_write_stan_file_dir" = test_dir)
-  file <- write_stan_file(stan_program)
-  expect_equal(repair_path(dirname(file)), repair_path(test_dir))
-  options("cmdstanr_write_stan_file_dir" = NULL)
+  test_dir <- withr::local_tempdir(pattern = "option_test")
+  local({
+    withr::local_options(list("cmdstanr_write_stan_file_dir" = test_dir))
+    file <- write_stan_file(stan_program)
+    expect_equal(repair_path(dirname(file)), repair_path(test_dir))
+  })
   file <- write_stan_file(stan_program)
   expect_equal(repair_path(dirname(file)), repair_path(base_dir))
-  if (!dir.exists(test_dir)) {
-    file.remove(test_dir)
-  }
 })
diff --git a/tests/testthat/test-failed-chains.R b/tests/testthat/test-failed-chains.R
index 7a67cf2c..b76e8700 100644
--- a/tests/testthat/test-failed-chains.R
+++ b/tests/testthat/test-failed-chains.R
@@ -1,4 +1,3 @@
-context("failed chains")
 set_cmdstan_path()
 stan_program <- testing_stan_file("chain_fails")
 stan_program_init_warnings <- testing_stan_file("init_warnings")
diff --git a/tests/testthat/test-fit-gq.R b/tests/testthat/test-fit-gq.R
index 33cd62be..6a47c3c3 100644
--- a/tests/testthat/test-fit-gq.R
+++ b/tests/testthat/test-fit-gq.R
@@ -1,5 +1,3 @@
-context("fitted-gq")
-
 set_cmdstan_path()
 fit <- testing_fit("bernoulli", method = "sample", seed = 123)
 fit_gq <- testing_fit("bernoulli_ppc", method = "generate_quantities", seed = 123, fitted_params = fit)
@@ -66,34 +64,41 @@ test_that("summary() method works after gq", {
 })
 
 test_that("print() method works after gq", {
-  expect_output(expect_s3_class(fit_gq$print(), "CmdStanGQ"), "variable")
-  expect_output(fit_gq$print(max_rows = 1), "# showing 1 of 11 rows")
-  expect_output(fit_gq$print(NULL, c("mad")), "mad")
-
-  expect_output(fit_gq$print(), "showing 10 of 11 rows")
-  expect_output(fit_gq$print(max_rows = 2), "showing 2 of 11 rows")
-  expect_output(fit_gq$print(max_rows = 11), "sum_y", fixed=TRUE) # last parameter
-  expect_output(fit_gq$print("y_rep", max_rows = 2), "showing 2 of 10 rows")
+  printed_fit_gq <- capture_print_snapshot(fit_gq$print())
+  expect_s3_class(printed_fit_gq$value, "CmdStanGQ")
+  expect_snapshot_output(cat(
+    paste(printed_fit_gq$lines, collapse = "\n"),
+    "\n",
+    paste(capture_print_snapshot(fit_gq$print(max_rows = 1))$lines, collapse = "\n"),
+    "\n",
+    paste(capture_print_snapshot(fit_gq$print(NULL, c("mad")))$lines, collapse = "\n"),
+    "\n",
+    sep = ""
+  ))
+  expect_snapshot_output(cat(
+    paste(capture_print_snapshot(fit_gq$print(max_rows = 2))$lines, collapse = "\n"),
+    "\n",
+    paste(capture_print_snapshot(fit_gq$print(max_rows = 11))$lines, collapse = "\n"),
+    "\n",
+    sep = ""
+  ))
+  expect_snapshot_output(cat(
+    paste(capture_print_snapshot(fit_gq$print("y_rep", max_rows = 2))$lines, collapse = "\n"),
+    "\n",
+    paste(capture_print_snapshot(fit_gq$print("y_rep"))$lines, collapse = "\n"),
+    "\n",
+    sep = ""
+  ))
+  expect_snapshot_output(cat(
+    paste(capture_print_snapshot(fit_gq$print(c("y_rep[1]", "sum_y", "y_rep[3]")))$lines, collapse = "\n"),
+    "\n",
+    sep = ""
+  ))
   expect_error(
     fit_gq$print(variable = "unknown", max_rows = 20),
     "Can't find the following variable(s): unknown",
     fixed = TRUE
-  ) # unknown parameter
-
-  out <- capture.output(fit_gq$print("y_rep"))
-  expect_length(out, 11) # columns names + 1 y_rep
-  expect_match(out[1], "variable")
-  expect_match(out[2], "y_rep[1]", fixed = TRUE)
-  expect_match(out[9], "y_rep[8]", fixed = TRUE)
-  expect_false(any(grepl("sum_y|theta", out)))
-
-  # make sure the row order is correct
-  out <- capture.output(fit_gq$print(c("y_rep[1]", "sum_y", "y_rep[3]")))
-  expect_length(out, 4)
-  expect_match(out[1], " variable")
-  expect_match(out[2], " y_rep[1]", fixed = TRUE)
-  expect_match(out[3], " sum_y")
-  expect_match(out[4], " y_rep[3]", fixed = TRUE)
+  )
 })
 
 test_that("output() method works after gq", {
diff --git a/tests/testthat/test-fit-init.R b/tests/testthat/test-fit-init.R
index 7b78c73a..e4aca759 100644
--- a/tests/testthat/test-fit-init.R
+++ b/tests/testthat/test-fit-init.R
@@ -1,4 +1,3 @@
-context("fitted-inits")
 set_cmdstan_path()
 
 data_list_schools <- testing_data("schools")
diff --git a/tests/testthat/test-fit-laplace.R b/tests/testthat/test-fit-laplace.R
index 8991a9e1..1acb8f8b 100644
--- a/tests/testthat/test-fit-laplace.R
+++ b/tests/testthat/test-fit-laplace.R
@@ -1,5 +1,3 @@
-context("fitted-laplace")
-
 set_cmdstan_path()
 fit_laplace <- testing_fit("logistic", method = "laplace", seed = 100)
 PARAM_NAMES <- c("alpha", "beta[1]", "beta[2]", "beta[3]")
@@ -15,8 +13,9 @@ test_that("summary() and print() methods works after laplace", {
   expect_equal(x$variable, c("lp__", "lp_approx__", PARAM_NAMES))
   expect_equal(colnames(x), c("variable", "mean", "sd"))
 
-  expect_output(expect_s3_class(fit_laplace$print(), "CmdStanLaplace"), "variable")
-  expect_output(fit_laplace$print(max_rows = 1), "# showing 1 of 6 rows")
+  expect_s3_class(fit_laplace$print(), "CmdStanLaplace")
+  expect_snapshot_output(fit_laplace$print())
+  expect_snapshot_output(fit_laplace$print(max_rows = 1))
 })
 
 test_that("draws() method returns posterior sample (reading csv works)", {
diff --git a/tests/testthat/test-fit-mcmc.R b/tests/testthat/test-fit-mcmc.R
index 6876fd5d..45fb7ee2 100644
--- a/tests/testthat/test-fit-mcmc.R
+++ b/tests/testthat/test-fit-mcmc.R
@@ -1,5 +1,3 @@
-context("fitted-mcmc")
-
 set_cmdstan_path()
 fit_mcmc <- testing_fit("logistic", method = "sample",
                         seed = 123, chains = 2)
@@ -125,37 +123,45 @@ test_that("summary() method works after mcmc", {
 })
 
 test_that("print() method works after mcmc", {
-  expect_output(expect_s3_class(fit_mcmc$print(), "CmdStanMCMC"), "variable")
-  expect_output(fit_mcmc$print(max_rows = 1), "# showing 1 of 5 rows")
-  expect_output(fit_mcmc$print(NULL, c("ess_sd")), "ess_sd")
-
+  printed_fit_mcmc <- capture_print_snapshot(fit_mcmc$print())
+  expect_s3_class(printed_fit_mcmc$value, "CmdStanMCMC")
+  expect_snapshot_output(cat(
+    paste(printed_fit_mcmc$lines, collapse = "\n"),
+    "\n",
+    paste(capture_print_snapshot(fit_mcmc$print(max_rows = 1))$lines, collapse = "\n"),
+    "\n",
+    paste(capture_print_snapshot(fit_mcmc$print(NULL, c("ess_sd")))$lines, collapse = "\n"),
+    "\n",
+    sep = ""
+  ))
   # test on model with more parameters
   fit <- cmdstanr_example("schools_ncp")
-  expect_output(fit$print(), "showing 10 of 19 rows")
-  expect_output(fit$print(max_rows = 2), "showing 2 of 19 rows")
-  expect_output(fit$print(max_rows = 19), "theta[8]", fixed=TRUE) # last parameter
-  expect_output(fit$print("theta", max_rows = 2), "showing 2 of 8 rows")
+  expect_snapshot_output(cat(
+    paste(capture_print_snapshot(fit$print())$lines, collapse = "\n"),
+    "\n",
+    paste(capture_print_snapshot(fit$print(max_rows = 2))$lines, collapse = "\n"),
+    "\n",
+    paste(capture_print_snapshot(fit$print(max_rows = 19))$lines, collapse = "\n"),
+    "\n",
+    sep = ""
+  ))
+  expect_snapshot_output(cat(
+    paste(capture_print_snapshot(fit$print("theta", max_rows = 2))$lines, collapse = "\n"),
+    "\n",
+    paste(capture_print_snapshot(fit$print("theta"))$lines, collapse = "\n"),
+    "\n",
+    sep = ""
+  ))
+  expect_snapshot_output(cat(
+    paste(capture_print_snapshot(fit$print(c("theta[1]", "tau", "mu", "theta_raw[3]")))$lines, collapse = "\n"),
+    "\n",
+    sep = ""
+  ))
   expect_error(
     fit$print(variable = "unknown", max_rows = 20),
     "Can't find the following variable(s): unknown",
     fixed = TRUE
-  ) # unknown parameter
-
-  out <- capture.output(fit$print("theta"))
-  expect_length(out, 9) # columns names + 8 thetas
-  expect_match(out[1], "variable")
-  expect_match(out[2], "theta[1]", fixed = TRUE)
-  expect_match(out[9], "theta[8]", fixed = TRUE)
-  expect_false(any(grepl("mu|tau|theta_raw", out)))
-
-  # make sure the row order is correct
-  out <- capture.output(fit$print(c("theta[1]", "tau", "mu", "theta_raw[3]")))
-  expect_length(out, 5)
-  expect_match(out[1], " variable")
-  expect_match(out[2], " theta[1]", fixed = TRUE)
-  expect_match(out[3], " tau")
-  expect_match(out[4], " mu")
-  expect_match(out[5], " theta_raw[3]", fixed = TRUE)
+  )
 })
 
 test_that("output() method works after mcmc", {
@@ -331,23 +337,20 @@ test_that("loo works for all draws storage formats", {
   skip_if_not_installed("loo")
   fit <- testing_fit("bernoulli_log_lik")
 
-  options(cmdstanr_draws_format = "draws_array")
+  withr::local_options(list(cmdstanr_draws_format = "draws_array"))
   expect_s3_class(suppressWarnings(fit$loo()), "loo")
 
-  options(cmdstanr_draws_format = "draws_df")
+  withr::local_options(list(cmdstanr_draws_format = "draws_df"))
   expect_s3_class(suppressWarnings(fit$loo()), "loo")
 
-  options(cmdstanr_draws_format = "draws_matrix")
+  withr::local_options(list(cmdstanr_draws_format = "draws_matrix"))
   expect_s3_class(suppressWarnings(fit$loo()), "loo")
 
-  options(cmdstanr_draws_format = "draws_list")
+  withr::local_options(list(cmdstanr_draws_format = "draws_list"))
   expect_s3_class(suppressWarnings(fit$loo()), "loo")
 
-  options(cmdstanr_draws_format = "draws_rvars")
+  withr::local_options(list(cmdstanr_draws_format = "draws_rvars"))
   expect_s3_class(suppressWarnings(fit$loo()), "loo")
-
-  # reset option
-  options(cmdstanr_draws_format = NULL)
 })
 
 test_that("draws() works for different formats", {
diff --git a/tests/testthat/test-fit-mle.R b/tests/testthat/test-fit-mle.R
index cd87a214..3f1d5de2 100644
--- a/tests/testthat/test-fit-mle.R
+++ b/tests/testthat/test-fit-mle.R
@@ -1,5 +1,3 @@
-context("fitted-mle")
-
 set_cmdstan_path()
 fit_mle <- testing_fit("logistic", method = "optimize", seed = 123)
 mod <- testing_model("bernoulli")
@@ -24,13 +22,10 @@ test_that("summary method works after optimization", {
 })
 
 test_that("print() method works after optimization", {
-  expect_output(expect_s3_class(fit_mle$print(), "CmdStanMLE"), "estimate")
-  expect_output(fit_mle$print(max_rows = 1), "# showing 1 of 5 rows")
-  expect_error(
-    fit_mle$print(variable = "unknown", max_rows = 20),
-    "Can't find the following variable(s): unknown",
-    fixed = TRUE
-  ) # unknown parameter
+  expect_s3_class(fit_mle$print(), "CmdStanMLE")
+  expect_snapshot_output(fit_mle$print())
+  expect_snapshot_output(fit_mle$print(max_rows = 1))
+  expect_snapshot_error(fit_mle$print(variable = "unknown", max_rows = 20))
 })
 
 test_that("time() method works after optimization", {
diff --git a/tests/testthat/test-fit-shared.R b/tests/testthat/test-fit-shared.R
index d422c47c..8f83be31 100644
--- a/tests/testthat/test-fit-shared.R
+++ b/tests/testthat/test-fit-shared.R
@@ -1,5 +1,3 @@
-context("fitted-shared-methods")
-
 set_cmdstan_path()
 fits <- list()
 fits[["sample"]] <- testing_fit("logistic", method = "sample",
diff --git a/tests/testthat/test-fit-vb.R b/tests/testthat/test-fit-vb.R
index e701e951..1763a629 100644
--- a/tests/testthat/test-fit-vb.R
+++ b/tests/testthat/test-fit-vb.R
@@ -1,5 +1,3 @@
-context("fitted-vb")
-
 set_cmdstan_path()
 fit_vb <- testing_fit("logistic", method = "variational", seed = 123)
 fit_vb_sci_not <- testing_fit("logistic", method = "variational", seed = 123, iter = 200000, adapt_iter = 100000)
@@ -20,29 +18,22 @@ test_that("summary() method works after vb", {
 })
 
 test_that("print() method works after vb", {
-  expect_output(expect_s3_class(fit_vb$print(), "CmdStanVB"), "variable")
-  expect_output(fit_vb$print(max_rows = 1), "# showing 1 of 6 rows")
+  expect_s3_class(fit_vb$print(), "CmdStanVB")
 
   # test on model with more parameters
   fit <- cmdstanr_example("schools_ncp", method = "variational", seed = 123)
-  expect_output(fit$print(), "lp_approx__")
-  expect_output(fit$print(), "showing 10 of 20 rows")
-  expect_output(fit$print(max_rows = 20), "theta[8]", fixed = TRUE) # last parameter
-  expect_error(
-    fit$print(variable = "unknown", max_rows = 20),
-    "Can't find the following variable(s): unknown",
-    fixed = TRUE
-  ) # unknown parameter
-
-  out <- capture.output(fit$print(c("theta", "tau", "lp__", "lp_approx__")))
-  expect_length(out, 13) # columns names + 8 thetas + tau + lp__ + lp_approx__ + empty + message
-  expect_match(out[1], " variable")
-  expect_match(out[2], " theta[1]", fixed = TRUE)
-  expect_match(out[9], " theta[8]", fixed = TRUE)
-  expect_match(out[10], " tau")
-  expect_match(out[11], " lp__")
-  expect_false(nzchar(out[12])) # empty line
-  expect_match(out[13], "10 of 11 rows")
+  expect_snapshot_output({
+    fit_vb$print()
+    fit_vb$print(max_rows = 1)
+  })
+  expect_snapshot_output({
+    fit$print()
+    fit$print(max_rows = 20)
+  })
+  expect_snapshot_output({
+    fit$print(c("theta", "tau", "lp__", "lp_approx__"))
+  })
+  expect_snapshot_error(fit$print(variable = "unknown", max_rows = 20))
 })
 
 test_that("draws() method returns posterior sample (reading csv works)", {
diff --git a/tests/testthat/test-install.R b/tests/testthat/test-install.R
index 0f2b9965..86011d5e 100644
--- a/tests/testthat/test-install.R
+++ b/tests/testthat/test-install.R
@@ -1,5 +1,3 @@
-context("install")
-
 # avoid parallel on Mac due to strange intermittent TBB errors on Github Actions
 CORES <- if (os_is_macos()) 1 else 2
 
@@ -69,18 +67,13 @@ test_that("install_cmdstan() errors if it times out", {
 })
 
 test_that("install_cmdstan() errors if invalid version or URL", {
-  expect_error(
-    install_cmdstan(version = "2.35.5", wsl = os_is_wsl()),
-    "Download of CmdStan failed with error: cannot open URL 'https://github.com/stan-dev/cmdstan/releases/download/v2.35.5/cmdstan-2.35.5.tar.gz'\nPlease check if the supplied version number is valid."
-  )
-  expect_error(
+  expect_snapshot_error(install_cmdstan(version = "2.35.5", wsl = os_is_wsl()))
+  expect_snapshot_error(
     install_cmdstan(release_url = "https://github.com/stan-dev/cmdstan/releases/download/v2.35.5/cmdstan-2.35.5.tar.gz",
-                    wsl = os_is_wsl()),
-    "Download of CmdStan failed with error: cannot open URL 'https://github.com/stan-dev/cmdstan/releases/download/v2.35.5/cmdstan-2.35.5.tar.gz'\nPlease check if the supplied release URL is valid."
+                    wsl = os_is_wsl())
   )
-  expect_error(
-    install_cmdstan(release_url = "https://github.com/stan-dev/cmdstan/releases/tag/v2.24.0", wsl = os_is_wsl()),
-    "cmdstanr supports installing from .tar.gz archives only"
+  expect_snapshot_error(
+    install_cmdstan(release_url = "https://github.com/stan-dev/cmdstan/releases/tag/v2.24.0", wsl = os_is_wsl())
   )
 })
 
@@ -124,26 +117,17 @@ test_that("install_cmdstan() works with version and release_url", {
 
 test_that("toolchain checks on Unix work", {
   skip_if(os_is_windows())
-  path_backup <- Sys.getenv("PATH")
-  Sys.setenv("PATH" = "")
-  if (os_is_macos()) {
-    err_msg_cpp <- "A suitable C++ compiler was not found. Please install the command line tools for Mac with 'xcode-select --install' or install Xcode from the app store. Then restart R and run cmdstanr::check_cmdstan_toolchain()."
-    err_msg_make <- "The 'make' tool was not found. Please install the command line tools for Mac with 'xcode-select --install' or install Xcode from the app store. Then restart R and run cmdstanr::check_cmdstan_toolchain()."
-  } else {
-    err_msg_cpp <- "A C++ compiler was not found. Please install the 'clang++' or 'g++' compiler, restart R, and run cmdstanr::check_cmdstan_toolchain()."
-    err_msg_make <- "The 'make' tool was not found. Please install 'make', restart R, and then run cmdstanr::check_cmdstan_toolchain()."
-  }
-  expect_error(
+  withr::local_envvar(c("PATH" = ""))
+  expect_snapshot(
     check_unix_cpp_compiler(),
-    err_msg_cpp,
-    fixed = TRUE
+    error = TRUE,
+    transform = transform_unix_toolchain_snapshot
   )
-  expect_error(
+  expect_snapshot(
     check_unix_make(),
-    err_msg_make,
-    fixed = TRUE
+    error = TRUE,
+    transform = transform_unix_toolchain_snapshot
   )
-  Sys.setenv("PATH" = path_backup)
 })
 
 test_that("clean and rebuild works", {
@@ -369,24 +353,20 @@ test_that("rtools4x_toolchain_path prefers ABI-compatible legacy fallback", {
   file.create(file.path(fake_rtools_home, "ucrt64", "bin", "g++.exe"))
 
   withr::with_envvar(setNames(fake_rtools_home, env_var), {
-    with_mocked_bindings(
-      {
-        expect_equal(
-          rtools4x_toolchain_path(),
-          repair_path(file.path(fake_rtools_home, "mingw64", "bin"))
-        )
-      },
-      is_ucrt_toolchain = function() FALSE
-    )
-    with_mocked_bindings(
-      {
-        expect_equal(
-          rtools4x_toolchain_path(),
-          repair_path(file.path(fake_rtools_home, "ucrt64", "bin"))
-        )
-      },
-      is_ucrt_toolchain = function() TRUE
-    )
+    local({
+      local_mocked_bindings(is_ucrt_toolchain = function() FALSE)
+      expect_equal(
+        rtools4x_toolchain_path(),
+        repair_path(file.path(fake_rtools_home, "mingw64", "bin"))
+      )
+    })
+    local({
+      local_mocked_bindings(is_ucrt_toolchain = function() TRUE)
+      expect_equal(
+        rtools4x_toolchain_path(),
+        repair_path(file.path(fake_rtools_home, "ucrt64", "bin"))
+      )
+    })
   })
 })
 
@@ -423,64 +403,70 @@ test_that("check_rtools4x_windows_toolchain reports checked toolchain paths", {
 })
 
 test_that("toolchain_PATH_env_var() handles missing and configured Rtools homes", {
-  with_mocked_bindings(
-    expect_null(toolchain_PATH_env_var()),
-    os_is_windows = function() FALSE
-  )
-  with_mocked_bindings(
-    expect_null(toolchain_PATH_env_var()),
-    os_is_windows = function() TRUE,
-    rtools4x_home_path = function() ""
-  )
-  with_mocked_bindings(
+  local({
+    local_mocked_bindings(os_is_windows = function() FALSE)
+    expect_null(toolchain_PATH_env_var())
+  })
+  local({
+    local_mocked_bindings(
+      os_is_windows = function() TRUE,
+      rtools4x_home_path = function() ""
+    )
+    expect_null(toolchain_PATH_env_var())
+  })
+  local({
+    local_mocked_bindings(
+      os_is_windows = function() TRUE,
+      rtools4x_home_path = function() "C:/rtools",
+      rtools4x_toolchain_path = function() "C:/rtools/ucrt64/bin",
+      repair_path = function(path) path
+    )
     expect_equal(
       toolchain_PATH_env_var(),
       "C:/rtools/usr/bin;C:/rtools/ucrt64/bin"
-    ),
-    os_is_windows = function() TRUE,
-    rtools4x_home_path = function() "C:/rtools",
-    rtools4x_toolchain_path = function() "C:/rtools/ucrt64/bin",
-    repair_path = function(path) path
-  )
+    )
+  })
 })
 
 test_that("check_rtools4x_windows_toolchain reports missing Rtools and make", {
   fake_rtools_home <- tempfile(pattern = "rtools-home-missing-", tmpdir = tempdir(check = TRUE))
   on.exit(unlink(fake_rtools_home, recursive = TRUE), add = TRUE)
 
-  with_mocked_bindings(
-    expect_error(
+  local({
+    local_mocked_bindings(
+      rtools4x_home_path = function() "",
+      rtools4x_version = function() "44"
+    )
+    expect_snapshot(
       check_rtools4x_windows_toolchain(),
-      "restart R, and then run cmdstanr::check_cmdstan_toolchain()",
-      fixed = TRUE
-    ),
-    rtools4x_home_path = function() "",
-    rtools4x_version = function() "44"
-  )
+      error = TRUE,
+      transform = transform_r_version_snapshot
+    )
+  })
 
   dir.create(file.path(fake_rtools_home, "usr", "bin"),
              recursive = TRUE, showWarnings = FALSE)
-  with_mocked_bindings(
+  local({
+    local_mocked_bindings(
+      rtools4x_home_path = function() fake_rtools_home,
+      rtools4x_version = function() "44"
+    )
     expect_error(
       check_rtools4x_windows_toolchain(),
       "restart R, and then run cmdstanr::check_cmdstan_toolchain()",
       fixed = TRUE
-    ),
-    rtools4x_home_path = function() fake_rtools_home,
-    rtools4x_version = function() "44"
-  )
+    )
+  })
 })
 
 test_that("check_rtools4x_windows_toolchain validates install path and empty candidates", {
-  with_mocked_bindings(
-    expect_error(
-      check_rtools4x_windows_toolchain(),
-      "Please reinstall the appropriate Rtools version for this R installation to a valid path",
-      fixed = TRUE
-    ),
-    rtools4x_home_path = function() "C:/Program Files/Rtools44",
-    rtools4x_version = function() "44"
-  )
+  local({
+    local_mocked_bindings(
+      rtools4x_home_path = function() "C:/Program Files/Rtools44",
+      rtools4x_version = function() "44"
+    )
+    expect_snapshot_error(check_rtools4x_windows_toolchain())
+  })
 
   fake_rtools_home <- tempfile(pattern = "rtools-home-empty-", tmpdir = tempdir(check = TRUE))
   on.exit(unlink(fake_rtools_home, recursive = TRUE), add = TRUE)
@@ -488,22 +474,20 @@ test_that("check_rtools4x_windows_toolchain validates install path and empty can
              recursive = TRUE, showWarnings = FALSE)
   file.create(file.path(fake_rtools_home, "usr", "bin", "make.exe"))
 
-  with_mocked_bindings(
+  local({
+    local_mocked_bindings(
+      rtools4x_home_path = function() fake_rtools_home,
+      rtools4x_version = function() "44",
+      rtools4x_toolchain_candidates = function() character()
+    )
     expect_error(
       check_rtools4x_windows_toolchain(),
       "restart R, and then run cmdstanr::check_cmdstan_toolchain()",
       fixed = TRUE
-    ),
-    rtools4x_home_path = function() fake_rtools_home,
-    rtools4x_version = function() "44",
-    rtools4x_toolchain_candidates = function() character()
-  )
+    )
+  })
 })
 
 test_that("check_cmdstan_toolchain(fix = TRUE) is deprecated", {
-  expect_warning(
-    check_cmdstan_toolchain(fix = TRUE),
-    "The 'fix' argument is deprecated and will be removed in a future release",
-    fixed = TRUE
-  )
+  expect_snapshot_warning(check_cmdstan_toolchain(fix = TRUE))
 })
diff --git a/tests/testthat/test-json.R b/tests/testthat/test-json.R
index c2c2657c..f7bd9d1b 100644
--- a/tests/testthat/test-json.R
+++ b/tests/testthat/test-json.R
@@ -1,30 +1,22 @@
-context("json")
-
 test_that("JSON output unboxing works", {
   temp_file <- tempfile()
   N <- 10
   write_stan_json(list(N = N), file = temp_file)
-  json_output <- readLines(temp_file)
-  expect_known_output(cat(json_output, sep = "\n"),
-                      file = test_path("answers", "json-unboxing.json"))
+  expect_snapshot_file(temp_file, "json-unboxing.json")
 })
 
 test_that("JSON output for boolean is correct", {
   temp_file <- tempfile()
   N <- c(TRUE, FALSE, TRUE)
   write_stan_json(list(N = N), file = temp_file)
-  json_output <- readLines(temp_file)
-  expect_known_output(cat(json_output, sep = "\n"),
-                      file = test_path("answers", "json-boolean.json"))
+  expect_snapshot_file(temp_file, "json-boolean.json")
 })
 
 test_that("JSON output for factors is correct", {
   temp_file <- tempfile()
   N <- factor(c(0,1,2,2,1,0), labels = c("c1", "c2", "c3"))
   write_stan_json(list(N = N), file = temp_file)
-  json_output <- readLines(temp_file)
-  expect_known_output(cat(json_output, sep = "\n"),
-                      file = test_path("answers", "json-factor.json"))
+  expect_snapshot_file(temp_file, "json-factor.json")
 })
 
 test_that("JSON output for integer vector is correct", {
@@ -32,9 +24,7 @@ test_that("JSON output for integer vector is correct", {
   N <- c(1.0, 2.0, 3, 4)
 
   write_stan_json(list(N = N), file = temp_file)
-  json_output <- readLines(temp_file)
-  expect_known_output(cat(json_output, sep = "\n"),
-                      file = test_path("answers", "json-integer.json"))
+  expect_snapshot_file(temp_file, "json-integer.json")
 })
 
 test_that("JSON output for data frame and matrix is correct", {
@@ -47,15 +37,12 @@ test_that("JSON output for data frame and matrix is correct", {
 
   write_stan_json(list(X = df), file = temp_file_df)
   write_stan_json(list(X = mat), file = temp_file_mat)
-  json_output_mat <- readLines(temp_file_df)
-  json_output_df <- readLines(temp_file_mat)
-  expect_identical(json_output_df, json_output_mat)
+  expect_identical(readLines(temp_file_df), readLines(temp_file_mat))
 
   # Floating-point error introduced in jsonlite 1.8.5
   # https://github.com/jeroen/jsonlite/issues/420
   if (packageVersion("jsonlite") != "1.8.5") {
-    expect_known_output(cat(json_output_df, sep = "\n"),
-                        file = test_path("answers", "json-df-matrix.json"))
+    expect_snapshot_file(temp_file_df, "json-df-matrix.json")
   }
 })
 
@@ -64,9 +51,7 @@ test_that("JSON output for list of vectors is correct", {
   N <- list(c(1,2,3), c(4,5,6))
 
   write_stan_json(list(N = N), file = temp_file)
-  json_output <- readLines(temp_file)
-  expect_known_output(cat(json_output, sep = "\n"),
-                      file = test_path("answers", "json-vector-lists.json"))
+  expect_snapshot_file(temp_file, "json-vector-lists.json")
 })
 
 test_that("JSON output for list of matrices is correct", {
@@ -76,9 +61,7 @@ test_that("JSON output for list of matrices is correct", {
     matrix(5:8, nrow = 2, byrow = TRUE)
   )
   write_stan_json(list(M = matrices), file = temp_file)
-  json_output <- readLines(temp_file)
-  expect_known_output(cat(json_output, sep = "\n"),
-                      file = test_path("answers", "json-matrix-lists.json"))
+  expect_snapshot_file(temp_file, "json-matrix-lists.json")
 })
 
 test_that("JSON output for table is correct", {
@@ -86,19 +69,13 @@ test_that("JSON output for table is correct", {
   f <- factor(rep(1:4, each = 5))
 
   write_stan_json(list(x = table(f)), file = temp_file)
-  json_output <- readLines(temp_file)
-  expect_known_output(cat(json_output, sep = "\n"),
-                      file = test_path("answers", "json-table-vector.json"))
+  expect_snapshot_file(temp_file, "json-table-vector.json")
 
   write_stan_json(list(x = table(f, f)), file = temp_file)
-  json_output <- readLines(temp_file)
-  expect_known_output(cat(json_output, sep = "\n"),
-                      file = test_path("answers", "json-table-matrix.json"))
+  expect_snapshot_file(temp_file, "json-table-matrix.json")
 
   write_stan_json(list(x = table(f, f, f)), file = temp_file)
-  json_output <- readLines(temp_file)
-  expect_known_output(cat(json_output, sep = "\n"),
-                      file = test_path("answers", "json-table-array.json"))
+  expect_snapshot_file(temp_file, "json-table-array.json")
 })
 
 test_that("write_stan_json errors if NAs", {
@@ -198,25 +175,7 @@ test_that("write_stan_json() errors if bad names", {
 test_that("write_stan_json() works with always_decimal = TRUE", {
   test_file <- tempfile(fileext = ".json")
   write_stan_json(list(a = 1L, b = 2), test_file, always_decimal = FALSE)
-  expect_match(
-    "  \"a\": 1,",
-    readLines(test_file)[2],
-    fixed = TRUE
-  )
-  expect_match(
-    "  \"b\": 2",
-    readLines(test_file)[3],
-    fixed = TRUE
-  )
+  expect_snapshot_file(test_file, "json-always-decimal-false.json")
   write_stan_json(list(a = 1L, b = 2), test_file, always_decimal = TRUE)
-  expect_match(
-    "  \"a\": 1,",
-    readLines(test_file)[2],
-    fixed = TRUE
-  )
-  expect_match(
-    "  \"b\": 2.0",
-    readLines(test_file)[3],
-    fixed = TRUE
-  )
+  expect_snapshot_file(test_file, "json-always-decimal-true.json")
 })
diff --git a/tests/testthat/test-knitr.R b/tests/testthat/test-knitr.R
index 7e4400c5..6141ac9c 100644
--- a/tests/testthat/test-knitr.R
+++ b/tests/testthat/test-knitr.R
@@ -1,5 +1,3 @@
-context("knitr engine")
-
 test_that("eng_cmdstan throws correct errors", {
   skip_if_not_installed("knitr")
   expect_error(eng_cmdstan(list(output.var = 1)), "must be a character string")
@@ -24,7 +22,7 @@ test_that("eng_cmdstan works", {
   ))
   expect_interactive_message(eng_cmdstan(opts), "Compiling Stan program")
   opts$eval <- FALSE
-  expect_silent(eng_cmdstan(opts))
+  expect_noninteractive_silent(eng_cmdstan(opts))
 })
 
 test_that("register_knitr_engine works with and without override", {
diff --git a/tests/testthat/test-model-code-print.R b/tests/testthat/test-model-code-print.R
index e3409ed9..2a7ce963 100644
--- a/tests/testthat/test-model-code-print.R
+++ b/tests/testthat/test-model-code-print.R
@@ -1,13 +1,10 @@
-context("model-code-print")
-
 set_cmdstan_path()
 stan_program <- testing_stan_file("bernoulli")
 mod <- testing_model("bernoulli")
 
-
 test_that("code() and print() methods work", {
-  expect_known_output(mod$print(), file = test_path("answers", "model-print-output.stan"))
-  expect_known_value(mod$code(), file = test_path("answers", "model-code-output.rds"))
+  expect_snapshot_output(mod$print())
+  expect_snapshot_output(cat(mod$code(), sep = "\n"))
 })
 
 test_that("code() and print() still work if file is removed", {
@@ -64,23 +61,12 @@ test_that("code() doesn't change when file changes (unless model is recreated)",
 
 test_that("code() warns and print() errors if only exe and no Stan file", {
   mod_exe <- cmdstan_model(exe_file = mod$exe_file())
-  expect_warning(
-    expect_null(mod_exe$code()),
-    "'$code()' will return NULL because the 'CmdStanModel' was not created with a Stan file",
-    fixed = TRUE
-  )
-  expect_error(
-    mod_exe$print(),
-    "'$print()' cannot be used because the 'CmdStanModel' was not created with a Stan file.",
-    fixed = TRUE
-  )
+  expect_snapshot_warning(code <- mod_exe$code())
+  expect_null(code)
+  expect_snapshot_error(mod_exe$print())
 })
 
 test_that("check_syntax() errors if only exe and no Stan file", {
   mod_exe <- cmdstan_model(exe_file = mod$exe_file())
-  expect_error(
-    mod_exe$check_syntax(),
-    "'$check_syntax()' cannot be used because the 'CmdStanModel' was not created with a Stan file.",
-    fixed = TRUE
-  )
+  expect_snapshot_error(mod_exe$check_syntax())
 })
diff --git a/tests/testthat/test-model-compile.R b/tests/testthat/test-model-compile.R
index 4bc60a95..2e20f7ee 100644
--- a/tests/testthat/test-model-compile.R
+++ b/tests/testthat/test-model-compile.R
@@ -1,10 +1,21 @@
-context("model-compile")
-
 set_cmdstan_path()
 stan_program <- cmdstan_example_file()
 mod <- cmdstan_model(stan_file = stan_program, compile = FALSE)
 cmdstan_make_local(cpp_options = list("PRECOMPILED_HEADERS"="false"))
 
+trim_stanc_invocations <- function(output) {
+  output <- gsub("\\\\", "/", output)
+  out <- grep("bin/stanc", output, value = TRUE, fixed = TRUE)
+  out <- sub("^.*?(bin/stanc(?:\\.exe)?)\\b", "\\1", out, perl = TRUE)
+  sub("( --o).*", "\\1", out)
+}
+
+stanc_snapshot_transform <- function(lines) {
+  lines <- gsub("\\\\", "/", lines)
+  lines <- sub("^.*?(bin/stanc)(?:\\.exe)?\\b", "\\1", lines, perl = TRUE)
+  gsub("--name=(['\"])?([^ '\"=]+)\\1", "--name='\\2'", lines, perl = TRUE)
+}
+
 test_that("object initialized correctly", {
   expect_equal(mod$stan_file(), stan_program)
   expect_equal(mod$exe_file(), character(0))
@@ -76,9 +87,10 @@ test_that("compile() method works with spaces in path", {
 
 test_that("compile() method overwrites binaries", {
   mod$compile(quiet = TRUE)
-  old_time = file.mtime(mod$exe_file())
+  old_time <- file.mtime(mod$exe_file()) - 10
+  Sys.setFileTime(mod$exe_file(), old_time)
+  old_time <- file.mtime(mod$exe_file())
   mod$compile(quiet = TRUE, force_recompile = TRUE)
-  new_time =
   expect_gt(file.mtime(mod$exe_file()), old_time)
 })
 
@@ -113,17 +125,24 @@ test_that("compilation works with include_paths", {
 })
 
 test_that("name in STANCFLAGS is set correctly", {
+  local_reproducible_output()
   out <- utils::capture.output(mod$compile(quiet = FALSE, force_recompile = TRUE))
-  if(os_is_windows() && !os_is_wsl()) {
-    out_no_name <- "bin/stanc.exe --name='bernoulli_model' --o"
-    out_name <- "bin/stanc.exe --name='bernoulli2_model' --o"
-  } else {
-    out_no_name <- "bin/stanc --name='bernoulli_model' --o"
-    out_name <- "bin/stanc --name='bernoulli2_model' --o"
-  }
-  expect_output(print(out), out_no_name)
-  out <- utils::capture.output(mod$compile(quiet = FALSE, force_recompile = TRUE, stanc_options = list(name = "bernoulli2_model")))
-  expect_output(print(out), out_name)
+  expect_snapshot(
+    cat(trim_stanc_invocations(out), sep = "\n"),
+    transform = stanc_snapshot_transform
+  )
+
+  out <- utils::capture.output(
+    mod$compile(
+      quiet = FALSE,
+      force_recompile = TRUE,
+      stanc_options = list(name = "bernoulli2_model")
+    )
+  )
+  expect_snapshot(
+    cat(trim_stanc_invocations(out), sep = "\n"),
+    transform = stanc_snapshot_transform
+  )
 })
 
 
@@ -299,11 +318,9 @@ test_that("check_syntax() works", {
 
   stan_file <- testing_stan_file("bernoulli")
   mod_ok <- cmdstan_model(stan_file, compile = FALSE)
-  expect_true(
-    expect_message(
-      mod_ok$check_syntax(),
-      "Stan program is syntactically correct"
-    )
+  expect_message(
+    mod_ok$check_syntax(),
+    "Stan program is syntactically correct"
   )
   expect_message(
     mod_ok$check_syntax(quiet = TRUE),
@@ -832,9 +849,9 @@ test_that("dirname of stan_file is used as include path if no other paths suppli
 })
 
 test_that("STANCFLAGS from get_cmdstan_flags() are included in compile output", {
+  local_reproducible_output()
   real_get_cmdstan_flags <- get_cmdstan_flags
-  out <- with_mocked_bindings(
-    utils::capture.output(mod$compile(quiet = FALSE, force_recompile = TRUE)),
+  local_mocked_bindings(
     get_cmdstan_flags = function(flag_name) {
       if (identical(flag_name, "STANCFLAGS")) {
         c("--O1", "--warn-pedantic")
@@ -843,12 +860,11 @@ test_that("STANCFLAGS from get_cmdstan_flags() are included in compile output",
       }
     }
   )
-  if(os_is_windows() && !os_is_wsl()) {
-    out_w_flags <- "bin/stanc.exe --name='bernoulli_model'[[:space:]]+--O1[[:space:]]+--warn-pedantic[[:space:]]+--o"
-  } else {
-    out_w_flags <- "bin/stanc --name='bernoulli_model'[[:space:]]+--O1[[:space:]]+--warn-pedantic[[:space:]]+--o"
-  }
-  expect_output(print(out), out_w_flags)
+  out <- utils::capture.output(mod$compile(quiet = FALSE, force_recompile = TRUE))
+  expect_snapshot(
+    cat(trim_stanc_invocations(out), sep = "\n"),
+    transform = stanc_snapshot_transform
+  )
 })
 
 test_that("compile() ignores directory chatter from MAKEFLAGS when reading STANCFLAGS", {
diff --git a/tests/testthat/test-model-data.R b/tests/testthat/test-model-data.R
index 2f91e575..1e1c2ed6 100644
--- a/tests/testthat/test-model-data.R
+++ b/tests/testthat/test-model-data.R
@@ -1,4 +1,3 @@
-context("model-data")
 # see separate test-json for testing writing data to JSON
 
 set_cmdstan_path()
@@ -33,5 +32,5 @@ test_that("empty data list doesn't error if no data block", {
   )
 
   # would error if fitting failed
-  expect_silent(fit$draws())
+  expect_no_error(fit$draws())
 })
diff --git a/tests/testthat/test-model-diagnose.R b/tests/testthat/test-model-diagnose.R
index 76ab4a47..0ef2db68 100644
--- a/tests/testthat/test-model-diagnose.R
+++ b/tests/testthat/test-model-diagnose.R
@@ -1,5 +1,3 @@
-context("model-diagnose")
-
 set_cmdstan_path()
 mod <- testing_model("bernoulli")
 data_list <- testing_data("bernoulli")
@@ -34,11 +32,11 @@ ok_arg_sci_nota_values <- list(
 test_that("diagnose() method runs when all arguments specified validly", {
   # specifying all arguments validly
   fit1 <- do.call(mod$diagnose, ok_arg_values)
-  expect_is(fit1, "CmdStanDiagnose")
+  expect_s3_class(fit1, "CmdStanDiagnose")
 
   # leaving all at default (except 'data' and 'seed')
   fit2 <- mod$diagnose(data = data_list, seed = 123)
-  expect_is(fit2, "CmdStanDiagnose")
+  expect_s3_class(fit2, "CmdStanDiagnose")
 })
 
 test_that("diagnose() method runs when arguments are specified in scientific notation", {
@@ -46,7 +44,7 @@ test_that("diagnose() method runs when arguments are specified in scientific not
 
   # specifying all arguments validly
   fit1 <- do.call(mod$diagnose, ok_arg_sci_nota_values)
-  expect_is(fit1, "CmdStanDiagnose")
+  expect_s3_class(fit1, "CmdStanDiagnose")
 })
 
 test_that("diagnose() method errors for any invalid argument before calling cmdstan", {
diff --git a/tests/testthat/test-model-expose-functions.R b/tests/testthat/test-model-expose-functions.R
index 9bcf1382..0cbc4f8b 100644
--- a/tests/testthat/test-model-expose-functions.R
+++ b/tests/testthat/test-model-expose-functions.R
@@ -1,5 +1,3 @@
-context("model-expose-functions")
-
 # Standalone functions not expected to work on WSL yet
 skip_if(os_is_wsl())
 
diff --git a/tests/testthat/test-model-generate_quantities.R b/tests/testthat/test-model-generate_quantities.R
index 7c641815..ec1924aa 100644
--- a/tests/testthat/test-model-generate_quantities.R
+++ b/tests/testthat/test-model-generate_quantities.R
@@ -1,5 +1,3 @@
-context("model-generate-quantities")
-
 set_cmdstan_path()
 fit <- testing_fit("bernoulli", method = "sample", seed = 123)
 mod_gq <- testing_model("bernoulli_ppc")
@@ -25,11 +23,11 @@ bad_arg_values <- list(
 test_that("generate_quantities() method runs when all arguments specified validly", {
   # specifying all arguments validly
   expect_gq_output(fit1 <- do.call(mod_gq$generate_quantities, ok_arg_values))
-  expect_is(fit1, "CmdStanGQ")
+  expect_s3_class(fit1, "CmdStanGQ")
 
   # leaving all at default (except 'data')
   expect_gq_output(fit2 <- mod_gq$generate_quantities(fitted_params = fit, data = data_list))
-  expect_is(fit2, "CmdStanGQ")
+  expect_s3_class(fit2, "CmdStanGQ")
 })
 
 test_that("generate_quantities() method errors for any invalid argument before calling cmdstan", {
diff --git a/tests/testthat/test-model-init.R b/tests/testthat/test-model-init.R
index ebe54677..e2e30dfc 100644
--- a/tests/testthat/test-model-init.R
+++ b/tests/testthat/test-model-init.R
@@ -1,5 +1,3 @@
-context("model-init")
-
 set_cmdstan_path()
 mod <- testing_model("bernoulli")
 data_list <- testing_data("bernoulli")
@@ -240,7 +238,7 @@ test_that("error if init function specified incorrectly", {
 })
 
 test_that("print message if not all parameters are initialized", {
-  options(cmdstanr_warn_inits = NULL) # should default to TRUE
+  withr::local_options(list(cmdstanr_warn_inits = NULL)) # should default to TRUE
   init_list <- list(
     list(
       alpha = 1
@@ -273,7 +271,7 @@ test_that("print message if not all parameters are initialized", {
 })
 
 test_that("No message printed if options(cmdstanr_warn_inits=FALSE)", {
-  options(cmdstanr_warn_inits = FALSE)
+  withr::local_options(list(cmdstanr_warn_inits = FALSE))
   expect_message(
     utils::capture.output(mod_logistic$optimize(data = data_list_logistic, init = list(list(a = 0)), seed = 123)),
     regexp = NA
@@ -286,7 +284,6 @@ test_that("No message printed if options(cmdstanr_warn_inits=FALSE)", {
     utils::capture.output(mod_logistic$sample(data = data_list_logistic, init = list(list(alpha = 1),list(alpha = 1)), chains = 2, seed = 123)),
     regexp = NA
   )
-  options(cmdstanr_warn_inits = TRUE)
 })
 
 test_that("Initial values for single-element containers treated correctly", {
diff --git a/tests/testthat/test-model-laplace.R b/tests/testthat/test-model-laplace.R
index 00dc78bb..961883c8 100644
--- a/tests/testthat/test-model-laplace.R
+++ b/tests/testthat/test-model-laplace.R
@@ -1,5 +1,3 @@
-context("model-laplace")
-
 set_cmdstan_path()
 mod <- testing_model("logistic")
 data_list <- testing_data("logistic")
@@ -48,7 +46,7 @@ test_that("laplace() method errors for any invalid argument before calling cmdst
 test_that("laplace() runs when all arguments specified validly", {
   # specifying all arguments validly
   expect_laplace_output(fit1 <- do.call(mod$laplace, ok_arg_values))
-  expect_is(fit1, "CmdStanLaplace")
+  expect_s3_class(fit1, "CmdStanLaplace")
 
   # check that correct arguments were indeed passed to CmdStan
   expect_equal(fit1$metadata()$refresh, ok_arg_values$refresh)
@@ -61,7 +59,7 @@ test_that("laplace() runs when all arguments specified validly", {
 
   # leaving all at default (except 'data')
   expect_laplace_output(fit2 <- mod$laplace(data = data_list, seed = 123))
-  expect_is(fit2, "CmdStanLaplace")
+  expect_s3_class(fit2, "CmdStanLaplace")
 })
 
 test_that("laplace() all valid 'mode' inputs give same results", {
@@ -72,12 +70,12 @@ test_that("laplace() all valid 'mode' inputs give same results", {
     fit3 <- mod$laplace(data = data_list, mode = NULL, seed = 100, refresh = 0)
   })
 
-  expect_is(fit1, "CmdStanLaplace")
-  expect_is(fit2, "CmdStanLaplace")
-  expect_is(fit3, "CmdStanLaplace")
-  expect_is(fit1$mode(), "CmdStanMLE")
-  expect_is(fit2$mode(), "CmdStanMLE")
-  expect_is(fit3$mode(), "CmdStanMLE")
+  expect_s3_class(fit1, "CmdStanLaplace")
+  expect_s3_class(fit2, "CmdStanLaplace")
+  expect_s3_class(fit3, "CmdStanLaplace")
+  expect_s3_class(fit1$mode(), "CmdStanMLE")
+  expect_s3_class(fit2$mode(), "CmdStanMLE")
+  expect_s3_class(fit3$mode(), "CmdStanMLE")
   expect_equal(fit1$mode()$mle(), fit2$mode()$mle())
   expect_equal(fit1$mode()$mle(), fit3$mode()$mle())
   expect_equal(fit1$lp(), fit2$lp())
diff --git a/tests/testthat/test-model-methods.R b/tests/testthat/test-model-methods.R
index 1e38ad25..c502af15 100644
--- a/tests/testthat/test-model-methods.R
+++ b/tests/testthat/test-model-methods.R
@@ -1,4 +1,3 @@
-context("model-methods")
 skip_if(os_is_wsl())
 
 set_cmdstan_path()
diff --git a/tests/testthat/test-model-optimize.R b/tests/testthat/test-model-optimize.R
index a3bd89c6..204a5b73 100644
--- a/tests/testthat/test-model-optimize.R
+++ b/tests/testthat/test-model-optimize.R
@@ -1,5 +1,3 @@
-context("model-optimize")
-
 set_cmdstan_path()
 mod <- testing_model("bernoulli")
 data_list <- testing_data("bernoulli")
@@ -45,17 +43,17 @@ ok_arg_sci_nota_values <- list(
 test_that("optimize() method runs when all arguments specified validly", {
   # specifying all arguments validly
   expect_optim_output(fit1 <- do.call(mod$optimize, ok_arg_values))
-  expect_is(fit1, "CmdStanMLE")
+  expect_s3_class(fit1, "CmdStanMLE")
 
   # leaving all at default (except 'data')
   expect_optim_output(fit2 <- mod$optimize(data = data_list, seed = 123))
-  expect_is(fit2, "CmdStanMLE")
+  expect_s3_class(fit2, "CmdStanMLE")
 })
 
 test_that("optimize() method runs when arguments are specified in scientific notation", {
   # specifying all arguments validly
   expect_optim_output(fit1 <- do.call(mod$optimize, ok_arg_sci_nota_values))
-  expect_is(fit1, "CmdStanMLE")
+  expect_s3_class(fit1, "CmdStanMLE")
 })
 
 test_that("optimize() warns if threads specified but not enabled", {
diff --git a/tests/testthat/test-model-output_dir.R b/tests/testthat/test-model-output_dir.R
index f5284827..0fceac4e 100644
--- a/tests/testthat/test-model-output_dir.R
+++ b/tests/testthat/test-model-output_dir.R
@@ -1,5 +1,3 @@
-context("model-output_dir-output-basename")
-
 set_cmdstan_path()
 sandbox <- file.path(tempdir(check = TRUE), "sandbox")
 if (!dir.exists(sandbox)) {
diff --git a/tests/testthat/test-model-pathfinder.R b/tests/testthat/test-model-pathfinder.R
index 58058010..7d4747f9 100644
--- a/tests/testthat/test-model-pathfinder.R
+++ b/tests/testthat/test-model-pathfinder.R
@@ -1,5 +1,3 @@
-context("model-pathfinder")
-
 set_cmdstan_path()
 stan_program <- testing_stan_file("bernoulli")
 mod <- testing_model("bernoulli")
@@ -103,15 +101,15 @@ expect_pathfinder_output <- function(object, num_chains = NULL) {
 
 test_that("Pathfinder Runs", {
   expect_pathfinder_output(fit <- mod$pathfinder(data=data_list, seed=1234, refresh = 0))
-  expect_is(fit, "CmdStanPathfinder")
+  expect_s3_class(fit, "CmdStanPathfinder")
 })
 
 test_that("pathfinder() method works with data files", {
   expect_pathfinder_output(fit_r <- mod$pathfinder(data = data_file_r))
-  expect_is(fit_r, "CmdStanPathfinder")
+  expect_s3_class(fit_r, "CmdStanPathfinder")
 
   expect_pathfinder_output(fit_json <- mod$pathfinder(data = data_file_json))
-  expect_is(fit_json, "CmdStanPathfinder")
+  expect_s3_class(fit_json, "CmdStanPathfinder")
 })
 
 test_that("pathfinder() method works with init file", {
@@ -132,7 +130,7 @@ test_that("pathfinder() method works with init function and default paths", {
 
 test_that("pathfinder() method runs when all arguments specified", {
   expect_pathfinder_output(fit <- do.call(mod$pathfinder, ok_arg_values))
-  expect_is(fit, "CmdStanPathfinder")
+  expect_s3_class(fit, "CmdStanPathfinder")
 })
 
 test_that("pathfinder() method runs when the stan file is removed", {
@@ -147,7 +145,7 @@ test_that("pathfinder() method runs when the stan file is removed", {
 
 test_that("no error when checking estimates after failure", {
   fit <- cmdstanr_example("schools", method = "pathfinder", seed = 123) # optim always fails for this
-  expect_silent(fit$summary()) # no error
+  expect_no_error(fit$summary())
 })
 
 test_that("no output with show_messages = FALSE", {
@@ -156,4 +154,3 @@ test_that("no output with show_messages = FALSE", {
   )
   expect_equal(length(output), 0)
 })
-
diff --git a/tests/testthat/test-model-sample-metric.R b/tests/testthat/test-model-sample-metric.R
index 422442fa..5789e583 100644
--- a/tests/testthat/test-model-sample-metric.R
+++ b/tests/testthat/test-model-sample-metric.R
@@ -1,5 +1,3 @@
-context("model-sample-metric")
-
 set_cmdstan_path()
 mod <- testing_model("bernoulli")
 data_list <- testing_data("bernoulli")
diff --git a/tests/testthat/test-model-sample.R b/tests/testthat/test-model-sample.R
index f9ce48b4..85c5dbb8 100644
--- a/tests/testthat/test-model-sample.R
+++ b/tests/testthat/test-model-sample.R
@@ -1,5 +1,3 @@
-context("model-sample")
-
 set_cmdstan_path()
 stan_program <- testing_stan_file("bernoulli")
 mod <- testing_model("bernoulli")
@@ -85,15 +83,15 @@ bad_arg_values_3 <- list(
 
 test_that("sample() method works with data list", {
   expect_sample_output(fit <- mod$sample(data = data_list, chains = 1), 1)
-  expect_is(fit, "CmdStanMCMC")
+  expect_s3_class(fit, "CmdStanMCMC")
 })
 
 test_that("sample() method works with data files", {
   expect_sample_output(fit_r <- mod$sample(data = data_file_r, chains = 1), 1)
-  expect_is(fit_r, "CmdStanMCMC")
+  expect_s3_class(fit_r, "CmdStanMCMC")
 
   expect_sample_output(fit_json <- mod$sample(data = data_file_json, chains = 1), 1)
-  expect_is(fit_json, "CmdStanMCMC")
+  expect_s3_class(fit_json, "CmdStanMCMC")
 })
 
 test_that("sample() method works with init file", {
@@ -109,7 +107,7 @@ test_that("sample() method works with init file", {
 
 test_that("sample() method runs when all arguments specified", {
   expect_sample_output(fit <- do.call(mod$sample, ok_arg_values), 2)
-  expect_is(fit, "CmdStanMCMC")
+  expect_s3_class(fit, "CmdStanMCMC")
 })
 
 test_that("sample() method runs when the stan file is removed", {
@@ -179,14 +177,14 @@ test_that("sampling in parallel works", {
 })
 
 test_that("mc.cores option detected", {
-  options(mc.cores = 3)
+  withr::local_options(list(mc.cores = 3))
   expect_output(
     mod$sample(data = data_list, chains = 3),
     "Running MCMC with 3 parallel chains",
     fixed = TRUE
   )
 
-  options(mc.cores = NULL)
+  withr::local_options(list(mc.cores = NULL))
   expect_output(
     mod$sample(data = data_list, chains = 2),
     "Running MCMC with 2 sequential chains",
@@ -198,7 +196,7 @@ test_that("sample() method runs when fixed_param = TRUE", {
   mod_fp$compile()
 
   expect_sample_output(fit_1000 <- mod_fp$sample(fixed_param = TRUE, iter_sampling = 1000), 4)
-  expect_is(fit_1000, "CmdStanMCMC")
+  expect_s3_class(fit_1000, "CmdStanMCMC")
   expect_equal(dim(fit_1000$draws()), c(1000,4,10))
 
   expect_sample_output(fit_500 <- mod_fp$sample(fixed_param = TRUE, iter_sampling = 500), 4)
@@ -221,15 +219,15 @@ test_that("sample() method runs when adapt_engaged = FALSE", {
 test_that("chain_ids work with sample()", {
   mod$compile()
   expect_sample_output(fit12 <- mod$sample(data = data_list, chains = 2, chain_ids = c(10,12)))
-  expect_is(fit12, "CmdStanMCMC")
+  expect_s3_class(fit12, "CmdStanMCMC")
   expect_equal(fit12$metadata()$id, c(10,12))
 
   expect_sample_output(fit12 <- mod$sample(data = data_list, chains = 2, chain_ids = c(100,7)))
-  expect_is(fit12, "CmdStanMCMC")
+  expect_s3_class(fit12, "CmdStanMCMC")
   expect_equal(fit12$metadata()$id, c(100,7))
 
   expect_sample_output(fit12 <- mod$sample(data = data_list, chains = 1, chain_ids = c(6)))
-  expect_is(fit12, "CmdStanMCMC")
+  expect_s3_class(fit12, "CmdStanMCMC")
   expect_equal(fit12$metadata()$id, c(6))
 
   expect_error(mod$sample(data = data_list, chains = 1, chain_ids = c(0)),
@@ -369,7 +367,7 @@ test_that("All output can be suppressed by show_messages", {
   stan_program <- testing_stan_file("bernoulli")
   data_list <- testing_data("bernoulli")
   mod <- cmdstan_model(stan_program, force_recompile = TRUE)
-  options("cmdstanr_verbose" = FALSE)
+  withr::local_options(list("cmdstanr_verbose" = FALSE))
   output <- capture.output(
     fit <- mod$sample(data = data_list, show_messages = FALSE)
   )
diff --git a/tests/testthat/test-model-sample_mpi.R b/tests/testthat/test-model-sample_mpi.R
index 1ea0f9a6..9463cbbb 100644
--- a/tests/testthat/test-model-sample_mpi.R
+++ b/tests/testthat/test-model-sample_mpi.R
@@ -1,5 +1,3 @@
-context("model-sample_mpi")
-
 test_that("sample_mpi() works", {
   skip_if(!mpi_toolchain_present())
   mpi_file <- write_stan_file("
@@ -35,9 +33,11 @@ test_that("sample_mpi() works", {
   if (os_is_wsl()) {
     # Default GHA WSL install runs as root, which MPI discourages
     # Specify that this is safe to ignore for this test
-    Sys.setenv("OMPI_ALLOW_RUN_AS_ROOT"=1)
-    Sys.setenv("OMPI_ALLOW_RUN_AS_ROOT_CONFIRM"=1)
-    Sys.setenv("WSLENV"="OMPI_ALLOW_RUN_AS_ROOT/u:OMPI_ALLOW_RUN_AS_ROOT_CONFIRM/u")
+    withr::local_envvar(c(
+      "OMPI_ALLOW_RUN_AS_ROOT" = "1",
+      "OMPI_ALLOW_RUN_AS_ROOT_CONFIRM" = "1",
+      "WSLENV" = "OMPI_ALLOW_RUN_AS_ROOT/u:OMPI_ALLOW_RUN_AS_ROOT_CONFIRM/u"
+    ))
   }
 
   utils::capture.output(
diff --git a/tests/testthat/test-model-variables.R b/tests/testthat/test-model-variables.R
index 5ca43ef2..ae71d29f 100644
--- a/tests/testthat/test-model-variables.R
+++ b/tests/testthat/test-model-variables.R
@@ -1,5 +1,3 @@
-context("model-variables")
-
 set_cmdstan_path()
 
 test_that("$variables() work correctly with example models", {
diff --git a/tests/testthat/test-model-variational.R b/tests/testthat/test-model-variational.R
index 6c917d1e..25b1fc77 100644
--- a/tests/testthat/test-model-variational.R
+++ b/tests/testthat/test-model-variational.R
@@ -1,5 +1,3 @@
-context("model-variational")
-
 set_cmdstan_path()
 mod <- testing_model("bernoulli")
 data_list <- testing_data("bernoulli")
@@ -46,11 +44,11 @@ bad_arg_values <- list(
 test_that("variational() method runs when all arguments specified validly", {
   # specifying all arguments validly
   expect_vb_output(fit1 <- do.call(mod$variational, ok_arg_values))
-  expect_is(fit1, "CmdStanVB")
+  expect_s3_class(fit1, "CmdStanVB")
 
   # leaving all at default (except data and seed)
   expect_vb_output(fit2 <- mod$variational(data = data_list, seed = 123))
-  expect_is(fit2, "CmdStanVB")
+  expect_s3_class(fit2, "CmdStanVB")
 })
 
 test_that("variational() warns if threads specified but not enabled", {
diff --git a/tests/testthat/test-opencl.R b/tests/testthat/test-opencl.R
index 92858141..44f5fc77 100644
--- a/tests/testthat/test-opencl.R
+++ b/tests/testthat/test-opencl.R
@@ -1,5 +1,3 @@
-context("opencl")
-
 set_cmdstan_path()
 fit <- testing_fit("bernoulli", method = "sample", seed = 123, chains = 1)
 
@@ -128,4 +126,3 @@ test_that("all methods run with valid opencl_ids", {
   expect_false(is.null(fit$metadata()$device))
   expect_false(is.null(fit$metadata()$platform))
 })
-
diff --git a/tests/testthat/test-path.R b/tests/testthat/test-path.R
index e9e7eda0..cc01acc5 100644
--- a/tests/testthat/test-path.R
+++ b/tests/testthat/test-path.R
@@ -1,5 +1,3 @@
-context("paths")
-
 Sys.unsetenv("CMDSTAN")
 PATH <- absolute_path(set_cmdstan_path())
 VERSION <- cmdstan_version()
@@ -19,34 +17,26 @@ test_that("Setting path works and confirms with message", {
 test_that("Setting bad path leads to warning (can't find directory)", {
   unset_cmdstan_path()
   expect_null(.cmdstanr$PATH)
-  expect_warning(
-    set_cmdstan_path("BAD_PATH"),
-    "Can't find directory"
-  )
+  expect_snapshot_warning(set_cmdstan_path("BAD_PATH"))
 })
 
 test_that("Setting bad path from env leads to warning (can't find directory)", {
   unset_cmdstan_path()
   .cmdstanr$WSL <- TRUE
-  Sys.setenv(CMDSTAN = "BAD_PATH")
-  expect_warning(
-    cmdstanr_initialize(),
-    "Can't find directory specified by environment variable"
-  )
+  withr::local_envvar(c(CMDSTAN = "BAD_PATH"))
+  expect_snapshot_warning(cmdstanr_initialize())
   expect_null(.cmdstanr$PATH)
   expect_null(.cmdstanr$VERSION)
   expect_false(isTRUE(.cmdstanr$WSL))
-  Sys.unsetenv("CMDSTAN")
 })
 
 test_that("Setting path from env var is detected", {
   unset_cmdstan_path()
   expect_true(is.null(.cmdstanr$VERSION))
-  Sys.setenv(CMDSTAN = PATH)
+  withr::local_envvar(c(CMDSTAN = PATH))
   expect_silent(cmdstanr_initialize())
   expect_equal(cmdstan_path(), PATH)
   expect_false(is.null(.cmdstanr$VERSION))
-  Sys.unsetenv("CMDSTAN")
 })
 
 test_that("Unsupported CmdStan path from env var is rejected", {
@@ -56,10 +46,9 @@ test_that("Unsupported CmdStan path from env var is rejected", {
   old_install <- file.path(parent_dir, "cmdstan-2.34.0")
   dir.create(old_install, recursive = TRUE, showWarnings = FALSE)
   on.exit(unlink(parent_dir, recursive = TRUE), add = TRUE)
-  on.exit(Sys.unsetenv("CMDSTAN"), add = TRUE)
   writeLines("CMDSTAN_VERSION := 2.34.0", con = file.path(old_install, "makefile"))
 
-  Sys.setenv(CMDSTAN = parent_dir)
+  withr::local_envvar(c(CMDSTAN = parent_dir))
   suppressWarnings(cmdstanr_initialize())
   expect_false(identical(.cmdstanr$PATH, absolute_path(old_install)))
   expect_false(identical(.cmdstanr$VERSION, "2.34.0"))
@@ -77,14 +66,9 @@ test_that("Existing CMDSTAN env path with no install resets cached state", {
   empty_parent <- file.path(tempdir(check = TRUE), "cmdstan-empty-parent")
   dir.create(empty_parent, recursive = TRUE, showWarnings = FALSE)
   on.exit(unlink(empty_parent, recursive = TRUE), add = TRUE)
-  on.exit(Sys.unsetenv("CMDSTAN"), add = TRUE)
 
-  Sys.setenv(CMDSTAN = empty_parent)
-  expect_warning(
-    cmdstanr_initialize(),
-    "No CmdStan installation found in the path specified by the environment variable 'CMDSTAN'.",
-    fixed = TRUE
-  )
+  withr::local_envvar(c(CMDSTAN = empty_parent))
+  expect_snapshot_warning(cmdstanr_initialize())
   expect_null(.cmdstanr$PATH)
   expect_null(.cmdstanr$VERSION)
   expect_false(isTRUE(.cmdstanr$WSL))
@@ -106,10 +90,7 @@ test_that("Getting a valid path works", {
 
 test_that("Getting missing path leads to error (path not set)", {
   unset_cmdstan_path()
-  expect_error(
-    cmdstan_path(),
-    "CmdStan path has not been set yet"
-  )
+  expect_snapshot_error(cmdstan_path())
   expect_null(.cmdstanr$PATH)
 })
 
@@ -124,16 +105,13 @@ test_that("cmdstan_version() behaves correctly when version is not set", {
   version <- .cmdstanr$VERSION
   on.exit(.cmdstanr$VERSION <- version)
   .cmdstanr$VERSION <- NULL
-  expect_error(cmdstan_version())
+  expect_snapshot_error(cmdstan_version())
   expect_null(cmdstan_version(error_on_NA = FALSE))
 })
 
 test_that("Warning message is thrown if can't detect version number", {
   path <- testthat::test_path("answers") # valid path but not cmdstan
-  expect_warning(
-    set_cmdstan_path(path),
-    "Can't find CmdStan makefile to detect version number"
-  )
+  expect_snapshot_warning(set_cmdstan_path(path))
 })
 
 test_that("Setting path rejects unsupported CmdStan versions", {
@@ -151,11 +129,7 @@ test_that("Setting path rejects unsupported CmdStan versions", {
   on.exit(unlink(path, recursive = TRUE), add = TRUE)
   writeLines("CMDSTAN_VERSION := 2.34.0", con = file.path(path, "makefile"))
 
-  expect_warning(
-    set_cmdstan_path(path),
-    "cmdstanr now requires CmdStan v2.35.0 or newer",
-    fixed = TRUE
-  )
+  expect_snapshot_warning(set_cmdstan_path(path))
   expect_null(.cmdstanr$PATH)
   expect_null(.cmdstanr$VERSION)
   expect_false(isTRUE(.cmdstanr$WSL))
diff --git a/tests/testthat/test-profiling.R b/tests/testthat/test-profiling.R
index 638f1fdc..dd45c45c 100644
--- a/tests/testthat/test-profiling.R
+++ b/tests/testthat/test-profiling.R
@@ -1,5 +1,3 @@
-context("profiling")
-
 set_cmdstan_path()
 
 
diff --git a/tests/testthat/test-threads.R b/tests/testthat/test-threads.R
index fb5eec61..5e9f04f7 100644
--- a/tests/testthat/test-threads.R
+++ b/tests/testthat/test-threads.R
@@ -1,5 +1,3 @@
-context("threading")
-
 set_cmdstan_path()
 stan_program <- testing_stan_file("bernoulli")
 stan_gq_program <- testing_stan_file("bernoulli_ppc")
diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R
index a2c39f15..56797f7f 100644
--- a/tests/testthat/test-utils.R
+++ b/tests/testthat/test-utils.R
@@ -1,5 +1,3 @@
-context("utils")
-
 set_cmdstan_path()
 fit_mcmc <- testing_fit("logistic", method = "sample",
                         seed = 123, chains = 2)
@@ -110,12 +108,12 @@ test_that("check_ebfmi and computing ebfmi works", {
   expect_equal(as.numeric(ebfmi(posterior::as_draws_list(energy_df))), check_val)
   expect_equal(as.numeric(ebfmi(posterior::as_draws_matrix(energy_df))), check_val)
   energy_df <- posterior::as_draws(data.frame("energy__" = 0))
-  expect_warning(check_ebfmi(energy_df), "E-BFMI not computed because it is undefined for posterior chains of length less than 3.")
-  expect_warning(ebfmi(energy_df), "E-BFMI not computed because it is undefined for posterior chains of length less than 3.")
+  expect_snapshot_warning(check_ebfmi(energy_df))
+  expect_snapshot_warning(ebfmi(energy_df))
 
   energy_df <- posterior::as_draws(data.frame("somethingelse" = 0))
-  expect_warning(check_ebfmi(energy_df), "E-BFMI not computed because the 'energy__' diagnostic could not be located.")
-  expect_warning(ebfmi(energy_df), "E-BFMI not computed because the 'energy__' diagnostic could not be located.")
+  expect_snapshot_warning(check_ebfmi(energy_df))
+  expect_snapshot_warning(ebfmi(energy_df))
 })