projpred: Projection predictive feature selection

Preliminary cv_varsel() run

To illustrate a preliminary cv_varsel() run with validate_search = FALSE, we set nterms_max to the number of predictor terms in the full model, i.e., nterms_max = 20. To speed up the building of the vignette (this is not recommended in general), we choose the "L1" search method and set refit_prj to FALSE.

# Preliminary cv_varsel() run:
cvvs_fast <- cv_varsel(
  refm_obj,
  validate_search = FALSE,
  ### Only for the sake of speed (not recommended in general):
  method = "L1",
  refit_prj = FALSE,
  ###
  nterms_max = 20,
  ### In interactive use, we recommend not to deactivate the verbose mode:
  verbose = FALSE
  ### 
)

In this case, we ignore the Pareto k̂ warnings due to the reduced values for chains and iter in the reference model fit.

To find a suitable value for nterms_max in subsequent cv_varsel() runs, we take a look at a plot of at least one predictive performance statistic in dependence of the submodel size. Here, we choose the mean log predictive density (MLPD; see the documentation for argument stats of summary.vsel() for details) as the only performance statistic. Since we will be using the following plot only to determine nterms_max for subsequent cv_varsel() runs, we can omit the predictor ranking from the plot by setting ranking_nterms_max to NA:

plot(cvvs_fast, stats = "mlpd", ranking_nterms_max = NA)

This plot suggests that the submodel MLPD levels off for the first time from submodel size 9 to 10 and that it increases only slightly afterwards. However, we used L1 search with refit_prj = FALSE, which means that the projections employed for the predictive performance evaluation are L1-penalized, which is usually undesired (Piironen, Paasiniemi, and Vehtari 2020, sec. 4). Thus, to investigate the impact of refit_prj = FALSE, we re-run cv_varsel(), but this time with the default of refit_prj = TRUE and re-using the search results (as well as the CV-related arguments such as validate_search; see section “Usage” in ?cv_varsel.vsel) from the former cv_varsel() call (this is done by applying cv_varsel() to cvvs_fast instead of refm_obj so that the cv_varsel() generic dispatches to the cv_varsel.vsel() method that was introduced in projpred 2.8.0⁶). To save time, we also set nclusters_pred to a comparatively low value of 20:

# Preliminary cv_varsel() run with `refit_prj = TRUE`:
cvvs_fast_refit <- cv_varsel(
  cvvs_fast,
  ### Only for the sake of speed (not recommended in general):
  nclusters_pred = 20,
  ###
  ### In interactive use, we recommend not to deactivate the verbose mode:
  verbose = FALSE
  ### 
)

As explained above, we ignore the Pareto k̂ warnings here. We also ignore the warning that SIS is used instead of PSIS (this is due to nclusters_pred = 20 which we used only to speed up the building of the vignette).

With the refit_prj = TRUE results, the predictive performance plot from above now looks as follows:

plot(cvvs_fast_refit, stats = "mlpd", ranking_nterms_max = NA)

This refined plot shows that the submodel MLPD does not change much after submodel size 8, so in our final cv_varsel() run, we set nterms_max to a value slightly higher than 8 (here: 9) to ensure that we see the MLPD leveling off. The search results from the initial cvvs_fast object could now be re-used again (via cv_varsel.vsel()) for investigating the sensitivity of the results to changes in nclusters_pred (or ndraws_pred). Here, we skip this for the sake of brevity and instead head over to the final cv_varsel() run.

Final cv_varsel() run

For this final cv_varsel() run (with validate_search = TRUE, as recommended), we use a K-fold CV with a small number of folds (K = 2) to make this vignette build faster. In practice, we recommend using either the default of cv_method = "LOO" or a larger value for K if this is possible in terms of computation time. Here, we also perform the K reference model refits outside of cv_varsel(). Although not strictly necessary here, this is helpful in practice because often, cv_varsel() needs to be re-run multiple times in order to try out different argument settings. We also illustrate how projpred’s CV (i.e., the CV comprising search and performance evaluation, after refitting the reference model K times) can be parallelized, even though this is of little use here (we have only K = 2 folds and the fold-wise searches and performance evaluations are quite fast, so the parallelization overhead eats up any runtime improvements). Note that for this final cv_varsel() run, we cannot make use of cv_varsel.vsel() (by applying cv_varsel() to cvvs_fast or cvvs_fast_refit instead of refm_obj) because of the change in nterms_max (here: 9, for cvvs_fast and cvvs_fast_refit: 20).

# Refit the reference model K times:
cv_fits <- run_cvfun(
  refm_obj,
  ### Only for the sake of speed (not recommended in general):
  K = 2
  ###
)
# For running projpred's CV in parallel (see cv_varsel()'s argument `parallel`):
doParallel::registerDoParallel(ncores)
# Final cv_varsel() run:
cvvs <- cv_varsel(
  refm_obj,
  cv_method = "kfold",
  cvfits = cv_fits,
  ### Only for the sake of speed (not recommended in general):
  method = "L1",
  nclusters_pred = 20,
  ###
  nterms_max = 9,
  parallel = TRUE,
  ### In interactive use, we recommend not to deactivate the verbose mode:
  verbose = FALSE
  ### 
)
# Tear down the CV parallelization setup:
doParallel::stopImplicitCluster()
foreach::registerDoSEQ()

We can now select a final submodel size by looking at a predictive performance plot similar to the one created for the preliminary cv_varsel() run above. By default, the performance statistics are plotted on their original scale, but with deltas = TRUE, they are plotted as differences⁷ from a baseline model (which is the reference model by default, at least in the most common cases). Since the differences⁸ and the (frequentist) uncertainty in their estimation are usually of more interest than the original-scale performance statistics (at least with regard to the decision for a final submodel size), we directly plot with deltas = TRUE here:

plot(cvvs, stats = "mlpd", deltas = TRUE)

Decision for final submodel size

Based on that final predictive performance plot, we decide for a submodel size. Usually, the aim is to find the smallest submodel size where the predictive performance of the submodels levels off and is close enough to the reference model’s predictive performance (the dashed red horizontal line). Sometimes (as here), the plot may be ambiguous because after reaching the reference model’s performance, the submodels’ performance may keep increasing (and hence become even better than the reference model’s performance⁹). In that case, one has to find a suitable trade-off between predictive performance (accuracy) and model size (sparsity) in the context of subject-matter knowledge. Here, we assume that the focus of our predictive model is sparsity (not accuracy). Hence, based on the plot, we decide for a submodel size of 6 because this is the smallest size where the submodel MLPD is close enough to the reference model MLPD:

size_decided <- 6

If the focus of our predictive model had been accuracy (not sparsity), size 7 (or even 8) would have been a natural choice because at size 7 (at the latest at size 8), the submodel MLPD can be considered to “level off”. Further down, the predictor ranking and the (CV) ranking proportions that are shown in the plot (below the submodel sizes on the x-axis) are explained in detail—and also how they could have been incorporated into our decision for a submodel size.

The suggest_size() function offered by projpred may help in the decision for a submodel size, but this is a rather heuristic method and needs to be interpreted with caution (see ?suggest_size):

suggest_size(cvvs, stat = "mlpd")

With this heuristic, we would get the same final submodel size (6) as by our manual (sparsity-based) decision.

summary() method for vsel objects

A tabular representation of the plot created by plot.vsel() can be achieved via summary.vsel(). For the output of summary.vsel(), there is a sophisticated print() method (print.vselsummary()) which is also called by the shortcut method print.vsel()¹⁰. Specifically, to create the summary table matching the predictive performance plot above as closely as possible (and to also adjust the minimum number of printed significant digits), we may call summary.vsel() and print.vselsummary() as follows:

smmry <- summary(cvvs, stats = "mlpd", type = c("mean", "lower", "upper"),
                 deltas = TRUE)
print(smmry, digits = 1)

Predictor ranking(s) and identification of the selected submodel

As indicated by its column name, the predictor ranking from column ranking_fulldata of the summary.vsel() output is based on the full-data search. In case of cv_varsel() with validate_search = TRUE, there is not only the full-data search, but also fold-wise searches, implying that there are also fold-wise predictor rankings. All of these predictor rankings (the full-data one and—if available—the fold-wise ones) can be retrieved via ranking():

rk <- ranking(cvvs)

In addition to inspecting the full-data predictor ranking, it usually makes sense to investigate the ranking proportions derived from the fold-wise predictor rankings (only available in case of cv_varsel() with validate_search = TRUE, which we have here) in order to get a sense for the variability in the ranking of the predictors. For a given predictor x and a given submodel size j, the ranking proportion is the proportion of CV folds which have predictor x at position j of their predictor ranking. To compute these ranking proportions, we use cv_proportions():

( pr_rk <- cv_proportions(rk) )

Here, the ranking proportions are of little use as we have used K = 2 (in the final cv_varsel() call above) for the sake of speed. Nevertheless, we can see that the two CV folds agree on the most relevant predictor term (X1). Since the column names of the matrix returned by cv_proportions() follow the full-data predictor ranking, we can infer that X1 is also the most relevant predictor term in the full-data predictor ranking. To see this more explicitly, we can access element fulldata of the ranking() output:

rk[["fulldata"]]

This is the same as column ranking_fulldata in the summary.vsel() output above (apart from the intercept). (As also stated by the message thrown by print.vselsummary(), column cv_proportions_diag of the summary.vsel() output contains the main diagonal of the matrix that we stored manually as pr_rk.)

Note that we have cut off the search at nterms_max = 9 (which is smaller than the number of predictor terms in the full model, 20 here), so the ranking proportions in the pr_rk matrix do not need to sum to 100 % (neither column-wise nor row-wise).

The transposed matrix of ranking proportions can be visualized via plot.cv_proportions():

plot(pr_rk)

Apart from visualizing the variability in the ranking of the predictors (here, this is of little use because of K = 2), this plot will be helpful later.

To retrieve the predictor terms of the final submodel (except for the intercept which is always included in the submodels), we combine the chosen submodel size of 6 with the full-data predictor ranking:

( predictors_final <- head(rk[["fulldata"]], size_decided) )

At this place, it is again helpful to take the ranking proportions into account, but now in a cumulated fashion:

plot(cv_proportions(rk, cumulate = TRUE))

This plot shows that the two fold-wise searches (as well as the full-data search, whose predictor ranking determines the order of the predictors on the y-axis) agree on the set of the 6 most relevant predictors: When looking at <=6 on the x-axis, all tiles above and including the 6^th main diagonal element are at 100 %. Similarly, the two CV folds agree on the set of the 7 most relevant predictors (and also on the set of the most relevant predictor, which is a singleton).

Although not demonstrated here, the cumulated ranking proportions could also have guided the decision for a submodel size (if we had not been willing to follow a strict rule based on accuracy or sparsity): From their plot, we can see that size 9 might have been an unfortunate choice because X10 (which—by cutting off the full-data predictor ranking at size 9—would then have been selected as the 9^th predictor in the final submodel) is not included among the first 9 terms by any CV fold. However, since K = 2 is too small for reliable statements regarding the variability of the predictor ranking, we did not take the cumulated ranking proportions into account when we made our decision for a submodel size above.

In a real-world application, we might also be able to incorporate the full-data predictor ranking into our decision for a submodel size (usually, this requires to also take into account the variability of the predictor ranking, as reflected by the—possibly cumulated—ranking proportions). For example, the predictors might be associated with different measurement costs, so that we might want to select a costly predictor only if the submodel size at which it would be selected (according to the full-data predictor ranking, but taking into account that there might be variability in the ranking of the predictors) comes with a considerable increase in predictive performance.

Marginals of the projected posterior

The posterior package provides a general way to deal with posterior distributions, so it can also be applied to our projected posterior. For example, to calculate summary statistics for the marginals of the projected posterior:

library(posterior)

prj_drws <- as_draws_matrix(prj_mat)
prj_smmry <- summarize_draws(
  prj_drws,
  "median", "mad", function(x) quantile(x, probs = c(0.025, 0.975))
)
# Coerce to a `data.frame` because pkgdown versions > 1.6.1 don't print the
# tibble correctly:
prj_smmry <- as.data.frame(prj_smmry)
print(prj_smmry, digits = 1)

A visualization of the projected posterior can be achieved with the bayesplot package, for example using its mcmc_intervals() function:

library(bayesplot)

bayesplot_theme_set(ggplot2::theme_bw())
mcmc_intervals(prj_mat) +
  ggplot2::coord_cartesian(xlim = c(-1.5, 1.6))

Note that we only visualize the 1-dimensional marginals of the projected posterior here. To gain a more complete picture, we would have to visualize at least some 2-dimensional marginals of the projected posterior (i.e., marginals for pairs of parameters).

For comparison, consider the marginal posteriors of the corresponding parameters in the reference model:

refm_mat <- as.matrix(refm_fit)
mcmc_intervals(refm_mat, pars = colnames(prj_mat)) +
  ggplot2::coord_cartesian(xlim = c(-1.5, 1.6))

Here, the reference model’s marginal posteriors differ only slightly from the marginals of the projected posterior. This does not necessarily have to be the case.

Predictions

Predictions from the final submodel can be made by proj_linpred() and proj_predict().

We start with proj_linpred(). For example, suppose we have the following new observations:

( dat_gauss_new <- setNames(
  as.data.frame(replicate(length(predictors_final), c(-1, 0, 1))),
  predictors_final
) )

Then proj_linpred() can calculate the linear predictors¹⁴ for all new observations from dat_gauss_new. Depending on argument integrated, these linear predictors can be averaged across the projected draws (within each new observation). For instance, the following computes the expected values of the new observations’ predictive distributions (beware that the following code refers to the Gaussian family with the identity link function; for other families—which usually come in combination with a different link function—one would typically have to use transform = TRUE in order to achieve such expected values):

prj_linpred <- proj_linpred(prj, newdata = dat_gauss_new, integrated = TRUE)
cbind(dat_gauss_new, linpred = as.vector(prj_linpred$pred))

If dat_gauss_new also contained response values (i.e., y values in this example), then proj_linpred() would also evaluate the log predictive density at these (conditional on each of the projected parameter draws if integrated = FALSE and integrated over the projected parameter draws—before taking the logarithm—if integrated = TRUE).

With proj_predict(), we can obtain draws from predictive distributions based on the final submodel. In contrast to proj_linpred(<...>, integrated = FALSE), this encompasses not only the uncertainty arising from parameter estimation, but also the uncertainty arising from the observation (or “sampling”) model for the response¹⁵. This is useful for what is usually termed a posterior predictive check (PPC), but would have to be termed something like a posterior-projection predictive check (PPPC) here:

prj_predict <- proj_predict(prj)
# Using the 'bayesplot' package:
ppc_dens_overlay(y = dat_gauss$y, yrep = prj_predict)

This PPPC shows that our final projection is able to generate predictions similar to the observed response values, which indicates that this model is reasonable, at least in this regard.

Non-convergence of predictive performance

Sometimes, the predictor ranking makes sense, but for an increasing submodel size, the predictive performance of the submodels does not approach the reference model’s predictive performance so that the submodels exhibit a predictive performance that stays worse than the reference model’s. There are different reasons that can explain this behavior (the following list might not be exhaustive, though):

1. The reference model’s posterior may be so wide that the default ndraws_pred could be too small. Usually, this comes in combination with a difference in predictive performance which is comparatively small. Increasing ndraws_pred should help, but it also increases the computational cost. Refitting the reference model and thereby ensuring a narrower posterior (usually by employing a stronger sparsifying prior) should have a similar effect.
2. For non-Gaussian models, the discrepancy may be due to the fact that the penalized iteratively reweighted least squares (PIRLS) algorithm might have convergence issues (Catalina, Bürkner, and Vehtari 2021). In this case, the latent-space approach by Catalina, Bürkner, and Vehtari (2021) might help, see also the latent-projection vignette linked above.

Overfitting

If you are using varsel(), then the lack of CV in varsel() may lead to overconfident and overfitted results. In this case, try running cv_varsel() instead of varsel() (which you should in any case for your final results).

Similarly, cv_varsel() with validate_search = FALSE is more prone to overfitting than cv_varsel() with validate_search = TRUE.

Issues with the traditional projection

For multilevel binomial models, the traditional projection may not work properly and give suboptimal results, see #353 on GitHub (the underlying issue is described in lme4 issue #682). With suboptimality of the results, we mean that in such cases, the relevance of the group-level terms can be underestimated. According to the simulation-based case study from #353, the latent projection should be considered as a currently available remedy.

For multilevel Poisson models, the traditional projection may take very long, see #353. According to the simulation-based case study from #353, the latent projection should be considered as a currently available remedy.

Finally, as illustrated in the Poisson example of the latent-projection vignette, the latent projection can be beneficial for non-multilevel models with a (non-Gaussian) family that is also supported by the traditional projection, at least in case of the Poisson family and L1 search.

Issues with the augmented-data projection

For multilevel models, the augmented-data projection seems to suffer from the same issue as the traditional projection for the binomial family (see above), i.e., it may not work properly and give suboptimal results, see #353 (the underlying issue is probably similar to the one described in lme4 issue #682). With suboptimality of the results, we mean that in such cases, the relevance of the group-level terms can be underestimated. According to the simulation-based case study from #353, the latent projection should be considered as a remedy in such cases.

Speed

There are many ways to speed up projpred, but in general, such speed-ups lead to results that are less accurate and hence should only be considered as preliminary results. Some speed-up possibilities are:

1. Using cv_varsel() with validate_search = FALSE instead of validate_search = TRUE. In case of cv_method = "LOO"²⁶, cv_varsel() with validate_search = FALSE has comparable runtime to varsel(), but accounts for some overfitting, namely that induced by varsel()’s in-sample predictions during the predictive performance evaluation. However, as explained in section “Variable selection” (see also section “Overfitting”), cv_varsel() with validate_search = FALSE is more prone to overfitting than cv_varsel() with validate_search = TRUE.

2. Using cv_varsel() with K-fold CV instead of PSIS-LOO CV. Whether this provides a speed improvement mainly depends on the number of observations and the complexity of the reference model. Note that PSIS-LOO CV is often more accurate than K-fold CV if argument K is (much) smaller than the number of observations.

3. Using a “custom” reference model object with a dimension reduction technique for the predictor data (e.g., by computing principal components from the original predictors, using these principal components as predictors when fitting the reference model, and then performing the variable selection in terms of the original predictor terms). Examples are given in Piironen, Paasiniemi, and Vehtari (2020) and Pavone et al. (2022). This approach makes sense if there is a large number of predictor variables, in which case this aims at improving the runtime required for fitting the reference model and hence improving the runtime of K-fold CV.

4. Using varsel() with its argument d_test for evaluating predictive performance on a hold-out dataset instead of doing this with cv_varsel()’s CV approach. Typically, the hold-out approach requires a large amount of data.

5. Reducing nterms_max in varsel() or cv_varsel(). The resulting predictive performance plot(s) should be inspected to check that the search is not terminated too early (i.e., before the submodel performance levels off), which would indicate that nterms_max has been reduced too much.

6. Reducing argument nclusters (of varsel() or cv_varsel()) below 20 and/or setting nclusters_pred to some non-NULL (and smaller than 400, the default for ndraws_pred) value. If setting nclusters_pred as low as nclusters (and using forward search), refit_prj can instead be set to FALSE, see below.

7. Using L1 search (see argument method of varsel() or cv_varsel()) instead of forward search. Note that L1 search implies nclusters = 1 and is not always supported. In general, forward search is more accurate than L1 search and hence more desirable (see section “Details” in ?varsel or ?cv_varsel for a more detailed comparison of the two). The issue demonstrated in the Poisson example from the latent-projection vignette is related to this.

8. Setting argument refit_prj (of varsel() or cv_varsel()) to FALSE, which basically means to set ndraws_pred = ndraws and nclusters_pred = nclusters, but in a more efficient (i.e., faster) way. In case of L1 search, this means that the L1-penalized projections of the regression coefficients are used for the predictive performance evaluation, which is usually undesired (Piironen, Paasiniemi, and Vehtari 2020, sec. 4). In case of forward search, this issue does not exist.

9. Parallelizing costly parts of the CV implied by cv_varsel() (this was demonstrated in the example above; see argument parallel of cv_varsel()). When using project(), parallelizing the projection might also help (see the general package documentation available online or by typing ?`projpred-package`).

10. Using varsel.vsel() or cv_varsel.vsel() to re-use previous search results for new performance evaluation(s). This is helpful if the performance evaluation part is run multiple times based on the same search results (e.g., when only arguments ndraws_pred or nclusters_pred of varsel() or cv_varsel() are changed). In the example above, this was illustrated when cv_varsel() was applied to cvvs_fast instead of refm_obj to yield cvvs_fast_refit. In that example, search and performance evaluation were effectively run separately by the cv_varsel() calls yielding cvvs_fast and cvvs_fast_refit, respectively. This is because cv_varsel() with refit_prj = FALSE (used for cvvs_fast) has almost only the computational cost of the search (the performance evaluation with refit_prj = FALSE has almost no computational cost) and cv_varsel.vsel() (used for cvvs_fast_refit) has almost only the computational cost of the performance evaluation (the search in cv_varsel.vsel() has no computational cost at all).

11. Using run_cvfun() in case of repeated K-fold CV with the same K reference model refits. The output of run_cvfun() is typically used as input for argument cvfits of cv_varsel.refmodel() (so in order to have a speed improvement, the output of run_cvfun() needs to be assigned to an object which is then re-used in multiple cv_varsel() calls).