From 19b6d0dad8cc00eba306bd14adc13fac47fb4c74 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Sat, 24 Oct 2026 14:23:29 +0300
Subject: [PATCH 1/6] Leave the left-out observation out of the fold
 projections in validated LOO

---
 R/cv_varsel.R                             | 13 ++++---
 tests/testthat/test_loo_fold_projection.R | 44 +++++++++++++++++++++++
 2 files changed, 53 insertions(+), 4 deletions(-)
 create mode 100644 tests/testthat/test_loo_fold_projection.R

diff --git a/R/cv_varsel.R b/R/cv_varsel.R
index fbb81b113..ef225c28d 100644
--- a/R/cv_varsel.R
+++ b/R/cv_varsel.R
@@ -1034,12 +1034,16 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws,
       # Run the search with the reweighted clusters (or thinned draws) (so the
       # *reweighted* fitted response values from the reference model act as
       # artifical response values in the projection (or L1-penalized
-      # projection)):
+      # projection)).  As in K-fold CV, the left-out observation is not used
+      # in the projections of its fold; it is used only for the evaluation
+      # (via `refmodel_fulldata`):
+      refmodel_fold <- refmodel
+      refmodel_fold$wobs[i] <- 0
       if (!search_out_rks_was_null) {
         search_path <- list(predictor_ranking = search_out_rks[[run_index]])
       } else {
         search_path <- .select(
-          refmodel = refmodel, ndraws = ndraws, nclusters = nclusters,
+          refmodel = refmodel_fold, ndraws = ndraws, nclusters = nclusters,
           reweighting_args = list(cl_ref = cl_sel, wdraws_ref = exp(lw[, i])),
           method = method, nterms_max = nterms_max, penalty = penalty,
           verbose = verbose_obs, verbose_line_length = 3,
@@ -1051,9 +1055,10 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws,
       # Run the performance evaluation for the submodels along the predictor
       # ranking:
       perf_eval_out <- perf_eval(
-        search_path = search_path, refmodel = refmodel, refit_prj = refit_prj,
-        ndraws = ndraws_pred, nclusters = nclusters_pred,
+        search_path = search_path, refmodel = refmodel_fold,
+        refit_prj = refit_prj, ndraws = ndraws_pred, nclusters = nclusters_pred,
         reweighting_args = list(cl_ref = cl_pred, wdraws_ref = exp(lw[, i])),
+        refmodel_fulldata = refmodel,
         indices_test = i, verbose = verbose_obs, verbose_line_length = 3,
         verbose_txt_add = vtxt_obs_i, ...
       )
diff --git a/tests/testthat/test_loo_fold_projection.R b/tests/testthat/test_loo_fold_projection.R
new file mode 100644
index 000000000..40b3926de
--- /dev/null
+++ b/tests/testthat/test_loo_fold_projection.R
@@ -0,0 +1,44 @@
+context("validated LOO: fold projections")
+
+# A Gaussian reference model with fixed draws and more predictors than
+# observations, so that submodels with at least n - 1 terms can interpolate
+# the reference predictions at the observations used in the projection.
+cre_ref_interp <- function(n = 8, p = 9, S = 20) {
+  set.seed(5309)
+  x <- matrix(rnorm(n * p), n, p, dimnames = list(NULL, paste0("x", seq_len(p))))
+  dat <- data.frame(y = x[, 1] + rnorm(n), x)
+  beta <- matrix(rnorm((p + 1) * S, sd = 0.3), p + 1, S)
+  beta[2, ] <- beta[2, ] + 1
+  ref_predfun <- function(fit, newdata = NULL) {
+    if (is.null(newdata)) newdata <- dat
+    cbind(1, as.matrix(newdata[, colnames(x), drop = FALSE])) %*% fit$beta
+  }
+  ref <- init_refmodel(
+    object = list(beta = beta), data = dat,
+    formula = reformulate(colnames(x), response = "y"), family = gaussian(),
+    ref_predfun = ref_predfun, dis = rep(1, S)
+  )
+  list(ref = ref, n = n, p = p, S = S)
+}
+
+test_that(paste(
+  "cv_varsel() with LOO and `validate_search = TRUE` does not use the",
+  "left-out observation in the projections of its fold"
+), {
+  skip_on_cran()
+  r <- cre_ref_interp()
+  cvvs <- suppressWarnings(cv_varsel(
+    r$ref, cv_method = "LOO", validate_search = TRUE, method = "forward",
+    nterms_max = r$p, nclusters = NULL, ndraws = 5, ndraws_pred = r$S,
+    refit_prj = TRUE, seed = 1, verbose = FALSE
+  ))
+  lppd <- sapply(cvvs$summaries$sub, `[[`, "lppd")
+  expect_equal(dim(lppd), c(r$n, r$p + 1))
+  # Columns `r$n` to `r$p + 1` are the sizes n - 1 to p. If the left-out
+  # observation were used in the fold projections, these submodels would
+  # interpolate the reference predictions at it, and their fold lppd would be
+  # the same for all these sizes.
+  sizes_interp <- seq(r$n, r$p + 1L)
+  spread <- apply(lppd[, sizes_interp, drop = FALSE], 1, function(z) diff(range(z)))
+  expect_gt(max(spread), 1e-2)
+})

From 9ee8f79e201f501fec84aee996bb5284eb3dd6e0 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Sun, 25 Oct 2026 20:25:05 +0200
Subject: [PATCH 2/6] Apply the observation weights in the augmented-data
 projection

---
 NEWS.md                                   |  2 +
 R/projfun.R                               | 11 +++-
 tests/testthat/test_loo_fold_projection.R | 68 +++++++++++++++++++++++
 3 files changed, 80 insertions(+), 1 deletion(-)

diff --git a/NEWS.md b/NEWS.md
index 94ffe21cd..8e667b0e6 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -11,6 +11,8 @@ If you read this from a place other than <https://mc-stan.org/projpred/news/inde
 
 ## Bug fixes and maintenance
 
+* 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: #<PR>)
+
 
 # projpred 2.11.0
 
diff --git a/R/projfun.R b/R/projfun.R
index aa79dbf31..00f852571 100644
--- a/R/projfun.R
+++ b/R/projfun.R
@@ -26,12 +26,21 @@ proj_to_submodl <- function(predictor_terms, p_ref, refmodel,
     verb_out("  Projecting onto `~ ", utils::tail(rhs_chr, 1), "`")
   }
 
+  projpred_ws_aug <- p_ref$mu
+  if (refmodel$family$for_augdat && any(refmodel$wobs != 1)) {
+    # The augmented-data submodel fitters use the probabilities of the response
+    # categories as weights and ignore `weights`, so apply the observation
+    # weights to these probabilities. The rows of the augmented-rows matrix are
+    # the observations within the response categories:
+    projpred_ws_aug <- projpred_ws_aug *
+      rep(refmodel$wobs, times = nrow(projpred_ws_aug) / length(refmodel$wobs))
+  }
   args_divmin <- list(formula = fml_divmin,
                       data = subset$data,
                       family = refmodel$family,
                       weights = refmodel$wobs,
                       projpred_var = p_ref$var,
-                      projpred_ws_aug = p_ref$mu,
+                      projpred_ws_aug = projpred_ws_aug,
                       verbose_divmin = verbose >= 3L)
   if (!is.null(search_control)) {
     args_divmin <- c(args_divmin, search_control)
diff --git a/tests/testthat/test_loo_fold_projection.R b/tests/testthat/test_loo_fold_projection.R
index 40b3926de..2ca90261a 100644
--- a/tests/testthat/test_loo_fold_projection.R
+++ b/tests/testthat/test_loo_fold_projection.R
@@ -42,3 +42,71 @@ test_that(paste(
   spread <- apply(lppd[, sizes_interp, drop = FALSE], 1, function(z) diff(range(z)))
   expect_gt(max(spread), 1e-2)
 })
+
+test_that(paste(
+  "an observation with weight 0 does not affect the projection (all model",
+  "types, families and projection types of the test setups)"
+), {
+  skip_if_not(run_prj)
+  # Changes the reference target of observation `i` in a way that changes the
+  # projection if the observation has a positive weight:
+  change_target <- function(p_ref, refm, i) {
+    fam <- refm$family
+    if (fam$for_augdat) {
+      # Reverse the order of the category probabilities of observation `i`:
+      n <- refm$nobs
+      idx <- i + n * (seq_len(nrow(p_ref$mu) / n) - 1L)
+      p_ref$mu[idx, ] <- p_ref$mu[rev(idx), ]
+    } else if (fam$for_latent || fam$family %in% c("gaussian", "Student_t")) {
+      p_ref$mu[i, ] <- p_ref$mu[i, ] + 5
+    } else if (fam$family == "binomial") {
+      p_ref$mu[i, ] <- 1 - p_ref$mu[i, ]
+    } else {
+      p_ref$mu[i, ] <- 3 * p_ref$mu[i, ] + 1
+    }
+    p_ref
+  }
+  # One linear term and, depending on the model type, one group-level term
+  # and one smooth term, so that each submodel fitter is used:
+  terms_case <- function(refm) {
+    trms <- setdiff(split_formula(refm$formula, data = refm$fetch_data()), "1")
+    trms <- trms[!grepl("\\+", trms)]
+    is_grp <- grepl("\\|", trms)
+    is_smooth <- grepl("^s\\(", trms)
+    c(trms[!is_grp & !is_smooth][1], trms[is_grp][1], trms[is_smooth][1])[
+      c(TRUE, any(is_grp), any(is_smooth))
+    ]
+  }
+  # One test setup per combination of model type, family, projection type and
+  # observation weights:
+  case_key <- function(tstsetup) {
+    refm <- refmods[[tstsetup]]
+    mod_nm <- sub("^[^.]+\\.([^.]+)\\..*$", "\\1", tstsetup)
+    paste(mod_nm, refm$family$family, refm$family$for_latent,
+          refm$family$for_augdat, any(refm$wobs != 1))
+  }
+  keys <- vapply(names(refmods), case_key, character(1))
+  tstsetups <- names(refmods)[!duplicated(keys)]
+  expect_gt(length(tstsetups), 5)
+  for (tstsetup in tstsetups) {
+    refm <- refmods[[tstsetup]]
+    trms <- terms_case(refm)
+    p_ref <- get_refdist(refm, ndraws = 3)
+    p_ref_chg <- change_target(p_ref, refm, i = 1)
+    refm_w0 <- refm
+    refm_w0$wobs[1] <- 0
+    # The predictions of the submodel fits for the training data:
+    prd <- function(p, r) {
+      outdmin <- suppressWarnings(suppressMessages(proj_to_submodl(
+        predictor_terms = trms, p_ref = p, refmodel = r, verbose = FALSE
+      )$outdmin))
+      unclass(refm$mu_fun(outdmin, offset = refm$offset))
+    }
+    # With weight 0, the changed target of observation 1 does not matter:
+    expect_equal(prd(p_ref_chg, refm_w0), prd(p_ref, refm_w0),
+                 tolerance = 1e-5, info = tstsetup)
+    # With its original weight, it does:
+    expect_gt(max(abs(prd(p_ref_chg, refm) - prd(p_ref, refm))), 1e-4,
+              label = tstsetup)
+  }
+})

From 4496f4f30b2f726539bc74c8d55a37436a007639 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Sun, 25 Oct 2026 23:17:57 +0200
Subject: [PATCH 3/6] Remove zero-weight observations before the multilevel and
 additive submodel fits

---
 R/divergence_minimizers.R                 | 32 +++++++++++++++++------
 tests/testthat/test_loo_fold_projection.R | 23 +++++++++++++---
 2 files changed, 43 insertions(+), 12 deletions(-)

diff --git a/R/divergence_minimizers.R b/R/divergence_minimizers.R
index 36345eda8..6b46b19ae 100644
--- a/R/divergence_minimizers.R
+++ b/R/divergence_minimizers.R
@@ -78,6 +78,23 @@ divmin <- function(
     )
   }
 
+  dot_args <- list(...)
+  # Observations with weight 0 (in the LOO CV with validated search, the
+  # left-out observation of a fold): remove them from the data for the
+  # multilevel and additive submodel fitters. Some of these fitters (for example
+  # lme4::lmer()) do not ignore the response of an observation with weight 0,
+  # and others (mgcv::gam()) use its predictor values for the smooth bases. The
+  # GLM fitters handle weight 0 exactly:
+  if ((has_grp || has_add) &&
+      is.data.frame(dot_args$data) && length(dot_args$weights) > 0 &&
+      any(dot_args$weights == 0) &&
+      length(dot_args$weights) == nrow(dot_args$data)) {
+    keep <- dot_args$weights != 0
+    dot_args$data <- dot_args$data[keep, , drop = FALSE]
+    dot_args$weights <- dot_args$weights[keep]
+    projpred_var <- projpred_var[keep, , drop = FALSE]
+  }
+
   if (is.null(getOption("projpred.parallel_proj_trigger")) &&
       !is.null(getOption("projpred.prll_prj_trigger"))) {
     warning(
@@ -104,13 +121,13 @@ divmin <- function(
         on.exit(utils::setTxtProgressBar(pb, s))
       }
       mssgs_warns_capt <- capt_mssgs_warns(
-        soutdmin <- sdivmin(
-          formula = formulas[[s]],
-          projpred_var = projpred_var[, s, drop = FALSE],
-          projpred_formula_no_random = projpred_formulas_no_random[[s]],
-          projpred_random = projpred_random,
-          ...
-        )
+        soutdmin <- do.call(sdivmin, c(
+          list(formula = formulas[[s]],
+               projpred_var = projpred_var[, s, drop = FALSE],
+               projpred_formula_no_random = projpred_formulas_no_random[[s]],
+               projpred_random = projpred_random),
+          dot_args
+        ))
       )
       return(nlist(soutdmin, mssgs_warns_capt))
     })
@@ -127,7 +144,6 @@ divmin <- function(
     } else {
       progressor_obj <- NULL
     }
-    dot_args <- list(...)
     `%do_projpred%` <- foreach::`%dopar%`
     outdmin <- foreach::foreach(
       formula_s = formulas,
diff --git a/tests/testthat/test_loo_fold_projection.R b/tests/testthat/test_loo_fold_projection.R
index 2ca90261a..d776d2087 100644
--- a/tests/testthat/test_loo_fold_projection.R
+++ b/tests/testthat/test_loo_fold_projection.R
@@ -85,8 +85,15 @@ test_that(paste(
     paste(mod_nm, refm$family$family, refm$family$for_latent,
           refm$family$for_augdat, any(refm$wobs != 1))
   }
+  # The augmented-data projection supports only these families:
+  is_projectable <- function(tstsetup) {
+    fam <- refmods[[tstsetup]]$family
+    !fam$for_augdat ||
+      fam$family %in% c("cumulative", "cumulative_rstanarm", "categorical")
+  }
   keys <- vapply(names(refmods), case_key, character(1))
   tstsetups <- names(refmods)[!duplicated(keys)]
+  tstsetups <- Filter(is_projectable, tstsetups)
   expect_gt(length(tstsetups), 5)
   for (tstsetup in tstsetups) {
     refm <- refmods[[tstsetup]]
@@ -102,11 +109,19 @@ test_that(paste(
       )$outdmin))
       unclass(refm$mu_fun(outdmin, offset = refm$offset))
     }
+    # An error in one setup is a failure for that setup only:
+    res <- tryCatch(
+      list(w0_chg = prd(p_ref_chg, refm_w0), w0 = prd(p_ref, refm_w0),
+           w1_chg = prd(p_ref_chg, refm), w1 = prd(p_ref, refm)),
+      error = function(e) conditionMessage(e)
+    )
+    if (is.character(res)) {
+      expect_true(FALSE, label = paste(tstsetup, "error:", res))
+      next
+    }
     # With weight 0, the changed target of observation 1 does not matter:
-    expect_equal(prd(p_ref_chg, refm_w0), prd(p_ref, refm_w0),
-                 tolerance = 1e-5, info = tstsetup)
+    expect_equal(res$w0_chg, res$w0, tolerance = 1e-5, info = tstsetup)
     # With its original weight, it does:
-    expect_gt(max(abs(prd(p_ref_chg, refm) - prd(p_ref, refm))), 1e-4,
-              label = tstsetup)
+    expect_gt(max(abs(res$w1_chg - res$w1)), 1e-4, label = tstsetup)
   }
 })

From 362fcc21f07872f1c1914b3f3cad7c6f27ad9ec8 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Mon, 26 Oct 2026 09:46:42 +0200
Subject: [PATCH 4/6] Add the PR number to the NEWS entry

---
 NEWS.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/NEWS.md b/NEWS.md
index 8e667b0e6..f3a97c0f7 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -11,7 +11,7 @@ If you read this from a place other than <https://mc-stan.org/projpred/news/inde
 
 ## Bug fixes and maintenance
 
-* 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: #<PR>)
+* 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)
 
 
 # projpred 2.11.0

From 3d3ca0c6d1cbf635673651e52f22648ba2332244 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Mon, 26 Oct 2026 19:26:31 +0200
Subject: [PATCH 5/6] Remove zero-weight observations before the augmented-data
 submodel fits, and predict lost factor levels as the reference level

---
 R/divergence_minimizers.R                 | 137 ++++++++++++++++++-
 R/refmodel.R                              |   3 +
 tests/testthat/test_loo_fold_projection.R | 152 +++++++++++++++++++++-
 3 files changed, 289 insertions(+), 3 deletions(-)

diff --git a/R/divergence_minimizers.R b/R/divergence_minimizers.R
index 6b46b19ae..1f175084a 100644
--- a/R/divergence_minimizers.R
+++ b/R/divergence_minimizers.R
@@ -84,13 +84,17 @@ divmin <- function(
   # multilevel and additive submodel fitters. Some of these fitters (for example
   # lme4::lmer()) do not ignore the response of an observation with weight 0,
   # and others (mgcv::gam()) use its predictor values for the smooth bases. The
-  # GLM fitters handle weight 0 exactly:
+  # GLM fitters handle weight 0 exactly. If `weights` and `data` differ in
+  # length (not in projpred's own calls), nothing is removed.
+  lost_lvls <- NULL
   if ((has_grp || has_add) &&
       is.data.frame(dot_args$data) && length(dot_args$weights) > 0 &&
       any(dot_args$weights == 0) &&
       length(dot_args$weights) == nrow(dot_args$data)) {
     keep <- dot_args$weights != 0
-    dot_args$data <- dot_args$data[keep, , drop = FALSE]
+    lost_lvls <- lost_levels(dot_args$data, keep = keep, formula = formula)
+    dot_args$data <- drop_lost_levels(dot_args$data[keep, , drop = FALSE],
+                                      lost_lvls)
     dot_args$weights <- dot_args$weights[keep]
     projpred_var <- projpred_var[keep, , drop = FALSE]
   }
@@ -182,9 +186,120 @@ divmin <- function(
   })
   warn_prj_drawwise(mssgs_warns_capts, throw_warn = throw_warn_sdivmin)
   check_conv(outdmin, lengths(mssgs_warns_capts), do_check = do_check_conv)
+  attr(outdmin, "projpred_lost_lvls") <- lost_lvls
   return(outdmin)
 }
 
+# Levels of factor-like fixed-effect predictors that occur only in the removed
+# rows (`!keep`), with the reference level that the fitter uses after it drops
+# these unused levels. Factor-like: factor, character and logical variables,
+# and numeric variables that the fixed-effects part of `formula` uses only
+# inside factor(), as.factor(), ordered() or as.ordered(). recode_lost_levels()
+# predicts these levels as the reference level. With treatment contrasts, this
+# is the prediction for a level without a coefficient. The GLM fitter keeps the
+# rows with weight 0 and gives the same prediction, because its ridge penalty
+# sets such a coefficient to almost zero. With other contrasts, it is the
+# prediction for the first remaining level. Grouping factors are excluded:
+# repair_re() treats a lost group as a new group. Returns `NULL` if no level is
+# lost.
+lost_levels <- function(data, keep, formula) {
+  fml_fixed <- reformulas::nobars(formula)
+  rhs_fixed <- fml_fixed[[length(fml_fixed)]]
+  grp_nms <- unlist(lapply(reformulas::findbars(formula), function(bar) {
+    all.vars(bar[[3]])
+  }))
+  # Variables inside factor() and similar calls, and variables elsewhere:
+  fac_fns <- c("factor", "as.factor", "ordered", "as.ordered")
+  find_vars <- function(e) {
+    if (is.call(e)) {
+      fn <- sub("^.*::", "", paste(deparse(e[[1]]), collapse = ""))
+      if (fn %in% fac_fns) {
+        return(list(inside = all.vars(e), outside = character()))
+      }
+      res <- lapply(as.list(e)[-1], find_vars)
+      return(list(inside = unlist(lapply(res, "[[", "inside")),
+                  outside = unlist(lapply(res, "[[", "outside"))))
+    }
+    if (is.name(e)) {
+      return(list(inside = character(), outside = as.character(e)))
+    }
+    return(list(inside = character(), outside = character()))
+  }
+  vars_rhs <- find_vars(rhs_fixed)
+  fac_only <- setdiff(vars_rhs$inside, vars_rhs$outside)
+  x_nms <- setdiff(intersect(all.vars(rhs_fixed), names(data)), grp_nms)
+  out <- lapply(setNames(nm = x_nms), function(x_nm) {
+    x <- data[[x_nm]]
+    if (!(is.factor(x) || is.character(x) || is.logical(x) ||
+          (is.numeric(x) && x_nm %in% fac_only))) {
+      return(NULL)
+    }
+    x_keep <- x[keep]
+    vals_keep <- unique(x_keep[!is.na(x_keep)])
+    vals_rm <- unique(x[!keep & !is.na(x)])
+    if (is.factor(x)) {
+      vals_keep <- as.character(vals_keep)
+      vals_rm <- as.character(vals_rm)
+    }
+    lost <- setdiff(vals_rm, vals_keep)
+    if (!length(lost) || !length(vals_keep)) {
+      return(NULL)
+    }
+    if (is.factor(x)) {
+      ref <- levels(droplevels(x_keep))[1]
+    } else {
+      ref <- sort(vals_keep)[1]
+    }
+    list(lost = lost, ref = ref)
+  })
+  out <- out[!vapply(out, is.null, logical(1))]
+  if (!length(out)) {
+    return(NULL)
+  }
+  return(out)
+}
+
+# Sets the levels that lost_levels() recorded (attribute `projpred_lost_lvls` of
+# the fits) to the reference level in `newdata`.
+recode_lost_levels <- function(newdata, lost_lvls) {
+  if (is.null(lost_lvls) || !is.data.frame(newdata)) {
+    return(newdata)
+  }
+  for (x_nm in intersect(names(lost_lvls), names(newdata))) {
+    x <- newdata[[x_nm]]
+    lost <- lost_lvls[[x_nm]]$lost
+    ref <- lost_lvls[[x_nm]]$ref
+    if (is.factor(x)) {
+      # Drop the lost levels, because the predict() method of
+      # mclogit::mblogit() uses all levels of `newdata`:
+      lvls <- union(setdiff(levels(x), lost), ref)
+      x_chr <- as.character(x)
+      x_chr[x_chr %in% lost] <- ref
+      x <- factor(x_chr, levels = lvls, ordered = is.ordered(x))
+    } else if (typeof(x) == typeof(ref) ||
+               (is.numeric(x) && is.numeric(ref))) {
+      x[x %in% lost] <- ref
+    }
+    newdata[[x_nm]] <- x
+  }
+  return(newdata)
+}
+
+# Drops the levels that lost_levels() recorded from the factors in `data` (the
+# data without the removed rows), so that fitters which keep unused levels in
+# the model matrix do not get a column of zeros.
+drop_lost_levels <- function(data, lost_lvls) {
+  for (x_nm in intersect(names(lost_lvls), names(data))) {
+    x <- data[[x_nm]]
+    if (is.factor(x)) {
+      data[[x_nm]] <- factor(x, levels = setdiff(levels(x),
+                                                 lost_lvls[[x_nm]]$lost),
+                             ordered = is.ordered(x))
+    }
+  }
+  return(data)
+}
+
 # Use projpred's own implementation to fit non-multilevel non-additive
 # submodels:
 fit_glm_ridge_callback <- function(formula, data,
@@ -710,6 +825,23 @@ divmin_augdat <- function(
   projpred_ws_aug <- unclass(projpred_ws_aug)
   attr(projpred_ws_aug, "ndiscrete") <- NULL
 
+  # Remove the rows of the observations with weight 0 (in validated LOO CV, the
+  # left-out observation) in all response categories before the submodel
+  # fitters. mclogit::mblogit() uses the predictor values of rows with weight
+  # 0, and MASS::polr() gives no coefficient for a factor level that only these
+  # rows have. The augmented data hold all observations for the first category,
+  # then all for the second, and so on. If the sizes do not match (not in
+  # projpred's own calls), nothing is removed.
+  lost_lvls <- NULL
+  if (is.data.frame(data) && length(weights) > 0 &&
+      any(weights == 0) && nrow(data) %% length(weights) == 0 &&
+      nrow(data) == nrow(projpred_ws_aug)) {
+    keep <- rep(weights != 0, times = nrow(data) / length(weights))
+    lost_lvls <- lost_levels(data, keep = keep, formula = formula)
+    data <- drop_lost_levels(data[keep, , drop = FALSE], lost_lvls)
+    projpred_ws_aug <- projpred_ws_aug[keep, , drop = FALSE]
+  }
+
   if (is.null(getOption("projpred.parallel_proj_trigger")) &&
       !is.null(getOption("projpred.prll_prj_trigger"))) {
     warning(
@@ -815,6 +947,7 @@ divmin_augdat <- function(
   })
   warn_prj_drawwise(mssgs_warns_capts, throw_warn = throw_warn_sdivmin)
   check_conv(outdmin, lengths(mssgs_warns_capts), do_check = do_check_conv)
+  attr(outdmin, "projpred_lost_lvls") <- lost_lvls
   return(outdmin)
 }
 
diff --git a/R/refmodel.R b/R/refmodel.R
index 51e4a55f5..648d58ffe 100644
--- a/R/refmodel.R
+++ b/R/refmodel.R
@@ -1404,6 +1404,9 @@ init_refmodel <- function(object, data, formula, family, ref_predfun = NULL,
     } else {
       stopifnot(length(offset) %in% c(1L, nrow(newdata)))
     }
+    # Predict factor levels that only the removed rows with weight 0 had as the
+    # reference level (see lost_levels()):
+    newdata <- recode_lost_levels(newdata, attr(fits, "projpred_lost_lvls"))
     pred_sub <- proj_predfun(fits, newdata = newdata)
     if (family$family %in% fams_neg_linpred()) {
       pred_sub <- pred_sub - offset
diff --git a/tests/testthat/test_loo_fold_projection.R b/tests/testthat/test_loo_fold_projection.R
index d776d2087..ecc3fbfca 100644
--- a/tests/testthat/test_loo_fold_projection.R
+++ b/tests/testthat/test_loo_fold_projection.R
@@ -102,6 +102,23 @@ test_that(paste(
     p_ref_chg <- change_target(p_ref, refm, i = 1)
     refm_w0 <- refm
     refm_w0$wobs[1] <- 0
+    # The same with other predictor values of observation 1 (in the training
+    # data of the projection only):
+    refm_w0_x <- refm_w0
+    refm_w0_x$fetch_data <- function(obs = NULL) {
+      dat <- refm_w0$fetch_data(obs = obs)
+      if (is.null(obs)) {
+        trms_ref <- delete.response(terms(refm$formula))
+        offs_vars <- unlist(lapply(attr(trms_ref, "offset"), function(i) {
+          all.vars(attr(trms_ref, "variables")[[i + 1L]])
+        }))
+        x_nms <- setdiff(intersect(all.vars(trms_ref), names(dat)), offs_vars)
+        x_nms <- x_nms[vapply(dat[x_nms], is.numeric, logical(1))]
+        dat[1, x_nms] <- dat[1, x_nms] + 3
+      }
+      dat
+    }
+    is_mlvl_add <- any(grepl("\\||^s\\(", trms))
     # The predictions of the submodel fits for the training data:
     prd <- function(p, r) {
       outdmin <- suppressWarnings(suppressMessages(proj_to_submodl(
@@ -112,7 +129,8 @@ test_that(paste(
     # An error in one setup is a failure for that setup only:
     res <- tryCatch(
       list(w0_chg = prd(p_ref_chg, refm_w0), w0 = prd(p_ref, refm_w0),
-           w1_chg = prd(p_ref_chg, refm), w1 = prd(p_ref, refm)),
+           w1_chg = prd(p_ref_chg, refm), w1 = prd(p_ref, refm),
+           w0_x = if (is_mlvl_add) prd(p_ref, refm_w0_x)),
       error = function(e) conditionMessage(e)
     )
     if (is.character(res)) {
@@ -121,7 +139,139 @@ test_that(paste(
     }
     # With weight 0, the changed target of observation 1 does not matter:
     expect_equal(res$w0_chg, res$w0, tolerance = 1e-5, info = tstsetup)
+    # The multilevel and additive fitters get the data without the rows with
+    # weight 0, so the predictor values of observation 1 must not matter
+    # either:
+    if (is_mlvl_add) {
+      expect_equal(res$w0_x, res$w0, tolerance = 1e-5, info = tstsetup)
+    }
     # With its original weight, it does:
     expect_gt(max(abs(res$w1_chg - res$w1)), 1e-4, label = tstsetup)
   }
 })
+
+test_that(paste(
+  "the submodel fits predict a factor level that only an observation with",
+  "weight 0 has as the reference level"
+), {
+  skip_on_cran()
+  set.seed(2)
+  n <- 40
+  dat <- data.frame(x = rnorm(n), g = factor(rep(paste0("g", 1:5), 8)),
+                    xc = factor(c(rep(c("a", "b"), length.out = n - 1), "r")))
+  for (s in 1:2) dat[[paste0(".y.", s)]] <- dat$x + rnorm(n)
+  # Observation n is the only one with level "r", and it has weight 0:
+  wobs <- c(rep(1, n - 1), 0)
+  dat_ref <- dat
+  dat_ref$xc[n] <- "a"
+  dat_ref$xc <- droplevels(dat_ref$xc)
+  divmin_w0 <- function(fml) {
+    suppressWarnings(divmin(
+      formula = fml, data = dat, family = extend_family(gaussian()),
+      weights = wobs, projpred_var = matrix(NA, n, 2),
+      projpred_ws_aug = matrix(NA, n, 2), verbose_divmin = FALSE
+    ))
+  }
+  fmls <- list(glmm = cbind(.y.1, .y.2) ~ x + xc + (1 | g),
+               gam = cbind(.y.1, .y.2) ~ s(x) + xc)
+  for (nm in names(fmls)) {
+    fits <- divmin_w0(fmls[[nm]])
+    expect_identical(attr(fits, "projpred_lost_lvls"),
+                     list(xc = list(lost = "r", ref = "a")), info = nm)
+    prd <- subprd(fits, recode_lost_levels(
+      dat, attr(fits, "projpred_lost_lvls")
+    ))
+    expect_identical(prd, subprd(fits, dat_ref), info = nm)
+  }
+  # Augmented-data projection (rows: the observations within the response
+  # categories):
+  ncat <- 3
+  lv <- paste0("c", seq_len(ncat))
+  dat_aug <- dat[rep(seq_len(n), ncat), ]
+  ws_aug <- matrix(runif(n * ncat * 2), ncol = 2)
+  for (fam_nm in c("cumulative", "categorical")) {
+    if (fam_nm == "cumulative") {
+      skip_if_not_installed("ordinal")
+    } else {
+      skip_if_not_installed("mclogit")
+    }
+    dat_aug$y <- factor(rep(lv, each = n), levels = lv,
+                        ordered = fam_nm == "cumulative")
+    fam <- structure(list(family = fam_nm, link = "logit"), class = "family")
+    # With and without group-level terms (clmm() or polr(), mblogit() or
+    # multinom()):
+    for (fml in list(y ~ x + xc + (1 | g), y ~ x + xc)) {
+      info_txt <- paste(fam_nm, deparse(fml))
+      fits <- suppressWarnings(divmin_augdat(
+        formula = fml, data = dat_aug, family = fam,
+        weights = wobs, projpred_var = matrix(NA, n * ncat, 2),
+        projpred_ws_aug = ws_aug, verbose_divmin = FALSE
+      ))
+      expect_identical(attr(fits, "projpred_lost_lvls"),
+                       list(xc = list(lost = "r", ref = "a")), info = info_txt)
+      prd <- subprd_augdat(fits, recode_lost_levels(
+        dat, attr(fits, "projpred_lost_lvls")
+      ))
+      expect_identical(prd, subprd_augdat(fits, dat_ref), info = info_txt)
+    }
+  }
+  # Other factor-like predictors: a lost first level, a character variable and
+  # a numeric variable inside factor():
+  dat2 <- dat
+  dat2$xf <- factor(c(rep(c("a", "b"), length.out = n - 1), "r"),
+                    levels = c("r", "a", "b"))
+  dat2$xs <- c(rep(c("u", "v"), length.out = n - 1), "w")
+  dat2$z <- c(rep(1:3, length.out = n - 1), 9)
+  fits2 <- suppressWarnings(divmin(
+    formula = cbind(.y.1, .y.2) ~ x + xf + xs + factor(z) + (1 | g),
+    data = dat2, family = extend_family(gaussian()), weights = wobs,
+    projpred_var = matrix(NA, n, 2), projpred_ws_aug = matrix(NA, n, 2),
+    verbose_divmin = FALSE
+  ))
+  expect_identical(attr(fits2, "projpred_lost_lvls"), list(
+    xf = list(lost = "r", ref = "a"), xs = list(lost = "w", ref = "u"),
+    z = list(lost = 9, ref = 1)
+  ))
+  dat2_ref <- dat2
+  dat2_ref$xf[n] <- "a"
+  dat2_ref$xf <- droplevels(dat2_ref$xf)
+  dat2_ref$xs[n] <- "u"
+  dat2_ref$z[n] <- 1
+  expect_identical(
+    subprd(fits2, recode_lost_levels(dat2, attr(fits2, "projpred_lost_lvls"))),
+    subprd(fits2, dat2_ref)
+  )
+  # A grouping factor is not recoded, also if it is a fixed-effect predictor:
+  # its lost level is a new group.
+  dat3 <- dat
+  dat3$g <- factor(c(as.character(dat$g[-n]), "g9"))
+  expect_null(lost_levels(dat3, keep = wobs != 0,
+                          formula = .y.1 ~ x + g + (1 | g)))
+  # Through proj_to_submodl() and the `mu_fun()` of a reference model:
+  n_drws <- 4
+  dat_r <- dat[, c("x", "g", "xc")]
+  dat_r$y <- dat$.y.1
+  refm <- suppressWarnings(init_refmodel(
+    object = list(beta = rnorm(n_drws, 1, 0.1)), data = dat_r,
+    formula = y ~ x + xc + (1 | g), family = gaussian(),
+    ref_predfun = function(fit, newdata = NULL) {
+      outer((newdata %||% dat_r)$x, fit$beta)
+    },
+    dis = rep(1, n_drws)
+  ))
+  refm$wobs[n] <- 0
+  submodl <- suppressWarnings(proj_to_submodl(
+    predictor_terms = c("x", "xc", "(1 | g)"),
+    p_ref = get_refdist(refm, ndraws = 2), refmodel = refm, verbose = FALSE
+  ))
+  expect_identical(attr(submodl$outdmin, "projpred_lost_lvls"),
+                   list(xc = list(lost = "r", ref = "a")))
+  dat_r_ref <- dat_r
+  dat_r_ref$xc[n] <- "a"
+  dat_r_ref$xc <- droplevels(dat_r_ref$xc)
+  expect_equal(refm$mu_fun(submodl$outdmin),
+               refm$mu_fun(submodl$outdmin, newdata = dat_r_ref))
+  # No level is lost if observation n has a positive weight:
+  wobs[n] <- 1
+  expect_null(attr(divmin_w0(fmls$glmm), "projpred_lost_lvls"))
+})

From 6782daf03d22724a7c7a5a99acd730f9c7298733 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Wed, 28 Oct 2026 20:40:18 +0200
Subject: [PATCH 6/6] State in a comment that data-dependent bases of GLM
 submodels still use the predictor values of weight-0 observations

---
 R/divergence_minimizers.R | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/R/divergence_minimizers.R b/R/divergence_minimizers.R
index 778b803b5..e7371464c 100644
--- a/R/divergence_minimizers.R
+++ b/R/divergence_minimizers.R
@@ -84,8 +84,10 @@ divmin <- function(
   # multilevel and additive submodel fitters. Some of these fitters (for example
   # lme4::lmer()) do not ignore the response of an observation with weight 0,
   # and others (mgcv::gam()) use its predictor values for the smooth bases. The
-  # GLM fitters handle weight 0 exactly. If `weights` and `data` differ in
-  # length (not in projpred's own calls), nothing is removed.
+  # GLM fitters handle weight 0 exactly, except that the data-dependent bases
+  # (for example, the knots of splines::ns()) still use its predictor values.
+  # If `weights` and `data` differ in length (not in projpred's own calls),
+  # nothing is removed.
   lost_lvls <- NULL
   if ((has_grp || has_add) &&
       is.data.frame(dot_args$data) && length(dot_args$weights) > 0 &&