\begin@twocolumnfalse

The concepts of Bayesian prediction, model comparison, and model selection have developed significantly over the last decade. As a result, the Bayesian community has witnessed a rapid growth in theoretical and applied contributions to building and selecting predictive models. Projection predictive inference in particular has shown promise to this end, finding application across a broad range of fields. It is less prone to over-fitting than naïve selection based purely on cross-validation or information criteria performance metrics, and has been known to out-perform other methods in terms of predictive performance. We survey the core concept and contemporary contributions to projection predictive inference, and present a safe, efficient, and modular workflow for prediction-oriented model selection therein. We also provide an interpretation of the projected posteriors achieved by projection predictive inference in terms of their limitations in causal settings.

\end@twocolumnfalse

§ INTRODUCTION

Bayesian inference has gained traction in recent years as computational availability has made it more broadly accessible. Alongside this surge in popularity and accessibility, the application of Bayesian models to predictive tasks has also become more widespread. In practice, we are often interested in identifying a model whose out-of-sample predictive performance is best while minimising measurement and inferential cost. This view of Bayesian inference was already discussed by 50 who developed a predictive model selection paradigm in which one constructs small models capable of approximately replicating the predictive performance of a so-called “reference model”. This is the best-performing model (in terms of its predictive performance) available which represents all uncertainty related to the task. Essentially, we fit smaller models to the fit of the reference model (we call this procedure the “projection”) and use these for their improved interpretability, or to decrease data collection cost [72]. The computational approximation proposed by 37 and 22 has made this technique feasible both computationally and in terms of a practical workflow. This procedure has become known as projection predictive model selection [101], or, more generally, projection predictive inference [72].

Due to its excellent performance, computational efficiency, and easy-to-use software, the projection predictive approach has gained popularity in many real-world problems (see Appendix <ref> for an enumeration of some of these). The projpred package  [71], based on projection predictive inference, is the most popular Bayesian variable selection package for R when comparing its download numbers to its alternatives (e.g. those recorded in Appendix <ref>) with the help of cranlogs [19] (last check: April 26, 2023).

The use of a reference model in model selection (i.e. substituting the observed data with the in-sample predictions made by the reference model) is not in itself a novel idea. In addition to 50, 48 previously motivated the use of the Kullback-Leibler (KL) divergence [46] to quantify the distance between models' predictive distributions, and 81 constructed a Bayesian model-averaged reference model. 6 later used a symmetric divergence measure, and 63 and 94 extended this to use different divergences. 97 discuss the interpretation implied by the different divergence measures used.

The idea of a reference model is also common in non-Bayesian literature, where, for example, 40 referred to them as a “gold standard” model one then seeks to approximate. Such approximation becomes increasingly important as the model space grows. 28 discussed this in the context of identifying risk groups within large neural networks. Indeed, as wide neural networks grow wider, it becomes more likely that one of their subnetworks is useful for prediction. This is closely related to the so-called “student-teacher” learning framework [103]. Under the terminology of “pre-conditioning,” 66 and 39 discuss fitting models to consistent estimates from the reference model. Across different fields and under different names, the notion of fitting to the fit of a reference model has become a popular approach to model selection.

Furthermore, projective inference is significantly cheaper than rerunning Markov chain Monte Carlo (MCMC) for submodel fitting (we will use the terms “submodel” and “restricted model” interchangeably throughout this paper), making it desirable also from a computational perspective. The projection has the advantage of being a deterministic transformation of the reference model's posterior and therefore requires a prior only for the reference model.

67 show that aside from the projection, other components of the projection predictive model selection framework contribute to the superior performance of projection predictive model selection. The use of cross-validation (CV) for selecting from the space of possible model sizes is one such component (see also Section <ref>). When we perform forward search through the model space (see Section <ref>), the objective optimised at each of its steps reduces the risk of over-fitting by operating in terms of the difference between predictive distributions and not the submodels' predictive performance directly. Within the projection, 67 argue, the use of the reference model's complete posterior predictive distribution (or some sensibly clustered or thinned approximation of it) allows the submodels to better incorporate uncertainty in their own posteriors.

74 presented projective inference in relation to other practical model selection techniques (including sparsifying priors). Later, 72 introduced the so-called “clustered projection”, combined projective model selection with Pareto-smoothed importance sampling leave-one-out cross-validation <cit.>, and investigated reasons for the good performance of projective inference. In more complex settings, 73 showed the potential for the procedure's use in Gaussian processes, and 1 demonstrated its use to create more explainable tree models. 53 motivated an application of the procedure to identify the order of auto-regressive moving-average models. Finally, 12, 11, and 104 demonstrated how the procedure can be elevated beyond the observation models used in generalised linear models (GLMs), allowing for even broader application.

§.§ Contributions and structure of this paper

An efficient projection predictive workflow. The nodes shown in [fill=paleaqua]blue relate to fitting and diagnosing a reference model, and are covered in Section <ref>. Once we have this reference model, we perform an initial full-data search (the node shown in [fill=palegoldenrod]yellow), which we deal with in Section <ref>. We then diagnose this initial search, and either cross-validate over multiple search paths if it does not pass our checks, or move directly to selecting a submodel size. These steps are shown in the [fill=palechestnut]red nodes, and will be discussed in Section <ref>. Finally we discuss in Section <ref> what can be done with the projected posterior (the final node shown in [fill=palegreen]green). Steps that are not always required are indicated by dashed arrows.

Unlike previous papers that have primarily sought to establish the theoretical foundations of projection predictive inference, we survey the contemporary landscape and detail an efficient workflow for practical projection predictive inference. We aim to bring structure and clarity to the model selection workflow of statisticians, especially those using Stan, and introduce the notion of projection predictive inference to the uninitiated.

Concretely, we:

* decompose the projective model selection workflow into clear and modular components (Figure <ref>);

* propose heuristics to achieve the most robust results possible with the least computational effort;

* provide an interpretation of the “projected posterior” (the distribution resulting from the projection of the reference model's posterior onto the parameter space of a given submodel), and its limitations in a causal setting.

Section <ref> recapitulates the role of a reference model in model selection, before we present the underlying theory of projection predictive inference specifically in Section <ref>; Section <ref> deals with the search component of our efficient workflow which is then followed by model validation and selection in Section <ref>; Section <ref> discusses how one might interpret the projected posterior we achieve through projection predictive inference. Section <ref> illustrates our workflow in various simulated and real-data case studies; finally Section <ref> completes this paper with further discussion on the role of projection predictive inference in a Bayesian workflow.

§.§ Motivating examples

To make the discussion of the properties of projection predictive inference easier to follow, we use a couple of examples throughout the paper. These examples use model selection to different ends, and have data that accentuate certain aspects of the workflow.

The first of these examples is a dataset used to predict lean body weight given weights, body measurements, and ages of $n=251$ individuals [69]. Since accurate measurements of body fat are costly, we wish to predict it using $p=13$ measurements. Some of these measurements require tape measurements of the patients, and we would like to minimise the human cost as much as possible. Thus our task with regards to model selection is purely predictive: identify the minimal set of these $13$ attributes required to accurately predict body fat. Namely, there is no causal interpretation to the results.

The second data provided by 18 relates Portuguese student performances in mathematics and Portuguese language to their demographics and various social and school factors ($n=395, p=30$). This is a purely observational study, making causal analysis difficult, but learning which predictors are relevant to prediction can help design an intervention study for causal analysis later. As such, model selection in this case is primarily used to identify which external factors best predict the requirement for teacher intervention.

§ THE ROLE OF A REFERENCE MODEL

The rich reference model in question is the best-performing (in terms of its predictive performance) model we have at our disposal, and is one we would be happy to use as-is. Model selection then becomes applicable when we

[label=(*), ref=*, nosep]

* have a predictive rich model, but would like to reduce computational burden,

* would like to use a model more robust to changes in the data-generating distribution, or

* would like to gain a better understanding of important correlation structures.

Given the importance of this model, it is reasonable to begin by describing how one might construct one, and the alternative uses of a reference model outwith projection predictive inference.

§.§ Building a reference model

When diagnosing a reference model, there are three primary dimensions we recommend the statistician to investigate:

* posterior sensitivity to the prior and likelihood;

* posterior predictive checks;

* cross-validation and the influence of data on the posterior.

Corresponding discussion on model diagnosis and visualisation therein has been covered by 43, 30, 34, and 99. The priorsense [42], bayesplot [29], and loo [98] packages are free resources in R to such ends, and ArviZ [47] provides an analogue in Python.

Our reference model should be fit with priors that allow for complexity without over-fitting the observed data. This can include sparsity-inducing priors such as the spike-and-slab [56], regularised horseshoe priors [75], or the $L_1$-ball prior [105]. More recently, we have witnessed the formation of joint priors over regression coefficients and model fit (namely $R^2$) such as the R2D2 [106, 112] and R2D2M2 priors [2]. Such constructions allow for the easy addition of predictors to the model in the knowledge that the joint prior behaves sensible with increasing dimensionality, which helps to avoid over-fitting. These priors can be verified with prior predictive checking as described by <cit.>, wherein we simulate from the prior distribution and investigate its implications for model predictions in the response space. We hope to see that simulating from our priors generates data which seem reasonable and inline with our intuition and observed data.

Once we have defined some reasonable priors, we look to diagnose the posterior predictions of our reference model. These checks are in many ways related to the prior predictive checking of before. Here, we simulate from the posterior predictive distribution of our fitted model to check that the simulated predictions resemble the observed data sufficiently well. This can mean sufficient coverage of the posterior predictive simulations, sufficient concentration of the posterior predictive simulations, or sufficiently close moments of the posterior predictive simulations to the observed data. For more discussion surrounding posterior predictive checks, the reader is advised to consult the work of <cit.> and <cit.>.

Finally, we can investigate the effect of individual data to the posterior of our reference model. Doing so allows us to better understand which regions of data are better represented by our model, and can also act as a sanity check for model misspecification. Since posterior predictive checks “double-dip” our data by using it both for fitting and diagnosis, it is liable to be overly optimistic in its conclusions. Cross-validation shows how the model would behave on unseen data. For example, investigating the leave-one-out (LOO) posterior predictive distributions by comparing the distribution of the LOO probability integral transform (LOO-PIT) values at the observed data to the uniform distribution which would hold under a calibrated model can reveal over- or under-dispersion of the LOO posterior predictive distributions [30, 34].

Further, PSIS-LOO-CV [102, 99] can provide insight into the behaviour of our cross-validation approximation. 102 and 99 claim that the LOO-CV estimator behaves well (meaning that with high probability, the estimation error is small) when the shape parameter of the Pareto distribution fit to the tail of the importance weights is $\hat k < 0.7$. As such, observing many Pareto $\hat k > 0.7$ values may be indicative of a potentially misspecified model. Such PSIS heuristics should be considered only as one tool in the broader arsenal of model diagnostics, but can nonetheless provide valuable information to the practitioner.

For the case examples discussed in Section <ref>, we find that it is sufficient to fit linear regressions to the response variables (respectively the body fat of an individual, and the Portuguese students' final mathematics grade on a scale of $0$ to $20$) including all available predictors. Since we understand that not all of these predictors will likely be relevant, we fit the model with an R2D2 [112] prior. More details of this prior and the specific hyperparameters used can be found in Appendix <ref>. We find that these reference models pass our checks and were it not for the need for model selection, we would be happy to use them as-is.

§.§ Using a reference model

While we interest ourselves primarily in the use of reference models as a predictive benchmark for projections, they can be used more generally as approximations of the underlying data-generating distribution. We presently discuss other use cases of reference models before returning to their role specifically in a projective model selection workflow.

Indeed, reference models can be combined with many other model selection procedures to increase selection stability, and 67 indicate that the more complex the models considered and the fewer data there are, the more valuable a reference model is in filtering some of the noise present in data. For instance, reference models can be used in minimal subset selection as part of a Bayesian step-wise selection strategy [67]. This might involve beginning with the reference model and at each step, excluding the predictor with associated regression coefficient $\theta$ having the highest Bayesian $p$-value. The selection continues for as long as the reduced model has a better estimate for the expected log (pointwise) predictive density (for a new dataset) <cit.> than the current model. 67 discuss how reference models can also be used in complete variable selection through either local false discovery rate minimisation [23, 24], empirical Bayesian median analysis [41], or based on posterior credible intervals [67].

§ PROJECTION PREDICTIVE INFERENCE

Given a reference model, we wish to achieve a parsimonious restricted model whose posterior predictive distribution most closely resembles that of the reference model.[Unless stated otherwise, the term “restricted model” will always refer to a candidate model that the reference model has been projected onto. Thus, the associated “posterior predictive distribution” is in fact a predictive distribution based on the projected posterior.] This resemblance is the objective of the projection.

§.§ Projection

To define the projection mathematically, we denote the number of observations by $n$ and the observed response vector by $y = (y_1, \dotsc, y_n)^T$.
Predictor data will always be conditioned on implicitly.
By $\tilde{y} = (\tilde{y}_1, \dotsc, \tilde{y}_n)^T$, we denote unobserved response values at the same observed predictor values. The KL divergence from a distribution $p(y)$ to a distribution $q(y)$ is written as $\mathbb{KL}\left\{p(y) \,\lVert \, q(y)\right\}$ and given by $\E{~p(y)}{\log \frac{p(y)}{q(y)}}$.

Letting $\Theta^\ast$ denote the reference model's parameter space, we assume the availability of $\mathcal{S}$ posterior draws $\theta^{\ast(s)} \in \Theta^\ast,\,s \in \{1, \dotsc, \mathcal{S}\}$.
To save computational resources, 72 suggest to cluster (or thin) these $\mathcal{S}$ posterior draws so that we may write $\{1, \dotsc, \mathcal{S}\} \supseteq \mathop{\bigcup}_{c = 1}^C \mathcal{I}^\ast_c$ with disjoint index sets $\mathcal{I}^\ast_c$.
For a given $c \in \{1, \dotsc, C\}$ and $i \in \{1, \dotsc, n\}$, we then slightly abuse notation in defining

\begin{equation}
p(\tilde{y}_i\mid\mathcal{I}^\ast_c) = \frac{1}{|\mathcal{I}^\ast_c|} \sum_{s \in \mathcal{I}^\ast_c} p(\tilde{y}_i\mid\theta^{\ast(s)}),
\end{equation}

where $p(\tilde{y}_i\mid\theta^{\ast(s)})$ denotes the reference model's predictive distribution for observation $i$, conditional on posterior draw $\theta^{\ast(s)}$.
Accordingly, expectations with respect to $p(\tilde{y}_i\mid\mathcal{I}^\ast_c)$ are denoted by $\E{\!}{\, \cdot \mid \mathcal{I}^\ast_c \,}$.

Given a restricted model with parameter space $\Theta_\perp$ (which may be completely unrelated to $\Theta^\ast$) and parameter-conditional predictive distribution $q(\tilde{y}_i\mid\theta_\perp)$ (for $\theta_\perp \in \Theta_\perp$ and $i \in \{1, \dotsc, n\}$), the projection problem for a given $c \in \{1, \dotsc, C\}$ may then be defined as

\begin{equation}
\theta_\perp^c = \argmin_{\theta_\perp \in \Theta_\perp}
\frac{1}{n} \sum_{i = 1}^n
\mathbb{KL}\left\{\, p(\tilde{y}_i\mid\mathcal{I}^\ast_c) \,\lVert\, q(\tilde{y}_i\mid\theta_\perp) \,\right\}, \label{eq:practical-projpred-orig}
\end{equation}

which simplifies to [72]

\begin{equation}
\theta_\perp^c = \argmax_{\theta_\perp \in \Theta_\perp}
\sum_{i = 1}^n
\E{\!}{\, \log q(\tilde{y}_i\mid\theta_\perp) \,\mid\, \mathcal{I}^\ast_c \,}. \label{eq:practical-projpred}
\end{equation}

By solving this projection problem for each $c \in \{1, \dotsc, C\}$, we obtain a collection of parameter “draws”[These “draws” are deterministic transformations of the reference model's posterior draws.] for the submodel that achieves the desired resemblance (in terms of predictive distribution) between the submodel and the reference model. Indeed, in the case of exponential family models, and in particular (hierarchical) GLMs, this objective can be efficiently computed by (marginal) maximum likelihood using, for example, the lme4 package [5] in R. For implementation details see the papers by 12, 13, 104, 72. In case of actual clustering (not thinning), the only difference compared to traditional MCMC draws is that in any post-projection inference, the cluster weights need to be taken into account.

§.§ An efficient projection predictive workflow

Having provided a theoretical motivation of projections, we present the primary contribution of this paper: an efficient workflow for projection predictive inference. The steps inherent to the workflow are visualised in Figure <ref>:

* begin with a reference model, which we will assume to be well-considered and appropriate to the task (see Section <ref> for more discussion);

* perform an initial search through the model space, or some heuristic subset of it, without cross-validation and using all available data to achieve a series of nested, increasingly complex submodels deemed the solution path (Section <ref>);

* diagnose this solution path to detect possible over-optimism using the cross-validated predictive performance of the submodels along this solution path (Section <ref>);

* if the initial solution path is deemed to be over-optimistic, then in a second stage include the search in the cross-validation, and then combine the resulting cross-validated performance evaluations of models along each fold-wise solution path (Section <ref>);

* select the minimal submodel size capable of producing similar predictive performance to our reference model (Section <ref>);

* finally, diagnose whether inference directly with our projected posterior (if such inference is desired) is reasonable with calibration tests (Section <ref>).

Before any of this, however, it is reasonable to ask whether model selection is required at all, and whether the projection predictive inference is the flavour best-suited to one's needs.

§.§ When model selection is not the answer

In a causal setting, it might be tempting to perform model selection in order to increase the precision of a treatment effect posterior of interest [16]. However, to achieve valid causal inference, one must account for the causal relationship between the response and measured predictors. This is done by including in the model what 68 calls the admissible set of predictors. When searching in the model space, all submodels must include this admissible set. If our main concern is to accurately infer a causally valid treatment effect, then no model selection is needed and we can use the reference model directly, provided it includes the admissible set.

It may also be unnecessary to involve projection predictive inference in one's model building workflow if one is dealing with relatively few (fewer than ten, say) models, especially if these models are not nested. In this case, it may be easier to compare the models either by estimating the difference in their predictive performance directly [99, 98], or by investigating their Bayesian model stacking or averaging weights [110, 109, 108].

For predictive tasks, the reference model can always be used as-is, so long as there is no need to reduce its complexity (e.g., due to high observation or “measurement” costs of the predictors).

Finally, if one is unable to achieve a suitable reference model, then projection predictive inference can be immediately discounted. And if the reference model can be constructed but is not very complex, then further simplification of it may not be worth the effort.

Crucially, we reiterate that we do not consider model selection to be a remedy to over-fitting, and we do not recommend it for such aims. Instead, we recommend that model iteration continue inline with the recommendations of 34 and that stronger joint priors over the predictive space be used.

If model selection is indeed performed, then it is important that any post-selection inference accounts for the model selection. Otherwise, post-selection inference would be prone to over-fitting. The projective part of projection predictive inference ensures that all predictive uncertainty inherent to the reference model is retained in the submodel, thus accounting for the model selection (apart from the decision for a submodel size, although based on the extensive experi ents that bias is negligible).

Returning then to our examples, we recall that the body fat analysis is purely predictive. As such, if measurement was zero cost then there would be no need for model selection, but as it stands we can benefit from reducing data collection cost. Likewise, we make no causal claims in the Portuguese student example.

§ SEARCH HEURISTICS

Decision tree diagram showing a forward search through (a heuristic subset of) the model space, starting with the empty model in the left-most node, and progressing through the path by minimising the KL divergence from the reference model's posterior predictive distribution (possibly coarsened by clustering or thinning the posterior draws) to the restricted model's posterior predictive distribution.

Once we have built our reference model and identified that the projection predictive inference fits our use case, we begin by producing some set of models to compare our reference model to. It is this search stage which produces our solution path, and there are two primary search methods we consider: KL divergence-based forward search and Lasso ($L_1$) regularisation search. These two have useful speed-accuracy balance and afford the possibility of an efficient stopping rule, but other algorithms can also be used.

§.§ Forward search

In forward search, we begin by projecting the reference model onto the “empty” (intercept-only) model, which acts as the root of our search tree. We then project the reference model onto all models with one predictor and the intercept term, and select the single-predictor model whose posterior predictive distribution is closest to the reference model's posterior predictive distribution (possibly coarsened by clustering or thinning the posterior draws) in KL divergence sense and as defined in Section <ref>. Denote this first predictor to be selected $x^{(1)}$. Following this, we fit all size-two models including the intercept and $x^{(1)}$ (“size-two” does not count the intercept here), and once more select the one closest to the reference model in terms of KL divergence of their posterior predictive distributions. Denote this second predictor to be selected $x^{(2)}$. This is repeated until either all predictors are selected, or some pre-defined limit on the model size is reached. A visualisation of this process is provided in Figure <ref>. Our solution path is then the list of predictors ordered by the stage at which they were selected: in our example $(\intercept, x^{(1)},x^{(2)},\dotsc)$.

We favour forward search over backward search for the ability to perform early-stopping along the path (i.e. only search up to some number of predictors fewer than the size of the reference model) starting with the least computationally expensive models (while backward search begins with the most computationally expensive calculations).

§.§ Lasso-type search

72 also discuss how in the case of single-cluster projections, one can reduce computational cost further by performing a Lasso-type search [93, 115], also called $L_1$ search. Concretely, this involves fitting a model with Lasso penalty to the in-sample predictions made by the reference model <cit.> and investigating the order in which the predictors enter the Lasso model. In our experience, Lasso-type search is computationally faster than forward search, but has higher variability in the search results (with that variability referring to variability across data realizations, as can be seen, e.g., by bootstrapping).

§.§ Output of the search heuristic

After our projection predictive workflow has identified a predictor ordering (the solution path), either through forward or $L_1$ search, we store this ordering so as to avoid repeated computations. The predictor ordering alone can also be interesting to the practitioner, who might identify some intuition about their problem from this alone. In a forward search, we can store (and later access) the KL divergence from the reference model to each submodel, so it might be the case that we can observe an elbow when plotting the KL divergence along the solution path, indicating some saturation of information being provided by predictors. In our experience, these elbows are usually very soft and difficult to infer from. In Section <ref>, predictive performance utilities are reasoned on instead of KL divergence for picking a submodel size.

§.§ Modularity of the search heuristic

Body fat example. The elpd point estimate and one standard error bars of submodels along the full-data solution path under forward search (in black) and $L_1$ search (in [fill=blue!20]blue) relative to the reference model (shown in [fill=red!20]red). The two methods identify different predictor orderings, leading to different elpd values along the solution path. In particular, forward search tends to produce smoother such elpd curves than $L_1$ search. We have omitted the submodel of size zero (the intercept-only model) as it is much worse than the submodel of size one, and is the same under both search methods trivially.

Body fat example. A comparison of the occurrence of predictors over multiple LOO-CV selection folds when using forward search and $L_1$ search heuristics. Shown are cumulative rates, i.e., the proportions of CV folds that include a given predictor from the $y$-axis at model sizes smaller than or equal to a given model size from the $x$-axis. In this case, the search paths are quite stable across the CV folds as most of the non-zero cumulative rates are 1. In case of collinearity of predictors, the used search method can have a large effect on the search order, but it is unlikely to have a big change in the predictive performance of the selected model.

Our workflow remains a modular initiative. Indeed, we do not require or assume that any specific search heuristic is implemented to continue.

Other search heuristics can be just as easily implemented, including exhaustive search [31], backward search [61], and stochastic search [35, 64]. Our general workflow can admit any of these without modification, although we have chosen the above two as candidates which remain efficient and stable in selection.

In particular, while Lasso search may achieve more stable paths over CV folds, forward search often achieves more predictive submodels and converges to the performance of the reference model faster. We have noted empirically that the Lasso path is also liable to fail in situations of high block-wise correlation. In general, it is possible for the two methods to identify different predictor orderings; in Figure <ref> we visualise how these different ordering are reflected by the different predictive performance of submodels along the two paths in our body fat example. We can investigate exactly which predictors differ in the ordering in Figure <ref>. Indeed, we find that in the body fat example, forward search is slightly less stable across cross-validation folds compared to lasso search. It is also clear from this figure that the order with which predictors enter the solution path is very different, hence why the elpds of smaller submodels differ between methods.

§.§ Computational efficiency of search heuristics

Since at each step in forward search, we project the reference model onto a subset of predictors, the complete forward search requires $\frac{p (p + 1)}{2} + 1$ projections if there are $p$ predictor terms (excluding the intercept).

We can perform the search using $C$ draws (or clusters), and then re-project the submodels along the solution path using more draws (or clusters), say $C^\prime$ (where $C' \gg C$) when evaluating the predictive performance of the submodels. If it takes $t_{\text{proj}}$ seconds to perform the projection for one model in the forward search using $C$ projected draws (or clusters of posterior draws), then a second projection using $C'$ projected draws will take approximately $\frac{C'}{C} \cdot t_{\text{proj}}$ seconds. In our experience, it is usually sufficient to perform the search with $20$ clusters, to perform the performance evaluation along the solution path with $400$ thinned posterior draws (not clusters), and finally to use all posterior draws when projecting onto the final selected model. This way, we optimise the balance between stability and efficiency.

We recommend one use forward search over Lasso search for final decision making. However, this is not to say that the latter is not a useful tool: when the number of predictors in the search space is large (greater than $40$, say), $L_1$ search is significantly faster to compute; one can use it to ascertain a heuristic upper bound on the submodel size, and then re-run forward search only up to this bound.

Naturally, performing the projections using more draws or clusters can improve the stability of the solution path at an increased computational cost. With extreme thinning, variability in model selection can increase. This can be reduced by smarter thinning <cit.> or ad hoc clustering. We note, however, that the cluster centroids may not be the best possible approximation to the posterior. As such, one may inadvertently introduce a small amount of bias into the model size selection by using clusters in cases of difficult (for example multi-modal, skewed, or thick-tailed) parameter posteriors. On the whole, this bias will likely remain small since we use the same draws over all projections (in a word, we repeat the same mistake at each step), so that the bias tends to be in the same direction intuitively.

§ SELECTION HEURISTICS AND DIAGNOSIS

Having achieved a predictor ordering in the search phase, we turn our attention to the predictive performance of submodels along the path. While in principle, this performance evaluation can be run on the same data that the search was run with, it is more sensible to run the performance evaluation on new data, typically leading to cross-validation. In the ideal case, the whole procedure (search and performance evaluation) is included in the cross-validation, meaning that the search is performed with the training data of each CV fold separately and the performance evaluation is performed with the test data of the respective CV fold. As this is computationally expensive, it makes sense to run a full-data search first, and approximately cross-validate only for the evaluation part using PSIS-LOO-CV. However, it may be that this is not sufficient, and in this case we include the search in the cross-validation.

In this section, we first detail two methods of cross-validating the predictive performance of individual models, before discussing how one might conduct cross-validation for the entire procedure should we over-fit the data in the search path.

§.§ K-fold cross-validation

In the $K$-fold-CV paradigm, we fit the reference model and repeat the model selection $K$ times [72]. We then estimate a pre-specified utility (in practice we interest ourselves in the elpd and thus it is used as an example here) of the model at position $k \in \{0, 1, \dotsc\}$ of the solution path and the reference model at observation $i$, denoted $u_k^{(i)}$ and $u_\ast^{(i)}$, respectively. The utility differential between a submodel of complexity $k$ and the reference model is then estimated with

\begin{equation}
\Delta\bar{U}_k = \sum_{i=1}^n(u_k^{(i)} - u_\ast^{(i)}).
\end{equation}

§.§ Leave-one-out cross-validation

The Bayesian LOO-CV estimate for a model's elpd on $n$ data observations is defined as

_    = ∑_i=1^nlogp(y_i|y_-i)

= ∑_i=1^n log∫p(y_i|θ)p(θ|y_-i) dθ,

where $y_{-i}$ denotes all observations omitting the $i^\text{th}$. In PSIS-LOO-CV, rather than naïvely computing this integral $n$ times (once for each fold, which requires refitting the reference model $n$ times), we instead estimate the posterior predictive distribution $p(y_i\mid y_{-i})$ based on the importance sampling estimate

\begin{equation}
p(y_i\mid y_{-i}) \approx \frac{\sum_{s=1}^\mathcal{S}p(y_i\mid\theta^{(s)})r(\theta^{(s)})}{\sum_{s=1}^\mathcal{S}r(\theta^{(s)})},
\end{equation}

wherein the weights takes the form

\begin{equation}
r(\theta^{(s)}) \propto \frac{1}{p(y_i\mid\theta^{(s)})}.
\end{equation}

In PSIS-LOO-CV, these importance weights are stabilised and diagnosed with Pareto smoothing as previously discussed in Section <ref> [102]. The primary advantage of PSIS-LOO-CV is that it does not require us to repeatedly refitting the reference model and repeating the search. When the number of observations or the complexity of the model grows, even this estimate can be expensive to compute and may occasionally fail. In this case, we fall back on $K$-fold-CV instead for the predictive utility.

§.§ Cross-validation over solution paths

Portuguese student example. We show point estimates for the elpd difference to the reference model of both the full-data search in [fill=gray!20]grey, and the cross-validated search in black, along with one standard error bars. When we include the search in the cross-validation, we find that instability in the predictor ordering leads to “jumpy” behaviour in the elpd differences. We smooth the elpd difference estimates with a monotonic spline shown in the [fill=blue!20]blue line, where the enveloping ribbon represents the one standard error interval around the point estimate. We do not include the null model in the spline, since we understand that the increase in predictive performance is so great that the intercept-only model is of little importance to later analysis.

Once we have cross-validated the predictive performance of the submodels along our initial full-data solution path, we can then immediately diagnose selection-induced over-optimism in the elpd estimates. This is apparent when we observe a “bulge” in the elpd of submodels along the solution path. Specifically, if a submodel smaller than the reference model has elpd much higher (better) than the latter, we can conclude that the cross-validation estimates for full data solution path are over-optimistic. Recall the predictive performance along the two solution paths from  Figure <ref>, where we see that both initial solution paths have over-optimistic performance estimates to the point where models of more than roughly five predictors achieve better elpd than the reference model of 13 predictors. Moreover, the over-optimistic submodel elpds are in some cases “significantly” better than that of the reference model.[By “significant”, we mean that the estimated elpd difference is at least one standard error larger than zero.]

In this case, understanding that we are witnessing over-optimism in the elpd estimates along the full data solution path, we should cross-validate the models along the solution path and cross-validate the model search. Due to the vastly increased computational requirements, we suggest that the search be re-run only up to the model complexity inducing the most over-optimistic elpd estimate (for example in Figure <ref>, we might re-run up to model size nine or ten to cut down some computational cost). Doing so achieves submodel performance evaluations that have smaller bias [74].

Consider now the Portuguese students motivating example. In Figure <ref> we see how the elpd estimates along full-data solution path are over-optimistic, inducing a bulge at model size roughly 12 which gradually diminishes to recover the elpd of the reference model. By including the search in the cross-validation, the bias in the elpd estimates is reduced. However, in doing so, we introduce some “jumpy” behaviour in the elpd differences at some model sizes. This is caused by instability across cross-validation folds in the predictor ordering.

§.§ Computational efficiency of cross-validation

The initial full-data search lends itself naturally to an estimate of the time required to perform $K$-fold-CV over the whole procedure (both search and evaluation): suppose each projection uses $C$ projected draws and takes $t_{\text{proj}}$ seconds. Then running the search for each of the $K$ CV folds with the same number $C$ of projected draws will take approximately

\begin{equation}
t_{\text{search}} = K \cdot \left(\frac{p (p + 1)}{2} + 1\right) \cdot t_{\text{proj}}
\end{equation}

seconds if not run in parallel and not taking into account the time required for refitting the reference model in $K$-fold-CV. Since a CV including the search is always a CV for predictive performance evaluation on the test data of each CV fold, it makes sense to combine this estimated search time with the estimated time for the corresponding performance evaluation based on $C'$ projected draws, which is

\begin{equation}
t_{\text{eval}} = K \cdot (p + 1) \cdot \frac{C'}{C} \cdot t_{\text{proj}}
\end{equation}

seconds, giving an estimate of

\begin{equation}
t_{\text{total}} = t_{\text{search}} + t_{\text{eval}}
\end{equation}

seconds in total. These heuristics assume that one wishes to perform an exhaustive search over the model space. That is, we search all model sizes up to the size of the full model, $p_{\text{max}} = p$.

As was previously discussed, we can use a smaller number of clusters in the search phase to achieve a solution path and then re-project the models along that path using more draws in order to achieve more stable cross-validation results. Typically, we have found that using approximately $20$ clusters in the search phase, and $400$ thinned posterior draws for cross-validation provides a good balance of stability and efficiency. We do not recommend using clustered projections for performance evaluation due to the risk of introducing bias as previously discussed in Section <ref> <cit.>.

Having then attained a cross-validated performance estimate for the model sizes along the solution path, we look to select the minimal set of predictors replicating the performance of the reference model to a sufficient degree.

§.§ Choosing a submodel size

We are now faced with the final decision of identifying the smallest submodel (minimal predictor set) whose predictive performance most closely resembles that of the reference model. For this task, we propose two primary heuristics: one based on differential utility intervals, and another on the mean elpd difference to the reference model.

The former is most common, and has been previously considered by 74, 72, 12, 71, and 104. In this case, we denote the utility of the reference model as $u_\ast$, and the utility of the submodel with $k$ predictors as $u_k$. Having cross-validated the predictive performance of the submodels, we then select the smallest submodel whose utility is one standard error less than the reference model's utility, with the standard error being that of the utility difference, denoted $s_k$. Formally, we select the smallest $k$ such that

\begin{equation}
\hat{u}_\ast \leq \hat{u}_k + s_k,\tag{SE}\label{eq:utility-ci}
\end{equation}

The second heuristic, motivated by 85, is to choose the smallest $k$ such that the estimated elpd of the submodel is at most four less than the estimated elpd of the reference model. That is, we choose the smallest $k$ such that

\begin{equation}
\hat{u}_k - \hat{u}_\ast \geq -4. \tag{$\Delta$utility}\label{eq:delta-elpd}
\end{equation}

We have found through our own experiments that the <ref> heuristic results in more stable selection, and selects better-performing submodels (results not shown in this paper). We have also observed that it may select larger models than the <ref> heuristic. In general, however, these have been noted to behave similarly.

In certain “large $n$, small $p$” regimes, it may be that even following a cross-validation including the search, the elpds at the model sizes along the path can be erratic and uncertain. This will be seen later in Figure <ref> in our second case study of Section <ref>. In order to stabilise selection of submodel size in these cases, and reduce the noise in the solution path, we smooth the submodels' estimated elpds. One well-established tool for smoothing a collection of points $\{(x_i,y_i)\}_{i=1}^n$ are B-splines [78, 27, 25, 26]. These are smoothly joined polynomial segments, knotted together at some (usually regular) intervals that allow us to model some functional form. For instance, the true elpd differences along model sizes are typically monotonically increasing and negative. In practice, we use the scam package [77] to compute these splines.

Usually the intercept-only model has so much worse predictive performance  than the model with just one additional predictor, that it is usually best to omit the intercept-only model (and even occasionally the size-one model) when fitting the spline. Further, we fit these splines to the normalised elpd difference values (elpd difference divided by its standard error), such that accuracy of each elpd difference estimate is taken into account, too. The smoothing of the elpd values is illustrated in Figure <ref>.

Once a minimal predictor set is identified, we may find that one or more predictors could be substituted with others to achieve very similarly performing submodels of the same complexity but at a lower measurement cost. Such instances can be diagnosed by inducing the distance metric between predictors with regression coefficients $\theta, \theta^\prime$,

θθ^' = 1 - |(    ỹ|ℳ_θ, ,

where $\E{}{\tilde{y}\mid\mathcal{M}_{\theta}, \data}$ are the in-sample predictions made by the single-predictor model, $\mathcal{M}_\theta$, containing only the intercept term and $\theta$ [65]. Computing then the pair-wise distances as defined in Equation <ref> between all predictors, we can achieve such dendrograms as in Figure <ref>. Specifically, we understand that should we identify a submodel admitting an individual's weekend alcohol consumption in the predictor set (named “walc” in the plot), but not their weekday alcohol consumption (“dalc”), that similar results might be achieved by including the latter instead of the former. This can be useful for building predicive models without explicit costs for the measurements. For instance, should it be found that similar predictions can be achieved in a medical setting through non-invasive measurements instead of invasive measurements, then the less intrusive measurement could be preferred.

Portuguese student example. Dendrogram of the predictors in the reference model, using the correlation of single-predictor models (Equation <ref>) as the distance metric. Such dendrograms can reveal predictive similarities between predictors and aid intuition in variable subset selection.

And so concludes projection predictive model selection: we have identified the minimal predictor set and identified correlation structures between them. What remains to be discussed is what can be done with this information.

§ INTERPRETING THE PROJECTED POSTERIOR

Simulated linear regression example. An example of overconfidence in the projected posterior of a linear regression on $95$ independent predictors of which only $15$ are relevant. The first three plots show a subset of marginal posteriors (namely $\beta_1, \beta_2$, and $\sigma$) from a reference model fit with priors that are liable to over-fit the data (in black), and a projected submodel including only the relevant predictors (in [fill=blue!20]blue). The dashed vertical lines indicate the true value of the parameters. The fourth plot depicts a kernel density estimate over the posterior predictive mean at each observation for each model, denoted $\tilde{y}$, with the shaded region indicating the observed data density. The $\beta_1$ and $\beta_2$ marginals of the projected posterior concentrate more sharply than in the reference model, and closer to the true parameter values. The additional variance in the $\beta_1$ and $\beta_2$ marginals of the reference model's posterior (compared to the $\beta_1$ and $\beta_2$ marginals of the projected posterior) is projected onto $\sigma$, inflating it.

The primary output of projection predictive inference is a collection of projected parameter draws for a submodel consisting of a minimal set of predictors that induces a predictive performance close to that of the reference model. However, it is interesting to understand under which circumstances valid post-selection inference can be done for these predictors using the marginals of the projected posterior directly.

First and foremost, should the aim of the analysis be to make inference on a treatment effect in a causal setting, we reiterate that the practitioners should use the reference model directly. If one chooses to use the projected posteriors, then they subject themselves to certain risks: the miscalibration of the projected posterior distributions due to model selection, and the projection of structured variance from the reference model to the unstructured dispersion parameter in the submodel.

Common folklore states that when using step-wise selection, credible intervals are not well-calibrated. However, using forward search method based on KL divergence, and cross-validation thereof, there is less selection-induced over-fitting. Furthermore, for the inference for the treatment effect to be well-calibrated, it is important to have a sensible reference model (which we show empirically in the case study of Section <ref>).

In certain instances, such as our motivating body fat example, we are only interested in a sparse model with good predictive performance. In the Portuguese student example however, we are also interested in the magnitude and sign of each coefficient, but can't make causal conclusions because of collinearity between measurements. In the presence of two collinear and truly relevant predictors, the removal of one will result in its effect being projected to the other. This will increase the effect of the remaining predictor to account for the removed correlated effect of the other. As such, we understand that our projection may distance parameter posteriors away from their theoretic “true” value in order to preserve the predictive performance of the larger reference model including collinear predictors.

As such, we discourage the reader from making causal decisions based on the projected posteriors in the presence of collinear predictors. Indeed, it is not clear that this issue disappears with the removal of collinearity in general. As such, we repeat that causal decisions should only be made based on the posterior parameter distributions of the reference model (assuming also the reference model is  causally sensible).

The other risk is that when the reference model is projected onto a submodel, any additional structured variation captured by the reference model is projected onto the unstructured dispersion parameter of the submodel, thereby inflating it [72]. Furthermore, in the case of an over-fitting reference model, we will demonstrate that the dispersion parameter may also absorb some of the posterior variance from the regression coefficients, leading to their overconfidence.

Consider a linear regression on $n=100$ observations of $95$ predictors of which only $15$ are relevant. We fit a reference model using independent Gaussian priors, which we know to be prone to over-fitting in these regimes. In Figure <ref>, we then project this reference model

onto only the set of truly relevant predictors. What we find is that the marginal posteriors of the regression coefficients (here we only show the first two) concentrate more sharply than in the reference model, and their means approach the true parameter values.[This convergence towards the true values is likely because we have projected onto the known true model, and is not the focus of our procedure.] The variance in the reference model's marginal posteriors of the regression coefficients is then transferred to the magnitude of $\sigma$. We note that this behaviour is to be expected, and is not a failure of the procedure. Indeed, this arises directly from the projection objective being focused on predictive inference with respect to the fit of the reference model.

As previously discussed in Section <ref>, when performing inference based on a projected submodel, we recommend that the user perform the draw-by-draw projection (i.e. $C = \mathcal{S}$).

§ CASE STUDIES

Presently we manifest the recommendations laid out in this paper through several simulated and real-data examples. The below case studies were carried out using the projpred package in R [71].[The code is freely available at <https://github.com/yannmclatchie/projpred-workflow/>.]

§.§ Highly correlated predictors

In the first case study, we employ projection predictive inference on simulated data where predictors are highly correlated. The setup is the same as in <cit.>, with data generated according to

x     ∼(0,R)

y     ∼(w^Tx,σ^2),

where the number of predictors is set to $p = 100$, $\normal(\mu, \Sigma)$ denotes a $p$-dimensional normal distribution with mean vector $\mu$ and covariance matrix $\Sigma$. The matrix $R \in \mathbb{R}^{p\times p}$ is block diagonal, each block being of dimension $5 \times 5$. Each predictor has mean zero and unit variance and is correlated with the other four predictors in its block with coefficient $\rho=0.9$, and uncorrelated with the predictors in all other blocks. Further, the weights $w$ are such that only the first $15$ predictors influence the target $y$ with weights $(w^{1:5},w^{6:10},w^{11:15}) = (\xi,0.5 \xi,0.25 \xi)$ and zero otherwise. We set $\xi = 0.59$ to fix $R^2 = 0.7$, and $\sigma^2=1$. We simulate $n=500$ data points according to this data-generating process (DGP).

Highly correlated predictor simulated case study. A dendrogram of the predictors (with associated regression coefficients $\theta$) present in the reference model, where the distance recorded on the $x$-axis is computed according to Equation <ref>, and where clusters containing chosen predictors have been coloured [fill=blue!20]blue.

Many weakly-revelant predictors case study. Over-optimstic predictive performance estimates based on a full-data search with many weakly-relevant predictors (in [fill=gray!20]grey), along with predictive performance estimates from PSIS-LOO-CV where the search has been included (in black) up to the lowest model size inducing the maximum degree of over-confidence (judged by elpd point estimate) in the original search. Since the cross-validated predictive performance estimates are very “jumpy”, we smooth them with a monotonic spline (shown in [fill=blue!20]blue) and perform our heuristic selection on that instead. Once more, we show the reference model's elpd by the horizontal dashed [fill=red!20]red line, and the selected model size in the vertical dashed [fill=gray!20]grey line.

Our reference model is the linear regression on all available predictors. To regularise our model, we choose an R2D2 prior on the predictor coefficients and the default priors in brms [9] for the intercept $\beta_0$ and the standard deviation $\sigma$ (see Appendix <ref> for more information on the construction of the R2D2 prior). As a result, the model passes basic computational diagnostics as well as basic model diagnostics, and we are comfortable using it as a reference model.[In this case, $\widehat{R}$ values are all smaller than 1.01 <cit.>, and bulk and tail effective sample sizes are high. All Pareto $\hat{k}$-values are satisfactory ($\hat{k}<0.7$) meaning  the second-order moment of the elpd estimate is likely finite, and the central limit theorem holds as a result [102].]

We begin by performing an initial forward search on the full data, where we find that there is an elpd "bulge" (not shown here, but similar to the full-data search path of Figure <ref>). That is, for some model sizes, the elpd point estimate and its one standard error interval lie completely above the point estimate of the reference model elpd, indicating potentially over-optimistic elpd estimates. We therefore include the search in the CV (and terminate the search at the model size inducing the highest over-opimistic elpd estimate, as wepreviously recommended), yielding more reliable elpd estimates and uncertainties for the different model sizes.

Since we have set up our experiment to have highly-correlated predictors, it is likely that our procedure was able to find a very small subset of predictors by choosing only one or two from each of the relevant correlated blocks, since they should repeat some predictive information. As such, we might ask whether any of the predictors in the minimal predictor set chosen could have been substituted with another predictor. To do so, we calculate the pairwise distances as previously described in Section <ref> and represent them with the dendrogram in Figure <ref>. We find that the predictors fall nicely into clusters of five, in line with the DGP. Further, it is likely that we can replace any of the selected predictors with any other predictor from the same cluster to achieve similar predictive performance.

This case study demonstrates that in the case of correlated predictors, projection predictive model selection is able to choose parsimonious submodels, and we can identify when predictors repeat or share predictive information with the distance metric we have previously motivated.

§.§ Many weakly-relevant predictors

Another situation not uncommon in statistical analyses is the presence of many predictors, each one contributing only a marginal amount of predictive information to the final model. Indeed, in situations where the data are already difficult to explain, having predictors that are only weakly-relevant can pose issues in model selection, and in particular in identifying a minimal predictor subset.

In this case study, we sample $n = 500$ observations from $p = 50$ predictors that are distributed according to a multivariate Gaussian, where all predictors are correlated with $\rho = 0.1$. As before, $y \sim \normal(w^Tx,\sigma^2)$, with $\sigma^2 = 1$ and $w$ computed and fixed such that $R^2 = 0.5$.

In these regimes, the predictor ordering in the solution path can be extremely erratic: since each predictor contributes so little predictive information, the search heuristic can be highly unstable and different runs may result in different orderings. As a remedy to this noise, we can smooth the difference of submodel elpds to the reference model elpd over the model sizes with a monotonic spline as motivated in Section <ref>.

This is seen in Figure <ref>, wherein the noise in the cross-validated predictive performance estimates for the submodels becomes immediately obvious. Indeed, we find that following the inclusion of the search in the cross-validation, we observe instability in the predictive performance estimates, likely caused by different predictor orderings, and instances where adding a new predictor can reduce the predictive performance estimate of the submodel. As previously discussed, this is merely an artifact of including the search in the cross-validation, and in general the addition of a new predictor should not worsen predictive performance on expectation. We find that smoothing these noisy elpd estimates produces more consistent interpretation, and we can select a model size by applying our <ref> heuristic to the spline point estimate – leading to model size $30$.

As such, even in regimes of weakly-relevant predictors, projection predictive inference is able to identify a saturation point in model size after which the addition of subsequent predictors do not add significant predictive benefit. Our smoothing technique also mitigates the noise induced by the predictors' predictive weakness.

§.§ Supervised principal component reference models

Simulated data case study. An example of miscalibration in the projected posterior in case of linear regression with independent predictors. The upper and lower rows compare the calibration of the treatment effect under the same model with two different priors: one using independent Gaussian priors for the regression coefficients, and a more predictive reference model using the R2D2 prior. The two columns then refer to the treatment effect in the reference model (on the left) and to the treatment effect after a projection of the reference model onto the selected minimal predictor set (on the right). A model is interpreted to be well-calibrated if the [fill=blue!20]navy line stays within the [fill=cyan!20]cyan envelope.

Although not often discussed, our reference model needs not be the largest of the set of nested models we search through. Indeed, it is entirely possible to search through a completely different model space than the one obtained by building increasingly complex models from the collection of the reference model's predictors. One such example is the use of principal components of the original predictor data as predictors in the reference model, whose fitted values are then used in the projection onto the original predictor space.

For example, we consider presently the task of predicting the presence of prostate cancer.[These data are freely available from <https://jundongl.github.io/scikit-feature>.] We observe $n = 102$ patients, each having $p = 5966$ highly collinear predictor measurements. To reduce the dimensionality of the data we compute $n-1 = 101$ iterative supervised principal components [76], and use them as predictors in our reference model. Having fit this reference model using the principal components and appropriate priors, we then project the predictive information back onto the original predictors so as to achieve an intuitive understanding of which measurements were most relevant in terms of predictive information. Using forward search in this case represents a large computational cost, so we first run $L_1$ search. In the case study example, this is sufficient as we did not observe a “bulge” in the predictive performance along the solution path (figure not illustrated here).

The use of principal components is then only to achieve a dimension reduction for the predictors in the reference model. This illustrates that we can project a reference model having a hard-to-interpret predictor space onto an interpretable predictor space. Naturally, one can equally use an additive Gaussian process [14], Bayesian additive regression trees [15], or any other non-parametric model to produce the reference model.

§.§ Understanding the effect of the reference model on projected calibration

We presently show how one might assess the calibration of a projected posterior after variable selection using simulation-based calibration <cit.> in the case of a linear regression. We generate data with $n=100$ observations on $p=70$ independent predictors from a standard normal distribution and a treatment predictor from a Bernoulli distribution. $10$ of the standard normal predictors and the treatment predictor were used to generate the noisy target, $y$, at each SBC iteration. We perform projection predictive variable selection using a $10$-fold cross-validation over the entire procedure (including forward search up to $p_{\text{max}} = 30$, and using the <ref> selection heuristic). We forced the treatment predictor to be included in the final submodel across all SBC iterations, and compared the calibration for its regression coefficient in the reference model posterior and in the projected posterior of the submodel selected by projection predictive inference. These calibration checks are visualised in Figure <ref>, where the top row belongs to an over-fitting reference model (using independent standard normal priors over the regression coefficients), and the bottom row belongs to a more predictive reference model (fitted using the R2D2 prior).

The calibration of the marginals in the projected posterior depends on the quality of the reference model: if the reference model over-fits the data, as is the case when independent standard normal priors are used, the projected posterior is over-confident. This is seen by the left tail of the PIT values exceeding the 95% simultaneous confidence band, and the right tail of the PIT values dropping far below it. In the case of the R2D2 sparsifying prior, the projected posterior is well calibrated.

As such, should we wish to perform post-selection inference using the projected posterior, we understand that SBC checks can appropriately diagnose over-confident projected posteriors.

Here, we presented results also for the “naïve” Gaussian prior reference model, but only for illustrative purposes. In a real-world application, we emphasise how important it is to build a sensible reference model in the first instance. In practice, different reference models are often compared in terms of their elpd, which can be interpreted as a proxy for parameter recovery in causally consistent GLMs [83]. In the present example, across 10 realisations of the data, the elpd of the R2D2 prior reference model was on average 33.7 higher than that of the Gaussian prior reference model (with an average standard error of 5.7). This indicates that the R2D2 prior reference model has considerably more predictive power than the Gaussian prior reference model for these data, and is therefore a more suitable reference model. Moreover, we can hypothesise that the gain in predictive power from using the R2D2 prior reference model positively impacts the calibration of the projected posterior, making it safer to use in post-selection inference.

§ DISCUSSION

It is worth reiterating the motivations for model selection. Namely, we do not consider model selection a remedy to over-fitting. Within a Bayesian framework, the statistician has the ability to mitigate over-fitting through the sensible use of priors. While we don't wish to distract from the core aim of this paper with a comprehensive discussion on sparsifying priors, we note that the spike-and-slab prior [56], regularised horseshoe prior [75], R2D2 prior [112, 106], and more recently the R2D2M2 [2] and the $L_1$-ball [105] priors have found success in explicating the statistician's desire for sparsity and consistent behavior with increasing dimensionality. Instead, we advertise model selection as a method to reduce measurement cost in future predictions, as a robustification technique against unrepresentative data, and to gain an understanding of predictor correlations.

Alongside projpred, the kulprit package [52] provides a projection predictive inference implementation in Python.

§.§ Future directions

While our proposed workflow mitigates many of the dangers associated with cross-validation and model selection in a Bayesian workflow, there remain some open pits one might fall into. For instance, in the case of forward search, adding a predictor to a model that does not improve its theoretical predictive performance can lead to overly-optimistic model selection in finite data regimes, resulting in selection-induced bias along the solution path. Future work may task to quantify this risk, and diagnose when the user is at risk with safe and robust stopping criteria.

The models we have considered thus far have remained relatively simple and all with univariate target. Indeed, an extension to the multivariate case presents itself as an interesting avenue to explore in time.

§.§.§ Acknowledgments

We acknowledge the computational resources provided by the Aalto Science-IT project. This paper was partially supported by the Research Council of Finland Flagship programme: Finnish Center for Artificial Intelligence, and Research Council of Finland project “Safe iterative model building” (340721).

[1]
Afrabandpey, H., Peltola, T., Piironen, J., Vehtari, A., and Kaski, S. (2020).
A decision-theoretic approach for model interpretability in
Bayesian framework.
Machine learning, 109(9):1855–1876.

[2]
Aguilar, J. E. and Bürkner, P.-C. (2022).
Intuitive joint priors for bayesian linear multilevel models: The
R2-D2-M2 prior.

[3]
Amirkhiz, R. G., Dixon, M. D., Palmer, J. S., and Swanson, D. L. (2021).
Investigating niches and distribution of a rare species in a
hierarchical framework: Virginia’s Warbler (Leiothlypis virginiae) at
its northeastern range limit.
Landscape Ecology, 36(4):1039–1054.

[4]
Bartonicek, A., Wickham, S. R., Pat, N., and Conner, T. S. (2021).
The value of Bayesian predictive projection for variable
selection: An example of selecting lifestyle predictors of young adult
BMC Public Health, 21(1):1–13.

[5]
Bates, D., Mächler, M., Bolker, B., and Walker, S. (2015).
Fitting linear mixed-effects models using lme4.
Journal of Statistical Software, 67(1).

[6]
Bernardo, J. M. and Rueda, R. (2002).
Bayesian hypothesis testing: a reference approach.
International Statistical Review, 70(3):351–372.

[7]
Bohn, M., Eckert, J., Hanus, D., Lugauer, B., Holtmann, J., and Haun, D. B. M.
Great ape cognition is structured by stable cognitive abilities and
predicted by developmental conditions.
Nature Ecology & Evolution.

[8]
Bratulic, S., Limeta, A., Dabestani, S., Birgisson, H., Enblad, G., Stålberg,
K., Hesselager, G., Häggman, M., Höglund, M., Simonson, O. E., Stålberg,
P., Lindman, H., Bång-Rudenstam, A., Ekstrand, M., Kumar, G., Cavarretta,
I., Alfano, M., Pellegrino, F., Mandel-Clausen, T., Salanti, A., Maccari, F.,
Galeotti, F., Volpi, N., Daugaard, M., Belting, M., Lundstam, S., Stierner,
U., Nyman, J., Bergman, B., Edqvist, P.-H., Levin, M., Salonia, A.,
Kjölhede, H., Jonasch, E., Nielsen, J., and Gatto, F. (2022).
Noninvasive detection of any-stage cancer using free
Proceedings of the National Academy of Sciences,

[9]
Bürkner, P.-C. (2017).
brms: An R package for Bayesian multilevel models using Stan.
Journal of Statistical Software, 80(1):1–28.

[10]
Carbonetto, P. and Stephens, M. (2012).
Scalable variational inference for Bayesian variable selection in
regression, and its accuracy in genetic association studies.
Bayesian Anal, 7(1):73–108.

[11]
Catalina, A., Bürkner, P.-C., and Vehtari, A. (2022).
Projection predictive inference for generalized linear and additive
multilevel models.
In Camps-Valls, G., Ruiz, F. J. R., and Valera, I., editors, Proceedings of The 25th International Conference on Artificial
Intelligence and Statistics, volume 151 of Proceedings of
Machine Learning Research, pages 4446–4461. PMLR.

[12]
Catalina, A., Bürkner, P., and Vehtari, A. (2021).
Latent space projection predictive inference.
arXiv:2109.04702 [stat].

[13]
Catalina, A., Bürkner, P.-C., and Vehtari, A. (2020).
Projection predictive inference for generalized linear and additive
multilevel models.
arXiv:2010.06994 [stat].

[14]
Cheng, L., Ramchandran, S., Vatanen, T., Lietzén, N., Lahesmaa, R., Vehtari,
A., and Lähdesmäki, H. (2019).
An additive Gaussian process regression model for interpretable
non-parametric analysis of longitudinal data.
Nature Communications, 10(1).

[15]
Chipman, H. A., George, E. I., and McCulloch, R. E. (2010).
BART: Bayesian additive regression trees.
The Annals of Applied Statistics, 4(1).

[16]
Cinelli, C., Forney, A., and Pearl, J. (2022).
A crash course in good and bad controls.
Sociological Methods & Research, page 00491241221099552.
Publisher: SAGE Publications Inc.

[17]
Clyde, M. (2022).
BAS: Bayesian variable selection and model averaging using
Bayesian adaptive sampling.
R package, version 1.6.4.

[18]
Cortez, P. and Silva, A. M. G. (2008).
Using data mining to predict secondary school student performance.

[19]
Csárdi, G. (2019).
cranlogs: Download logs from the 'RStudio' 'CRAN'
R package, version 2.1.1.

[20]
Dean, R. T. and Taylor, J. R. (2023).
Continuous affect responses to a large diverse set of unfamiliar
music: Bayesian time-series and cluster analyses.
Psychomusicology: Music, Mind, and Brain.

[21]
Digby, A., Eason, D., Catalina, A., Lierz, M., Galla, S., Urban, L., Le Lec,
M. F., Guhlin, J., Steeves, T. E., Dearden, P. K., Joustra, T., Lees, C.,
Davis, T., Vercoe, D., and Kākāpō Recovery Team (2023).
Hidden impacts of conservation management on fertility of the
critically endangered kākāpō.
PeerJ, 11:e14675.

[22]
Dupuis, J. and Robert, C. (2003).
Variable selection in qualitative models via an entropic explanatory
Journal of Statistical Planning and Inference, 111:77–94.

[23]
Efron, B. (2008).
Microarrays, empirical bayes and the two-groups model.
Statistical Science, 23(1).
arXiv:0808.0572 [stat].

[24]
Efron, B. (2010).
Large-Scale Inference: Empirical Bayes Methods for Estimation,
Testing, and Prediction.
Cambridge University Press, 1 edition.

[25]
Eilers, P. H. C. (2003).
A perfect smoother.
Analytical Chemistry, 75(14):3631–3636.

[26]
Eilers, P. H. C. (2005).
Unimodal smoothing.
Journal of Chemometrics, 19(5-7):317–328.

[27]
Eilers, P. H. C. and Marx, B. D. (1996).
Flexible smoothing with B-splines and penalties.
Statistical Science, 11(2).

[28]
Faraggi, D., LeBlanc, M., and Crowley, J. (2001).
Understanding neural networks using regression trees: an application
to multiple myeloma survival data.
Statistics in Medicine, 20(19):2965–2976.

[29]
Gabry, J. and Mahr, T. (2022).
bayesplot: Plotting for bayesian models.
R package version 1.9.0.

[30]
Gabry, J., Simpson, D., Vehtari, A., Betancourt, M., and Gelman, A. (2019).
Visualization in bayesian workflow.
J. R. Stat. Soc. A, 182:389–402.

[31]
Galatenko, V. V., Shkurnikov, M. Y., Samatov, T. R., Galatenko, A. V.,
Mityakina, I. A., Kaprin, A. D., Schumacher, U., and Tonevitsky, A. G.
Highly informative marker sets consisting of genes with low
individual degree of differential expression.
Scientific Reports, 5(1):14967.

[32]
Garcia-Donato, G. and Forte, A. (2018).
Bayesian testing, variable selection and model averaging in
linear models using R with BayesVarSel.
R J, 10(1):155–174.

[33]
Gelman, A. and Rubin, D. B. (1992).
Inference from iterative simulation using multiple sequences.
Statistical Science, 7(4).

[34]
Gelman, A., Vehtari, A., Simpson, D., Margossian, C. C., Carpenter, B., Yao,
Y., Kennedy, L., Gabry, J., Bürkner, P.-C., and Modrák, M. (2020).
Bayesian workflow.

[35]
George, E. I. and McCulloch, R. E. (1993).
Variable selection via gibbs sampling.
Journal of the American Statistical Association,

[36]
Goto, S., Mori, H., Uchiyama, K., Ishizuka, W., Taneda, H., Kono, M.,
Kajiya-Kanegae, H., and Iwata, H. (2021).
Genetic dissection of growth and eco-physiological traits associated
with altitudinal adaptation in sakhalin fir (abies sachalinensis) based on
qtl mapping.
Genes, 12(8):1110.

[37]
Goutis, C. and Robert, C. (1998).
Model choice in generalised linear models: A Bayesian approach via
Kullback-Leibler projections.
Biometrika, 85:29–37.

[38]
Greenop, A., Cook, S. M., Wilby, A., Pywell, R. F., and Woodcock, B. A. (2020).
Invertebrate community structure predicts natural pest control
resilience to insecticide exposure.
Journal of Applied Ecology, 57(12):2441–2453.

[39]
Hahn, P. R. and Carvalho, C. M. (2014).
Decoupling shrinkage and selection in Bayesian linear models: a
posterior summary perspective.
arXiv:1408.0464 [stat].

[40]
Harrell, F. E. (2001).
Regression Modeling Strategies: With Applications to Linear
Models, Logistic Regression, and Survival Analysis.
Springer Series in Statistics. Springer New York, New York, NY.

[41]
Johnstone, I. M. and Silverman, B. W. (2004).
Needles and straw in haystacks: Empirical Bayes estimates of
possibly sparse sequences.
The Annals of Statistics, 32(4).

[42]
Kallioinen, N. (2022).
priorsense: an R package for prior diagnostics and sensitivity.

[43]
Kallioinen, N., Paananen, T., Bürkner, P.-C., and Vehtari, A. (2022).
Detecting and diagnosing prior and likelihood sensitivity with

[44]
Karsdorp, F. and Fonteyn, L. (2019).
Cultural entrenchment of folktales is encoded in language.
Palgrave Communications, 5(1):25.

[45]
Kuck, K., Naik, B. I., Domino, K. B., Posner, K. L., Saager, L., Stuart, A. R.,
Johnson, K. B., Alpert, S. B., Durieux, M. E., Sinha, A. K., Brummett, C. M.,
Aziz, M. F., Cummings, K. C., Gaudet, J. G., Kurz, A., Rijsdijk, M.,
Wanderer, J. P., Pace, N. L., the Multicenter Perioperative Outcomes Group
Enhanced Observation Study Investigator Group for the Multicenter
Perioperative Outcomes Group Enhanced Observation Study Collaborator Group,
Berman, M. F., Corradini, B. R., Mincer, S. L., Rose, S. E., van Klei, W. A.,
Edwards, D. A., Lamers, O. O., and Vaughn, M. T. (2023).
Prolonged opioid use and pain outcome and associated factors after
surgery under general anesthesia: A prospective cohort association
multicenter study.
Anesthesiology, 138(5):462–476.

[46]
Kullback, S. and Leibler, R. A. (1951).
On information and sufficiency.
The Annals of Mathematical Statistics, 22(1):79 – 86.

[47]
Kumar, R., Carroll, C., Hartikainen, A., and Martin, O. (2019).
Arviz a unified library for exploratory analysis of bayesian models
in python.
Journal of Open Source Software, 4(33):1143.

[48]
Leamer, E. E. (1979).
Information criteria for choice of regression models: A comment.
Econometrica, 47(2):507–510.

[49]
Leung Soo, C., Bhatnagar, S., Bartlett, S. J., Esmail, A., Dheda, K., and
Pant Pai, N. (2023).
Development and evaluation of a digital hiv risk assessment tool
incorporated within an app-based self-testing program.
JAIDS Journal of Acquired Immune Deficiency Syndromes, Publish
Ahead of Print.

[50]
Lindley, D. V. (1968).
The choice of variables in multiple regression.
Journal of the Royal Statistical Society. Series B
(Methodological), 30(1):31–66.

[51]
Liquet, B. and Sutton, M. (2017).
MBSGS: Multivariate Bayesian sparse group selection with
spike and slab.
R package, version 1.1.0.

[52]
McLatchie, Y., Martin, O., and Capretto, T. (2022a).
Kulprit: Kullback-Leibler projections for Bayesian model
selection in Python.

[53]
McLatchie, Y., Matamoros, A. A., Kohns, D., and Vehtari, A. (2022b).
Bayesian order identification of ARMA models with projection
predictive inference.

[54]
Mercer, J. J., Liefert, D. T., and Williams, D. G. (2020).
Atmospheric vapour and precipitation are not in isotopic equilibrium
in a continental mountain environment.
Hydrological Processes, 34(14):3078–3101.

[55]
Milne, A. J., Dean, R., and Bulger, D. (2021).
The effects of rhythmic structure on tapping accuracy.
preprint, PsyArXiv.

[56]
Mitchell, T. J. and Beauchamp, J. J. (1988).
Bayesian variable selection in linear regression.
Journal of the American Statistical Association,

[57]
Modrák, M., Moon, A. H., Kim, S., Bürkner, P., Huurre, N., Faltejsková, K.,
Gelman, A., and Vehtari, A. (2022).
Simulation-based calibration checking for bayesian computation: The
choice of test quantities shapes sensitivity.
arXiv:2211.02383 [stat].

[58]
Mori, H., Ueno, S., Ujino-Ihara, T., Fujiwara, T., Yamashita, K., Kanetani, S.,
Endo, R., Matsumoto, A., Uchiyama, K., Yoshida, T., Sakai, Y., Moriguchi, Y.,
Kusano, R., and Tsumura, Y. (2022).
Genotype-by-environment interaction and genetic dissection of
heartwood color in Cryptomeria japonica based on multiple common gardens
and quantitative trait loci mapping.
PLOS ONE, 17(7):e0270522.

[59]
Ng, A. M. B., Pontius, M. L., De Ruiter, S. L., and Proppe, D. S. (2020).
Noise, avian abundance, and productivity at banding stations across
the Continental United States.
Avian Conservation and Ecology, 15(2):art4.

[60]
Nikooienejad, A. and Johnson, V. E. (2020).
BVSNLP: Bayesian variable selection in high dimensional
settings using nonlocal priors.
R package, version 1.1.9.

[61]
Nilsson, N. J. (1998).
Artificial Intelligence: A New Synthesis.

[62]
Nordström, H. and Laukka, P. (2019).
The time course of emotion recognition in speech and music.
The Journal of the Acoustical Society of America,

[63]
Nott, D. J. and Leng, C. (2010).
Bayesian projection approaches to variable selection in generalized
linear models.
Computational Statistics & Data Analysis, 54(12):3227–3241.

[64]
Ntzoufras, I., Forster, J. J., and Dellaportas, P. (2000).
Stochastic search variable selection for log-linear models.
Journal of Statistical Computation and Simulation,

[65]
Paasiniemi, M. (2018).
Methods and Tools for Interpretable Bayesian Variable
Master's thesis.

[66]
Paul, D., Bair, E., Hastie, T., and Tibshirani, R. (2008).
“Preconditioning” for feature selection and regression in
high-dimensional problems.
The Annals of Statistics, 36(4).

[67]
Pavone, F., Piironen, J., Bürkner, P.-C., and Vehtari, A. (2022).
Using reference models in variable selection.
Computational Statistics.

[68]
Pearl, J. (2009).
Causal inference in statistics: An overview.
Statistics Surveys, 3:96–146.
Publisher: Amer. Statist. Assoc., the Bernoulli Soc., the Inst. Math.
Statist., and the Statist. Soc. Canada.

[69]
Penrose, K. W., Nelson, A. G., and Fisher, A. G. (1985).
Generalized body composition prediction equation for men using simple
measurement techniques.
Medicine & Science in Sports & Exercise, 17(2):189.

[70]
Perez, P. and de los Campos, G. (2014).
Genome-wide regression and prediction with the BGLR statistical
Genetics, 198(2):483–495.

[71]
Piironen, J., Paasiniemi, M., Catalina, A., Weber, F., and Vehtari, A. (2023).
projpred: Projection predictive feature selection.
R package, version 2.5.0.

[72]
Piironen, J., Paasiniemi, M., and Vehtari, A. (2020).
Projective inference in high-dimensional problems: Prediction and
feature selection.
Electronic Journal of Statistics, 14(1):2155–2197.

[73]
Piironen, J. and Vehtari, A. (2016).
Projection predictive model selection for Gaussian processes.
In 2016 IEEE 26th International Workshop on Machine Learning
for Signal Processing (MLSP), pages 1–6. IEEE.

[74]
Piironen, J. and Vehtari, A. (2017a).
Comparison of Bayesian predictive methods for model selection.
Statistics and Computing, 27(3):711–735.

[75]
Piironen, J. and Vehtari, A. (2017b).
Sparsity information and regularization in the horseshoe and other
shrinkage priors.
Electronic Journal of Statistics, 11(2).

[76]
Piironen, J. and Vehtari, A. (2018).
Iterative supervised principal components.
In Storkey, A. and Perez-Cruz, F., editors, Proceedings of the
Twenty-First International Conference on Artificial Intelligence and
Statistics, volume 84 of Proceedings of Machine Learning Research,
pages 106–114. PMLR.

[77]
Pya, N. (2022).
scam: Shape constrained additive models.

[78]
Ramsay, J. O. (1988).
Monotone Regression Splines in Action.
Statistical Science, 3(4).

[79]
Rossell, D., Cook, J. D., Telesca, D., Roebuck, P., Abril, O., and Torrens, M.
mombf: Model selection with Bayesian methods and
information criteria.
R package, version 3.3.1.

[80]
Rosso, V., Albrecht, J., Roberts, F., and Janey, J. M. (2019).
Uniting laboratory automation, DoE data, and modeling techniques to
accelerate chemical process development.
Reaction Chemistry & Engineering, 4(9):1646–1657.

[81]
San Martini, A. and Spezzaferri, F. (1984).
A predictive model selection criterion.
Journal of the Royal Statistical Society: Series B
(Methodological), 46(2):296–303.

[82]
Scheipl, F. (2011).
spikeSlabGAM: Bayesian variable selection, model choice and
regularization for generalized additive mixed models in R.
J Stat Softw, 43(14):1–24.

[83]
Scholz, M. and Bürkner, P.-C. (2022).
Prediction can be safely used as a proxy for explanation in causally
consistent Bayesian generalized linear models.
arXiv:2210.06927 [stat].

[84]
Schrittenlocher, S., Bachmann, B., Tiurbe, A. M., Tuac, O., Velten, K.,
Schmidt, D., and Cursiefen, C. (2019).
Impact of preoperative visual acuity on descemet membrane endothelial
keratoplasty (DMEK) outcome.
Graefe's Archive for Clinical and Experimental Ophthalmology,

[85]
Sivula, T., Magnusson, M., Matamoros, A. A., and Vehtari, A. (2022).
Uncertainty in bayesian leave-one-out cross-validation based model
arXiv:2008.10296 [stat].

[86]
Soliman, S. A., Haque, A., Vanarsa, K., Zhang, T., Ismail, F., Lee, K. H.,
Pedroza, C., Greenbaum, L. A., Mason, S., Hicks, M. J., Wenderfer, S. E., and
Mohan, C. (2022).
Urine ALCAM, PF4 and VCAM-1 surpass conventional metrics in
identifying nephritis disease activity in childhood-onset systemic lupus
Frontiers in Immunology, 13:885307.

[87]
South, L. F., Riabiz, M., Teymur, O., and Oates, C. J. (2022).
Postprocessing of MCMC.
Annual Review of Statistics and Its Application, 9(1):529–555.

[88]
Stell, L. and Sabatti, C. (2015).
ptycho: Bayesian variable selection with hierarchical priors.
R package, version 1.1-4.

[89]
Sullivan, B., Barker, E., Williams, P., MacGregor, L., Bhamber, R., Thomas, M.,
Gurney, S., Hyams, C., Whiteway, A., Cooper, J. A., McWilliams, C., Turner,
K., Dowsey, A. W., and Albur, M. (2022).
Conventional and Bayesian workflows for clinical prediction
modelling of severe Covid-19 outcomes based on clinical biomarker test
results: LabMarCS: Laboratory Markers of COVID-19 Severity -
Bristol Cohort.
preprint, Infectious Diseases (except HIV/AIDS).

[90]
Säilynoja, T., Bürkner, P.-C., and Vehtari, A. (2021).
Graphical test for discrete uniformity and its applications in
goodness of fit evaluation and multiple sample comparison.
arXiv:2103.10522 [stat].

[91]
Talts, S., Betancourt, M., Simpson, D., Vehtari, A., and Gelman, A. (2020).
Validating bayesian inference algorithms with simulation-based
arXiv:1804.06788 [stat].
arXiv: 1804.06788.

[92]
Teng, K. T.-y., Martinez Avilés, M., Ugarte-Ruiz, M., Barcena, C., de la
Torre, A., Lopez, G., Moreno, M. A., Dominguez, L., and Alvarez, J. (2020).
Spatial trends in salmonella infection in pigs in spain.
Frontiers in Veterinary Science, 7:345.

[93]
Tibshirani, R. (2011).
Regression shrinkage selection via the lasso.
Journal of the Royal Statistical Society Series B, 73:273–282.

[94]
Tran, M.-N., Nott, D., and Leng, C. (2010).
The predictive lasso.
arXiv:1009.2302 [stat].

[95]
Trier, H., O'Reilly, J., Spiering, L., Ma, S., Kolling, N., Rushworth, M., and
Scholl, J. (2023).
Emotions and individual differences shape foraging under threat.
preprint, PsyArXiv.

[96]
Vanarsa, K., Castillo, J., Wang, L., Lee, K. H., Pedroza, C., Lotan, Y., and
Mohan, C. (2023).
Comprehensive proteomics and platform validation of urinary
biomarkers for bladder cancer diagnosis and staging.
BMC Medicine, 21(1):133.

[97]
Vehtari, A. (2021).
On reparameterization invariant Bayesian point estimates and
credible regions.
arXiv preprint arXiv:2109.10843.

[98]
Vehtari, A., Gabry, J., Magnusson, M., Yao, Y., Bürkner, P.-C., Paananen, T.,
and Gelman, A. (2022).
loo: Efficient leave-one-out cross-validation and waic for Bayesian
R package version 2.5.0.

[99]
Vehtari, A., Gelman, A., and Gabry, J. (2016).
Practical Bayesian model evaluation using leave-one-out
cross-validation and WAIC.
Statistics and Computing, 27(5):1413–1432.

[100]
Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Bürkner, P.-C.
Rank-normalization, folding, and localization: An improved $\widehat
R$ for assessing convergence of MCMC (with discussion).
Bayesian Analysis, 16(2):667–718.

[101]
Vehtari, A. and Ojanen, J. (2012).
A survey of Bayesian predictive methods for model assessment,
selection and comparison.
Statistics Surveys, 6(none):142 – 228.

[102]
Vehtari, A., Simpson, D., Gelman, A., Yao, Y., and Gabry, J. (2015).
Pareto smoothed importance sampling.

[103]
Wang, L. and Yoon, K.-J. (2022).
Knowledge distillation and student-teacher learning for visual
intelligence: A review and new outlooks.
IEEE Transactions on Pattern Analysis and Machine Intelligence,
arXiv:2004.05937 [cs].

[104]
Weber, F. and Vehtari, A. (2023).
Projection predictive variable selection for discrete response
families with finite support.

[105]
Xu, M. and Duan, L. L. (2020).
Bayesian inference with the l1-ball prior: Solving combinatorial
problems with exact zeros.

[106]
Yanchenko, E., Bondell, H. D., and Reich, B. J. (2021).
The R2D2 prior for generalized linear mixed models.

[107]
Yanchenko, E., Bondell, H. D., and Reich, B. J. (2023).
R2D2 goes to space! a principled approach to setting prior
distributions on spatial parameters.
arXiv:2301.09951 [stat].

[108]
Yao, Y., Carvalho, L. M., Mesquita, D., and McLatchie, Y. (2023).
Locking and Quacking: Stacking Bayesian model predictions by
log-pooling and superposition.
arXiv:2305.07334 [cs, stat].

[109]
Yao, Y., Pirš, G., Vehtari, A., and Gelman, A. (2022).
Bayesian hierarchical stacking: Some models are (somewhere) useful.
Bayesian Analysis, 17(4).

[110]
Yao, Y., Vehtari, A., Simpson, D., and Gelman, A. (2018).
Using stacking to average Bayesian predictive distributions (with
Bayesian Analysis, 13(3):917–1007.

[111]
Zhang, Y. and Peleato, N. M. (2021).
Predicting cyanobacteria abundance with Bayesian zero-inflated
negative binomial models.
SSRN Electronic Journal.

[112]
Zhang, Y. D., Naughton, B. P., Bondell, H. D., and Reich, B. J. (2022).
Bayesian regression using a prior on the model fit: The r2-d2
shrinkage prior.
Journal of the American Statistical Association,

[113]
Zhao, Z., Banterle, M., Bottolo, L., Richardson, S., Lewin, A., and Zucknick,
M. (2021).
BayesSUR: An R package for high-dimensional multivariate
Bayesian variable and covariance selection in linear regression.
J Stat Softw, 100(11):1–32.

[114]
Zhou, G., Lee, M.-C., Wang, X., Zhong, D., and Yan, G. (2022).
Mapping potential malaria vector larval habitats for larval source
management: Introduction to multi-model ensembling approaches.
preprint, Ecology.

[115]
Zou, H. (2006).
The adaptive lasso and its oracle properties.
Journal of the American Statistical Association,

§ APPLICATIONS OF PROJECTION PREDICTIVE INFERENCE

Projection predictive inference has found application in many fields, including medicine, psychology, and life sciences <cit.>, chemistry <cit.>, hydrology <cit.>, and cultural studies and linguistics <cit.>.

Alongside projpred, there exist other model selection packages in R: BayesVarSel [32]; BAS [17]; varbvs [10]; spikeSlabGAM [82]; BVSNLP [60], ptycho [88]; BayesSUR [113]; BGLR [70]; MBSGS [51]; and mombf [79]. A check of cranlogs [19] revealed that projpred is the most popular of these in terms of number of downloads (last check: April 26, 2023).

§ THE R2D2 PRIOR

Case study     Dataset     $n$     $p$     $\xi$     $\mu_{R^2}$     $\varphi_{R^2}$
Section <ref>     Highly-correlated     $500$     $100$     $1$     $0.3$     $5$

Section <ref>     Weakly-relevant     $500$     $50$     $10$     $0.3$     $5$

Section <ref>     Prostate cancer     $102$     $5966$     $1$     $0.3$     $5$

Section <ref>     Linear regression     $100$     $60$     $0.2$     $0.3$     $20$

R2D2 prior hyperparameter values used across case studies.

Fitting high-dimensional models on relatively sparse data is a difficult task in general. Joint shrinkage priors implore the practitioner to encode a prior jointly over the predictor-related parameter space and the predictive space. The R2D2 prior [112, 106, 2, 107] allows the statistician to encode some prior belief on the model's $R^2$ and the number of predictors necessary to achieve it. In doing so, one mitigates the risk of over-fitting the data present when using more simple independent Gaussians over the regression coefficients. Formally, the R2D2 model as we use it is defined as:

y_i     ∼(β_0 + ∑_k=1^px_k,iβ_k,σ^2)

β_0     ∼_3(0, 2.5)

σ    ∼^+_3(0, 2.5)

β_k     ∼(0, σ^2 τ^2 ϕ_k)

R^2     ∼(μ_R^2, φ_R^2)

ϕ    ∼(ξ,…,ξ)

τ^2     = R^2/1-R^2.

By our notation, $\studentt_\nu(\mu, \sigma)$ denotes Student's $t$-distribution with $\nu$ degrees of freedom, location $\mu$, and scale $\sigma$. The super-script $+$ indicates a half-Student-$t$ distribution. Note that we use the mean and a pseudo-precision to parameterise the beta distribution for $R^2$ as opposed to the conventional shape parameters $a > 0$ and $b > 0$. The relationship between the two parameterisations is as follows:

μ_R^2     = a/a+b,

φ_R^2     = a + b.

The hyperparameters used throughout the three case studies are tabulated in Table <ref>. The Dirichlet concentration parameter $\xi$ (which we choose to be constant across dimensions) controls sparsity in the regression coefficients: larger $\xi$ imply more uniform regression coefficients (in magnitude), typically with more predictors explaining the variance of the response. In Section <ref>, when assessing the calibration, the R2D2 prior was applied to all predictors except the one for which the calibration was being assessed. This predictor under inspection was given a standard normal prior.