From 2424a48c46eeeaa557a0986f7282cd1a3e9e01c5 Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Tue, 12 Jul 2022 17:14:09 +0200
Subject: [PATCH 01/16] Simplify the creation of `d_test` in case of
 `is.null(d_test)`.

---
 R/varsel.R | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/R/varsel.R b/R/varsel.R
index 27c990006..dcec4efde 100644
--- a/R/varsel.R
+++ b/R/varsel.R
@@ -202,10 +202,9 @@ varsel.refmodel <- function(object, d_test = NULL, method = NULL,
 
   if (is.null(d_test)) {
     d_type <- "train"
-    test_points <- seq_len(NROW(refmodel$y))
-    d_test <- nlist(
-      y = refmodel$y, test_points, data = NULL, weights = refmodel$wobs,
-      type = d_type, offset = refmodel$offset
+    d_test <- list(
+      y = refmodel$y, test_points = seq_len(NROW(refmodel$y)), data = NULL,
+      weights = refmodel$wobs, type = d_type, offset = refmodel$offset
     )
   } else {
     d_type <- d_test$type

From 5596d7601c92aa33d308e1ff6c2f4603d1286b5e Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Tue, 12 Jul 2022 17:16:12 +0200
Subject: [PATCH 02/16] Further simplify `d_test`-related code (now `d_type`).

---
 R/varsel.R | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/R/varsel.R b/R/varsel.R
index dcec4efde..d54a83cd3 100644
--- a/R/varsel.R
+++ b/R/varsel.R
@@ -201,13 +201,10 @@ varsel.refmodel <- function(object, d_test = NULL, method = NULL,
   search_terms <- args$search_terms
 
   if (is.null(d_test)) {
-    d_type <- "train"
     d_test <- list(
       y = refmodel$y, test_points = seq_len(NROW(refmodel$y)), data = NULL,
-      weights = refmodel$wobs, type = d_type, offset = refmodel$offset
+      weights = refmodel$wobs, type = "train", offset = refmodel$offset
     )
-  } else {
-    d_type <- d_test$type
   }
 
   ## reference distributions for selection and prediction after selection
@@ -243,7 +240,7 @@ varsel.refmodel <- function(object, d_test = NULL, method = NULL,
     ntest <- NROW(refmodel$y)
     ref <- list(mu = rep(NA, ntest), lppd = rep(NA, ntest))
   } else {
-    if (d_type == "train") {
+    if (d_test$type == "train") {
       mu_test <- refmodel$mu
       if (!all(refmodel$offset == 0)) {
         mu_test <- refmodel$family$linkinv(

From efb67f225e7136041af4c14ff7f98e3aaf0ff6ad Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Tue, 12 Jul 2022 21:07:33 +0200
Subject: [PATCH 03/16] Fix argument `d_test` of `varsel()`: The predictive
 performance of the

submodels also needs to be evaluated on the test data (not only the predictive
performance of the reference model).
---
 R/summary_funs.R | 14 ++++++++------
 R/varsel.R       |  8 +++++---
 man/varsel.Rd    |  3 ++-
 3 files changed, 15 insertions(+), 10 deletions(-)

diff --git a/R/summary_funs.R b/R/summary_funs.R
index 87f698b73..7aa5e05e0 100644
--- a/R/summary_funs.R
+++ b/R/summary_funs.R
@@ -1,13 +1,15 @@
-.get_sub_summaries <- function(submodels, test_points, refmodel) {
+.get_sub_summaries <- function(submodels, test_points, refmodel,
+                               y = refmodel$y[test_points],
+                               wobs = refmodel$wobs[test_points],
+                               newdata = NULL,
+                               offset = refmodel$offset[test_points]) {
   lapply(submodels, function(model) {
     .weighted_summary_means(
-      y_test = list(y = refmodel$y[test_points],
-                    weights = refmodel$wobs[test_points]),
+      y_test = list(y = y, weights = wobs),
       family = refmodel$family,
       wsample = model$weights,
-      mu = refmodel$family$mu_fun(model$submodl,
-                                  obs = test_points,
-                                  offset = refmodel$offset[test_points]),
+      mu = refmodel$family$mu_fun(model$submodl, obs = test_points,
+                                  newdata = newdata, offset = offset),
       dis = model$dis
     )
   })
diff --git a/R/varsel.R b/R/varsel.R
index d54a83cd3..ca022730b 100644
--- a/R/varsel.R
+++ b/R/varsel.R
@@ -12,7 +12,8 @@
 #'   [get_refmodel()].
 #' @param d_test For internal use only. A `list` providing information about the
 #'   test set which is used for evaluating the predictive performance of the
-#'   reference model. If not provided, the training set is used.
+#'   submodels as well as of the reference model. If `NULL`, the training set is
+#'   used.
 #' @param method The method for the search part. Possible options are `"L1"` for
 #'   L1 search and `"forward"` for forward search. If `NULL`, then internally,
 #'   `"L1"` is used, except if the reference model has multilevel or additive
@@ -227,8 +228,9 @@ varsel.refmodel <- function(object, d_test = NULL, method = NULL,
     ...
   )
   sub <- .get_sub_summaries(
-    submodels = submodels, test_points = seq_along(refmodel$y),
-    refmodel = refmodel
+    submodels = submodels, test_points = d_test$test_points,
+    refmodel = refmodel, y = d_test$y, wobs = d_test$weights,
+    newdata = d_test$data, offset = d_test$offset
   )
 
   ## predictive statistics of the reference model on test data. if no test data
diff --git a/man/varsel.Rd b/man/varsel.Rd
index 8ed681910..a01c2400f 100644
--- a/man/varsel.Rd
+++ b/man/varsel.Rd
@@ -42,7 +42,8 @@ the latter only if \code{refit_prj} is \code{TRUE}).}
 
 \item{d_test}{For internal use only. A \code{list} providing information about the
 test set which is used for evaluating the predictive performance of the
-reference model. If not provided, the training set is used.}
+submodels as well as of the reference model. If \code{NULL}, the training set is
+used.}
 
 \item{method}{The method for the search part. Possible options are \code{"L1"} for
 L1 search and \code{"forward"} for forward search. If \code{NULL}, then internally,

From 667725c639ede3ee8192ec25d4a47c8d4cabe3b5 Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Tue, 12 Jul 2022 21:33:48 +0200
Subject: [PATCH 04/16] Element `test_points` of `d_test` is actually not
 needed (and could

even cause confusion and be used incorrectly, given the new
`.get_sub_summaries()` arguments).
---
 R/cv_varsel.R | 9 ++-------
 R/varsel.R    | 8 +++-----
 2 files changed, 5 insertions(+), 12 deletions(-)

diff --git a/R/cv_varsel.R b/R/cv_varsel.R
index ffde602f6..dee4c9800 100644
--- a/R/cv_varsel.R
+++ b/R/cv_varsel.R
@@ -516,12 +516,8 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws,
   summ_ref <- list(lppd = loo_ref, mu = mu_ref)
   summaries <- list(sub = summ_sub, ref = summ_ref)
 
-  d_test <- list(
-    y = refmodel$y, type = "LOO",
-    test_points = seq_along(refmodel$y),
-    weights = refmodel$wobs,
-    data = NULL, offset = refmodel$offset
-  )
+  d_test <- list(y = refmodel$y, type = "LOO", weights = refmodel$wobs,
+                 data = NULL, offset = refmodel$offset)
 
   out_list <- nlist(solution_terms_cv = solution_terms_mat, summaries, d_test)
   if (!validate_search) {
@@ -643,7 +639,6 @@ kfold_varsel <- function(refmodel, method, nterms_max, ndraws,
   d_cv <- rbind2list(lapply(list_cv, function(fold) {
     list(y = fold$d_test$y,
          weights = fold$d_test$weights,
-         test_points = fold$d_test$omitted,
          offset = fold$d_test$offset)
   }))
 
diff --git a/R/varsel.R b/R/varsel.R
index ca022730b..a0a0f8eed 100644
--- a/R/varsel.R
+++ b/R/varsel.R
@@ -202,10 +202,8 @@ varsel.refmodel <- function(object, d_test = NULL, method = NULL,
   search_terms <- args$search_terms
 
   if (is.null(d_test)) {
-    d_test <- list(
-      y = refmodel$y, test_points = seq_len(NROW(refmodel$y)), data = NULL,
-      weights = refmodel$wobs, type = "train", offset = refmodel$offset
-    )
+    d_test <- list(y = refmodel$y, data = NULL, weights = refmodel$wobs,
+                   type = "train", offset = refmodel$offset)
   }
 
   ## reference distributions for selection and prediction after selection
@@ -228,7 +226,7 @@ varsel.refmodel <- function(object, d_test = NULL, method = NULL,
     ...
   )
   sub <- .get_sub_summaries(
-    submodels = submodels, test_points = d_test$test_points,
+    submodels = submodels, test_points = seq_along(d_test$y),
     refmodel = refmodel, y = d_test$y, wobs = d_test$weights,
     newdata = d_test$data, offset = d_test$offset
   )

From cc4a378d0408e2b1f0bd7d17352178ffa5190585 Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Tue, 12 Jul 2022 21:43:02 +0200
Subject: [PATCH 05/16] Use a consistent order of arguments like `newdata`,
 `offset`, `wobs`, etc.

---
 R/cv_varsel.R    | 13 ++++++-------
 R/summary_funs.R |  7 +++----
 R/varsel.R       | 16 +++++++++-------
 3 files changed, 18 insertions(+), 18 deletions(-)

diff --git a/R/cv_varsel.R b/R/cv_varsel.R
index dee4c9800..db7118e17 100644
--- a/R/cv_varsel.R
+++ b/R/cv_varsel.R
@@ -481,9 +481,9 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws,
         p_ref = p_pred, refmodel = refmodel, regul = opt$regul,
         refit_prj = refit_prj, ...
       )
-      summaries_sub <- .get_sub_summaries(
-        submodels = submodels, test_points = c(i), refmodel = refmodel
-      )
+      summaries_sub <- .get_sub_summaries(submodels = submodels,
+                                          refmodel = refmodel,
+                                          test_points = i)
       for (k in seq_along(summaries_sub)) {
         loo_sub[[k]][i] <- summaries_sub[[k]]$lppd
         mu_sub[[k]][i] <- summaries_sub[[k]]$mu
@@ -602,10 +602,9 @@ kfold_varsel <- function(refmodel, method, nterms_max, ndraws,
   # Perform the evaluation of the submodels for each fold (and make sure to
   # combine the results from the K folds into a single results list):
   get_summaries_submodel_cv <- function(submodels, fold) {
-    .get_sub_summaries(
-      submodels = submodels, test_points = fold$d_test$omitted,
-      refmodel = refmodel
-    )
+    .get_sub_summaries(submodels = submodels,
+                       refmodel = refmodel,
+                       test_points = fold$d_test$omitted)
   }
   sub_cv_summaries <- mapply(get_summaries_submodel_cv, submodels_cv, list_cv)
   if (is.null(dim(sub_cv_summaries))) {
diff --git a/R/summary_funs.R b/R/summary_funs.R
index 7aa5e05e0..d541c675d 100644
--- a/R/summary_funs.R
+++ b/R/summary_funs.R
@@ -1,8 +1,7 @@
-.get_sub_summaries <- function(submodels, test_points, refmodel,
-                               y = refmodel$y[test_points],
+.get_sub_summaries <- function(submodels, refmodel, test_points, newdata = NULL,
+                               offset = refmodel$offset[test_points],
                                wobs = refmodel$wobs[test_points],
-                               newdata = NULL,
-                               offset = refmodel$offset[test_points]) {
+                               y = refmodel$y[test_points]) {
   lapply(submodels, function(model) {
     .weighted_summary_means(
       y_test = list(y = y, weights = wobs),
diff --git a/R/varsel.R b/R/varsel.R
index a0a0f8eed..8b5e120d7 100644
--- a/R/varsel.R
+++ b/R/varsel.R
@@ -202,8 +202,8 @@ varsel.refmodel <- function(object, d_test = NULL, method = NULL,
   search_terms <- args$search_terms
 
   if (is.null(d_test)) {
-    d_test <- list(y = refmodel$y, data = NULL, weights = refmodel$wobs,
-                   type = "train", offset = refmodel$offset)
+    d_test <- list(type = "train", data = NULL, offset = refmodel$offset,
+                   weights = refmodel$wobs, y = refmodel$y)
   }
 
   ## reference distributions for selection and prediction after selection
@@ -225,11 +225,13 @@ varsel.refmodel <- function(object, d_test = NULL, method = NULL,
     p_ref = p_pred, refmodel = refmodel, regul = regul, refit_prj = refit_prj,
     ...
   )
-  sub <- .get_sub_summaries(
-    submodels = submodels, test_points = seq_along(d_test$y),
-    refmodel = refmodel, y = d_test$y, wobs = d_test$weights,
-    newdata = d_test$data, offset = d_test$offset
-  )
+  sub <- .get_sub_summaries(submodels = submodels,
+                            refmodel = refmodel,
+                            test_points = seq_along(d_test$y),
+                            newdata = d_test$data,
+                            offset = d_test$offset,
+                            wobs = d_test$weights,
+                            y = d_test$y)
 
   ## predictive statistics of the reference model on test data. if no test data
   ## are provided,

From ec44c09bfc6be9064fdff532ff5cd467183787c5 Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Tue, 12 Jul 2022 21:47:17 +0200
Subject: [PATCH 06/16] `varsel()`: Pick an even more straightforward value for

`test_points`, namely `NULL` (see argument `obs` of `fetch_data()`).
---
 R/varsel.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/varsel.R b/R/varsel.R
index 8b5e120d7..9edf62a26 100644
--- a/R/varsel.R
+++ b/R/varsel.R
@@ -227,7 +227,7 @@ varsel.refmodel <- function(object, d_test = NULL, method = NULL,
   )
   sub <- .get_sub_summaries(submodels = submodels,
                             refmodel = refmodel,
-                            test_points = seq_along(d_test$y),
+                            test_points = NULL,
                             newdata = d_test$data,
                             offset = d_test$offset,
                             wobs = d_test$weights,

From 246e949c41636776b78153395419280239857f19 Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Tue, 12 Jul 2022 21:55:18 +0200
Subject: [PATCH 07/16] Define `d_test` consistently (see also `varsel()`).

---
 R/cv_varsel.R | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/R/cv_varsel.R b/R/cv_varsel.R
index db7118e17..f251b4bbf 100644
--- a/R/cv_varsel.R
+++ b/R/cv_varsel.R
@@ -516,8 +516,8 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws,
   summ_ref <- list(lppd = loo_ref, mu = mu_ref)
   summaries <- list(sub = summ_sub, ref = summ_ref)
 
-  d_test <- list(y = refmodel$y, type = "LOO", weights = refmodel$wobs,
-                 data = NULL, offset = refmodel$offset)
+  d_test <- list(type = "LOO", data = NULL, offset = refmodel$offset,
+                 weights = refmodel$wobs, y = refmodel$y)
 
   out_list <- nlist(solution_terms_cv = solution_terms_mat, summaries, d_test)
   if (!validate_search) {
@@ -636,14 +636,14 @@ kfold_varsel <- function(refmodel, method, nterms_max, ndraws,
   # Combine the K separate test "datasets" (rather "information objects") into a
   # single list:
   d_cv <- rbind2list(lapply(list_cv, function(fold) {
-    list(y = fold$d_test$y,
+    list(offset = fold$d_test$offset,
          weights = fold$d_test$weights,
-         offset = fold$d_test$offset)
+         y = fold$d_test$y)
   }))
 
   return(nlist(solution_terms_cv,
                summaries = nlist(sub, ref),
-               d_test = c(d_cv, type = "kfold")))
+               d_test = c(list(type = "kfold", data = NULL), d_cv)))
 }
 
 # Re-fit the reference model K times (once for each fold; `cvfun` case) or fetch

From 1de3c2285ad7d2d1bfc88900559ca3445a826deb Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Wed, 13 Jul 2022 09:39:54 +0200
Subject: [PATCH 08/16] Requiring the user to always specify `d_test$type =
 "test"` is

tedious and not necessary. So add that element internally.
---
 R/misc.R   | 4 ++++
 R/varsel.R | 3 +++
 2 files changed, 7 insertions(+)

diff --git a/R/misc.R b/R/misc.R
index e5a46582a..05cc763ce 100644
--- a/R/misc.R
+++ b/R/misc.R
@@ -4,6 +4,10 @@
   packageStartupMessage(msg)
 }
 
+nms_d_test <- function() {
+  c("type", "data", "offset", "weights", "y")
+}
+
 weighted.sd <- function(x, w, na.rm = FALSE) {
   if (na.rm) {
     ind <- !is.na(w) & !is.na(x)
diff --git a/R/varsel.R b/R/varsel.R
index 9edf62a26..4d4dfdc2a 100644
--- a/R/varsel.R
+++ b/R/varsel.R
@@ -204,6 +204,9 @@ varsel.refmodel <- function(object, d_test = NULL, method = NULL,
   if (is.null(d_test)) {
     d_test <- list(type = "train", data = NULL, offset = refmodel$offset,
                    weights = refmodel$wobs, y = refmodel$y)
+  } else {
+    d_test$type <- "test"
+    d_test <- d_test[nms_d_test()]
   }
 
   ## reference distributions for selection and prediction after selection

From 24c6d0bc0ad5e23ffddf69d96695d43c33725ce2 Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Wed, 13 Jul 2022 11:07:06 +0200
Subject: [PATCH 09/16] Docs: Declare `d_test` as available for everyone (not
 for internal

use only). Also update `NEWS.md`.
---
 NEWS.md       |  2 ++
 R/varsel.R    | 21 +++++++++++++++++----
 man/varsel.Rd | 21 +++++++++++++++++----
 3 files changed, 36 insertions(+), 8 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 4d8dbf775..3400116b4 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,6 +5,7 @@
 * Several improvements in the documentation (especially in the explanation of the `suggest_size()` heuristic).
 * At multiple places throughout the package: Improvement of the numerical stability for some link functions, achieved by avoiding unnecessary back-and-forth transformations between latent space and response space. (GitHub: #337, #338)
 * All arguments `seed` and `.seed` are now allowed to be `NA` for not calling `set.seed()` internally at all.
+* Argument `d_test` of `varsel()` is not considered as an internal feature anymore. This was possible after fixing a bug for `d_test` (see below).
 
 ## Bug fixes
 
@@ -16,6 +17,7 @@
 * Fix GitHub issue #331.
 * `plot.vsel()` now draws the dashed red horizontal line for the reference model (and---if present---the dotted black horizontal line for the baseline model) first (i.e., before the submodel-specific graphical elements), to avoid overplotting.
 * Fix GitHub issue #339. (GitHub: #340)
+* Fix argument `d_test` of `varsel()`: Not only the predictive performance of the *reference model* needs to be evaluated on the test data, but also the predictive performance of the *submodels*.
 
 # projpred 2.1.2
 
diff --git a/R/varsel.R b/R/varsel.R
index 4d4dfdc2a..c82738249 100644
--- a/R/varsel.R
+++ b/R/varsel.R
@@ -10,10 +10,10 @@
 #' @param object An object of class `refmodel` (returned by [get_refmodel()] or
 #'   [init_refmodel()]) or an object that can be passed to argument `object` of
 #'   [get_refmodel()].
-#' @param d_test For internal use only. A `list` providing information about the
-#'   test set which is used for evaluating the predictive performance of the
-#'   submodels as well as of the reference model. If `NULL`, the training set is
-#'   used.
+#' @param d_test A `list` of the structure outlined in section "Argument
+#'   `d_test`" below, providing test data for evaluating the predictive
+#'   performance of the submodels as well as of the reference model. If `NULL`,
+#'   the training data is used.
 #' @param method The method for the search part. Possible options are `"L1"` for
 #'   L1 search and `"forward"` for forward search. If `NULL`, then internally,
 #'   `"L1"` is used, except if the reference model has multilevel or additive
@@ -83,6 +83,19 @@
 #'   minimizer (during a forward search and also during the evaluation part, but
 #'   the latter only if `refit_prj` is `TRUE`).
 #'
+#' @details
+#'
+#' # Argument `d_test`
+#'
+#' If not `NULL`, then `d_test` needs to be a `list` with the following
+#' elements:
+#' * `data`: a `data.frame` containing the predictor variables for the test set.
+#' * `offset`: a numeric vector containing the offset values for the test set
+#' (if there is no offset, use a vector of zeros).
+#' * `weights`: a numeric vector containing the observation weights for the test
+#' set (if there are no observation weights, use a vector of ones).
+#' * `y`: a numeric vector containing the response values for the test set.
+#'
 #' @details Arguments `ndraws`, `nclusters`, `nclusters_pred`, and `ndraws_pred`
 #'   are automatically truncated at the number of posterior draws in the
 #'   reference model (which is `1` for `datafit`s). Using less draws or clusters
diff --git a/man/varsel.Rd b/man/varsel.Rd
index a01c2400f..2853866ba 100644
--- a/man/varsel.Rd
+++ b/man/varsel.Rd
@@ -40,10 +40,10 @@ varsel(object, ...)
 minimizer (during a forward search and also during the evaluation part, but
 the latter only if \code{refit_prj} is \code{TRUE}).}
 
-\item{d_test}{For internal use only. A \code{list} providing information about the
-test set which is used for evaluating the predictive performance of the
-submodels as well as of the reference model. If \code{NULL}, the training set is
-used.}
+\item{d_test}{A \code{list} of the structure outlined in section "Argument
+\code{d_test}" below, providing test data for evaluating the predictive
+performance of the submodels as well as of the reference model. If \code{NULL},
+the training data is used.}
 
 \item{method}{The method for the search part. Possible options are \code{"L1"} for
 L1 search and \code{"forward"} for forward search. If \code{NULL}, then internally,
@@ -184,6 +184,19 @@ two candidates. At the last model size of 4, \code{y ~ x1 + x2 + x3} would be
 the only candidate. As another example, to exclude \code{x1} from the search,
 specify \code{search_terms = c("x2", "x3", "x2 + x3")}.
 }
+\section{Argument \code{d_test}}{
+If not \code{NULL}, then \code{d_test} needs to be a \code{list} with the following
+elements:
+\itemize{
+\item \code{data}: a \code{data.frame} containing the predictor variables for the test set.
+\item \code{offset}: a numeric vector containing the offset values for the test set
+(if there is no offset, use a vector of zeros).
+\item \code{weights}: a numeric vector containing the observation weights for the test
+set (if there are no observation weights, use a vector of ones).
+\item \code{y}: a numeric vector containing the response values for the test set.
+}
+}
+
 \examples{
 if (requireNamespace("rstanarm", quietly = TRUE)) {
   # Data:

From 97f5244a3f7e5293e974e85f2b97ad20c7db59f7 Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Wed, 13 Jul 2022 11:25:55 +0200
Subject: [PATCH 10/16] Also in case of K-fold CV: Order the subelements of

`<vsel_object>$summaries` and `<vsel_object>$d_test` as in the original dataset (or test dataset, in case of a non-`NULL` argument `d_test` of `varsel()`, but that doesn't concern K-fold CV).
---
 NEWS.md       |  1 +
 R/cv_varsel.R | 12 ++++++++++++
 2 files changed, 13 insertions(+)

diff --git a/NEWS.md b/NEWS.md
index 3400116b4..27ffa6de5 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -6,6 +6,7 @@
 * At multiple places throughout the package: Improvement of the numerical stability for some link functions, achieved by avoiding unnecessary back-and-forth transformations between latent space and response space. (GitHub: #337, #338)
 * All arguments `seed` and `.seed` are now allowed to be `NA` for not calling `set.seed()` internally at all.
 * Argument `d_test` of `varsel()` is not considered as an internal feature anymore. This was possible after fixing a bug for `d_test` (see below).
+* The order of the observations in the subelements of `<vsel_object>$summaries` and `<vsel_object>$d_test` now corresponds to the order of the observations in the original dataset if `<vsel_object>` was created by a call to `cv_varsel([...], cv_method = "kfold")` (formerly, in that case, the observations in those subelements were ordered by fold). Thereby, the order of the observations in those subelements now always corresponds to the order of the observations in the original dataset, except if `<vsel_object>` was created by a call to `varsel([...], d_test = <non-NULL_d_test_object>)`, in which case the order of the observations in those subelements corresponds to the order of the observations in `<non-NULL_d_test_object>`.
 
 ## Bug fixes
 
diff --git a/R/cv_varsel.R b/R/cv_varsel.R
index f251b4bbf..ec99bb5fb 100644
--- a/R/cv_varsel.R
+++ b/R/cv_varsel.R
@@ -613,7 +613,14 @@ kfold_varsel <- function(refmodel, method, nterms_max, ndraws,
     dim(sub_cv_summaries) <- rev(summ_dim)
   }
   sub <- apply(sub_cv_summaries, 1, rbind2list)
+  idxs_sorted_by_fold <- unlist(lapply(list_cv, function(fold) {
+    fold$d_test$omitted
+  }))
   sub <- lapply(sub, function(summ) {
+    summ$mu <- summ$mu[order(idxs_sorted_by_fold)]
+    summ$lppd <- summ$lppd[order(idxs_sorted_by_fold)]
+
+    # Add weights (see GitHub issue #330 for why this needs to be clarified):
     summ$w <- rep(1, length(summ$mu))
     summ$w <- summ$w / sum(summ$w)
     return(summ)
@@ -632,6 +639,8 @@ kfold_varsel <- function(refmodel, method, nterms_max, ndraws,
       dis = fold$refmodel$dis
     )
   }))
+  ref$mu <- ref$mu[order(idxs_sorted_by_fold)]
+  ref$lppd <- ref$lppd[order(idxs_sorted_by_fold)]
 
   # Combine the K separate test "datasets" (rather "information objects") into a
   # single list:
@@ -640,6 +649,9 @@ kfold_varsel <- function(refmodel, method, nterms_max, ndraws,
          weights = fold$d_test$weights,
          y = fold$d_test$y)
   }))
+  d_cv <- as.list(
+    as.data.frame(d_cv)[order(idxs_sorted_by_fold), , drop = FALSE]
+  )
 
   return(nlist(solution_terms_cv,
                summaries = nlist(sub, ref),

From c0540b83096136cfef899470729cf4495d74fae2 Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Wed, 13 Jul 2022 11:20:16 +0200
Subject: [PATCH 11/16] Update the tests.

---
 tests/testthat/helpers/testers.R | 35 +++++++++++---------------------
 tests/testthat/setup.R           |  6 ++----
 tests/testthat/test_varsel.R     | 19 +++++++++--------
 3 files changed, 24 insertions(+), 36 deletions(-)

diff --git a/tests/testthat/helpers/testers.R b/tests/testthat/helpers/testers.R
index ee324d9be..47f6c971c 100644
--- a/tests/testthat/helpers/testers.R
+++ b/tests/testthat/helpers/testers.R
@@ -1076,7 +1076,9 @@ pp_tester <- function(pp,
 # @param with_cv A single logical value indicating whether `vs` was created by
 #   cv_varsel() (`TRUE`) or not (`FALSE`).
 # @param refmod_expected The expected `vs$refmodel` object.
-# @param dtest_expected The expected `vs$d_test` object.
+# @param dtest_expected If `vs` was created with a non-`NULL` argument `d_test`
+#   (which is only possible for varsel()), then this needs to be the expected
+#   `vs$d_test` object. Otherwise, this needs to be `NULL`.
 # @param solterms_len_expected A single numeric value giving the expected number
 #   of solution terms (not counting the intercept, even for the intercept-only
 #   model).
@@ -1122,8 +1124,6 @@ vsel_tester <- function(
     info_str = ""
 ) {
   # Preparations:
-  dtest_type <- "train"
-  dtest_nms <- c("y", "test_points", "data", "weights", "type", "offset")
   if (with_cv) {
     vsel_nms <- vsel_nms_cv
     vsel_smmrs_sub_nms <- c("lppd", "mu", "w")
@@ -1136,13 +1136,7 @@ vsel_tester <- function(
       valsearch_expected <- TRUE
     }
 
-    dtest_type <- cv_method_expected
-    if (cv_method_expected == "LOO") {
-      # Re-order:
-      dtest_nms <- dtest_nms[c(1, 5, 2, 4, 3, 6)]
-    } else if (cv_method_expected == "kfold") {
-      # Re-order and remove `"data"`:
-      dtest_nms <- dtest_nms[c(1, 4, 2, 6, 5)]
+    if (cv_method_expected == "kfold") {
       # Re-order:
       vsel_smmrs_sub_nms[1:2] <- vsel_smmrs_sub_nms[2:1]
       vsel_smmrs_ref_nms[1:2] <- vsel_smmrs_ref_nms[2:1]
@@ -1263,21 +1257,16 @@ vsel_tester <- function(
   # d_test
   if (is.null(dtest_expected)) {
     expect_type(vs$d_test, "list")
-    expect_named(vs$d_test, dtest_nms, info = info_str)
-    if (identical(cv_method_expected, "kfold")) {
-      expect_identical(vs$d_test$y[order(vs$d_test$test_points)],
-                       vs$refmodel$y, info = info_str)
-      expect_identical(vs$d_test$test_points[order(vs$d_test$test_points)],
-                       seq_len(nobsv), info = info_str)
-      expect_identical(vs$d_test$weights[order(vs$d_test$test_points)],
-                       vs$refmodel$wobs, info = info_str)
-    } else {
-      expect_identical(vs$d_test$y, vs$refmodel$y, info = info_str)
-      expect_identical(vs$d_test$test_points, seq_len(nobsv), info = info_str)
-      expect_identical(vs$d_test$weights, vs$refmodel$wobs, info = info_str)
+    expect_named(vs$d_test, nms_d_test(), info = info_str)
+    dtest_type <- cv_method_expected
+    if (length(dtest_type) == 0) {
+      dtest_type <- "train"
     }
-    expect_null(vs$d_test$data, info = info_str)
     expect_identical(vs$d_test$type, dtest_type, info = info_str)
+    expect_null(vs$d_test$data, info = info_str)
+    expect_identical(vs$d_test$offset, vs$refmodel$offset, info = info_str)
+    expect_identical(vs$d_test$weights, vs$refmodel$wobs, info = info_str)
+    expect_identical(vs$d_test$y, vs$refmodel$y, info = info_str)
   } else {
     expect_identical(vs$d_test, dtest_expected, info = info_str)
   }
diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R
index 6c0318a99..bc403be9a 100644
--- a/tests/testthat/setup.R
+++ b/tests/testthat/setup.R
@@ -1227,8 +1227,6 @@ vsel_nms_cv <- c(
 vsel_nms_pred <- c("summaries", "solution_terms", "kl", "suggested_size",
                    "summary")
 vsel_nms_pred_opt <- c("solution_terms", "suggested_size")
-# Related to `d_test`:
-vsel_nms_dtest <- c("d_test", setdiff(vsel_nms_pred, c("solution_terms", "kl")))
 # Related to `nloo`:
 vsel_nms_cv_nloo <- c("summaries", "pct_solution_terms_cv", "suggested_size",
                       "summary")
@@ -1239,8 +1237,8 @@ vsel_nms_cv_valsearch <- c("validate_search", "summaries",
                            "summary")
 vsel_nms_cv_valsearch_opt <- c("suggested_size")
 # Related to `cvfits`:
-vsel_nms_cv_cvfits <- c("refmodel", "d_test", "summaries",
-                        "pct_solution_terms_cv", "summary", "suggested_size")
+vsel_nms_cv_cvfits <- c("refmodel", "summaries", "pct_solution_terms_cv",
+                        "summary", "suggested_size")
 vsel_nms_cv_cvfits_opt <- c("pct_solution_terms_cv", "suggested_size")
 vsel_smmrs_sub_nms <- vsel_smmrs_ref_nms <- c("mu", "lppd")
 
diff --git a/tests/testthat/test_varsel.R b/tests/testthat/test_varsel.R
index 743aec03e..09e2fddff 100644
--- a/tests/testthat/test_varsel.R
+++ b/tests/testthat/test_varsel.R
@@ -87,12 +87,10 @@ test_that("`d_test` works", {
     fam_crr <- args_vs_i$fam_nm
     refmod_crr <- refmods[[tstsetup_ref]]
     d_test_crr <- list(
-      y = refmod_crr$y,
-      test_points = seq_along(refmod_crr$y),
       data = dat,
+      offset = refmod_crr$offset,
       weights = refmod_crr$wobs,
-      type = "test",
-      offset = refmod_crr$offset
+      y = refmod_crr$y
     )
     vs_repr <- do.call(varsel, c(
       list(object = refmod_crr, d_test = d_test_crr),
@@ -105,7 +103,7 @@ test_that("`d_test` works", {
     vsel_tester(
       vs_repr,
       refmod_expected = refmod_crr,
-      dtest_expected = d_test_crr,
+      dtest_expected = c(list(type = "test"), d_test_crr),
       solterms_len_expected = args_vs_i$nterms_max,
       method_expected = meth_exp_crr,
       nprjdraws_search_expected = args_vs_i$nclusters,
@@ -115,10 +113,13 @@ test_that("`d_test` works", {
         all(grepl("\\+", args_vs_i$search_terms)),
       info_str = tstsetup
     )
-    expect_identical(vs_repr$d_test, d_test_crr, info = tstsetup)
-    expect_equal(vs_repr[setdiff(names(vs_repr), vsel_nms_dtest)],
-                 vss[[tstsetup]][setdiff(names(vss[[tstsetup]]),
-                                         vsel_nms_dtest)],
+    expect_equal(vs_repr[setdiff(names(vs_repr), "d_test")],
+                 vss[[tstsetup]][setdiff(names(vss[[tstsetup]]), "d_test")],
+                 info = tstsetup)
+    expect_equal(vs_repr$d_test[setdiff(names(vs_repr$d_test),
+                                        c("type", "data"))],
+                 vss[[tstsetup]]$d_test[setdiff(names(vss[[tstsetup]]$d_test),
+                                                c("type", "data"))],
                  info = tstsetup)
   }
 })

From a80156986a6e7756fe88e3c3b4d5ecfff9cb2d2f Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Wed, 13 Jul 2022 16:56:12 +0200
Subject: [PATCH 12/16] Improve the structure of `test_varsel.R`.

---
 tests/testthat/test_varsel.R | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/tests/testthat/test_varsel.R b/tests/testthat/test_varsel.R
index 09e2fddff..22ff11fd9 100644
--- a/tests/testthat/test_varsel.R
+++ b/tests/testthat/test_varsel.R
@@ -73,7 +73,11 @@ test_that("`seed` works (and restores the RNG state afterwards)", {
   }
 })
 
-test_that("`d_test` works", {
+## d_test -----------------------------------------------------------------
+
+test_that(paste(
+  "`d_test` set to the training data gives the same results as its default"
+), {
   skip_if_not(run_vs)
   tstsetups <- names(vss)
   ### Alternative with less test setups:
@@ -421,7 +425,7 @@ test_that("for L1 search, `penalty` has an expected effect", {
   }
 })
 
-# search_terms ------------------------------------------------------------
+## search_terms -----------------------------------------------------------
 
 test_that(paste(
   "including all terms in `search_terms` gives the same results as the default",
@@ -564,6 +568,8 @@ test_that("`seed` works (and restores the RNG state afterwards)", {
   }
 })
 
+## nloo -------------------------------------------------------------------
+
 test_that("invalid `nloo` fails", {
   for (tstsetup in names(refmods)) {
     # Use suppressWarnings() because of occasional warnings concerning Pareto k
@@ -647,6 +653,8 @@ test_that("setting `nloo` smaller than the number of observations works", {
   }
 })
 
+## validate_search --------------------------------------------------------
+
 test_that("`validate_search` works", {
   skip_if_not(run_cvvs)
   tstsetups <- grep("\\.default_cvmeth", names(cvvss), value = TRUE)
@@ -742,6 +750,8 @@ test_that("`validate_search` works", {
   expect_true(sum(!suggsize_cond, na.rm = TRUE) <= sum_as_unexpected)
 })
 
+## Arguments specific to K-fold CV ----------------------------------------
+
 test_that("invalid `K` fails", {
   expect_error(cv_varsel(refmods[[1]], cv_method = "kfold", K = 1),
                "^`K` must be at least 2\\.$")

From 759136749be5ff0c748a94c228a679b8efe40726 Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Thu, 14 Jul 2022 09:20:02 +0200
Subject: [PATCH 13/16] Improve the `d_test` test.

---
 tests/testthat/test_varsel.R | 30 ++++++++++++++++++++++++------
 1 file changed, 24 insertions(+), 6 deletions(-)

diff --git a/tests/testthat/test_varsel.R b/tests/testthat/test_varsel.R
index 22ff11fd9..d02c6d1dc 100644
--- a/tests/testthat/test_varsel.R
+++ b/tests/testthat/test_varsel.R
@@ -89,15 +89,33 @@ test_that(paste(
     pkg_crr <- args_vs_i$pkg_nm
     mod_crr <- args_vs_i$mod_nm
     fam_crr <- args_vs_i$fam_nm
-    refmod_crr <- refmods[[tstsetup_ref]]
+    if (!all(refmods[[tstsetup_ref]]$offset == 0)) {
+      offs_crr <- offs_tst
+    } else {
+      offs_crr <- rep(0, nobsv)
+    }
+    if (!all(refmods[[tstsetup_ref]]$wobs == 1)) {
+      wobs_crr <- wobs_tst
+    } else {
+      wobs_crr <- rep(1, nobsv)
+    }
+    formul_fit_crr <- args_fit[[args_vs_i$tstsetup_fit]]$formula
+    dat_crr <- get_dat_formul(formul_crr = formul_fit_crr,
+                              needs_adj = grepl("\\.spclformul", tstsetup))
     d_test_crr <- list(
       data = dat,
-      offset = refmod_crr$offset,
-      weights = refmod_crr$wobs,
-      y = refmod_crr$y
+      offset = offs_crr,
+      weights = wobs_crr,
+      y = dat_crr[[gsub(
+        "\\(|\\)",
+        "",
+        as.character(
+          rm_addresp(rm_cbind(formul_fit_crr))
+        )[2]
+      )]]
     )
     vs_repr <- do.call(varsel, c(
-      list(object = refmod_crr, d_test = d_test_crr),
+      list(object = refmods[[tstsetup_ref]], d_test = d_test_crr),
       excl_nonargs(args_vs_i)
     ))
     meth_exp_crr <- args_vs_i$method
@@ -106,7 +124,7 @@ test_that(paste(
     }
     vsel_tester(
       vs_repr,
-      refmod_expected = refmod_crr,
+      refmod_expected = refmods[[tstsetup_ref]],
       dtest_expected = c(list(type = "test"), d_test_crr),
       solterms_len_expected = args_vs_i$nterms_max,
       method_expected = meth_exp_crr,

From b7ac09e55002e7a6fea51afac527f4ef8889936e Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Wed, 13 Jul 2022 17:35:14 +0200
Subject: [PATCH 14/16] Generate independent test data for the tests.

---
 tests/testthat/setup.R | 43 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 43 insertions(+)

diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R
index bc403be9a..2483126a6 100644
--- a/tests/testthat/setup.R
+++ b/tests/testthat/setup.R
@@ -305,6 +305,49 @@ dat_offs_new <- within(dat, {
   offs_col_new <- seq(-2, 2, length.out = nobsv)
 })
 
+nobsv_indep <- tail(nobsv_tst, 1)
+dis_indep <- runif(1L, 1, 2)
+offs_indep <- rnorm(nobsv_indep)
+wobs_indep <- sample(1:4, nobsv_indep, replace = TRUE)
+idxs_indep <- sample.int(nobsv, size = nobsv_indep, replace = TRUE)
+dat_indep <- lapply(mod_nms, function(mod_nm) {
+  lapply(fam_nms, function(fam_nm) {
+    pred_link <- get(paste0("eta_", mod_nm))
+    pred_link <- pred_link[idxs_indep, , drop = FALSE]
+    if (fam_nm != "brnll" && !mod_nm %in% c("gam", "gamm")) {
+      # For the "brnll" `fam_nm`, offsets are simply not added to have some
+      # scenarios without offsets.
+      # For GAMs, offsets are not added because of rstanarm issue #546 (see
+      # also further below).
+      # For GAMMs, offsets are not added because of rstanarm issue #253 (see
+      # also further below).
+      pred_link <- pred_link + offs_indep
+    }
+    pred_resp <- get(paste0("f_", fam_nm))$linkinv(pred_link)
+    if (fam_nm == "gauss") {
+      return(rnorm(nobsv_indep, mean = pred_resp, sd = dis_indep))
+    } else if (fam_nm == "brnll") {
+      return(rbinom(nobsv_indep, 1, pred_resp))
+    } else if (fam_nm == "binom") {
+      return(rbinom(nobsv_indep, wobs_indep, pred_resp))
+    } else if (fam_nm == "poiss") {
+      return(rpois(nobsv_indep, pred_resp))
+    } else {
+      stop("Unknown `fam_nm`.")
+    }
+  })
+})
+dat_indep <- unlist(dat_indep, recursive = FALSE)
+names(dat_indep) <- paste("y", gsub("\\.", "_", names(dat_indep)), sep = "_")
+dat_indep <- cbind(
+  as.data.frame(dat_indep),
+  dat[idxs_indep,
+      grep("^y_", names(dat), value = TRUE, invert = TRUE),
+      drop = FALSE]
+)
+dat_indep$wobs_col <- wobs_indep
+dat_indep$offs_col <- offs_indep
+
 # Fits --------------------------------------------------------------------
 
 ## Setup ------------------------------------------------------------------

From fc39bba0d195eeadb191218246bc9c3ed7827496 Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Thu, 14 Jul 2022 07:57:11 +0200
Subject: [PATCH 15/16] Use the independent test data for testing equivalence
 between

`<vsel_object>summaries$sub` and `proj_linpred()` as well as between
`<vsel_object>summaries$ref` and `posterior_epred()` / `log_lik()` in case of
a non-`NULL` object `d_test`.
---
 tests/testthat/helpers/testers.R |  12 ++-
 tests/testthat/test_varsel.R     | 146 ++++++++++++++++++++++++++++++-
 2 files changed, 151 insertions(+), 7 deletions(-)

diff --git a/tests/testthat/helpers/testers.R b/tests/testthat/helpers/testers.R
index 47f6c971c..3dacf0329 100644
--- a/tests/testthat/helpers/testers.R
+++ b/tests/testthat/helpers/testers.R
@@ -1281,10 +1281,14 @@ vsel_tester <- function(
       nloo_expected <- nobsv
     }
   }
+  nobsv_summ <- nobsv
+  if (!is.null(dtest_expected)) {
+    nobsv_summ <- nrow(dtest_expected$data)
+  }
   for (j in seq_along(vs$summaries$sub)) {
     expect_named(vs$summaries$sub[[!!j]], vsel_smmrs_sub_nms, info = info_str)
     expect_type(vs$summaries$sub[[!!j]]$mu, "double")
-    expect_length(vs$summaries$sub[[!!j]]$mu, nobsv)
+    expect_length(vs$summaries$sub[[!!j]]$mu, nobsv_summ)
     if (with_cv) {
       expect_identical(sum(!is.na(vs$summaries$sub[[!!j]]$mu)),
                        nloo_expected, info = info_str)
@@ -1292,7 +1296,7 @@ vsel_tester <- function(
       expect_true(all(!is.na(vs$summaries$sub[[!!j]]$mu)), info = info_str)
     }
     expect_type(vs$summaries$sub[[!!j]]$lppd, "double")
-    expect_length(vs$summaries$sub[[!!j]]$lppd, nobsv)
+    expect_length(vs$summaries$sub[[!!j]]$lppd, nobsv_summ)
     if (with_cv) {
       expect_identical(sum(!is.na(vs$summaries$sub[[!!j]]$lppd)),
                        nloo_expected, info = info_str)
@@ -1314,13 +1318,13 @@ vsel_tester <- function(
   }
   expect_type(vs$summaries$ref, "list")
   expect_named(vs$summaries$ref, vsel_smmrs_ref_nms, info = info_str)
-  expect_length(vs$summaries$ref$mu, nobsv)
+  expect_length(vs$summaries$ref$mu, nobsv_summ)
   if (!from_datafit) {
     expect_true(all(!is.na(vs$summaries$ref$mu)), info = info_str)
   } else {
     expect_true(all(is.na(vs$summaries$ref$mu)), info = info_str)
   }
-  expect_length(vs$summaries$ref$lppd, nobsv)
+  expect_length(vs$summaries$ref$lppd, nobsv_summ)
   if (!from_datafit) {
     expect_true(all(!is.na(vs$summaries$ref$lppd)), info = info_str)
   } else {
diff --git a/tests/testthat/test_varsel.R b/tests/testthat/test_varsel.R
index d02c6d1dc..f9a9344c6 100644
--- a/tests/testthat/test_varsel.R
+++ b/tests/testthat/test_varsel.R
@@ -80,9 +80,6 @@ test_that(paste(
 ), {
   skip_if_not(run_vs)
   tstsetups <- names(vss)
-  ### Alternative with less test setups:
-  # tstsetups <- grep("\\.glm\\.", names(vss), value = TRUE)
-  ###
   for (tstsetup in tstsetups) {
     args_vs_i <- args_vs[[tstsetup]]
     tstsetup_ref <- args_vs_i$tstsetup_ref
@@ -146,6 +143,149 @@ test_that(paste(
   }
 })
 
+test_that(paste(
+  "`d_test` set to actual test data gives `<vsel_object>summaries$sub` results",
+  "that can be reproduced by proj_linpred() and `<vsel_object>summaries$ref`",
+  "results that can be reproduced by posterior_epred() and log_lik()"
+), {
+  skip_if_not(run_vs)
+  if (exists(".Random.seed", envir = .GlobalEnv)) {
+    rng_old <- get(".Random.seed", envir = .GlobalEnv)
+  }
+  tstsetups <- names(vss)
+  ### TODO (GAMMs): Currently, the following test setup leads to the error
+  ### ```
+  ### Error in t(as.matrix(b$reTrms$Zt[ii, ])) %*%
+  ### as.matrix(c(as.matrix(ranef[[i]]))) :
+  ###   non-conformable arguments
+  ### ```
+  ### thrown by predict.gamm4(). This needs to be fixed. For now, exclude this
+  ### test setup:
+  tstsetups <- grep("brms\\.gamm\\.binom", tstsetups, value = TRUE,
+                    invert = TRUE)
+  ###
+  for (tstsetup in tstsetups) {
+    args_vs_i <- args_vs[[tstsetup]]
+    tstsetup_ref <- args_vs_i$tstsetup_ref
+    pkg_crr <- args_vs_i$pkg_nm
+    mod_crr <- args_vs_i$mod_nm
+    fam_crr <- args_vs_i$fam_nm
+    if (!all(refmods[[tstsetup_ref]]$offset == 0)) {
+      offs_crr <- offs_indep
+    } else {
+      offs_crr <- rep(0, nobsv_indep)
+    }
+    if (!all(refmods[[tstsetup_ref]]$wobs == 1)) {
+      wobs_crr <- wobs_indep
+    } else {
+      wobs_crr <- rep(1, nobsv_indep)
+    }
+    formul_fit_crr <- args_fit[[args_vs_i$tstsetup_fit]]$formula
+    dat_indep_crr <- get_dat_formul(
+      formul_crr = formul_fit_crr,
+      needs_adj = grepl("\\.spclformul", tstsetup),
+      dat_crr = dat_indep
+    )
+    d_test_crr <- list(
+      data = dat_indep,
+      offset = offs_crr,
+      weights = wobs_crr,
+      y = dat_indep_crr[[gsub(
+        "\\(|\\)",
+        "",
+        as.character(
+          rm_addresp(rm_cbind(formul_fit_crr))
+        )[2]
+      )]]
+    )
+    vs_indep <- do.call(varsel, c(
+      list(object = refmods[[tstsetup_ref]], d_test = d_test_crr),
+      excl_nonargs(args_vs_i)
+    ))
+    meth_exp_crr <- args_vs_i$method
+    if (is.null(meth_exp_crr)) {
+      meth_exp_crr <- ifelse(mod_crr == "glm", "L1", "forward")
+    }
+    vsel_tester(
+      vs_indep,
+      refmod_expected = refmods[[tstsetup_ref]],
+      dtest_expected = c(list(type = "test"), d_test_crr),
+      solterms_len_expected = args_vs_i$nterms_max,
+      method_expected = meth_exp_crr,
+      nprjdraws_search_expected = args_vs_i$nclusters,
+      nprjdraws_eval_expected = args_vs_i$nclusters_pred,
+      search_trms_empty_size =
+        length(args_vs_i$search_terms) &&
+        all(grepl("\\+", args_vs_i$search_terms)),
+      info_str = tstsetup
+    )
+
+    ### Summaries for the submodels -------------------------------------------
+
+    # For getting the correct seed in proj_linpred():
+    set.seed(args_vs_i$seed)
+    p_sel_dummy <- .get_refdist(refmods[[tstsetup_ref]],
+                                nclusters = vs_indep$nprjdraws_search)
+    # As soon as GitHub issues #168 and #211 are fixed, we can use `refit_prj =
+    # FALSE` here:
+    pl_indep <- proj_linpred(vs_indep,
+                             newdata = dat_indep_crr,
+                             offsetnew = d_test_crr$offset,
+                             weightsnew = d_test_crr$weights,
+                             transform = TRUE,
+                             integrated = TRUE,
+                             .seed = NA,
+                             nterms = c(0L, seq_along(vs_indep$solution_terms)),
+                             nclusters = args_vs_i$nclusters_pred,
+                             seed = NA)
+    summ_sub_ch <- lapply(pl_indep, function(pl_indep_k) {
+      names(pl_indep_k)[names(pl_indep_k) == "pred"] <- "mu"
+      names(pl_indep_k)[names(pl_indep_k) == "lpd"] <- "lppd"
+      pl_indep_k$mu <- unname(drop(pl_indep_k$mu))
+      pl_indep_k$lppd <- drop(pl_indep_k$lppd)
+      return(pl_indep_k)
+    })
+    names(summ_sub_ch) <- NULL
+    expect_equal(vs_indep$summaries$sub, summ_sub_ch,
+                 tolerance = .Machine$double.eps, info = tstsetup)
+
+    ### Summaries for the reference model -------------------------------------
+
+    if (pkg_crr == "rstanarm") {
+      mu_new <- rstantools::posterior_epred(refmods[[tstsetup_ref]]$fit,
+                                            newdata = dat_indep,
+                                            offset = d_test_crr$offset)
+      if (grepl("\\.without_wobs", tstsetup)) {
+        lppd_new <- rstantools::log_lik(refmods[[tstsetup_ref]]$fit,
+                                        newdata = dat_indep,
+                                        offset = d_test_crr$offset)
+      } else {
+        ### Currently, rstanarm issue #567 causes an error to be thrown when
+        ### calling log_lik(). Therefore, use the following dummy which
+        ### guarantees test success:
+        lppd_new <- matrix(vs_indep$summaries$ref$lppd,
+                           nrow = nrefdraws, ncol = nobsv_indep, byrow = TRUE)
+        ###
+      }
+    } else if (pkg_crr == "brms") {
+      mu_new <- rstantools::posterior_epred(refmods[[tstsetup_ref]]$fit,
+                                            newdata = dat_indep)
+      lppd_new <- rstantools::log_lik(refmods[[tstsetup_ref]]$fit,
+                                      newdata = dat_indep)
+    }
+    summ_ref_ch <- list(
+      mu = unname(colMeans(mu_new)),
+      lppd = unname(apply(lppd_new, 2, log_sum_exp) - log(nrefdraws))
+    )
+    expect_equal(vs_indep$summaries$ref, summ_ref_ch,
+                 tolerance = 1e2 * .Machine$double.eps, info = tstsetup)
+    lppd_ref_ch2 <- unname(loo::elpd(lppd_new)$pointwise[, "elpd"])
+    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)
+})
+
 ## Regularization ---------------------------------------------------------
 
 # In fact, `regul` is already checked in `test_project.R`, so the `regul` tests

From 3aba93790425df996511a2ea933164dbd1af41fa Mon Sep 17 00:00:00 2001
From: fweber144 <fweber144@protonmail.com>
Date: Thu, 14 Jul 2022 22:24:34 +0200
Subject: [PATCH 16/16] Simplify small parts of the `d_test` tests.

---
 tests/testthat/test_varsel.R | 16 ++--------------
 1 file changed, 2 insertions(+), 14 deletions(-)

diff --git a/tests/testthat/test_varsel.R b/tests/testthat/test_varsel.R
index f9a9344c6..c7a98b8ca 100644
--- a/tests/testthat/test_varsel.R
+++ b/tests/testthat/test_varsel.R
@@ -103,13 +103,7 @@ test_that(paste(
       data = dat,
       offset = offs_crr,
       weights = wobs_crr,
-      y = dat_crr[[gsub(
-        "\\(|\\)",
-        "",
-        as.character(
-          rm_addresp(rm_cbind(formul_fit_crr))
-        )[2]
-      )]]
+      y = dat_crr[[stdize_lhs(formul_fit_crr)$y_nm]]
     )
     vs_repr <- do.call(varsel, c(
       list(object = refmods[[tstsetup_ref]], d_test = d_test_crr),
@@ -190,13 +184,7 @@ test_that(paste(
       data = dat_indep,
       offset = offs_crr,
       weights = wobs_crr,
-      y = dat_indep_crr[[gsub(
-        "\\(|\\)",
-        "",
-        as.character(
-          rm_addresp(rm_cbind(formul_fit_crr))
-        )[2]
-      )]]
+      y = dat_indep_crr[[stdize_lhs(formul_fit_crr)$y_nm]]
     )
     vs_indep <- do.call(varsel, c(
       list(object = refmods[[tstsetup_ref]], d_test = d_test_crr),