Leave the left-out observation out of the fold projections in validated LOO CV

NEWS.md

@@ -13,6 +13,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 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

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, ...
       )

R/divergence_minimizers.R

@@ -78,6 +78,29 @@ 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, 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 &&
+      any(dot_args$weights == 0) &&
+      length(dot_args$weights) == nrow(dot_args$data)) {
+    keep <- dot_args$weights != 0
+    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]
+  }
+
   if (is.null(getOption("projpred.parallel_proj_trigger")) &&
       !is.null(getOption("projpred.prll_prj_trigger"))) {
     warning(
@@ -104,13 +127,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 +150,6 @@ divmin <- function(
     } else {
       progressor_obj <- NULL
     }
-    dot_args <- list(...)
     `%do_projpred%` <- foreach::`%dopar%`
     outdmin <- foreach::foreach(
       formula_s = formulas,
@@ -166,9 +188,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,
@@ -696,6 +829,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(
@@ -801,6 +951,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)
 }
 

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)

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