diff --git a/DESCRIPTION b/DESCRIPTION
index a90548193..948779c53 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -72,6 +72,7 @@ Suggests:
     doRNG,
     unix,
     testthat,
+    withr (>= 3.0.3),
     vdiffr,
     knitr,
     rmarkdown,
diff --git a/NEWS.md b/NEWS.md
index 94ffe21cd..bd31bcf95 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -11,6 +11,8 @@ If you read this from a place other than <https://mc-stan.org/projpred/news/inde
 
 ## Bug fixes and maintenance
 
+* Tests: Leave the options, the random seed and the global environment of the R session unchanged after the tests, also when a test fails, and remove the manual cleanup of the test workspace. (GitHub: #586)
+
 
 # projpred 2.11.0
 
diff --git a/tests/testthat/helpers/skips.R b/tests/testthat/helpers/skips.R
new file mode 100644
index 000000000..3eb8ee1f7
--- /dev/null
+++ b/tests/testthat/helpers/skips.R
@@ -0,0 +1,9 @@
+# Skips -----------------------------------------------------------------------
+
+# Skip a test that only loops over the test setups in `x` (for example
+# `refmods`) if `x` is empty. On CRAN, the setup fits no reference models, so
+# such a test would run no expectation; the skip reports the reason.
+skip_if_no_setups <- function(x) {
+  skip_if(length(x) == 0,
+          paste0("`", deparse(substitute(x)), "` is empty (no test setups)"))
+}
diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R
index 74bcef4b9..616daf734 100755
--- a/tests/testthat/setup.R
+++ b/tests/testthat/setup.R
@@ -4,6 +4,17 @@
 
 # General setup -----------------------------------------------------------
 
+# When the tests end, restore the global options (and remove the ones that the
+# tests added). testthat restores its own options:
+opts_orig <- options()
+opts_orig <- opts_orig[!startsWith(names(opts_orig), "testthat")]
+withr::defer({
+  opts_new <- setdiff(names(options()), names(opts_orig))
+  opts_new <- opts_new[!startsWith(opts_new, "testthat")]
+  options(stats::setNames(vector("list", length(opts_new)), opts_new))
+  options(opts_orig)
+}, envir = testthat::teardown_env())
+
 options(warn = 1)
 
 # These switches may be set to `FALSE` to save time (e.g., when debugging
@@ -156,6 +167,9 @@ use_fac <- TRUE
 # Use polym() instead of poly()?:
 use_polym <- FALSE
 
+source(testthat::test_path("helpers", "skips.R"), local = TRUE)
+# When the tests end, restore the state of the random number generator:
+withr::local_preserve_seed(.local_envir = testthat::teardown_env())
 source(testthat::test_path("helpers", "unlist_cust.R"), local = TRUE)
 source(testthat::test_path("helpers", "testers.R"), local = TRUE)
 source(testthat::test_path("helpers", "args.R"), local = TRUE)
@@ -617,12 +631,17 @@ offss_tst <- list(with_offs = list(), # with_offs = list(offset = offs_tst),
 # rstanarm:::kfold.stanreg() seems to assume these objects to lie in the
 # global environment. Since `testthat` uses a new environment for running
 # the tests (see `?testthat::test_env`), we need the following code to be
-# able to run devtools::test():
+# able to run devtools::test(). The withr::defer() below removes them again
+# when the tests end:
+f_objs_global <- character()
 for (obj_symb_chr in c(paste0("f_", fam_nms_trad))) {
   if (!exists(obj_symb_chr, envir = .GlobalEnv)) {
     assign(obj_symb_chr, get(obj_symb_chr), envir = .GlobalEnv)
+    f_objs_global <- c(f_objs_global, obj_symb_chr)
   }
 }
+withr::defer(rm(list = f_objs_global, envir = .GlobalEnv),
+             envir = testthat::teardown_env())
 
 args_fit <- lapply(pkg_nms, function(pkg_nm) {
   # Depending on the package used for fitting the reference model, `mod_nms`
diff --git a/tests/testthat/test_as_matrix.R b/tests/testthat/test_as_matrix.R
index 484f3a0a2..3b768f399 100644
--- a/tests/testthat/test_as_matrix.R
+++ b/tests/testthat/test_as_matrix.R
@@ -1,5 +1,3 @@
-context("as.matrix.projection()")
-
 test_that("as.matrix.projection() works", {
   skip_if_not(run_prj)
   for (tstsetup in names(prjs)) {
@@ -221,12 +219,12 @@ test_that("as.matrix.projection() works", {
     ###
     if (run_snaps) {
       if (testthat_ed_max2) local_edition(3)
-      width_orig <- options(width = 145)
-      expect_snapshot({
-        print(tstsetup)
-        print(rlang::hash(m)) # message(m)
+      withr::with_options(list(width = 145), {
+        expect_snapshot({
+          print(tstsetup)
+          print(rlang::hash(m)) # message(m)
+        })
       })
-      options(width_orig)
       if (testthat_ed_max2) local_edition(2)
     }
   }
@@ -306,11 +304,11 @@ test_that(paste(
 
 if (run_snaps) {
   if (testthat_ed_max2) local_edition(3)
-  width_orig <- options(width = 145)
 
   test_that(paste(
     "as.matrix.projection() works for projections based on varsel() output"
   ), {
+    withr::local_options(list(width = 145))
     skip_if_not(run_vs)
     for (tstsetup in names(prjs_vs)) {
       if (args_prj_vs[[tstsetup]]$mod_nm == "gam") {
@@ -350,6 +348,7 @@ if (run_snaps) {
   test_that(paste(
     "as.matrix.projection() works for projections based on cv_varsel() output"
   ), {
+    withr::local_options(list(width = 145))
     skip_if_not(run_cvvs)
     for (tstsetup in names(prjs_cvvs)) {
       if (args_prj_cvvs[[tstsetup]]$mod_nm == "gam") {
@@ -386,6 +385,5 @@ if (run_snaps) {
     }
   })
 
-  options(width_orig)
   if (testthat_ed_max2) local_edition(2)
 }
diff --git a/tests/testthat/test_augdat.R b/tests/testthat/test_augdat.R
index 35dcbc8ba..0968b3d33 100644
--- a/tests/testthat/test_augdat.R
+++ b/tests/testthat/test_augdat.R
@@ -1,10 +1,3 @@
-context("augmented-data projection")
-
-# Setup -------------------------------------------------------------------
-
-# Needed to clean up the workspace afterwards (i.e, after this test file):
-ls_bu <- ls()
-
 # augdat.R ----------------------------------------------------------------
 
 nobs_orig_tst <- 4L
@@ -188,9 +181,7 @@ test_that(paste(
   "the traditional projection"
 ), {
   skip_if_not(run_prj)
-  if (exists(".Random.seed", envir = .GlobalEnv)) {
-    rng_old <- get(".Random.seed", envir = .GlobalEnv)
-  }
+  withr::local_preserve_seed()
   tstsetups <- grep("\\.brnll\\..*\\.augdat\\.", names(prjs), value = TRUE)
   for (tstsetup in tstsetups) {
     tstsetup_trad <- sub("\\.augdat\\.", ".trad_compare.", tstsetup)
@@ -214,7 +205,6 @@ test_that(paste(
     expect_equal(structure(unclass(eta_aug), ndiscrete = NULL), eta_trad,
                      info = tstsetup)
   }
-  if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
 })
 
 ## Projection -------------------------------------------------------------
@@ -454,15 +444,16 @@ test_that(paste(
     # Sometimes, there are seemingly large differences on logit scale which are
     # probably due to numerical overflow of probabilities (i.e., on probability
     # scale) towards 1 or underflow towards zero. Thus, compare on probability
-    # scale if the comparison on logit scale fails:
-    tryCatch(expect_equal(summs_sub_mu, summs_sub_mu_trad,
-                          tolerance = 1e-6, info = tstsetup),
-             error = function(e) {
-               stopifnot(is_kfold)
-               expect_equal(ilinkfun_raw(summs_sub_mu, link_nm = "logit"),
-                            ilinkfun_raw(summs_sub_mu_trad, link_nm = "logit"),
-                            tolerance = 1e-5, info = tstsetup)
-             })
+    # scale if the comparison on logit scale fails (only for K-fold CV):
+    if (is_kfold &&
+        !isTRUE(all.equal(summs_sub_mu, summs_sub_mu_trad, tolerance = 1e-6))) {
+      expect_equal(ilinkfun_raw(summs_sub_mu, link_nm = "logit"),
+                   ilinkfun_raw(summs_sub_mu_trad, link_nm = "logit"),
+                   tolerance = 1e-5, info = tstsetup)
+    } else {
+      expect_equal(summs_sub_mu, summs_sub_mu_trad, tolerance = 1e-6,
+                   info = tstsetup)
+    }
     summs_sub_lppd <- do.call(cbind, lapply(summs_sub, "[[", "lppd"))
     summs_sub_lppd_trad <- unname(
       do.call(cbind, lapply(summs_sub_trad, "[[", "lppd"))
@@ -470,14 +461,16 @@ test_that(paste(
     # Sometimes, there are seemingly large differences on log scale which are
     # probably due to numerical overflow of probabilities (i.e., on exp scale)
     # towards 1 or underflow towards zero. Thus, compare on exp scale if the
-    # comparison on log scale fails:
-    tryCatch(expect_equal(summs_sub_lppd, summs_sub_lppd_trad,
-                          tolerance = 1e-6, info = tstsetup),
-             error = function(e) {
-               stopifnot(is_kfold)
-               expect_equal(exp(summs_sub_lppd), exp(summs_sub_lppd_trad),
-                            tolerance = 1e-5, info = tstsetup)
-             })
+    # comparison on log scale fails (only for K-fold CV):
+    if (is_kfold &&
+        !isTRUE(all.equal(summs_sub_lppd, summs_sub_lppd_trad,
+                          tolerance = 1e-6))) {
+      expect_equal(exp(summs_sub_lppd), exp(summs_sub_lppd_trad),
+                   tolerance = 1e-5, info = tstsetup)
+    } else {
+      expect_equal(summs_sub_lppd, summs_sub_lppd_trad, tolerance = 1e-6,
+                   info = tstsetup)
+    }
     summs_ref <- cvvs$summaries$ref
     summs_ref$mu <- structure(unclass(summs_ref$mu), ndiscrete = NULL)
     summs_ref$mu <- summs_ref$mu[(nobsv + 1):(2 * nobsv)]
@@ -526,7 +519,7 @@ test_that(paste(
     } else {
       tol_smmry <- 1e-6
     }
-    compare_exp <- function(e) {
+    compare_exp <- function() {
       stopifnot(is_kfold)
       # Check that we have the default `stats` in this case, meaning only the
       # ELPD:
@@ -548,9 +541,16 @@ test_that(paste(
       expect_equal(smmry_pd$elpd.se, smmry_pd_trad$elpd.se, tolerance = 1e-3,
                    info = tstsetup)
     }
-    tryCatch(expect_equal(smmry_cvvs$perf_sub, smmry_cvvs_trad$perf_sub,
-                          tolerance = tol_smmry, info = tstsetup),
-             error = compare_exp)
+    # For K-fold CV, compare on exp scale if the comparison on log scale fails
+    # (see compare_exp()):
+    if (is_kfold &&
+        !isTRUE(all.equal(smmry_cvvs$perf_sub, smmry_cvvs_trad$perf_sub,
+                          tolerance = tol_smmry))) {
+      compare_exp()
+    } else {
+      expect_equal(smmry_cvvs$perf_sub, smmry_cvvs_trad$perf_sub,
+                   tolerance = tol_smmry, info = tstsetup)
+    }
   }
 })
 
@@ -606,7 +606,7 @@ test_that(paste(
     } else {
       tol_smmry <- 1e-6
     }
-    compare_exp <- function(e) {
+    compare_exp <- function() {
       stopifnot(is_kfold)
 
       smmry_pd <- smmry_cvvs$perf_sub
@@ -638,15 +638,17 @@ test_that(paste(
                    smmry_pd_trad[se_cols_trad],
                    tolerance = 1e-3, info = tstsetup)
     }
-    tryCatch(expect_equal(smmry_cvvs$perf_sub, smmry_cvvs_trad$perf_sub,
-                          tolerance = tol_smmry, info = tstsetup),
-             error = compare_exp)
+    # For K-fold CV, compare on exp scale if the comparison on log scale fails
+    # (see compare_exp()):
+    if (is_kfold &&
+        !isTRUE(all.equal(smmry_cvvs$perf_sub, smmry_cvvs_trad$perf_sub,
+                          tolerance = tol_smmry))) {
+      compare_exp()
+    } else {
+      expect_equal(smmry_cvvs$perf_sub, smmry_cvvs_trad$perf_sub,
+                   tolerance = tol_smmry, info = tstsetup)
+    }
     expect_equal(smmry_cvvs$perf_ref, smmry_cvvs_trad$perf_ref,
                  tolerance = tol_smmry, info = tstsetup)
   }
 })
-
-# Teardown ----------------------------------------------------------------
-
-# Clean up the workspace:
-rm(list = setdiff(ls(), ls_bu))
diff --git a/tests/testthat/test_cvindices.R b/tests/testthat/test_cvindices.R
index 4e49b8444..a9474cfe7 100644
--- a/tests/testthat/test_cvindices.R
+++ b/tests/testthat/test_cvindices.R
@@ -1,5 +1,3 @@
-context("cv-indices")
-
 # tests for cv_folds and cv_ids
 
 test_that("k is valid", {
diff --git a/tests/testthat/test_datafit.R b/tests/testthat/test_datafit.R
index 4f56d499f..8e0b769d1 100644
--- a/tests/testthat/test_datafit.R
+++ b/tests/testthat/test_datafit.R
@@ -1,5 +1,3 @@
-context("datafit")
-
 # Setup -------------------------------------------------------------------
 
 # Note: Since PR #351, offsets are not supported anymore for `datafit`s. Here,
@@ -194,6 +192,7 @@ if (run_vs) {
 ## Reference model --------------------------------------------------------
 
 test_that("init_refmodel(): `object` of class `datafit` works", {
+  skip_if_no_setups(datafits)
   for (tstsetup in names(datafits)) {
     tstsetup_fit <- args_datafit[[tstsetup]]$tstsetup_fit
     with_spclformul_crr <- grepl("\\.spclformul", tstsetup)
@@ -222,6 +221,7 @@ test_that("init_refmodel(): `object` of class `datafit` works", {
 })
 
 test_that("predict.refmodel(): `object` of class `datafit` fails", {
+  skip_if_no_setups(datafits)
   for (tstsetup in names(datafits)) {
     expect_error(
       predict(datafits[[tstsetup]], newdata = dat),
@@ -430,12 +430,12 @@ test_that(paste(
               info_str = paste(tstsetup, "with_args", sep = "__"))
     if (run_snaps) {
       if (testthat_ed_max2) local_edition(3)
-      width_orig <- options(width = 145)
-      expect_snapshot({
-        print(tstsetup)
-        print(rlang::hash(pl_with_args))
+      withr::with_options(list(width = 145), {
+        expect_snapshot({
+          print(tstsetup)
+          print(rlang::hash(pl_with_args))
+        })
       })
-      options(width_orig)
       if (testthat_ed_max2) local_edition(2)
     }
   }
@@ -471,12 +471,12 @@ test_that(paste(
               info_str = paste(tstsetup, "with_args", sep = "__"))
     if (run_snaps) {
       if (testthat_ed_max2) local_edition(3)
-      width_orig <- options(width = 145)
-      expect_snapshot({
-        print(tstsetup)
-        print(rlang::hash(pp_with_args))
+      withr::with_options(list(width = 145), {
+        expect_snapshot({
+          print(tstsetup)
+          print(rlang::hash(pp_with_args))
+        })
       })
-      options(width_orig)
       if (testthat_ed_max2) local_edition(2)
     }
   }
@@ -485,6 +485,7 @@ test_that(paste(
 ## summary.vsel() ---------------------------------------------------------
 
 test_that("summary.vsel(): `object` of class `datafit` fails", {
+  skip_if_no_setups(datafits)
   for (tstsetup in names(datafits)) {
     expect_error(
       summary.vsel(datafits[[tstsetup]]),
@@ -537,12 +538,12 @@ test_that(paste(
     )
     if (run_snaps) {
       if (testthat_ed_max2) local_edition(3)
-      width_orig <- options(width = 145)
-      expect_snapshot({
-        print(tstsetup)
-        print(smmry, digits = 6)
+      withr::with_options(list(width = 145), {
+        expect_snapshot({
+          print(tstsetup)
+          print(smmry, digits = 6)
+        })
       })
-      options(width_orig)
       if (testthat_ed_max2) local_edition(2)
     }
   }
@@ -580,12 +581,12 @@ test_that(paste(
     )
     if (run_snaps) {
       if (testthat_ed_max2) local_edition(3)
-      width_orig <- options(width = 145)
-      expect_snapshot({
-        print(tstsetup)
-        print(smmry, digits = 6)
+      withr::with_options(list(width = 145), {
+        expect_snapshot({
+          print(tstsetup)
+          print(smmry, digits = 6)
+        })
       })
-      options(width_orig)
       if (testthat_ed_max2) local_edition(2)
     }
   }
@@ -606,9 +607,7 @@ test_that(paste(
   # This test sometimes behaves inpredictably when run in `R CMD check`, so skip
   # it on CRAN:
   skip_on_cran()
-  if (exists(".Random.seed", envir = .GlobalEnv)) {
-    rng_old <- get(".Random.seed", envir = .GlobalEnv)
-  }
+  withr::local_preserve_seed()
   suppressWarnings(RNGversion("3.5.0"))
   set.seed(1235)
   n <- 100
@@ -754,5 +753,4 @@ test_that(paste(
     expect_identical(ind, predictor_ranking_lasso, info = as.character(i))
   }
   RNGversion(getRversion())
-  if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
 })
diff --git a/tests/testthat/test_div_minimizer.R b/tests/testthat/test_div_minimizer.R
index ad1bb076f..98a218980 100644
--- a/tests/testthat/test_div_minimizer.R
+++ b/tests/testthat/test_div_minimizer.R
@@ -1,9 +1,6 @@
-context("div_minimizer")
-
-# Setup -------------------------------------------------------------------
-
-# Needed to clean up the workspace afterwards (i.e, after this test file):
-ls_bu <- ls()
+# The divmin() and the divmin_augdat() test store their `brnll` results here for
+# the comparison test at the end of this file:
+outdmin_brnll <- new.env()
 
 # divmin() ----------------------------------------------------------------
 
@@ -88,7 +85,7 @@ test_that("divmin() works", {
                    with_offs = grepl("\\.with_offs", tstsetup),
                    info_str = tstsetup)
   }
-  assign("outdmin_brnll_trad", outdmin_brnll_tmp, envir = .GlobalEnv)
+  outdmin_brnll$trad <- outdmin_brnll_tmp
 })
 
 # divmin_augdat() ---------------------------------------------------------
@@ -214,7 +211,7 @@ test_that("divmin_augdat() works", {
                    check_y_from_resp = FALSE,
                    info_str = tstsetup)
   }
-  assign("outdmin_brnll_augdat", outdmin_brnll_tmp, envir = .GlobalEnv)
+  outdmin_brnll$augdat <- outdmin_brnll_tmp
 })
 
 # Comparison of divmin() and divmin_augdat() ------------------------------
@@ -224,10 +221,13 @@ test_that(paste(
   "`brnll` family)"
 ), {
   skip_on_cran()
-  for (tstsetup in names(outdmin_brnll_augdat)) {
+  skip_if(is.null(outdmin_brnll$trad) || is.null(outdmin_brnll$augdat),
+          "the divmin() or the divmin_augdat() test did not finish")
+  skip_if_no_setups(outdmin_brnll$augdat)
+  for (tstsetup in names(outdmin_brnll$augdat)) {
     args_ref_i <- args_ref[[tstsetup]]
-    outdmin_augdat <- outdmin_brnll_augdat[[tstsetup]]
-    outdmin_trad <- outdmin_brnll_trad[[args_ref_i$tstsetup_fit]]
+    outdmin_augdat <- outdmin_brnll$augdat[[tstsetup]]
+    outdmin_trad <- outdmin_brnll$trad[[args_ref_i$tstsetup_fit]]
     expect_length(outdmin_augdat, 1)
     expect_length(outdmin_trad, 1)
     tol_coefs <- ifelse(args_ref_i$mod_nm == "glmm", 1e-4, 1e-5)
@@ -245,8 +245,3 @@ test_that(paste(
                  info = tstsetup)
   }
 })
-
-# Teardown ----------------------------------------------------------------
-
-# Clean up the workspace:
-rm(list = setdiff(ls(), ls_bu))
diff --git a/tests/testthat/test_extend_family.R b/tests/testthat/test_extend_family.R
index 500d5e2b6..047f3a770 100644
--- a/tests/testthat/test_extend_family.R
+++ b/tests/testthat/test_extend_family.R
@@ -1,5 +1,3 @@
-context("extend_family()")
-
 test_that("extend_family() works for the traditional projection", {
   fam_nms_tst <- c("gauss", "binom", "poiss")
   ### TODO: Also test other families extensible by projpred:
diff --git a/tests/testthat/test_formula.R b/tests/testthat/test_formula.R
index 5cdec246b..e5d1a8edb 100644
--- a/tests/testthat/test_formula.R
+++ b/tests/testthat/test_formula.R
@@ -1,5 +1,3 @@
-context("formula")
-
 test_that(paste(
   "check that we recover the correct terms for a simple linear",
   "model without interactions or group terms"
diff --git a/tests/testthat/test_glm_elnet.R b/tests/testthat/test_glm_elnet.R
index f720fcec9..3aed36ff3 100644
--- a/tests/testthat/test_glm_elnet.R
+++ b/tests/testthat/test_glm_elnet.R
@@ -1,14 +1,6 @@
-context("elnet")
-
 # tests for glm_elnet
 
-# Needed to clean up the workspace afterwards (i.e, after this test file):
-ls_bu <- ls()
-
-if (exists(".Random.seed", envir = .GlobalEnv)) {
-  rng_old <- get(".Random.seed", envir = .GlobalEnv)
-}
-set.seed(1235)
+withr::local_seed(1235)
 n <- 40
 nterms <- 10
 nterms_fit <- nterms - 5
@@ -221,7 +213,3 @@ test_that("glm_elnet with alpha=0 and glm_ridge give the same result.", {
     }
   }
 })
-
-if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
-# Clean up the workspace:
-rm(list = setdiff(ls(), ls_bu))
diff --git a/tests/testthat/test_glm_ridge.R b/tests/testthat/test_glm_ridge.R
index c7ccca48c..d82aee67d 100644
--- a/tests/testthat/test_glm_ridge.R
+++ b/tests/testthat/test_glm_ridge.R
@@ -1,5 +1,3 @@
-context("ridge")
-
 # tests for ridge regression, currently untested combinations
 # - gaussian with inverse-link
 # - binomial with log or cloglog-link
@@ -7,13 +5,7 @@ context("ridge")
 # - Gamma with inverse or id-link
 # - everything except gaussian with id-link for ridge penalty
 
-# Needed to clean up the workspace afterwards (i.e, after this test file):
-ls_bu <- ls()
-
-if (exists(".Random.seed", envir = .GlobalEnv)) {
-  rng_old <- get(".Random.seed", envir = .GlobalEnv)
-}
-set.seed(1235)
+withr::local_seed(1235)
 n <- 40
 nterms <- 10
 nterms_fit <- nterms - 5
@@ -270,7 +262,3 @@ test_that("glm_ridge: poisson, log-link, no intercept, lambda = 0", {
 #   expect_equal(unname(coef(glmfit)), c(ridgefit$beta0, ridgefit$beta),
 #                tolerance = tol)
 # })
-
-if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
-# Clean up the workspace:
-rm(list = setdiff(ls(), ls_bu))
diff --git a/tests/testthat/test_helpers_tst.R b/tests/testthat/test_helpers_tst.R
new file mode 100644
index 000000000..29b529fc6
--- /dev/null
+++ b/tests/testthat/test_helpers_tst.R
@@ -0,0 +1,6 @@
+# Helpers of the tests (tests/testthat/helpers/) -----------------------------
+
+test_that("skip_if_no_setups() skips only for an empty list", {
+  expect_condition(skip_if_no_setups(list()), class = "skip")
+  expect_no_condition(skip_if_no_setups(list(a = 1)), class = "skip")
+})
diff --git a/tests/testthat/test_latent.R b/tests/testthat/test_latent.R
index 242ded393..9c9dc31dd 100644
--- a/tests/testthat/test_latent.R
+++ b/tests/testthat/test_latent.R
@@ -1,13 +1,6 @@
-context("latent projection")
-
 # Setup -------------------------------------------------------------------
 
-# Needed to clean up the workspace afterwards (i.e, after this test file):
-ls_bu <- ls()
-
-if (exists(".Random.seed", envir = .GlobalEnv)) {
-  rng_old <- get(".Random.seed", envir = .GlobalEnv)
-}
+withr::local_preserve_seed()
 
 # latent.R ----------------------------------------------------------------
 
@@ -101,9 +94,7 @@ test_that(paste(
   "traditional projection"
 ), {
   skip_if_not(run_prj)
-  if (exists(".Random.seed", envir = .GlobalEnv)) {
-    rng_old <- get(".Random.seed", envir = .GlobalEnv)
-  }
+  withr::local_preserve_seed()
   tstsetups <- grep("\\.brnll\\..*\\.latent\\.", names(prjs), value = TRUE)
   for (tstsetup in tstsetups) {
     tstsetup_trad <- sub("\\.latent\\.", ".trad_compare.", tstsetup)
@@ -154,7 +145,6 @@ test_that(paste(
     expect_equal(mu_offs_lat_oscale_cl, pref_trad$mu_offs, tolerance = 1e-10,
                  info = tstsetup)
   }
-  if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
 })
 
 ## Gaussian family --------------------------------------------------------
@@ -232,12 +222,6 @@ test_that(paste(
   }
 })
 
-# Teardown ----------------------------------------------------------------
-
-if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
-# Clean up the workspace:
-rm(list = setdiff(ls(), ls_bu))
-
 # Censored observations ---------------------------------------------------
 
 test_that(paste0(
diff --git a/tests/testthat/test_methods_vsel.R b/tests/testthat/test_methods_vsel.R
index a7f12f14c..632f0cd6e 100644
--- a/tests/testthat/test_methods_vsel.R
+++ b/tests/testthat/test_methods_vsel.R
@@ -1,7 +1,5 @@
 # Common tests ------------------------------------------------------------
 
-context("summary(), plot(), suggest_size()")
-
 test_that("invalid `object` fails", {
   objs_invalid <- nlist(
     NULL,
@@ -61,8 +59,6 @@ test_that("invalid `stats` fails", {
 
 # summary() ---------------------------------------------------------------
 
-context("summary()")
-
 test_that("invalid `baseline` fails", {
   skip_if_not(run_vs)
   for (tstsetup in head(names(vss), 1)) {
@@ -124,8 +120,6 @@ test_that(paste(
 
 # performances() ----------------------------------------------------------
 
-context("performances()")
-
 test_that("`object` of class `vsel` (created by varsel()) works", {
   skip_if_not(run_vs)
   for (tstsetup in names(perfs_vs)) {
@@ -185,8 +179,6 @@ test_that("performances.vsel() is a shortcut", {
 
 # print() -----------------------------------------------------------------
 
-context("print()")
-
 test_that("`x` of class `vselsummary` (based on varsel()) works", {
   skip_if_not(run_vs)
   for (tstsetup in names(smmrys_vs)) {
@@ -198,12 +190,12 @@ test_that("`x` of class `vselsummary` (based on varsel()) works", {
     expect_identical(print_obj, smmrys_vs[[tstsetup]], info = tstsetup)
     if (run_snaps) {
       if (testthat_ed_max2) local_edition(3)
-      width_orig <- options(width = 145)
-      expect_snapshot({
-        print(tstsetup)
-        print(smmrys_vs[[tstsetup]], digits = 6)
+      withr::with_options(list(width = 145), {
+        expect_snapshot({
+          print(tstsetup)
+          print(smmrys_vs[[tstsetup]], digits = 6)
+        })
       })
-      options(width_orig)
       if (testthat_ed_max2) local_edition(2)
     }
   }
@@ -221,12 +213,12 @@ test_that("`x` of class `vselsummary` (based on cv_varsel())  works", {
     expect_identical(print_obj, smmrys_cvvs[[tstsetup]], info = tstsetup)
     if (run_snaps) {
       if (testthat_ed_max2) local_edition(3)
-      width_orig <- options(width = 145)
-      expect_snapshot({
-        print(tstsetup)
-        print(smmrys_cvvs[[tstsetup]], digits = 6)
+      withr::with_options(list(width = 145), {
+        expect_snapshot({
+          print(tstsetup)
+          print(smmrys_cvvs[[tstsetup]], digits = 6)
+        })
       })
-      options(width_orig)
       if (testthat_ed_max2) local_edition(2)
     }
   }
@@ -284,8 +276,6 @@ test_that(paste(
 
 # plot() ------------------------------------------------------------------
 
-context("plot()")
-
 args_common_for_rk_NA <- c(
   "tstsetup_vsel", "deltas", "text_angle", "nterms_max", "ranking_nterms_max",
   "size_position", "show_cv_proportions"
@@ -538,8 +528,6 @@ test_that(paste(
 
 # suggest_size() ----------------------------------------------------------
 
-context("suggest_size()")
-
 test_that("`stat` of invalid length fails", {
   stopifnot(length(stats_common) > 1)
   skip_if_not(run_vs)
@@ -609,8 +597,6 @@ test_that("returns nterms_max when no criterion is met", {
 
 # ranking() ---------------------------------------------------------------
 
-context("ranking()")
-
 test_that("`object` of class `vsel` (created by varsel()) works", {
   skip_if_not(run_vs)
   for (tstsetup in names(rks_vs)) {
@@ -643,8 +629,6 @@ test_that("`object` of class `vsel` (created by cv_varsel()) works", {
 
 # cv_proportions() --------------------------------------------------------
 
-context("cv_proportions()")
-
 test_that("`object` of class `ranking` (based on varsel() output) fails", {
   skip_if_not(run_vs)
   expect_length(prs_vs, 0)
@@ -691,9 +675,6 @@ test_that("cv_proportions.vsel() is a shortcut", {
   }
 })
 
-# Needed to clean up the workspace afterwards:
-ls_bu <- ls()
-
 ntrms <- 9L
 rk_fdata <- rev(head(letters, ntrms))
 rk_fwise <- do.call(rbind, lapply(seq_len(ntrms), function(idx_trm) {
@@ -740,13 +721,8 @@ test_that("ranking proportions are computed correctly", {
   expect_true(has_1_colwise, info = "cumulate = TRUE")
 })
 
-# Clean up the workspace:
-rm(list = setdiff(ls(), ls_bu))
-
 # plot.cv_proportions() ---------------------------------------------------
 
-context("plot.cv_proportions()")
-
 test_that("`x` of class `cv_proportions` works", {
   skip_if_not(run_cvvs)
   for (tstsetup in names(plotprs)) {
@@ -773,9 +749,6 @@ test_that("plot.ranking() is a shortcut", {
   }
 })
 
-# Needed to clean up the workspace afterwards:
-ls_bu <- ls()
-
 ntrms <- 20L
 pr_dummy <- matrix(seq(0, 1, length.out = ntrms^2), nrow = ntrms, ncol = ntrms,
                    dimnames = list("size" = as.character(seq_len(ntrms)),
@@ -818,6 +791,3 @@ test_that("the ggplot can be modified", {
     vdiffr::expect_doppelganger("plotprc_dummy_mod", plotprc_dummy_mod)
   }
 })
-
-# Clean up the workspace:
-rm(list = setdiff(ls(), ls_bu))
diff --git a/tests/testthat/test_misc.R b/tests/testthat/test_misc.R
index 942d517aa..1d46cd3e3 100644
--- a/tests/testthat/test_misc.R
+++ b/tests/testthat/test_misc.R
@@ -1,5 +1,3 @@
-context("miscellaneous")
-
 # get_refdist() ----------------------------------------------------------
 
 ## ndraws and nclusters ---------------------------------------------------
@@ -8,10 +6,8 @@ test_that(paste(
   "get_refdist(): `ndraws = NULL` and `nclusters = NULL` leads to",
   "`ndraws = S` (and `nclusters = NULL`)"
 ), {
-  if (exists(".Random.seed", envir = .GlobalEnv)) {
-    rng_old <- get(".Random.seed", envir = .GlobalEnv)
-  }
-  set.seed(seed2_tst)
+  skip_if_no_setups(refmods)
+  withr::local_seed(seed2_tst)
   for (tstsetup in names(refmods)) {
     refdist <- get_refdist(refmods[[tstsetup]])
     if (args_ref[[tstsetup]]$prj_nm == "augdat") {
@@ -31,17 +27,14 @@ test_that(paste(
       info_str = tstsetup
     )
   }
-  if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
 })
 
 test_that(paste(
   "`ndraws` and/or `nclusters` too big causes them to be cut off at the number",
   "of posterior draws in the reference model"
 ), {
-  if (exists(".Random.seed", envir = .GlobalEnv)) {
-    rng_old <- get(".Random.seed", envir = .GlobalEnv)
-  }
-  set.seed(seed2_tst)
+  skip_if_no_setups(refmods)
+  withr::local_seed(seed2_tst)
   for (tstsetup in names(refmods)) {
     for (ndraws_crr in list(nrefdraws + 1L)) {
       for (nclusters_crr in list(NULL, nrefdraws + 1L)) {
@@ -64,7 +57,6 @@ test_that(paste(
       }
     }
   }
-  if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
 })
 
 # force_search_terms() ----------------------------------------------------
@@ -193,10 +185,7 @@ test_that(paste(
 # <https://github.com/stan-dev/projpred/pull/496#issuecomment-2481313265>), then
 # changed function name and added seed reset:
 test_that(".srs_diff_est_w() works as expected", {
-  if (exists(".Random.seed", envir = .GlobalEnv)) {
-    rng_old <- get(".Random.seed", envir = .GlobalEnv)
-  }
-  set.seed(1234)
+  withr::local_seed(1234)
   N <- 1000
   y_true <- 1:N
   sigma_hat_true <- sqrt(N * sum((y_true - mean(y_true))^2) / length(y_true))
@@ -226,15 +215,11 @@ test_that(".srs_diff_est_w() works as expected", {
   expect_equal(res$y_hat, 500500, tol = 0.0001)
   expect_equal(res$v_y_hat, 0, tol = 0.0001)
   expect_equal(sqrt(res$hat_v_y), sigma_hat_true, tol = 0.1)
-  if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
 })
 
 # Adapted from the test above:
 test_that(".srs_diff_est_w() works as expected with `wobs` not full of ones", {
-  if (exists(".Random.seed", envir = .GlobalEnv)) {
-    rng_old <- get(".Random.seed", envir = .GlobalEnv)
-  }
-  set.seed(seed3_tst)
+  withr::local_seed(seed3_tst)
   N <- 1e3
   y_true <- seq_len(N)
   wobs <- sample(2:7, size = N, replace = TRUE)
@@ -277,13 +262,12 @@ test_that(".srs_diff_est_w() works as expected with `wobs` not full of ones", {
   expect_equal(res$y_hat, 500845, tol = 1e-4)
   expect_equal(res$v_y_hat, 0, tol = 1e-4)
   expect_equal(sqrt(res$hat_v_y), sigma_hat_true, tol = 1e-1)
-  if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
 })
 
 test_that(".srs_diff_est_w() propagates input `NA`s to its output", {
   nloo_tst <- nobsv %/% 5L
   loo_inds_tst <- ceiling(seq(1L, nobsv, length.out = nloo_tst))
-  skip_if_not(length(unique(loo_inds_tst)) == nloo_tst)
+  stopifnot(length(unique(loo_inds_tst)) == nloo_tst)
   srs_diff_res_NAs <- list(
     m = nloo_tst,
     N = nobsv,
diff --git a/tests/testthat/test_parallel.R b/tests/testthat/test_parallel.R
index 5de734c51..e3ff83412 100644
--- a/tests/testthat/test_parallel.R
+++ b/tests/testthat/test_parallel.R
@@ -1,15 +1,17 @@
-context("parallel")
-
 # Setup -------------------------------------------------------------------
 
+# withr resets the options and the parallel backend at the end of this file,
+# also if the file stops with an error.
 if (run_prll) {
-  trigger_default <- options(projpred.parallel_proj_trigger = 0L)
+  withr::local_options(list(projpred.parallel_proj_trigger = 0L))
 
   if (dopar_backend == "doParallel") {
     doParallel::registerDoParallel(ncores)
+    withr::defer(doParallel::stopImplicitCluster())
   } else if (dopar_backend == "doFuture") {
     doFuture::registerDoFuture()
-    export_default <- options(doFuture.foreach.export = ".export")
+    withr::local_options(list(doFuture.foreach.export = ".export"))
+    withr::defer(future::plan(future::sequential))
     if (future_plan == "multicore") {
       future::plan(future::multicore, workers = ncores)
     } else if (future_plan == "multisession") {
@@ -79,20 +81,3 @@ test_that("cv_varsel() in parallel gives the same results as sequentially", {
     expect_equal(cvvs_repr, cvvss[[tstsetup]], info = tstsetup)
   }
 })
-
-# Teardown ----------------------------------------------------------------
-
-if (run_prll) {
-  if (dopar_backend == "doParallel") {
-    doParallel::stopImplicitCluster()
-  } else if (dopar_backend == "doFuture") {
-    future::plan(future::sequential)
-    options(export_default)
-    rm(export_default)
-  } else {
-    stop("Unrecognized `dopar_backend`.")
-  }
-
-  options(trigger_default)
-  rm(trigger_default)
-}
diff --git a/tests/testthat/test_proj_epred.R b/tests/testthat/test_proj_epred.R
index eb77be9ab..d03f6c4ba 100644
--- a/tests/testthat/test_proj_epred.R
+++ b/tests/testthat/test_proj_epred.R
@@ -1,7 +1,5 @@
 # proj_epred() -------------------------------------------------------------
 
-context("proj_epred()")
-
 ## object ------------------------------------------------------------------
 
 test_that("pe: `object` of class `projection` works", {
diff --git a/tests/testthat/test_proj_pred.R b/tests/testthat/test_proj_pred.R
index 2edd4eb6a..5fde3f6f0 100644
--- a/tests/testthat/test_proj_pred.R
+++ b/tests/testthat/test_proj_pred.R
@@ -1,7 +1,5 @@
 # proj_linpred() ----------------------------------------------------------
 
-context("proj_linpred()")
-
 ## object -----------------------------------------------------------------
 
 test_that("pl: `object` of class `projection` works", {
@@ -27,12 +25,12 @@ test_that("pl: `object` of class `projection` works", {
               info_str = tstsetup)
     if (run_snaps) {
       if (testthat_ed_max2) local_edition(3)
-      width_orig <- options(width = 145)
-      expect_snapshot({
-        print(tstsetup)
-        print(rlang::hash(pls[[tstsetup]]))
+      withr::with_options(list(width = 145), {
+        expect_snapshot({
+          print(tstsetup)
+          print(rlang::hash(pls[[tstsetup]]))
+        })
       })
-      options(width_orig)
       if (testthat_ed_max2) local_edition(2)
     }
   }
@@ -76,12 +74,12 @@ test_that(paste(
               info_str = tstsetup)
     if (run_snaps) {
       if (testthat_ed_max2) local_edition(3)
-      width_orig <- options(width = 145)
-      expect_snapshot({
-        print(tstsetup)
-        print(rlang::hash(pls_vs[[tstsetup]]))
+      withr::with_options(list(width = 145), {
+        expect_snapshot({
+          print(tstsetup)
+          print(rlang::hash(pls_vs[[tstsetup]]))
+        })
       })
-      options(width_orig)
       if (testthat_ed_max2) local_edition(2)
     }
   }
@@ -125,12 +123,12 @@ test_that(paste(
               info_str = tstsetup)
     if (run_snaps) {
       if (testthat_ed_max2) local_edition(3)
-      width_orig <- options(width = 145)
-      expect_snapshot({
-        print(tstsetup)
-        print(rlang::hash(pls_cvvs[[tstsetup]]))
+      withr::with_options(list(width = 145), {
+        expect_snapshot({
+          print(tstsetup)
+          print(rlang::hash(pls_cvvs[[tstsetup]]))
+        })
       })
-      options(width_orig)
       if (testthat_ed_max2) local_edition(2)
     }
   }
@@ -797,13 +795,11 @@ test_that("`transform` works", {
     if (args_prj[[tstsetup]]$prj_nm != "latent") {
       pred_false2true <- prjs[[tstsetup]]$refmodel$family$linkinv(pred_false)
     } else {
-      if (exists(".Random.seed", envir = .GlobalEnv)) {
-        rng_old <- get(".Random.seed", envir = .GlobalEnv)
-      }
-      set.seed(args_prj[[tstsetup]]$seed)
-      clust_ref <- get_refdist(prjs[[tstsetup]]$refmodel,
-                               ndraws = args_prj[[tstsetup]]$ndraws,
-                               nclusters = args_prj[[tstsetup]]$nclusters)
+      clust_ref <- withr::with_seed(args_prj[[tstsetup]]$seed, {
+        get_refdist(prjs[[tstsetup]]$refmodel,
+                    ndraws = args_prj[[tstsetup]]$ndraws,
+                    nclusters = args_prj[[tstsetup]]$nclusters)
+      })
       pred_false2true <- structure(
         prjs[[tstsetup]]$refmodel$family$latent_ilink(
           t(pred_false), cl_ref = clust_ref$cl
@@ -812,7 +808,6 @@ test_that("`transform` works", {
       )
     }
     expect_equal(pred_false2true, pred_true, info = tstsetup)
-    if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
   }
 })
 
@@ -1067,8 +1062,6 @@ test_that(paste(
 
 # proj_predict() ----------------------------------------------------------
 
-context("proj_predict()")
-
 ## seed -------------------------------------------------------------------
 
 test_that("`.seed` works (and restores the RNG state afterwards)", {
@@ -1115,12 +1108,12 @@ test_that("pp: `object` of class `projection` works", {
               info_str = tstsetup)
     if (run_snaps) {
       if (testthat_ed_max2) local_edition(3)
-      width_orig <- options(width = 145)
-      expect_snapshot({
-        print(tstsetup)
-        print(rlang::hash(pps[[tstsetup]]))
+      withr::with_options(list(width = 145), {
+        expect_snapshot({
+          print(tstsetup)
+          print(rlang::hash(pps[[tstsetup]]))
+        })
       })
-      options(width_orig)
       if (testthat_ed_max2) local_edition(2)
     }
   }
@@ -1149,12 +1142,12 @@ test_that(paste(
               info_str = tstsetup)
     if (run_snaps) {
       if (testthat_ed_max2) local_edition(3)
-      width_orig <- options(width = 145)
-      expect_snapshot({
-        print(tstsetup)
-        print(rlang::hash(pps_vs[[tstsetup]]))
+      withr::with_options(list(width = 145), {
+        expect_snapshot({
+          print(tstsetup)
+          print(rlang::hash(pps_vs[[tstsetup]]))
+        })
       })
-      options(width_orig)
       if (testthat_ed_max2) local_edition(2)
     }
   }
@@ -1183,12 +1176,12 @@ test_that(paste(
               info_str = tstsetup)
     if (run_snaps) {
       if (testthat_ed_max2) local_edition(3)
-      width_orig <- options(width = 145)
-      expect_snapshot({
-        print(tstsetup)
-        print(rlang::hash(pps_cvvs[[tstsetup]]))
+      withr::with_options(list(width = 145), {
+        expect_snapshot({
+          print(tstsetup)
+          print(rlang::hash(pps_cvvs[[tstsetup]]))
+        })
       })
-      options(width_orig)
       if (testthat_ed_max2) local_edition(2)
     }
   }
diff --git a/tests/testthat/test_proj_predfun.R b/tests/testthat/test_proj_predfun.R
index caa083a4b..8bf4df8a0 100644
--- a/tests/testthat/test_proj_predfun.R
+++ b/tests/testthat/test_proj_predfun.R
@@ -1,5 +1,3 @@
-context("proj_predfun")
-
 test_that("repair_re() works for GLMMs", {
   # Setup -------------------------------------------------------------------
 
@@ -9,10 +7,7 @@ test_that("repair_re() works for GLMMs", {
   # stan-dev/rstanarm#610, so don't run this test on CRAN:
   skip_on_cran()
 
-  if (exists(".Random.seed", envir = .GlobalEnv)) {
-    rng_old <- get(".Random.seed", envir = .GlobalEnv)
-  }
-  set.seed(seed2_tst)
+  withr::local_seed(seed2_tst)
 
   # Data --------------------------------------------------------------------
 
@@ -169,9 +164,6 @@ test_that("repair_re() works for GLMMs", {
     })
   expect_equal(lpreds, lpreds_man, tolerance = 1e-15)
 
-  # Teardown ----------------------------------------------------------------
-
-  if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
 })
 
 test_that("repair_re() works for multilevel cumulative() models", {
@@ -179,10 +171,7 @@ test_that("repair_re() works for multilevel cumulative() models", {
 
   # Setup -------------------------------------------------------------------
 
-  if (exists(".Random.seed", envir = .GlobalEnv)) {
-    rng_old <- get(".Random.seed", envir = .GlobalEnv)
-  }
-  set.seed(seed2_tst)
+  withr::local_seed(seed2_tst)
 
   # Data --------------------------------------------------------------------
 
@@ -331,23 +320,14 @@ test_that("repair_re() works for multilevel cumulative() models", {
        }))
   expect_equal(unname(lpreds), lpreds_man, tolerance = 1e-14)
 
-  # Teardown ----------------------------------------------------------------
-
-  if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
 })
 
 # Multilevel categorical models -------------------------------------------
 
 ## Setup ------------------------------------------------------------------
 
-# Needed to clean up the workspace afterwards (i.e, after this test file):
-ls_bu <- ls()
-
 # Set seed:
-if (exists(".Random.seed", envir = .GlobalEnv)) {
-  rng_old <- get(".Random.seed", envir = .GlobalEnv)
-}
-set.seed(seed2_tst)
+withr::local_seed(seed2_tst)
 
 ## Data -------------------------------------------------------------------
 
@@ -910,10 +890,3 @@ test_that(paste(
   expect_equal(as.numeric(lpreds), as.numeric(lpreds_man),
                tolerance = 1e-14)
 })
-
-# Teardown ----------------------------------------------------------------
-
-if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
-
-# Clean up the workspace:
-rm(list = setdiff(ls(), ls_bu))
diff --git a/tests/testthat/test_project.R b/tests/testthat/test_project.R
index 7c2162606..7ce8ad3d3 100644
--- a/tests/testthat/test_project.R
+++ b/tests/testthat/test_project.R
@@ -1,5 +1,3 @@
-context("project()")
-
 # object, predictor_terms, ndraws / nclusters, nterms ---------------------
 
 test_that(paste(
diff --git a/tests/testthat/test_refmodel.R b/tests/testthat/test_refmodel.R
index 6a87083b0..810bc098c 100644
--- a/tests/testthat/test_refmodel.R
+++ b/tests/testthat/test_refmodel.R
@@ -1,8 +1,7 @@
 # get_refmodel() ----------------------------------------------------------
 
-context("get_refmodel()")
-
 test_that("`object` of class `stanreg` or `brmsfit` works", {
+  skip_if_no_setups(refmods)
   for (tstsetup in names(refmods)) {
     tstsetup_fit <- args_ref[[tstsetup]]$tstsetup_fit
     with_spclformul_crr <- grepl("\\.spclformul", tstsetup)
@@ -199,6 +198,7 @@ test_that("extra arguments in s() or t2() terms fail", {
 })
 
 test_that("get_refmodel() is idempotent", {
+  skip_if_no_setups(refmods)
   for (tstsetup in names(refmods)) {
     expect_identical(get_refmodel(refmods[[tstsetup]]),
                      refmods[[tstsetup]],
@@ -208,8 +208,6 @@ test_that("get_refmodel() is idempotent", {
 
 # predict.refmodel() ------------------------------------------------------
 
-context("predict.refmodel()")
-
 test_that("invalid `type` fails", {
   skip_if_not(length(fits) > 0)
   expect_error(predict(refmods[[1]], dat, type = "zzz"),
@@ -226,6 +224,7 @@ test_that(paste(
   "`object` of class `refmodel` and arguments `newdata`, `ynew`, and `type`",
   "work"
 ), {
+  skip_if_no_setups(refmods)
   for (tstsetup in names(refmods)) {
     pkg_crr <- args_ref[[tstsetup]]$pkg_nm
     mod_crr <- args_ref[[tstsetup]]$mod_nm
@@ -355,15 +354,15 @@ test_that(paste(
     # Snapshots:
     if (run_snaps) {
       if (testthat_ed_max2) local_edition(3)
-      width_orig <- options(width = 145)
-      expect_snapshot({
-        print(tstsetup)
-        print(rlang::hash(predref_resp))
-        print(rlang::hash(predref_link))
-        print(rlang::hash(predref_ynew_resp))
-        print(rlang::hash(predref_ynew_link))
+      withr::with_options(list(width = 145), {
+        expect_snapshot({
+          print(tstsetup)
+          print(rlang::hash(predref_resp))
+          print(rlang::hash(predref_link))
+          print(rlang::hash(predref_ynew_resp))
+          print(rlang::hash(predref_ynew_link))
+        })
       })
-      options(width_orig)
       if (testthat_ed_max2) local_edition(2)
     }
   }
diff --git a/tests/testthat/test_varsel.R b/tests/testthat/test_varsel.R
index 9f5dda555..e1f8808dc 100644
--- a/tests/testthat/test_varsel.R
+++ b/tests/testthat/test_varsel.R
@@ -1,7 +1,5 @@
 # varsel() ----------------------------------------------------------------
 
-context("varsel()")
-
 test_that(paste(
   "`object` of class `refmodel` and arguments `method`, `nterms_max`,",
   "`nclusters`, and `nclusters_pred` work"
@@ -34,6 +32,7 @@ test_that("invalid `object` fails", {
 })
 
 test_that("invalid `method` fails", {
+  skip_if_no_setups(refmods)
   for (tstsetup in names(refmods)) {
     expect_error(varsel(refmods[[tstsetup]], method = "k-fold"),
                  "^Unexpected value for argument `method`\\.$",
@@ -181,9 +180,7 @@ test_that(paste(
   "posterior_epred() and log_lik()"
 ), {
   skip_if_not(run_vs)
-  if (exists(".Random.seed", envir = .GlobalEnv)) {
-    rng_old <- get(".Random.seed", envir = .GlobalEnv)
-  }
+  withr::local_preserve_seed()
   tstsetups <- names(vss)
   ### TODO (GAMMs): Currently, the following test setups (can) lead to the error
   ### ```
@@ -495,7 +492,6 @@ test_that(paste(
     expect_equal(vs_indep$summaries$ref$lppd, lppd_ref_ch2,
                  tolerance = 1e2 * .Machine$double.eps, info = tstsetup)
   }
-  if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
 })
 
 ## refit_prj --------------------------------------------------------------
@@ -553,7 +549,7 @@ if (run_more) {
     "the warning for full Gaussian multilevel models is thrown correctly"
   ), {
     skip_if_not(run_vs)
-    warn_instable_orig <- options(projpred.warn_instable_projections = NULL)
+    withr::local_options(list(projpred.warn_instable_projections = NULL))
     tstsetups <- names(vss)
     pick_these <- sapply(tstsetups, function(tstsetup) {
       refmod_i <- refmods[[args_vs[[tstsetup]]$tstsetup_ref]]
@@ -584,7 +580,6 @@ if (run_more) {
         info = tstsetup
       )
     }
-    options(warn_instable_orig)
   })
 }
 
@@ -594,7 +589,7 @@ if (run_more) {
   test_that("the message when cutting off the search is thrown correctly", {
     skip_if_not(run_vs)
     skip_if_not_installed("rstanarm")
-    mssg_cut_search_orig <- options(projpred.mssg_cut_search = NULL)
+    withr::local_options(list(projpred.mssg_cut_search = NULL))
     dat_gauss <- data.frame(y = df_gaussian$y, df_gaussian$x)
     stopifnot(sum(grepl("^X", names(dat_gauss))) > 19)
     fit_exceed <- suppressWarnings(rstanarm::stan_glm(
@@ -628,7 +623,6 @@ if (run_more) {
         info = paste("fit_not_exceed, nterms_max =", nterms_max_i %||% "NULL")
       )
     }
-    options(mssg_cut_search_orig)
   })
 }
 
@@ -959,7 +953,7 @@ test_that("for L1 search, `penalty` has an expected effect", {
 test_that("L1 search handles three-way (second-order) interactions correctly", {
   skip_if_not(run_vs)
   skip_if_not_installed("rstanarm")
-  warn_L1_ia_orig <- options(projpred.warn_L1_interactions = TRUE)
+  withr::local_options(list(projpred.warn_L1_interactions = TRUE))
   main_terms_in_ia <- c("xca.2", "xco.3", "xco.1")
   all_ias_split <- lapply(seq_along(main_terms_in_ia), combn,
                           x = main_terms_in_ia, simplify = FALSE)
@@ -967,6 +961,7 @@ test_that("L1 search handles three-way (second-order) interactions correctly", {
     lapply(ia_split, all_ia_perms, is_split = TRUE)
   }))
   trms_universe_split_bu <- trms_universe_split
+  withr::defer(trms_universe_split <<- trms_universe_split_bu)
   trms_universe_split <<- union(trms_universe_split, all_ias)
   tstsetup <- head(grep("^rstanarm\\.glm\\.", names(fits), value = TRUE), 1)
   args_fit_i <- args_fit[[tstsetup]]
@@ -1020,8 +1015,6 @@ test_that("L1 search handles three-way (second-order) interactions correctly", {
     })),
     info = tstsetup
   )
-  trms_universe_split <<- trms_universe_split_bu
-  options(warn_L1_ia_orig)
 })
 
 ## search_terms -----------------------------------------------------------
@@ -1170,8 +1163,6 @@ test_that("varsel.vsel() works for `vsel` objects from cv_varsel()", {
 
 # cv_varsel() -------------------------------------------------------------
 
-context("cv_varsel()")
-
 test_that(paste(
   "`object` of class `refmodel` and arguments `method`, `cv_method`,",
   "`nterms_max`, `nclusters`, and `nclusters_pred` work"
@@ -1213,6 +1204,7 @@ test_that("invalid `object` fails", {
 })
 
 test_that("invalid `method` fails", {
+  skip_if_no_setups(refmods)
   for (tstsetup in names(refmods)) {
     expect_error(cv_varsel(refmods[[tstsetup]], method = "k-fold"),
                  "^Unexpected value for argument `method`\\.$",
@@ -1235,6 +1227,7 @@ test_that("invalid `method` fails", {
 })
 
 test_that("invalid `cv_method` fails", {
+  skip_if_no_setups(refmods)
   for (tstsetup in names(refmods)) {
     expect_error(cv_varsel(refmods[[tstsetup]], cv_method = "k-fold"),
                  "^Unknown `cv_method`\\.$",
@@ -1340,7 +1333,7 @@ if (run_more) {
   test_that("the message when cutting off the search is thrown correctly", {
     skip_if_not(run_vs)
     skip_if_not_installed("rstanarm")
-    mssg_cut_search_orig <- options(projpred.mssg_cut_search = NULL)
+    withr::local_options(list(projpred.mssg_cut_search = NULL))
     dat_gauss <- data.frame(y = df_gaussian$y, df_gaussian$x)
     stopifnot(sum(grepl("^X", names(dat_gauss))) > 19)
     fit_exceed <- suppressWarnings(rstanarm::stan_glm(
@@ -1376,7 +1369,6 @@ if (run_more) {
         info = paste("fit_not_exceed, nterms_max =", nterms_max_i %||% "NULL")
       )
     }
-    options(mssg_cut_search_orig)
   })
 }
 
@@ -1619,6 +1611,9 @@ test_that(paste(
     }
     tstsetups <- tstsetups_tmp
   }
+  # Test setups for which rstanarm::kfold() fails (the skip at the end of this
+  # test lists them):
+  kfold_failed <- character()
   for (tstsetup in tstsetups) {
     args_cvvs_i <- args_cvvs[[tstsetup]]
     tstsetup_fit <- args_cvvs_i$tstsetup_fit
@@ -1638,13 +1633,10 @@ test_that(paste(
       # belonging to the same level of a variable with group-level effects are
       # in the same fold, so prediction is performed for new levels (see, e.g.,
       # brms's GitHub issue #1286):
-      if (exists(".Random.seed", envir = .GlobalEnv)) {
-        rng_old <- get(".Random.seed", envir = .GlobalEnv)
-      }
       # Make the construction of the CV folds reproducible:
-      set.seed(seed2_tst * 3L)
-      folds_vec <- loo::kfold_split_grouped(K = K_crr, x = dat$z.1)
-      if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
+      folds_vec <- withr::with_seed(seed2_tst * 3L, {
+        loo::kfold_split_grouped(K = K_crr, x = dat$z.1)
+      })
     } else {
       folds_vec <- cv_folds(nobsv, K = K_crr, seed = seed2_tst)
     }
@@ -1659,9 +1651,9 @@ test_that(paste(
                                             cores = 1)),
                      silent = TRUE)
     if (inherits(kfold_obj, "try-error")) {
-      message("Could not test `tstsetup = \"", tstsetup, "\"` in the rstanarm ",
-              "`cvfits` test. Error message: \"",
-              attr(kfold_obj, "condition")$message, "\"")
+      kfold_failed <- c(kfold_failed, paste0(
+        tstsetup, " (", attr(kfold_obj, "condition")$message, ")"
+      ))
       next
     }
     kfold_obj <- structure(kfold_obj$fits[, "fit"], folds = folds_vec)
@@ -1718,6 +1710,9 @@ test_that(paste(
       )
     }
   }
+  skip_if(length(kfold_failed) > 0,
+          paste("rstanarm::kfold() failed for:",
+                paste(kfold_failed, collapse = "; ")))
 })
 
 test_that(paste(
@@ -1752,13 +1747,10 @@ test_that(paste(
       # belonging to the same level of a variable with group-level effects are
       # in the same fold, so prediction is performed for new levels (see, e.g.,
       # brms's GitHub issue #1286):
-      if (exists(".Random.seed", envir = .GlobalEnv)) {
-        rng_old <- get(".Random.seed", envir = .GlobalEnv)
-      }
       # Make the construction of the CV folds reproducible:
-      set.seed(seed2_tst + 10L)
-      folds_vec <- loo::kfold_split_grouped(K = K_crr, x = dat$z.1)
-      if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
+      folds_vec <- withr::with_seed(seed2_tst + 10L, {
+        loo::kfold_split_grouped(K = K_crr, x = dat$z.1)
+      })
     } else if (grepl("\\.gam\\.", tstsetup)) {
       folds_vec <- cv_folds(nobsv, K = K_crr, seed = seed2_tst + 10L)
     } else {
@@ -2638,9 +2630,7 @@ test_that("argument `folds` of run_cvfun() works", {
   if (!run_more) {
     tstsetups <- head(tstsetups, 1)
   }
-  if (exists(".Random.seed", envir = .GlobalEnv)) {
-    rng_old <- get(".Random.seed", envir = .GlobalEnv)
-  }
+  withr::local_preserve_seed()
   for (tstsetup in tstsetups) {
     set.seed(seed3_tst)
     folds_sep <- cv_folds(nobsv, K = K_tst)
@@ -2648,5 +2638,4 @@ test_that("argument `folds` of run_cvfun() works", {
     expect_identical(lapply(cvfits_sep, as.matrix),
                      lapply(cvfitss[[tstsetup]], as.matrix), info = tstsetup)
   }
-  if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)
 })