Variable selection
Now, projpred comes into play.
In projpred, the projection predictive variable selection relies on a so-called search part and a so-called evaluation part. The search part determines the predictor ranking (also known as solution path), i.e., the best submodel for each submodel size (the size is given by the number of predictor terms). The evaluation part determines the predictive performance of the increasingly complex submodels along the predictor ranking.
There are two functions for running the combination of search and
evaluation: varsel() and cv_varsel(). In
contrast to varsel(), cv_varsel() performs a
cross-validation (CV). With cv_method = "LOO" (the
default), cv_varsel() runs a Pareto-smoothed importance
sampling leave-one-out CV (PSIS-LOO CV, see Vehtari, Gelman, and Gabry 2017;
Vehtari et al. 2022). With
cv_method = "kfold", cv_varsel() runs a \(K\)-fold CV. The extent of the CV depends
on cv_varsel()’s argument validate_search: If
validate_search = TRUE (the default), the search part is
run with the training data of each CV fold separately and the evaluation
part is run with the corresponding test data of each CV fold. If
validate_search = FALSE (which is currently only available
for cv_method = "LOO"), the search is excluded from the CV
so that only a single full-data search is run. Because of its most
thorough protection against overfitting4,
cv_varsel() with validate_search = TRUE is
recommended over varsel() and cv_varsel() with
validate_search = FALSE. Nonetheless, a preliminary and
comparatively fast run of varsel() or
cv_varsel() with validate_search = FALSE can
give a rough idea of the performance of the submodels and can be used
for finding a suitable value for argument nterms_max in
subsequent runs (argument nterms_max imposes a limit on the
submodel size up to which the search is continued and is thus able to
reduce the runtime considerably).
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 nclusters_pred to a comparatively low value of
20.
# Preliminary cv_varsel() run:
cvvs_fast <- cv_varsel(
refm_fit,
validate_search = FALSE,
### Only for the sake of speed (not recommended in general):
method = "L1",
nclusters_pred = 20,
###
nterms_max = 20,
### In interactive use, we recommend not to deactivate the verbose mode:
verbose = FALSE
###
)Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.Warning in warn_pareto(n07 = sum(pareto_k > 0.7), n05 = sum(0.7 >= pareto_k & :
In the calculation of the reference model's PSIS-LOO CV weights, 6 (out of 100)
Pareto k-values are in the interval (0.5, 0.7]. Moment matching (see the loo
package), mixture importance sampling (see the loo package), and `reloo`-ing
(see the brms package) are not supported by projpred. If these techniques (run
outside of projpred, i.e., for the reference model only; note that `reloo`-ing
may be computationally costly) result in a markedly different reference model
ELPD estimate than ordinary PSIS-LOO CV does, we recommend to use K-fold CV
within projpred.Warning in loo_varsel(refmodel = refmodel, method = method, nterms_max =
nterms_max, : The projected draws used for the performance evaluation have
different (i.e., nonconstant) weights, so using standard importance sampling
(SIS) instead of Pareto-smoothed importance sampling (PSIS). In general, PSIS
is recommended over SIS.In this case, we ignore the Pareto \(\hat{k}\) warnings due to the reduced
values for chains and iter in the reference
model fit. 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).
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:
This plot (see ?plot.vsel for a description) 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: nterms_max = 9) to
ensure that we see the MLPD leveling off.
For this final cv_varsel() run, 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"
(with validate_search = TRUE) or a larger value for
K if this is possible in terms of computation time. We also
illustrate how projpred’s CV 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).
# For the CV parallelization (cv_varsel()'s argument `parallel`):
doParallel::registerDoParallel(ncores)
# Final cv_varsel() run:
cvvs <- cv_varsel(
refm_fit,
cv_method = "kfold",
### Only for the sake of speed (not recommended in general):
K = 2,
method = "L1",
nclusters_pred = 20,
###
nterms_max = 9,
parallel = TRUE,
### In interactive use, we recommend not to deactivate the verbose mode:
verbose = FALSE
###
)Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess
Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess# Tear down the CV parallelization setup:
doParallel::stopImplicitCluster()
foreach::registerDoSEQ()Again, we ignore the bulk-ESS warnings due to the reduced values for
chains and iter in the reference model
fit.
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:
Based on that 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 performance5). 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:
If the focus of our predictive model had been accuracy (not sparsity), size 8 would have been a natural choice because at size 8, the submodel MLPD “levels off” (in fact, size 8 here even comes with the maximum submodel MLPD among all plotted sizes, but in general, this does not need to be the case). Further below, 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):
[1] 6With this heuristic, we would get the same final submodel size
(6) as by our manual (sparsity-based) decision.
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()6.
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)
Family: gaussian
Link function: identity
Formula: y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 + X10 + X11 +
X12 + X13 + X14 + X15 + X16 + X17 + X18 + X19 + X20
Observations: 100
Projection method: traditional
CV method: K-fold CV with K = 2 and search included (i.e., fold-wise searches)
Search method: L1
Maximum submodel size for the search: 9
Number of projected draws in the search: 1 (from clustered projection)
Number of projected draws in the performance evaluation: 20 (from clustered projection)
Performance evaluation summary with `deltas = TRUE` and `cumulate = FALSE`:
size solution_terms cv_proportions_diag mlpd.kfold lower upper
0 <NA> NA -0.92 -1.01 -0.83
1 X1 1.0 -0.70 -0.79 -0.61
2 X14 1.0 -0.57 -0.66 -0.48
3 X5 0.5 -0.52 -0.63 -0.41
4 X20 1.0 -0.33 -0.42 -0.24
5 X6 0.0 -0.16 -0.25 -0.08
6 X3 0.0 0.05 -0.01 0.12
7 X8 0.5 0.08 0.04 0.12
8 X11 0.0 0.11 0.07 0.14
9 X10 0.0 0.09 0.06 0.13Column `solution_terms` contains the full-data predictor ranking. To retrieve the fold-wise predictor rankings, use the ranking() function, possibly followed by cv_proportions() for computing the ranking proportions (which can be visualized by plot.cv_proportions()). The main diagonal of the matrix returned by cv_proportions() (with `cumulate = FALSE`) is contained in column `cv_proportions_diag`.As highlighted by the message above, the predictor ranking from
column solution_terms 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():
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():
predictor
size X1 X14 X5 X20 X6 X3 X8 X11 X10
1 1 0 0.0 0 0.0 0 0.0 0.0 0
2 0 1 0.0 0 0.0 0 0.0 0.0 0
3 0 0 0.5 0 0.5 0 0.0 0.0 0
4 0 0 0.0 1 0.0 0 0.0 0.0 0
5 0 0 0.0 0 0.0 1 0.0 0.0 0
6 0 0 0.5 0 0.5 0 0.0 0.0 0
7 0 0 0.0 0 0.0 0 0.5 0.5 0
8 0 0 0.0 0 0.0 0 0.5 0.0 0
9 0 0 0.0 0 0.0 0 0.0 0.5 0
attr(,"class")
[1] "cv_proportions"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 set of the two most relevant predictor terms
(X1 and X14) as well as on their order. Since
the column names of the matrix returned by cv_proportions()
follow the full-data predictor ranking, we can infer that
X1 and X14 are also the two most relevant
predictor terms in the full-data predictor ranking. To see this more
explicitly, we can access element fulldata of the
ranking() output:
[1] "X1" "X14" "X5" "X20" "X6" "X3" "X8" "X11" "X10"This is the same as column solution_terms in the
summary.vsel() output above. (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():
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:
[1] "X1" "X14" "X5" "X20" "X6" "X3" At this place, it is again helpful to take the ranking proportions into account, but now in a cumulated fashion:
This plot confirms (from a slightly different perspective) 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 two most relevant predictors (X1
and X14): When looking at <=2 on the
x-axis, all tiles above and including the second main diagonal element
are at 100 %. Similarly, the two CV folds agree on the sets of
the six and nine 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 ninth
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.
