diff --git a/R/methods.R b/R/methods.R
index 6138136b7..cdd776cfd 100644
--- a/R/methods.R
+++ b/R/methods.R
@@ -287,9 +287,6 @@ proj_predict_aux <- function(proj, mu, weights, ...) {
 #'
 #' @inheritParams summary.vsel
 #' @param x An object of class `vsel` (returned by [varsel()] or [cv_varsel()]).
-#' @param baseline Either `"ref"` or `"best"` indicating whether the baseline is
-#'   the reference model or the best submodel (in terms of `stats[1]`),
-#'   respectively.
 #'
 #' @examples
 #' if (requireNamespace("rstanarm", quietly = TRUE)) {
@@ -329,10 +326,9 @@ plot.vsel <- function(
   ## compute all the statistics and fetch only those that were asked
   nfeat_baseline <- .get_nfeat_baseline(object, baseline, stats[1])
   tab <- rbind(
-    .tabulate_stats(object, stats,
-                    alpha = alpha,
-                    nfeat_baseline = nfeat_baseline),
-    .tabulate_stats(object, stats, alpha = alpha)
+    .tabulate_stats(object, stats, alpha = alpha,
+                    nfeat_baseline = nfeat_baseline, ...),
+    .tabulate_stats(object, stats, alpha = alpha, ...)
   )
   stats_table <- subset(tab, tab$delta == deltas)
   stats_ref <- subset(stats_table, stats_table$size == Inf)
@@ -430,10 +426,12 @@ plot.vsel <- function(
 #'   * `"mlpd"`: mean log predictive density, that is, `"elpd"` divided by the
 #'   number of observations.
 #'   * `"mse"`: mean squared error.
-#'   * `"rmse"`: root mean squared error.
+#'   * `"rmse"`: root mean squared error. For the corresponding standard error,
+#'   bootstrapping is used.
 #'   * `"acc"` (or its alias, `"pctcorr"`): classification accuracy
 #'   ([binomial()] family only).
-#'   * `"auc"`: area under the ROC curve ([binomial()] family only).
+#'   * `"auc"`: area under the ROC curve ([binomial()] family only). For the
+#'   corresponding standard error, bootstrapping is used.
 #' @param type One or more items from `"mean"`, `"se"`, `"lower"`, `"upper"`,
 #'   `"diff"`, and `"diff.se"` indicating which of these to compute for each
 #'   item from `stats` (mean, standard error, lower and upper confidence
@@ -449,10 +447,14 @@ plot.vsel <- function(
 #'   normal-approximation confidence intervals. For example, `alpha = 0.32`
 #'   corresponds to a coverage of 68%, i.e., one-standard-error intervals
 #'   (because of the normal approximation).
-#' @param baseline Only relevant if `deltas` is `TRUE`. Either `"ref"` or
-#'   `"best"` indicating whether the baseline is the reference model or the best
-#'   submodel (in terms of `stats[1]`), respectively.
-#' @param ... Currently ignored.
+#' @param baseline For [summary.vsel()]: Only relevant if `deltas` is `TRUE`.
+#'   For [plot.vsel()]: Always relevant. Either `"ref"` or `"best"`, indicating
+#'   whether the baseline is the reference model or the best submodel found (in
+#'   terms of `stats[1]`), respectively.
+#' @param ... Arguments passed to the internal function which is used for
+#'   bootstrapping (if applicable; see argument `stats`). Currently, relevant
+#'   arguments are `b` (the number of bootstrap samples, defaulting to `2000`)
+#'   and `seed` (see [set.seed()], defaulting to `NULL`).
 #'
 #' @examples
 #' if (requireNamespace("rstanarm", quietly = TRUE)) {
@@ -513,9 +515,9 @@ summary.vsel <- function(
   if (deltas) {
     nfeat_baseline <- .get_nfeat_baseline(object, baseline, stats[1])
     tab <- .tabulate_stats(object, stats, alpha = alpha,
-                           nfeat_baseline = nfeat_baseline)
+                           nfeat_baseline = nfeat_baseline, ...)
   } else {
-    tab <- .tabulate_stats(object, stats, alpha = alpha)
+    tab <- .tabulate_stats(object, stats, alpha = alpha, ...)
   }
   stats_table <- subset(tab, tab$size != Inf) %>%
     dplyr::group_by(.data$statistic) %>%
@@ -642,12 +644,15 @@ print.vselsummary <- function(x, digits = 1, ...) {
 #' @param ... Further arguments passed to [summary.vsel()] (apart from
 #'   argument `digits` which is passed to [print.vselsummary()]).
 #'
-#' @return The `data.frame` returned by [summary.vsel()] (invisible).
+#' @return The output of [summary.vsel()] (invisible).
 #'
 #' @export
 print.vsel <- function(x, ...) {
-  stats <- summary.vsel(x, ...)
-  print(stats, ...)
+  dot_args <- list(...)
+  stats <- do.call(summary.vsel, c(list(object = x),
+                                   dot_args[names(dot_args) != "digits"]))
+  do.call(print, c(list(x = stats),
+                   dot_args[names(dot_args) == "digits"]))
   return(invisible(stats))
 }
 
@@ -663,35 +668,30 @@ print.vsel <- function(x, ...) {
 #'   [cv_varsel()]).
 #' @param stat Statistic used for the decision. See [summary.vsel()] for
 #'   possible choices.
-#' @param alpha A number determining the (nominal) coverage `1 - alpha` of the
-#'   normal-approximation confidence intervals based on which the decision is
-#'   made. For example, `alpha = 0.32` corresponds to a coverage of 68%, i.e.,
-#'   one-standard-error intervals (because of the normal approximation). See
-#'   section "Details" below for more information.
 #' @param pct A number giving the relative proportion (*not* percents) between
 #'   baseline model and null model utilities one is willing to sacrifice. See
 #'   section "Details" below for more information.
 #' @param type Either `"upper"` or `"lower"` determining whether the decision is
 #'   based on the upper or lower confidence interval bound, respectively. See
 #'   section "Details" below for more information.
-#' @param baseline Either `"ref"` or `"best"` indicating whether the baseline is
-#'   the reference model or the best submodel (in terms of `stat[1]`),
-#'   respectively.
 #' @param warnings Mainly for internal use. A single logical value indicating
 #'   whether to throw warnings if automatic suggestion fails. Usually there is
 #'   no reason to set this to `FALSE`.
-#' @param ... Currently ignored.
+#' @param ... Arguments passed to [summary.vsel()], except for `object`, `stats`
+#'   (which is set to `stat`), `type`, and `deltas` (which is set to `TRUE`).
+#'   See section "Details" below for some important arguments which may be
+#'   passed here.
 #'
 #' @details The suggested model size is the smallest model size for which either
 #'   the lower or upper bound (depending on argument `type`) of the
 #'   normal-approximation confidence interval (with nominal coverage `1 -
-#'   alpha`) for \eqn{u_k - u_{\mbox{base}}}{u_k - u_base} (with \eqn{u_k}
-#'   denoting the \eqn{k}-th submodel's utility and
-#'   \eqn{u_{\mbox{base}}}{u_base} denoting the baseline model's utility) falls
-#'   above (or is equal to) \deqn{\mbox{pct} * (u_0 - u_{\mbox{base}})}{pct *
-#'   (u_0 - u_base)} where \eqn{u_0} denotes the null model utility. The
-#'   baseline is either the reference model or the best submodel found (see
-#'   argument `baseline`).
+#'   alpha`, see argument `alpha` of [summary.vsel()]) for \eqn{u_k -
+#'   u_{\mbox{base}}}{u_k - u_base} (with \eqn{u_k} denoting the \eqn{k}-th
+#'   submodel's utility and \eqn{u_{\mbox{base}}}{u_base} denoting the baseline
+#'   model's utility) falls above (or is equal to) \deqn{\mbox{pct} * (u_0 -
+#'   u_{\mbox{base}})}{pct * (u_0 - u_base)} where \eqn{u_0} denotes the null
+#'   model utility. The baseline is either the reference model or the best
+#'   submodel found (see argument `baseline` of [summary.vsel()]).
 #'
 #'   For example, `alpha = 0.32`, `pct = 0`, and `type = "upper"` means that we
 #'   select the smallest model size for which the upper bound of the confidence
@@ -742,10 +742,8 @@ suggest_size <- function(object, ...) {
 suggest_size.vsel <- function(
   object,
   stat = "elpd",
-  alpha = 0.32,
   pct = 0,
   type = "upper",
-  baseline = if (!inherits(object$refmodel, "datafit")) "ref" else "best",
   warnings = TRUE,
   ...
 ) {
@@ -771,9 +769,10 @@ suggest_size.vsel <- function(
   }
   bound <- type
   stats <- summary.vsel(object,
-                        stats = stat, alpha = alpha,
+                        stats = stat,
                         type = c("mean", "upper", "lower"),
-                        baseline = baseline, deltas = TRUE)$selection
+                        deltas = TRUE,
+                        ...)$selection
   util_null <- sgn * unlist(unname(subset(
     stats, stats$size == 0,
     paste0(stat, suffix)
diff --git a/R/misc.R b/R/misc.R
index 9235701cf..f86dd5ec6 100644
--- a/R/misc.R
+++ b/R/misc.R
@@ -56,7 +56,7 @@ auc <- function(x) {
 # Bootstrap an arbitrary quantity `fun` that takes the sample `x` as the first
 # input. Other arguments of `fun` can be passed by `...`. Example:
 # `boostrap(x, mean)`.
-bootstrap <- function(x, fun = mean, b = 1000, seed = NULL, ...) {
+bootstrap <- function(x, fun = mean, b = 2000, seed = NULL, ...) {
   # set random seed but ensure the old RNG state is restored on exit
   if (exists(".Random.seed")) {
     rng_state_old <- .Random.seed
diff --git a/R/summary_funs.R b/R/summary_funs.R
index cb7e49177..87f698b73 100644
--- a/R/summary_funs.R
+++ b/R/summary_funs.R
@@ -46,7 +46,7 @@
 # statistics relative to the baseline model of that size (`nfeat_baseline = Inf`
 # means that the baseline model is the reference model).
 .tabulate_stats <- function(varsel, stats, alpha = 0.05,
-                            nfeat_baseline = NULL) {
+                            nfeat_baseline = NULL, ...) {
   stat_tab <- data.frame()
   summ_ref <- varsel$summaries$ref
   summ_sub <- varsel$summaries$sub
@@ -75,7 +75,7 @@
     ## reference model statistics
     summ <- summ_ref
     res <- get_stat(summ$mu, summ$lppd, varsel$d_test, stat, mu.bs = mu.bs,
-                    lppd.bs = lppd.bs, weights = summ$w, alpha = alpha)
+                    lppd.bs = lppd.bs, weights = summ$w, alpha = alpha, ...)
     row <- data.frame(
       data = varsel$d_test$type, size = Inf, delta = delta, statistic = stat,
       value = res$value, lq = res$lq, uq = res$uq, se = res$se, diff = NA,
@@ -93,10 +93,10 @@
         ## scale
         res_ref <- get_stat(summ_ref$mu, summ_ref$lppd, varsel$d_test,
                             stat, mu.bs = NULL, lppd.bs = NULL,
-                            weights = summ_ref$w, alpha = alpha)
+                            weights = summ_ref$w, alpha = alpha, ...)
         res_diff <- get_stat(summ$mu, summ$lppd, varsel$d_test, stat,
                              mu.bs = summ_ref$mu, lppd.bs = summ_ref$lppd,
-                             weights = summ$w, alpha = alpha)
+                             weights = summ$w, alpha = alpha, ...)
         val <- res_ref$value + res_diff$value
         val.se <- sqrt(res_ref$se^2 + res_diff$se^2)
         lq <- qnorm(alpha / 2, mean = val, sd = val.se)
@@ -109,10 +109,10 @@
       } else {
         ## normal case
         res <- get_stat(summ$mu, summ$lppd, varsel$d_test, stat, mu.bs = mu.bs,
-                        lppd.bs = lppd.bs, weights = summ$w, alpha = alpha)
+                        lppd.bs = lppd.bs, weights = summ$w, alpha = alpha, ...)
         diff <- get_stat(summ$mu, summ$lppd, varsel$d_test, stat,
                          mu.bs = summ_ref$mu, lppd.bs = summ_ref$lppd,
-                         weights = summ$w, alpha = alpha)
+                         weights = summ$w, alpha = alpha, ...)
         row <- data.frame(
           data = varsel$d_test$type, size = k - 1, delta = delta,
           statistic = stat, value = res$value, lq = res$lq, uq = res$uq,
@@ -127,7 +127,7 @@
 }
 
 get_stat <- function(mu, lppd, d_test, stat, mu.bs = NULL, lppd.bs = NULL,
-                     weights = NULL, alpha = 0.1, seed = 1208499, B = 2000) {
+                     weights = NULL, alpha = 0.1, ...) {
   ##
   ## Calculates given statistic stat with standard error and confidence bounds.
   ## mu.bs and lppd.bs are the pointwise mu and lppd for another model that is
@@ -195,16 +195,14 @@ get_stat <- function(mu, lppd, d_test, stat, mu.bs = NULL, lppd.bs = NULL,
         function(resid2) {
           sqrt(mean(weights * resid2, na.rm = TRUE))
         },
-        b = B,
-        seed = seed
+        ...
       )
       value.bootstrap2 <- bootstrap(
         (mu.bs - y)^2,
         function(resid2) {
           sqrt(mean(weights * resid2, na.rm = TRUE))
         },
-        b = B,
-        seed = seed
+        ...
       )
       value.se <- sd(value.bootstrap1 - value.bootstrap2)
     } else {
@@ -214,8 +212,7 @@ get_stat <- function(mu, lppd, d_test, stat, mu.bs = NULL, lppd.bs = NULL,
         function(resid2) {
           sqrt(mean(weights * resid2, na.rm = TRUE))
         },
-        b = B,
-        seed = seed
+        ...
       )
       value.se <- sd(value.bootstrap)
     }
@@ -240,12 +237,12 @@ get_stat <- function(mu, lppd, d_test, stat, mu.bs = NULL, lppd.bs = NULL,
       mu[is.na(mu.bs)] <- NA # for which both mu and mu.bs are non-NA
       auc.data.bs <- cbind(y, mu.bs, weights)
       value <- auc(auc.data) - auc(auc.data.bs)
-      value.bootstrap1 <- bootstrap(auc.data, auc, b = B, seed = seed)
-      value.bootstrap2 <- bootstrap(auc.data.bs, auc, b = B, seed = seed)
+      value.bootstrap1 <- bootstrap(auc.data, auc, ...)
+      value.bootstrap2 <- bootstrap(auc.data.bs, auc, ...)
       value.se <- sd(value.bootstrap1 - value.bootstrap2, na.rm = TRUE)
     } else {
       value <- auc(auc.data)
-      value.bootstrap <- bootstrap(auc.data, auc, b = B, seed = seed)
+      value.bootstrap <- bootstrap(auc.data, auc, ...)
       value.se <- sd(value.bootstrap, na.rm = TRUE)
     }
   }
diff --git a/man/plot.vsel.Rd b/man/plot.vsel.Rd
index 2fd16ee91..c160dbc55 100644
--- a/man/plot.vsel.Rd
+++ b/man/plot.vsel.Rd
@@ -29,10 +29,12 @@ calculate. Available statistics are:
 \item \code{"mlpd"}: mean log predictive density, that is, \code{"elpd"} divided by the
 number of observations.
 \item \code{"mse"}: mean squared error.
-\item \code{"rmse"}: root mean squared error.
+\item \code{"rmse"}: root mean squared error. For the corresponding standard error,
+bootstrapping is used.
 \item \code{"acc"} (or its alias, \code{"pctcorr"}): classification accuracy
 (\code{\link[=binomial]{binomial()}} family only).
-\item \code{"auc"}: area under the ROC curve (\code{\link[=binomial]{binomial()}} family only).
+\item \code{"auc"}: area under the ROC curve (\code{\link[=binomial]{binomial()}} family only). For the
+corresponding standard error, bootstrapping is used.
 }}
 
 \item{deltas}{If \code{TRUE}, the submodel statistics are estimated as differences
@@ -44,11 +46,15 @@ normal-approximation confidence intervals. For example, \code{alpha = 0.32}
 corresponds to a coverage of 68\%, i.e., one-standard-error intervals
 (because of the normal approximation).}
 
-\item{baseline}{Either \code{"ref"} or \code{"best"} indicating whether the baseline is
-the reference model or the best submodel (in terms of \code{stats[1]}),
-respectively.}
+\item{baseline}{For \code{\link[=summary.vsel]{summary.vsel()}}: Only relevant if \code{deltas} is \code{TRUE}.
+For \code{\link[=plot.vsel]{plot.vsel()}}: Always relevant. Either \code{"ref"} or \code{"best"}, indicating
+whether the baseline is the reference model or the best submodel found (in
+terms of \code{stats[1]}), respectively.}
 
-\item{...}{Currently ignored.}
+\item{...}{Arguments passed to the internal function which is used for
+bootstrapping (if applicable; see argument \code{stats}). Currently, relevant
+arguments are \code{b} (the number of bootstrap samples, defaulting to \code{2000})
+and \code{seed} (see \code{\link[=set.seed]{set.seed()}}, defaulting to \code{NULL}).}
 }
 \description{
 This is the \code{\link[=plot]{plot()}} method for \code{vsel} objects (returned by \code{\link[=varsel]{varsel()}} or
diff --git a/man/print.vsel.Rd b/man/print.vsel.Rd
index e957f7750..a1884c00b 100644
--- a/man/print.vsel.Rd
+++ b/man/print.vsel.Rd
@@ -13,7 +13,7 @@
 argument \code{digits} which is passed to \code{\link[=print.vselsummary]{print.vselsummary()}}).}
 }
 \value{
-The \code{data.frame} returned by \code{\link[=summary.vsel]{summary.vsel()}} (invisible).
+The output of \code{\link[=summary.vsel]{summary.vsel()}} (invisible).
 }
 \description{
 This is the \code{\link[=print]{print()}} method for \code{vsel} objects (returned by \code{\link[=varsel]{varsel()}} or
diff --git a/man/suggest_size.Rd b/man/suggest_size.Rd
index 2ae66f878..0ca84f8ba 100644
--- a/man/suggest_size.Rd
+++ b/man/suggest_size.Rd
@@ -10,10 +10,8 @@ suggest_size(object, ...)
 \method{suggest_size}{vsel}(
   object,
   stat = "elpd",
-  alpha = 0.32,
   pct = 0,
   type = "upper",
-  baseline = if (!inherits(object$refmodel, "datafit")) "ref" else "best",
   warnings = TRUE,
   ...
 )
@@ -22,17 +20,14 @@ suggest_size(object, ...)
 \item{object}{An object of class \code{vsel} (returned by \code{\link[=varsel]{varsel()}} or
 \code{\link[=cv_varsel]{cv_varsel()}}).}
 
-\item{...}{Currently ignored.}
+\item{...}{Arguments passed to \code{\link[=summary.vsel]{summary.vsel()}}, except for \code{object}, \code{stats}
+(which is set to \code{stat}), \code{type}, and \code{deltas} (which is set to \code{TRUE}).
+See section "Details" below for some important arguments which may be
+passed here.}
 
 \item{stat}{Statistic used for the decision. See \code{\link[=summary.vsel]{summary.vsel()}} for
 possible choices.}
 
-\item{alpha}{A number determining the (nominal) coverage \code{1 - alpha} of the
-normal-approximation confidence intervals based on which the decision is
-made. For example, \code{alpha = 0.32} corresponds to a coverage of 68\%, i.e.,
-one-standard-error intervals (because of the normal approximation). See
-section "Details" below for more information.}
-
 \item{pct}{A number giving the relative proportion (\emph{not} percents) between
 baseline model and null model utilities one is willing to sacrifice. See
 section "Details" below for more information.}
@@ -41,10 +36,6 @@ section "Details" below for more information.}
 based on the upper or lower confidence interval bound, respectively. See
 section "Details" below for more information.}
 
-\item{baseline}{Either \code{"ref"} or \code{"best"} indicating whether the baseline is
-the reference model or the best submodel (in terms of \code{stat[1]}),
-respectively.}
-
 \item{warnings}{Mainly for internal use. A single logical value indicating
 whether to throw warnings if automatic suggestion fails. Usually there is
 no reason to set this to \code{FALSE}.}
@@ -59,13 +50,13 @@ final decision based on what is most appropriate for the problem at hand.
 \details{
 The suggested model size is the smallest model size for which either
 the lower or upper bound (depending on argument \code{type}) of the
-normal-approximation confidence interval (with nominal coverage \code{1 - alpha}) for \eqn{u_k - u_{\mbox{base}}}{u_k - u_base} (with \eqn{u_k}
-denoting the \eqn{k}-th submodel's utility and
-\eqn{u_{\mbox{base}}}{u_base} denoting the baseline model's utility) falls
-above (or is equal to) \deqn{\mbox{pct} * (u_0 - u_{\mbox{base}})}{pct *
-  (u_0 - u_base)} where \eqn{u_0} denotes the null model utility. The
-baseline is either the reference model or the best submodel found (see
-argument \code{baseline}).
+normal-approximation confidence interval (with nominal coverage \code{1 - alpha}, see argument \code{alpha} of \code{\link[=summary.vsel]{summary.vsel()}}) for \eqn{u_k -
+  u_{\mbox{base}}}{u_k - u_base} (with \eqn{u_k} denoting the \eqn{k}-th
+submodel's utility and \eqn{u_{\mbox{base}}}{u_base} denoting the baseline
+model's utility) falls above (or is equal to) \deqn{\mbox{pct} * (u_0 -
+  u_{\mbox{base}})}{pct * (u_0 - u_base)} where \eqn{u_0} denotes the null
+model utility. The baseline is either the reference model or the best
+submodel found (see argument \code{baseline} of \code{\link[=summary.vsel]{summary.vsel()}}).
 
 For example, \code{alpha = 0.32}, \code{pct = 0}, and \code{type = "upper"} means that we
 select the smallest model size for which the upper bound of the confidence
diff --git a/man/summary.vsel.Rd b/man/summary.vsel.Rd
index 60469a1b5..87bb00c6a 100644
--- a/man/summary.vsel.Rd
+++ b/man/summary.vsel.Rd
@@ -31,10 +31,12 @@ calculate. Available statistics are:
 \item \code{"mlpd"}: mean log predictive density, that is, \code{"elpd"} divided by the
 number of observations.
 \item \code{"mse"}: mean squared error.
-\item \code{"rmse"}: root mean squared error.
+\item \code{"rmse"}: root mean squared error. For the corresponding standard error,
+bootstrapping is used.
 \item \code{"acc"} (or its alias, \code{"pctcorr"}): classification accuracy
 (\code{\link[=binomial]{binomial()}} family only).
-\item \code{"auc"}: area under the ROC curve (\code{\link[=binomial]{binomial()}} family only).
+\item \code{"auc"}: area under the ROC curve (\code{\link[=binomial]{binomial()}} family only). For the
+corresponding standard error, bootstrapping is used.
 }}
 
 \item{type}{One or more items from \code{"mean"}, \code{"se"}, \code{"lower"}, \code{"upper"},
@@ -55,11 +57,15 @@ normal-approximation confidence intervals. For example, \code{alpha = 0.32}
 corresponds to a coverage of 68\%, i.e., one-standard-error intervals
 (because of the normal approximation).}
 
-\item{baseline}{Only relevant if \code{deltas} is \code{TRUE}. Either \code{"ref"} or
-\code{"best"} indicating whether the baseline is the reference model or the best
-submodel (in terms of \code{stats[1]}), respectively.}
+\item{baseline}{For \code{\link[=summary.vsel]{summary.vsel()}}: Only relevant if \code{deltas} is \code{TRUE}.
+For \code{\link[=plot.vsel]{plot.vsel()}}: Always relevant. Either \code{"ref"} or \code{"best"}, indicating
+whether the baseline is the reference model or the best submodel found (in
+terms of \code{stats[1]}), respectively.}
 
-\item{...}{Currently ignored.}
+\item{...}{Arguments passed to the internal function which is used for
+bootstrapping (if applicable; see argument \code{stats}). Currently, relevant
+arguments are \code{b} (the number of bootstrap samples, defaulting to \code{2000})
+and \code{seed} (see \code{\link[=set.seed]{set.seed()}}, defaulting to \code{NULL}).}
 }
 \description{
 This is the \code{\link[=summary]{summary()}} method for \code{vsel} objects (returned by \code{\link[=varsel]{varsel()}} or
diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R
index afc8755e2..7360b3166 100644
--- a/tests/testthat/setup.R
+++ b/tests/testthat/setup.R
@@ -621,6 +621,7 @@ type_tst <- c("mean", "lower", "upper", "se")
 
 seed_tst <- 74345
 seed2_tst <- 866028
+seed3_tst <- 1208499
 
 ## Reference model --------------------------------------------------------
 
@@ -965,9 +966,15 @@ if (run_vs) {
   args_smmry_vs <- unlist_cust(args_smmry_vs)
 
   smmrys_vs <- lapply(args_smmry_vs, function(args_smmry_vs_i) {
+    if (any(c("rmse", "auc") %in% args_smmry_vs_i$stats)) {
+      smmry_seed <- list(seed = seed3_tst)
+    } else {
+      smmry_seed <- list()
+    }
     do.call(summary, c(
       list(object = vss[[args_smmry_vs_i$tstsetup_vsel]]),
-      excl_nonargs(args_smmry_vs_i)
+      excl_nonargs(args_smmry_vs_i),
+      smmry_seed
     ))
   })
 }
@@ -989,9 +996,15 @@ if (run_cvvs) {
   args_smmry_cvvs <- unlist_cust(args_smmry_cvvs)
 
   smmrys_cvvs <- lapply(args_smmry_cvvs, function(args_smmry_cvvs_i) {
+    if (any(c("rmse", "auc") %in% args_smmry_cvvs_i$stats)) {
+      smmry_seed <- list(seed = seed3_tst)
+    } else {
+      smmry_seed <- list()
+    }
     do.call(summary, c(
       list(object = cvvss[[args_smmry_cvvs_i$tstsetup_vsel]]),
-      excl_nonargs(args_smmry_cvvs_i)
+      excl_nonargs(args_smmry_cvvs_i),
+      smmry_seed
     ))
   })
 }
diff --git a/tests/testthat/test_datafit.R b/tests/testthat/test_datafit.R
index 1f2bc229d..3de2207d2 100644
--- a/tests/testthat/test_datafit.R
+++ b/tests/testthat/test_datafit.R
@@ -519,7 +519,8 @@ test_that(paste(
   for (tstsetup in tstsetups) {
     smmry <- summary(vss_datafit[[tstsetup]],
                      stats = stats_common,
-                     type = type_tst)
+                     type = type_tst,
+                     seed = seed3_tst)
     smmry_tester(
       smmry,
       vsel_expected = vss_datafit[[tstsetup]],
@@ -555,7 +556,8 @@ test_that(paste(
   for (tstsetup in tstsetups) {
     smmry <- summary(cvvss_datafit[[tstsetup]],
                      stats = stats_common,
-                     type = type_tst)
+                     type = type_tst,
+                     seed = seed3_tst)
     smmry_tester(
       smmry,
       vsel_expected = cvvss_datafit[[tstsetup]],
diff --git a/tests/testthat/test_methods_vsel.R b/tests/testthat/test_methods_vsel.R
index 5fb9f10a5..c2d8611b4 100644
--- a/tests/testthat/test_methods_vsel.R
+++ b/tests/testthat/test_methods_vsel.R
@@ -154,10 +154,16 @@ test_that(paste(
   skip_if_not(run_vs)
   for (tstsetup in head(names(smmrys_vs), 1)) {
     args_smmry_vs_i <- args_smmry_vs[[tstsetup]]
+    if (any(c("rmse", "auc") %in% args_smmry_vs_i$stats)) {
+      smmry_seed <- list(seed = seed3_tst)
+    } else {
+      smmry_seed <- list()
+    }
     expect_output(
       print_obj <- do.call(print, c(
         list(x = vss[[args_smmry_vs_i$tstsetup_vsel]]),
-        excl_nonargs(args_smmry_vs_i)
+        excl_nonargs(args_smmry_vs_i),
+        smmry_seed
       )),
       "Family:.*Link function:.*Formula:.*Observations:",
       info = tstsetup
@@ -173,10 +179,16 @@ test_that(paste(
   skip_if_not(run_cvvs)
   for (tstsetup in head(names(smmrys_cvvs), 1)) {
     args_smmry_cvvs_i <- args_smmry_cvvs[[tstsetup]]
+    if (any(c("rmse", "auc") %in% args_smmry_cvvs_i$stats)) {
+      smmry_seed <- list(seed = seed3_tst)
+    } else {
+      smmry_seed <- list()
+    }
     expect_output(
       print_obj <- do.call(print, c(
         list(x = cvvss[[args_smmry_cvvs_i$tstsetup_vsel]]),
-        excl_nonargs(args_smmry_cvvs_i)
+        excl_nonargs(args_smmry_cvvs_i),
+        smmry_seed
       )),
       "Family:.*Link function:.*Formula:.*Observations:",
       info = tstsetup
@@ -273,10 +285,15 @@ test_that("`stat` works", {
                           "common_stats")
     stat_vec <- stats_tst[[stat_crr_nm]]$stats
     for (stat_crr in stat_vec) {
+      if (stat_crr %in% c("rmse", "auc")) {
+        suggsize_seed <- seed3_tst
+      } else {
+        suggsize_seed <- NULL
+      }
       # Warnings are suppressed, but a suggested size of `NA` (because of a
       # search which was terminated too early) is tested below:
       suggsize <- suppressWarnings(
-        suggest_size(vss[[tstsetup_vs]], stat = stat_crr)
+        suggest_size(vss[[tstsetup_vs]], stat = stat_crr, seed = suggsize_seed)
       )
       expect_type(suggsize, "double")
       expect_length(suggsize, 1)