Variable selection
Now, projpred comes into play.
From the reference model fit (called refm_fit here), we
create a reference model object (i.e., an object of class
refmodel) since this avoids redundant calculations in the
remainder of this vignette3:
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 mainly
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, the
search is excluded from the CV so that only a single full-data search is
run. Because of its more 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).
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 = 0
###
)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 (which requires
size_position = "primary_x_bottom", which we set here via
its global option, just like argument text_angle):
options(projpred.plot_vsel_size_position = "primary_x_bottom")
options(projpred.plot_vsel_text_angle = 0)
plot(cvvs_fast, stats = "mlpd", ranking_nterms_max = NA)
This plot suggests that the submodel MLPD levels off from submodel size
8 on. 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.05). 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 = 0
###
)Using standard importance sampling (SIS) due to a small number of clusters.Here, we 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:
This refined plot confirms 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"
(possibly subsampled, see argument nloo of
cv_varsel()) 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 = 0
###
)
# Tear down the CV parallelization setup:
doParallel::stopImplicitCluster()
foreach::registerDoSEQ()Predictive performance plot from final cv_varsel()
run
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 actual scale and the uncertainty bars
match this scale, but argument deltas of
plot.vsel() offers two more options:
- With
deltas = TRUE, the performance statistics are plotted on difference scale, i.e., as differences6 from the baseline model7 and the uncertainty bars match this scale, - With
deltas = "mixed", the performance statistics (i.e., their point estimates) are plotted on the actual scale, but the uncertainty bars visualize the difference-scale uncertainty.
Since the difference-scale uncertainty is usually more helpful than
the actual-scale uncertainty (at least with regard to the decision for a
final submodel size), we plot with deltas = TRUE here
(deltas = "mixed" would be another good choice). We also
set show_cv_proportions = TRUE (via its global option) for
illustrative purposes, which becomes clearer in section “Predictor ranking(s) from final cv_varsel()
run and identification of the selected submodel”:
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, 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 performance8). 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 decide for a submodel size of 7 because it seems to provide the best trade-off between sparsity and accuracy (size 7 is the smallest size where the submodel MLPD is close enough to the reference model MLPD and from size 7 on, the submodel MLPD levels off).
In section “Predictor ranking(s) from final
cv_varsel() run and identification of the selected
submodel”, 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] 7In this case, that heuristic gives the same final submodel size
(7) as our manual decision.
Predictive performance table from final cv_varsel()
run
A tabular representation of the plot created by
plot.vsel() can be achieved via summary.vsel()
or performances(). In contrast to
performances(), the output of summary.vsel()
contains more information than just the predictive performance results,
which is also why there is a sophisticated print() method
for objects of class vselsummary (the output of
summary.vsel()). This method
print.vselsummary() is also called by the shortcut method
print.vsel() which can be applied to an object resulting
from varsel() or cv_varsel().
Specifically, to create the 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)
Argument `refit_prj`: TRUE
Submodel performance evaluation summary with `deltas = TRUE` and `cumulate = FALSE`:
size ranking_fulldata cv_proportions_diag mlpd mlpd.lower mlpd.upper
0 (Intercept) NA -0.976 -1.05 -0.903
1 X1 1.0 -0.786 -0.87 -0.704
2 X14 1.0 -0.555 -0.62 -0.487
3 X5 0.5 -0.553 -0.64 -0.469
4 X20 0.0 -0.519 -0.61 -0.429
5 X6 0.0 -0.208 -0.27 -0.151
6 X3 0.0 -0.214 -0.26 -0.169
7 X8 0.5 -0.003 -0.02 0.016
8 X11 0.0 -0.010 -0.03 0.009
9 X7 0.0 0.003 -0.01 0.021
Reference model performance evaluation summary with `deltas = TRUE`:
mlpd mlpd.lower mlpd.upper
0 0 0 The generic function performances() (with main method
performances.vselsummary() and the shortcut method
performances.vsel()) essentially extracts the predictive
performance results from the output of summary.vsel():
List of 2
$ submodels :'data.frame': 10 obs. of 4 variables:
..$ size : num [1:10] 0 1 2 3 4 5 6 7 8 9
..$ mlpd : num [1:10] -0.976 -0.786 -0.555 -0.553 -0.519 ...
..$ mlpd.lower: num [1:10] -1.049 -0.867 -0.623 -0.637 -0.609 ...
..$ mlpd.upper: num [1:10] -0.903 -0.704 -0.487 -0.469 -0.429 ...
$ reference_model: Named num [1:3] 0 0 0
..- attr(*, "names")= chr [1:3] "mlpd" "mlpd.lower" "mlpd.upper"
- attr(*, "class")= chr "performances"Predictor ranking(s) from final cv_varsel() run 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. This full-data predictor ranking is
also what is shown in the second line of the x-axis tick labels of the
predictive performance plot from section “Predictive performance plot from final
cv_varsel() run”.
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 often
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 X7
1 1 0 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 0
3 0 0 0.5 0 0.5 0.0 0.0 0.0 0
4 0 0 0.5 0 0.0 0.5 0.0 0.0 0
5 0 0 0.0 1 0.0 0.0 0.0 0.0 0
6 0 0 0.0 0 0.5 0.0 0.5 0.0 0
7 0 0 0.0 0 0.0 0.5 0.5 0.0 0
8 0 0 0.0 0 0.0 0.0 0.0 0.0 0
9 0 0 0.0 0 0.0 0.0 0.0 0.5 0
attr(,"class")
[1] "cv_proportions"The main diagonal of this matrix is contained in column
cv_proportions_diag of the summary.vsel()
output and is also shown in the third line of the x-axis tick labels of
the predictive performance plot from section “Predictive performance plot from final
cv_varsel() run” (due to
show_cv_proportions = TRUE set via its global option).
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) and the
second most relevant predictor term (X14). 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 most relevant predictor terms (in this
order) 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" "X7" This is the same as column ranking_fulldata in the
summary.vsel() output above (apart from the intercept).
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 7 with the full-data predictor ranking:
[1] "X1" "X14" "X5" "X20" "X6" "X3" "X8" At this place, it is again helpful to take the ranking proportions into account, but now in a cumulated fashion:
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 7 most
relevant predictors: When looking at <=7 on the x-axis,
all tiles above and including the 7th main diagonal element
are at 100 %. (Similarly, the two CV folds also agree on the set of the
most relevant predictor—which is a singleton—and the set of the two most
relevant predictors, but we already observed this above.)
Although not demonstrated here, the cumulated ranking proportions
also could have guided the decision for a submodel size: From the plot
of the cumulated ranking proportions, we can see that size 8 might have
been an unfortunate choice because X11 (which—by cutting
off the full-data predictor ranking at size 8—would then have been
selected as the 8th predictor in the final submodel) is not
included among the first 8 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.
