diff --git a/NEWS.md b/NEWS.md
index baf486604..a93f4f838 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -11,6 +11,7 @@ If you read this from a place other than <https://mc-stan.org/projpred/news/inde
 * The standard error for performance statistic `"rmse"` is now computed via the delta method (instead of bootstrapping). The confidence interval for `"rmse"` is now based on a log-normal approximation (instead of bootstrapping) if argument `deltas` of `summary.vsel()` or `plot.vsel()` is `FALSE` and based on a normal approximation (instead of bootstrapping) if `deltas` is `TRUE`. (GitHub: #496)
 * Performance statistic `"R2"` (R-squared) has been added, see argument `stats` of `summary.vsel()` and `plot.vsel()`; argument `stat` of `suggest_size()` supports it as well. (GitHub: #483, #496)
 * The performance evaluation part of `cv_varsel()` with `cv_method = "LOO"` and `validate_search = FALSE` now always applies Pareto smoothing when computing the importance sampling weights (as long as the number of importance ratios in the tail is large enough; otherwise, no Pareto smoothing is applied). Previously, in case of projected draws with nonconstant weights (i.e., in case of clustering), no Pareto smoothing had been applied. (GitHub: #496, #507)
+* Argument `deltas` of `plot.vsel()` has gained the option `"mixed"` which combines the point estimates from `deltas = FALSE` with the uncertainty bars from `deltas = TRUE`. (GitHub: #511)
 
 ## Minor changes
 
diff --git a/R/methods.R b/R/methods.R
index a60d640b2..3f3ce17b7 100644
--- a/R/methods.R
+++ b/R/methods.R
@@ -594,6 +594,17 @@ proj_predict_aux <- function(proj, newdata, offsetnew, weightsnew,
 #'
 #' @inheritParams summary.vsel
 #' @param x An object of class `vsel` (returned by [varsel()] or [cv_varsel()]).
+#' @param deltas May be set to `FALSE`, `TRUE`, or `"mixed"`. If `FALSE`, the
+#'   submodel performance statistics are plotted on their actual scale and the
+#'   uncertainty bars match this scale. If `TRUE`, the submodel statistics are
+#'   plotted relatively to the baseline model (see argument `baseline`) and the
+#'   uncertainty bars match this scale. For the GMPD, the term "relatively"
+#'   refers to the *ratio* vs. the baseline model (i.e., the submodel statistic
+#'   divided by the baseline model statistic). For all other `stats`,
+#'   "relatively" refers to the *difference* from the baseline model (i.e., the
+#'   submodel statistic minus the baseline model statistic). If set to
+#'   `"mixed"`, the `deltas = FALSE` point estimates are combined with the
+#'   uncertainty bars from the `deltas = TRUE` plot.
 #' @param thres_elpd Only relevant if `any(stats %in% c("elpd", "mlpd",
 #'   "gmpd"))`. The threshold for the ELPD difference (taking the submodel's
 #'   ELPD minus the baseline model's ELPD) above which the submodel's ELPD is
@@ -673,9 +684,10 @@ proj_predict_aux <- function(proj, newdata, offsetnew, weightsnew,
 #' the baseline model's performance is shown as a dotted black horizontal line.
 #' If `!is.na(thres_elpd)` and `any(stats %in% c("elpd", "mlpd", "gmpd"))`, the
 #' value supplied to `thres_elpd` (which is automatically adapted internally in
-#' case of the MLPD or the GMPD or `deltas = FALSE`) is shown as a dot-dashed
-#' gray horizontal line for the reference model and, if `baseline = "best"`, as
-#' a long-dashed green horizontal line for the baseline model.
+#' case of the MLPD or the GMPD or `deltas = FALSE` or `deltas = "mixed"`) is
+#' shown as a dot-dashed gray horizontal line for the reference model and, if
+#' `baseline = "best"`, as a long-dashed green horizontal line for the baseline
+#' model.
 #'
 #' @examplesIf requireNamespace("rstanarm", quietly = TRUE)
 #' # Data:
@@ -732,21 +744,57 @@ plot.vsel <- function(
   } else if (!is.null(ranking_repel)) {
     stopifnot(isTRUE(ranking_repel %in% c("text", "label")))
   }
+  stopifnot(isTRUE(deltas) || isFALSE(deltas) || identical(deltas, "mixed"))
 
   # Define `nfeat_baseline` and a slightly modified variant that can be used for
   # .tabulate_stats()'s argument `nfeat_baseline`:
   nfeat_baseline <- get_nfeat_baseline(object, baseline, stats[1],
                                        resp_oscale = resp_oscale)
-  if (deltas) {
+  nfeat_baseline_mixed <- NULL
+  if (isTRUE(deltas)) {
     nfeat_baseline_for_tab <- nfeat_baseline
   } else {
     nfeat_baseline_for_tab <- NULL
+    if (identical(deltas, "mixed")) {
+      nfeat_baseline_mixed <- nfeat_baseline
+    }
   }
 
   # Compute the predictive performance statistics:
   stats_table_all <- .tabulate_stats(object, stats, alpha = alpha,
                                      nfeat_baseline = nfeat_baseline_for_tab,
+                                     nfeat_baseline_diff = nfeat_baseline_mixed,
                                      resp_oscale = resp_oscale, ...)
+  if (identical(deltas, "mixed")) {
+    stats_bs_init <- subset(stats_table_all,
+                            stats_table_all$size == nfeat_baseline)
+    names(stats_bs_init)[names(stats_bs_init) == "value"] <- "value_bs"
+    stats_table_all <- merge(stats_table_all,
+                             stats_bs_init[, c("statistic", "value_bs")],
+                             by = "statistic", all.x = TRUE, all.y = FALSE,
+                             sort = FALSE)
+    is_gmpd_entry <- stats_table_all[["statistic"]] == "gmpd"
+    stats_table_all[["diff.lq"]][is_gmpd_entry] <-
+      log(stats_table_all[["diff.lq"]][is_gmpd_entry])
+    stats_table_all[["diff.uq"]][is_gmpd_entry] <-
+      log(stats_table_all[["diff.uq"]][is_gmpd_entry])
+    stats_table_all[["value_bs"]][is_gmpd_entry] <-
+      log(stats_table_all[["value_bs"]][is_gmpd_entry])
+    stats_table_all[["lq"]] <- stats_table_all[["diff.lq"]] +
+      stats_table_all[["value_bs"]]
+    stats_table_all[["uq"]] <- stats_table_all[["diff.uq"]] +
+      stats_table_all[["value_bs"]]
+    stats_table_all[["diff.lq"]][is_gmpd_entry] <-
+      exp(stats_table_all[["diff.lq"]][is_gmpd_entry])
+    stats_table_all[["diff.uq"]][is_gmpd_entry] <-
+      exp(stats_table_all[["diff.uq"]][is_gmpd_entry])
+    stats_table_all[["value_bs"]][is_gmpd_entry] <-
+      exp(stats_table_all[["value_bs"]][is_gmpd_entry])
+    stats_table_all[["lq"]][is_gmpd_entry] <-
+      exp(stats_table_all[["lq"]][is_gmpd_entry])
+    stats_table_all[["uq"]][is_gmpd_entry] <-
+      exp(stats_table_all[["uq"]][is_gmpd_entry])
+  }
   stats_ref <- subset(stats_table_all, stats_table_all$size == Inf)
   stats_sub <- subset(stats_table_all, stats_table_all$size != Inf)
   stats_bs <- subset(stats_table_all, stats_table_all$size == nfeat_baseline)
@@ -783,7 +831,7 @@ plot.vsel <- function(
   } else {
     baseline_pretty <- "best submodel"
   }
-  if (deltas) {
+  if (isTRUE(deltas)) {
     if (all(stats != "gmpd")) {
       ylab <- paste0("Difference vs. ", baseline_pretty)
     } else if (all(stats == "gmpd")) {
@@ -1070,14 +1118,31 @@ plot.vsel <- function(
     ci_type <- "bootstrap "
   } else if (all(stats %in% c("gmpd"))) {
     ci_type <- "exponentiated normal-approximation "
-  } else if (all(stats %in% c("mse", "rmse")) && !deltas) {
+  } else if (all(stats %in% c("mse", "rmse")) && isFALSE(deltas)) {
     ci_type <- "log-normal-approximation "
   } else if (all(!stats %in% c("auc", "gmpd", "mse", "rmse")) ||
-             (all(!stats %in% c("auc", "gmpd")) && deltas)) {
+             (all(!stats %in% c("auc", "gmpd")) &&
+              (identical(deltas, "mixed") || isTRUE(deltas)))) {
     ci_type <- "normal-approximation "
   } else {
     ci_type <- ""
   }
+  interval_description <- paste0("Vertical bars indicate ",
+                                 round(100 * (1 - alpha), 1), "% ", ci_type,
+                                 "intervals")
+  if (identical(deltas, "mixed")) {
+    interval_description <- paste0(interval_description,
+                                   ", \nshowing uncertainty for ")
+    if (all(stats != "gmpd")) {
+      diff_pretty <- "difference"
+    } else if (all(stats == "gmpd")) {
+      diff_pretty <- "ratio"
+    } else {
+      diff_pretty <- "difference (ratio for GMPD)"
+    }
+    interval_description <- paste0(interval_description, diff_pretty, " vs. ",
+                                   baseline_pretty)
+  }
   if (identical(size_position, "secondary_x")) {
     tick_labs_x_sec <- as.character(rk_dfr[
       order(match(rk_dfr[["size"]], breaks), na.last = NA),
@@ -1108,9 +1173,7 @@ plot.vsel <- function(
                        labels = tick_labs_x,
                        sec.axis = x_axis_sec) +
     labs(x = xlab, y = ylab, title = "Predictive performance",
-         subtitle = paste0("Vertical bars indicate ",
-                           round(100 * (1 - alpha), 1), "% ", ci_type,
-                           "intervals")) +
+         subtitle = interval_description) +
     theme(axis.text.x.bottom = element_text(angle = text_angle,
                                             hjust = hjust_val,
                                             vjust = vjust_val)) +
@@ -1179,12 +1242,14 @@ plot.vsel <- function(
 #'   * `"mse"`: mean squared error (only available in the situations mentioned
 #'   in section "Details" below). For the corresponding confidence interval, a
 #'   log-normal approximation is used if `deltas` is `FALSE` and a normal
-#'   approximation is used if `deltas` is `TRUE`.
+#'   approximation is used if `deltas` is `TRUE` (or `"mixed"`, in case of
+#'   [plot.vsel()]).
 #'   * `"rmse"`: root mean squared error (only available in the situations
 #'   mentioned in section "Details" below). For the corresponding standard
 #'   error, the delta method is used. For the corresponding confidence interval,
 #'   a log-normal approximation is used if `deltas` is `FALSE` and a normal
-#'   approximation is used if `deltas` is `TRUE`.
+#'   approximation is used if `deltas` is `TRUE` (or `"mixed"`, in case of
+#'   [plot.vsel()]).
 #'   * `"R2"`: R-squared, i.e., coefficient of determination (only available in
 #'   the situations mentioned in section "Details" below). For the corresponding
 #'   standard error, the delta method is used. For the corresponding confidence
@@ -1202,20 +1267,24 @@ plot.vsel <- function(
 #'   supported in case of subsampled LOO-CV (see argument `nloo` of
 #'   [cv_varsel()]).
 #' @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
-#'   interval bounds, mean difference to the corresponding statistic of the
-#'   reference model, and standard error of this difference, respectively; note
-#'   that for the GMPD, `"diff"`, and `"diff.se"` actually refer to the ratio
-#'   vs. the reference model, not the difference). The confidence interval
-#'   bounds belong to confidence intervals with (nominal) coverage `1 - alpha`.
-#'   Items `"diff"` and `"diff.se"` are only supported if `deltas` is `FALSE`.
-#' @param deltas If `TRUE`, the submodel statistics are estimated relatively to
-#'   the baseline model (see argument `baseline`). For the GMPD, the term
-#'   "relatively" refers to the ratio vs. the baseline model (i.e., the submodel
-#'   statistic divided by the baseline model statistic). For all other `stats`,
-#'   "relatively" refers to the difference from the baseline model (i.e., the
-#'   submodel statistic minus the baseline model statistic).
+#'   `"diff"`, `"diff.lower"`, `"diff.upper"`, and `"diff.se"` indicating which
+#'   of these to compute for each item from `stats` (mean, standard error, lower
+#'   and upper confidence interval bounds, mean difference to the corresponding
+#'   statistic of the reference model, lower and upper confidence interval bound
+#'   for this difference, and standard error of this difference, respectively;
+#'   note that for the GMPD, `"diff"`, `"diff.lower"`, `"diff.upper"`, and
+#'   `"diff.se"` actually refer to the ratio vs. the reference model, not the
+#'   difference). The confidence interval bounds belong to confidence intervals
+#'   with (nominal) coverage `1 - alpha`. Items `"diff"`, `"diff.lower"`,
+#'   `"diff.upper"`, and `"diff.se"` are only supported if `deltas` is `FALSE`.
+#' @param deltas May be set to `FALSE` or `TRUE`. If `FALSE`, the submodel
+#'   performance statistics are estimated on their actual scale. If `TRUE`, the
+#'   submodel statistics are estimated relatively to the baseline model (see
+#'   argument `baseline`). For the GMPD, the term "relatively" refers to the
+#'   *ratio* vs. the baseline model (i.e., the submodel statistic divided by the
+#'   baseline model statistic). For all other `stats`, "relatively" refers to
+#'   the *difference* from the baseline model (i.e., the submodel statistic
+#'   minus the baseline model statistic).
 #' @param alpha A number determining the (nominal) coverage `1 - alpha` of the
 #'   confidence intervals. For example, in case of a normal-approximation
 #'   confidence interval, `alpha = 2 * pnorm(-1)` corresponds to a confidence
@@ -1407,8 +1476,8 @@ summary.vsel <- function(
 # reference model performance and one table for the submodel performance):
 mk_colnms_smmry <- function(type, stats, deltas) {
   # Pre-process `type`:
-  if (is.null(deltas) || deltas) {
-    type <- setdiff(type, c("diff", "diff.se"))
+  if (is.null(deltas) || isTRUE(deltas)) {
+    type <- setdiff(type, c("diff", "diff.lower", "diff.upper", "diff.se"))
   }
   type_dot <- paste0(".", type)
   type_dot[type_dot == ".mean"] <- ""
@@ -1418,6 +1487,8 @@ mk_colnms_smmry <- function(type, stats, deltas) {
   nms_old[nms_old == "mean"] <- "value"
   nms_old[nms_old == "upper"] <- "uq"
   nms_old[nms_old == "lower"] <- "lq"
+  nms_old[nms_old == "diff.upper"] <- "diff.uq"
+  nms_old[nms_old == "diff.lower"] <- "diff.lq"
   # The clean column names that should be used in the output table:
   nms_new <- lapply(stats, paste0, type_dot)
   return(nlist(nms_old, nms_new))
diff --git a/R/misc.R b/R/misc.R
index b02b226fe..8e8052e60 100644
--- a/R/misc.R
+++ b/R/misc.R
@@ -204,12 +204,12 @@ validate_baseline <- function(vsel_obj, baseline, deltas) {
   if (!(baseline %in% c("ref", "best"))) {
     stop("Argument 'baseline' must be either 'ref' or 'best'.")
   }
-  if (baseline == "ref" && deltas == TRUE &&
+  if (baseline == "ref" && (identical(deltas, "mixed") || isTRUE(deltas)) &&
       inherits(vsel_obj$refmodel, "datafit")) {
     # no reference model (or the results missing for some other reason),
     # so cannot compute differences (or ratios) vs. the reference model
-    stop("Cannot use deltas = TRUE and baseline = 'ref' when there is no ",
-         "reference model.")
+    stop("Cannot use `deltas = TRUE` (or `deltas = \"mixed\"`) and ",
+         "`baseline = \"ref\"` when there is no reference model.")
   }
   if (baseline == "best" && vsel_obj$cv_method == "LOO" &&
       isTRUE(vsel_obj$nloo < vsel_obj$refmodel$nobs)) {
diff --git a/R/summary_funs.R b/R/summary_funs.R
index e6f45ec36..9de6fe11a 100644
--- a/R/summary_funs.R
+++ b/R/summary_funs.R
@@ -86,9 +86,14 @@ weighted_summary_means <- function(y_wobs_test, family, wdraws, mu, dis, cl_ref,
 # errors, and confidence intervals with coverage `1 - alpha` based on the
 # variable selection output. If `nfeat_baseline` is given, then compute the
 # statistics relative to the baseline model of that size (`nfeat_baseline = Inf`
-# means that the baseline model is the reference model).
+# means that the baseline model is the reference model). Argument
+# `nfeat_baseline_diff` is currently only used (i.e., non-`NULL`) in
+# plot.vsel()'s `deltas = "mixed"` case. In that case, `nfeat_baseline_diff`
+# gives the size of the baseline model analogously to `nfeat_baseline` (in
+# particular, `nfeat_baseline_diff = Inf` designates the reference model).
 .tabulate_stats <- function(varsel, stats, alpha = 0.05,
-                            nfeat_baseline = NULL, resp_oscale = TRUE, ...) {
+                            nfeat_baseline = NULL, nfeat_baseline_diff = NULL,
+                            resp_oscale = TRUE, ...) {
   stat_tab <- data.frame()
   summaries_ref <- varsel$summaries$ref
   summaries_sub <- varsel$summaries$sub
@@ -230,6 +235,15 @@ weighted_summary_means <- function(y_wobs_test, family, wdraws, mu, dis, cl_ref,
     }
     delta <- TRUE
   }
+  if (is.null(nfeat_baseline_diff)) {
+    summaries_baseline_diff <- summaries_ref
+  } else {
+    if (nfeat_baseline_diff == Inf) {
+      summaries_baseline_diff <- summaries_ref
+    } else {
+      summaries_baseline_diff <- summaries_sub[[nfeat_baseline_diff + 1]]
+    }
+  }
 
   for (s in seq_along(stats)) {
     stat <- stats[s]
@@ -242,8 +256,8 @@ weighted_summary_means <- function(y_wobs_test, family, wdraws, mu, dis, cl_ref,
                     varsel$y_wobs_test, stat, alpha = alpha, ...)
     row <- data.frame(
       data = varsel$type_test, size = Inf, delta = delta, statistic = stat,
-      value = res$value, lq = res$lq, uq = res$uq, se = res$se, diff = NA,
-      diff.se = NA
+      value = res$value, lq = res$lq, uq = res$uq, se = res$se,
+      diff = NA, diff.lq = NA, diff.uq = NA, diff.se = NA
     )
     stat_tab <- rbind(stat_tab, row)
 
@@ -254,15 +268,18 @@ weighted_summary_means <- function(y_wobs_test, family, wdraws, mu, dis, cl_ref,
                       summaries_fast = summaries_fast_sub[[k]],
                       loo_inds = varsel$loo_inds,
                       varsel$y_wobs_test, stat, alpha = alpha, ...)
+      # TODO: If `!is.null(nfeat_baseline) && nfeat_baseline == Inf` and
+      # `is.null(nfeat_baseline_diff)`, we should be able to set `diff <- res`
+      # and hence avoid a second call to get_stat().
       diff <- get_stat(summaries = summaries_sub[[k]],
-                       summaries_baseline = summaries_ref,
+                       summaries_baseline = summaries_baseline_diff,
                        summaries_fast = summaries_fast_sub[[k]],
                        loo_inds = varsel$loo_inds,
                        varsel$y_wobs_test, stat, alpha = alpha, ...)
       row <- data.frame(
         data = varsel$type_test, size = k - 1, delta = delta, statistic = stat,
         value = res$value, lq = res$lq, uq = res$uq, se = res$se,
-        diff = diff$value, diff.se = diff$se
+        diff = diff$value, diff.lq = diff$lq, diff.uq = diff$uq, diff.se = diff$se
       )
       stat_tab <- rbind(stat_tab, row)
     }
diff --git a/man/plot.vsel.Rd b/man/plot.vsel.Rd
index ee30aab2d..625112582 100644
--- a/man/plot.vsel.Rd
+++ b/man/plot.vsel.Rd
@@ -63,12 +63,14 @@ of the MLPD.
 \item \code{"mse"}: mean squared error (only available in the situations mentioned
 in section "Details" below). For the corresponding confidence interval, a
 log-normal approximation is used if \code{deltas} is \code{FALSE} and a normal
-approximation is used if \code{deltas} is \code{TRUE}.
+approximation is used if \code{deltas} is \code{TRUE} (or \code{"mixed"}, in case of
+\code{\link[=plot.vsel]{plot.vsel()}}).
 \item \code{"rmse"}: root mean squared error (only available in the situations
 mentioned in section "Details" below). For the corresponding standard
 error, the delta method is used. For the corresponding confidence interval,
 a log-normal approximation is used if \code{deltas} is \code{FALSE} and a normal
-approximation is used if \code{deltas} is \code{TRUE}.
+approximation is used if \code{deltas} is \code{TRUE} (or \code{"mixed"}, in case of
+\code{\link[=plot.vsel]{plot.vsel()}}).
 \item \code{"R2"}: R-squared, i.e., coefficient of determination (only available in
 the situations mentioned in section "Details" below). For the corresponding
 standard error, the delta method is used. For the corresponding confidence
@@ -87,12 +89,17 @@ supported in case of subsampled LOO-CV (see argument \code{nloo} of
 \code{\link[=cv_varsel]{cv_varsel()}}).
 }}
 
-\item{deltas}{If \code{TRUE}, the submodel statistics are estimated relatively to
-the baseline model (see argument \code{baseline}). For the GMPD, the term
-"relatively" refers to the ratio vs. the baseline model (i.e., the submodel
-statistic divided by the baseline model statistic). For all other \code{stats},
-"relatively" refers to the difference from the baseline model (i.e., the
-submodel statistic minus the baseline model statistic).}
+\item{deltas}{May be set to \code{FALSE}, \code{TRUE}, or \code{"mixed"}. If \code{FALSE}, the
+submodel performance statistics are plotted on their actual scale and the
+uncertainty bars match this scale. If \code{TRUE}, the submodel statistics are
+plotted relatively to the baseline model (see argument \code{baseline}) and the
+uncertainty bars match this scale. For the GMPD, the term "relatively"
+refers to the \emph{ratio} vs. the baseline model (i.e., the submodel statistic
+divided by the baseline model statistic). For all other \code{stats},
+"relatively" refers to the \emph{difference} from the baseline model (i.e., the
+submodel statistic minus the baseline model statistic). If set to
+\code{"mixed"}, the \code{deltas = FALSE} point estimates are combined with the
+uncertainty bars from the \code{deltas = TRUE} plot.}
 
 \item{alpha}{A number determining the (nominal) coverage \code{1 - alpha} of the
 confidence intervals. For example, in case of a normal-approximation
@@ -241,9 +248,10 @@ shown in the plot as a dashed red horizontal line. If \code{baseline = "best"},
 the baseline model's performance is shown as a dotted black horizontal line.
 If \code{!is.na(thres_elpd)} and \code{any(stats \%in\% c("elpd", "mlpd", "gmpd"))}, the
 value supplied to \code{thres_elpd} (which is automatically adapted internally in
-case of the MLPD or the GMPD or \code{deltas = FALSE}) is shown as a dot-dashed
-gray horizontal line for the reference model and, if \code{baseline = "best"}, as
-a long-dashed green horizontal line for the baseline model.
+case of the MLPD or the GMPD or \code{deltas = FALSE} or \code{deltas = "mixed"}) is
+shown as a dot-dashed gray horizontal line for the reference model and, if
+\code{baseline = "best"}, as a long-dashed green horizontal line for the baseline
+model.
 }
 
 \examples{
diff --git a/man/summary.vsel.Rd b/man/summary.vsel.Rd
index cc131e62d..bc2172d50 100644
--- a/man/summary.vsel.Rd
+++ b/man/summary.vsel.Rd
@@ -53,12 +53,14 @@ of the MLPD.
 \item \code{"mse"}: mean squared error (only available in the situations mentioned
 in section "Details" below). For the corresponding confidence interval, a
 log-normal approximation is used if \code{deltas} is \code{FALSE} and a normal
-approximation is used if \code{deltas} is \code{TRUE}.
+approximation is used if \code{deltas} is \code{TRUE} (or \code{"mixed"}, in case of
+\code{\link[=plot.vsel]{plot.vsel()}}).
 \item \code{"rmse"}: root mean squared error (only available in the situations
 mentioned in section "Details" below). For the corresponding standard
 error, the delta method is used. For the corresponding confidence interval,
 a log-normal approximation is used if \code{deltas} is \code{FALSE} and a normal
-approximation is used if \code{deltas} is \code{TRUE}.
+approximation is used if \code{deltas} is \code{TRUE} (or \code{"mixed"}, in case of
+\code{\link[=plot.vsel]{plot.vsel()}}).
 \item \code{"R2"}: R-squared, i.e., coefficient of determination (only available in
 the situations mentioned in section "Details" below). For the corresponding
 standard error, the delta method is used. For the corresponding confidence
@@ -78,21 +80,25 @@ supported in case of subsampled LOO-CV (see argument \code{nloo} of
 }}
 
 \item{type}{One or more items from \code{"mean"}, \code{"se"}, \code{"lower"}, \code{"upper"},
-\code{"diff"}, and \code{"diff.se"} indicating which of these to compute for each
-item from \code{stats} (mean, standard error, lower and upper confidence
-interval bounds, mean difference to the corresponding statistic of the
-reference model, and standard error of this difference, respectively; note
-that for the GMPD, \code{"diff"}, and \code{"diff.se"} actually refer to the ratio
-vs. the reference model, not the difference). The confidence interval
-bounds belong to confidence intervals with (nominal) coverage \code{1 - alpha}.
-Items \code{"diff"} and \code{"diff.se"} are only supported if \code{deltas} is \code{FALSE}.}
+\code{"diff"}, \code{"diff.lower"}, \code{"diff.upper"}, and \code{"diff.se"} indicating which
+of these to compute for each item from \code{stats} (mean, standard error, lower
+and upper confidence interval bounds, mean difference to the corresponding
+statistic of the reference model, lower and upper confidence interval bound
+for this difference, and standard error of this difference, respectively;
+note that for the GMPD, \code{"diff"}, \code{"diff.lower"}, \code{"diff.upper"}, and
+\code{"diff.se"} actually refer to the ratio vs. the reference model, not the
+difference). The confidence interval bounds belong to confidence intervals
+with (nominal) coverage \code{1 - alpha}. Items \code{"diff"}, \code{"diff.lower"},
+\code{"diff.upper"}, and \code{"diff.se"} are only supported if \code{deltas} is \code{FALSE}.}
 
-\item{deltas}{If \code{TRUE}, the submodel statistics are estimated relatively to
-the baseline model (see argument \code{baseline}). For the GMPD, the term
-"relatively" refers to the ratio vs. the baseline model (i.e., the submodel
-statistic divided by the baseline model statistic). For all other \code{stats},
-"relatively" refers to the difference from the baseline model (i.e., the
-submodel statistic minus the baseline model statistic).}
+\item{deltas}{May be set to \code{FALSE} or \code{TRUE}. If \code{FALSE}, the submodel
+performance statistics are estimated on their actual scale. If \code{TRUE}, the
+submodel statistics are estimated relatively to the baseline model (see
+argument \code{baseline}). For the GMPD, the term "relatively" refers to the
+\emph{ratio} vs. the baseline model (i.e., the submodel statistic divided by the
+baseline model statistic). For all other \code{stats}, "relatively" refers to
+the \emph{difference} from the baseline model (i.e., the submodel statistic
+minus the baseline model statistic).}
 
 \item{alpha}{A number determining the (nominal) coverage \code{1 - alpha} of the
 confidence intervals. For example, in case of a normal-approximation
diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R
index 8bfd0257f..cc6e69b60 100755
--- a/tests/testthat/setup.R
+++ b/tests/testthat/setup.R
@@ -991,7 +991,8 @@ angle_tst <- list(
 
 deltas_tst_plot <- list(
   default_deltas = list(),
-  deltas_TRUE = list(deltas = TRUE)
+  deltas_TRUE = list(deltas = TRUE),
+  deltas_mixed = list(deltas = "mixed")
 )
 
 ### nterms ----------------------------------------------------------------
diff --git a/tests/testthat/test_datafit.R b/tests/testthat/test_datafit.R
index 42d056bbe..d6e32df7c 100644
--- a/tests/testthat/test_datafit.R
+++ b/tests/testthat/test_datafit.R
@@ -500,8 +500,8 @@ test_that("summary.vsel(): `baseline = \"ref\"` and `deltas = TRUE` fails", {
   for (tstsetup in head(names(vss_datafit), 1)) {
     expect_error(
       summary(vss_datafit[[tstsetup]], baseline = "ref", deltas = TRUE),
-      paste("^Cannot use deltas = TRUE and baseline = 'ref' when there is no",
-            "reference model\\.$"),
+      paste0("^Cannot use `deltas = TRUE` .+or `deltas = .+mixed.+`.+ and ",
+             "`baseline = .+ref.+ when there is no reference model\\.$"),
       info = tstsetup
     )
   }
diff --git a/vignettes/latent.Rmd b/vignettes/latent.Rmd
index d7aeaf200..8d9c6b555 100644
--- a/vignettes/latent.Rmd
+++ b/vignettes/latent.Rmd
@@ -197,7 +197,7 @@ The message telling that `<refmodel>$dis` consists of only `NA`s will not concer
 
 In order to decide for a submodel size, we first inspect the `plot()` results:
 ```{r plot_vsel_lat}
-( gg_lat <- plot(vs_lat, stats = "mlpd", deltas = TRUE) )
+( gg_lat <- plot(vs_lat, stats = "mlpd", deltas = "mixed") )
 ```
 
 Although the submodels' MLPDs seem to be very close to the reference model's MLPD from a submodel size of 6 on, a zoomed plot reveals that there is still some discrepancy at sizes 6 to 11 and that size 12 would be a better choice (further down below in the `summary()` output, we will also see that on absolute scale, the discrepancy at sizes 6 to 11 is not negligible):
@@ -269,7 +269,7 @@ rm(warn_instable_orig)
 ```
 ```{r post_vs_trad}
 print(time_trad)
-( gg_trad <- plot(vs_trad, stats = "mlpd", deltas = TRUE) )
+( gg_trad <- plot(vs_trad, stats = "mlpd", deltas = "mixed") )
 smmry_trad <- summary(vs_trad, stats = "mlpd",
                       type = c("mean", "lower", "upper", "diff"))
 print(smmry_trad, digits = 2)
@@ -365,7 +365,7 @@ The message concerning `latent_ilink` can be safely ignored here (the internal d
 
 Again, we first inspect the `plot()` results to decide for a submodel size:
 ```{r plot_vsel_nebin}
-( gg_nebin <- plot(vs_nebin, stats = "mlpd", deltas = TRUE) )
+( gg_nebin <- plot(vs_nebin, stats = "mlpd", deltas = "mixed") )
 ```
 
 Again, a zoomed plot is more helpful:
diff --git a/vignettes/projpred.Rmd b/vignettes/projpred.Rmd
index 73dced63a..e46310ec4 100755
--- a/vignettes/projpred.Rmd
+++ b/vignettes/projpred.Rmd
@@ -228,10 +228,14 @@ foreach::registerDoSEQ()
 <!-- Again, we ignore the tail-ESS warnings due to the reduced values for `chains` and `iter` in the reference model fit. -->
 
 We can now select a final submodel size by looking at a predictive performance plot similar to the one created for the preliminary `cv_varsel()` run above.
-By default, the performance statistics are plotted on their original scale, but with `deltas = TRUE`, they are plotted as differences^[For the geometric mean predictive density (GMPD, see argument `stats` of `summary.vsel()`), `deltas = TRUE` means to calculate the GMPD *ratio* (not difference) vs. the baseline model.] from a baseline model (which is the reference model by default, at least in the most common cases).
-Since the differences^[For the GMPD, this is again a ratio.] and the (frequentist) uncertainty in their estimation are usually of more interest than the original-scale performance statistics (at least with regard to the decision for a final submodel size), we directly plot with `deltas = TRUE` here:
+By default, the performance statistics are plotted on their actual scale and the uncertainty bars match this scale, but argument `deltas` offers two more options:
+
+*   With `deltas = TRUE`, the performance statistics are plotted as differences^[For the geometric mean predictive density (GMPD, see argument `stats` of `summary.vsel()` and `plot.vsel()`), `deltas = TRUE` plots the GMPD *ratio* (not difference) vs. the baseline model. Despite this special case, we will call `deltas = TRUE` the "difference scale" for simplicity.] from the baseline model^[For the definition of the baseline model, see argument `baseline` of `summary.vsel()` and `plot.vsel()`; in the most common cases, the default baseline model is the reference model.] and the uncertainty bars match this scale,
+*   With `deltas = "mixed"`, the performance statistics (i.e., their point estimates) are plotted on the actual scale, but the uncertainty bars visualize the difference-scale uncertainty.
+
+Since the difference-scale uncertainty is usually of more interest than the actual-scale uncertainty (at least with regard to the decision for a final submodel size) and the actual-scale point estimates are often of more interest than the difference-scale point estimates, we plot with `deltas = "mixed"` here:
 ```{r plot_vsel}
-plot(cvvs, stats = "mlpd", deltas = TRUE)
+plot(cvvs, stats = "mlpd", deltas = "mixed")
 ```
 
 ### Decision for final submodel size
@@ -260,8 +264,10 @@ A tabular representation of the plot created by `plot.vsel()` can be achieved vi
 For the output of `summary.vsel()`, there is a sophisticated `print()` method (`print.vselsummary()`) which is also called by the shortcut method `print.vsel()`^[`print.vsel()` is the method that is called when simply printing an object resulting from `varsel()` or `cv_varsel()`.].
 Specifically, to create the summary table matching the predictive performance plot above as closely as possible (and to also adjust the minimum number of printed significant digits), we may call `summary.vsel()` and `print.vselsummary()` as follows:
 ```{r smmry_vsel}
-smmry <- summary(cvvs, stats = "mlpd", type = c("mean", "lower", "upper"),
-                 deltas = TRUE)
+smmry <- summary(cvvs,
+                 stats = "mlpd",
+                 type = c("mean", "diff.lower", "diff.upper"),
+                 deltas = FALSE)
 print(smmry, digits = 1)
 ```