Fix full-data dis in K-fold reference models of get_refmodel.stanreg()

NEWS.md

@@ -27,6 +27,7 @@ If you read this from a place other than <https://mc-stan.org/projpred/news/inde
 
 * Fixed bug in the predictions of submodels for new data: Terms with a data-dependent basis, such as `poly()` or `splines::ns()`, used a basis computed from the new data instead of the basis of the training data. This affected `proj_linpred()` and `proj_predict()` with `newdata`, `varsel()` with `d_test`, and K-fold CV, and caused `cv_varsel()` with `validate_search = TRUE` and LOO CV to fail with `poly()` terms. (GitHub: #573)
 * Tests: Adapt to version 2.1-0 of package lme4. (GitHub: #577)
+* Fixed bug in `cv_varsel()` with `cv_method = "kfold"` for **rstanarm** reference models with a dispersion parameter (e.g., the `gaussian()` family). (GitHub: #587)
 * Docs: Explain how custom latent-projection functions can use fold-wise draws in K-fold CV (`?extend_family`, latent-projection vignette). (GitHub: #588)
 * Fixed bias in `cv_varsel()` with `cv_method = "LOO"` and `validate_search = TRUE`: the projections of each fold used the left-out observation, so the submodels, mostly the larger ones, looked too good. The elpd estimates change, and `suggest_size()` can select larger submodels if many Pareto k values are high. (GitHub: #580)
 * Fixed the error `solve(): solution not found` (and possible wrong fits) when a predictor of a GLM submodel is constant over the observations with nonzero weight, e.g., a rare binary predictor in a CV fold. Such a predictor now gets the coefficient 0. (GitHub: #589)

R/refmodel.R

@@ -999,8 +999,12 @@ get_refmodel.stanreg <- function(object, latent = FALSE, dis = NULL, ...) {
     )$fits[, "fit"]
   }
 
+  # Keep the user's `dis` (possibly `NULL`) for the folds. The default set below
+  # uses the full-data `sigma`; with `NULL`, each fold takes `sigma` from its
+  # own refit.
+  dis_user <- dis
   cvrefbuilder <- function(cvfit) {
-    get_refmodel(cvfit, latent = latent, dis = dis,
+    get_refmodel(cvfit, latent = latent, dis = dis_user,
                  called_from_cvrefbuilder = TRUE, ...)
   }
 

tests/testthat/test_refmodel.R

@@ -206,6 +206,28 @@ test_that("get_refmodel() is idempotent", {
   }
 })
 
+test_that(paste(
+  "the `cvrefbuilder` of get_refmodel.stanreg() takes the default `dis` from",
+  "the fold-wise fit and passes a user-supplied `dis` unchanged"
+), {
+  tstsetups <- grep("^rstanarm\\..*\\.gauss\\.", names(fits), value = TRUE)
+  skip_if_not(length(tstsetups) >= 2)
+  # The second fit stands in for a fold refit:
+  fit_full <- fits[[tstsetups[1]]]
+  fit_fold <- fits[[tstsetups[2]]]
+  sigma_full <- data.frame(fit_full)[, "sigma"]
+  sigma_fold <- data.frame(fit_fold)[, "sigma"]
+  skip_if(isTRUE(all.equal(sigma_full, sigma_fold)))
+
+  refmod <- get_refmodel(fit_full)
+  expect_identical(refmod$dis, sigma_full)
+  expect_identical(refmod$cvrefbuilder(fit_fold)$dis, sigma_fold)
+
+  dis_user <- rep(0.7, length(sigma_full))
+  refmod_user <- get_refmodel(fit_full, dis = dis_user)
+  expect_identical(refmod_user$cvrefbuilder(fit_fold)$dis, dis_user)
+})
+
 # predict.refmodel() ------------------------------------------------------
 
 test_that("invalid `type` fails", {